| blob | 789ce9afc07e49e7518eec0c6605af3643bfdea9 |
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
1379 itf=MIN(ite,ide-1)
1380 jtf=MIN(jte,jde-1)
1381 ktf=MIN(kte,kde-1)
1383 ! Notes on map scale factors (WCS, 2 march 2008)
1384 ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1385 ! -(1/my) my mu v d/dy(phi)
1386 ! - omega d/d_eta(phi)
1387 ! +mu g w/my
1388 !
1389 ! A little further explanation...
1390 ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1391 ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1392 ! solution is computed in subourine advance_w. The formulation dates from the early
1393 ! days of the mass coordinate model when we were testing both a flux and an advective formulation
1394 ! for the geopotential equation and different forms of the vertical momentum equation and the
1395 ! vertically implicit solver.
1397 ! advective form for the geopotential equation
1399 DO j = jts, jtf
1401 DO k = 2, kte
1402 DO i = its, itf
1403 wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
1404 *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1405 ENDDO
1406 ENDDO
1408 ! RHS term 3 is: - omega partial d/dnu(phi)
1410 DO k = 2, kte-1
1411 DO i = its, itf
1412 ph_tend(i,k,j) = ph_tend(i,k,j) &
1413 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1414 ENDDO
1415 ENDDO
1417 ENDDO
1419 IF (non_hydrostatic) THEN ! add in "gw" term.
1420 DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1421 ! after the timestep to give us "w"
1422 DO i = its, itf
1423 ph_tend(i,kde,j) = 0.
1424 ENDDO
1426 DO k = 2, kte
1427 DO i = its, itf
1428 ! phi equation RHS term 4
1429 ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1430 ENDDO
1431 ENDDO
1433 ENDDO
1435 END IF
1437 ! Notes on map scale factors:
1438 ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1439 ! -(1/my) my v mu partial d/dy(phi)
1441 IF (advective_order <= 2) THEN
1443 ! y (v) advection
1445 i_start = its
1446 j_start = jts
1447 itf=MIN(ite,ide-1)
1448 jtf=MIN(jte,jde-1)
1450 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1451 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
1453 DO j = j_start, jtf
1455 DO k = 2, kte-1
1456 DO i = i_start, itf
1457 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1458 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1459 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1460 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1461 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1462 ENDDO
1463 ENDDO
1465 k = kte
1466 DO i = i_start, itf
1467 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1468 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1469 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1470 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1471 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1472 ENDDO
1474 ENDDO
1476 ! x (u) advection
1478 i_start = its
1479 j_start = jts
1480 itf=MIN(ite,ide-1)
1481 jtf=MIN(jte,jde-1)
1483 IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1484 IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
1486 DO j = j_start, jtf
1488 DO k = 2, kte-1
1489 DO i = i_start, itf
1490 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1491 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1492 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1493 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1494 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1495 ENDDO
1496 ENDDO
1498 k = kte
1499 DO i = i_start, itf
1500 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1501 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1502 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1503 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1504 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1505 ENDDO
1507 ENDDO
1509 ELSE IF (advective_order <= 4) THEN
1511 ! y (v) advection
1513 i_start = its
1514 j_start = jts
1515 itf=MIN(ite,ide-1)
1516 jtf=MIN(jte,jde-1)
1518 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1519 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
1521 DO j = j_start, jtf
1523 DO k = 2, kte-1
1524 DO i = i_start, itf
1525 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*( &
1526 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1527 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
1528 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1529 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1530 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1531 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1534 ENDDO
1535 ENDDO
1537 k = kte
1538 DO i = i_start, itf
1539 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
1540 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1541 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
1542 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1543 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1544 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1545 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1547 ENDDO
1549 ENDDO
1551 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1553 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
1555 j = jds+1
1556 DO k = 2, kte-1
1557 DO i = i_start, itf
1558 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1559 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1560 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1561 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1562 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1563 ENDDO
1564 ENDDO
1566 k = kte
1567 DO i = i_start, itf
1568 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1569 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1570 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1571 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1572 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1573 ENDDO
1575 END IF
1577 IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
1579 j = jde-2
1580 DO k = 2, kte-1
1581 DO i = i_start, itf
1582 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1583 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1584 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1585 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1586 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1587 ENDDO
1588 ENDDO
1590 k = kte
1591 DO i = i_start, itf
1592 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1593 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1594 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1595 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1596 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1597 ENDDO
1599 END IF
1601 ! x (u) advection
1603 i_start = its
1604 j_start = jts
1605 itf=MIN(ite,ide-1)
1606 jtf=MIN(jte,jde-1)
1608 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1609 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
1611 DO j = j_start, jtf
1613 DO k = 2, kte-1
1614 DO i = i_start, itf
1615 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1616 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1617 +muu(i,j )*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
1618 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1619 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1620 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1621 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1622 ENDDO
1623 ENDDO
1625 k = kte
1626 DO i = i_start, itf
1627 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1628 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1629 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
1630 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1631 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1632 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1633 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1634 ENDDO
1636 ENDDO
1638 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1640 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1642 i = ids + 1
1644 DO j = j_start, jtf
1645 DO k = 2, kte-1
1646 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1647 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1648 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1649 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1650 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1651 ENDDO
1652 ENDDO
1654 k = kte
1655 DO j = j_start, jtf
1656 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1657 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1658 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1659 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1660 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1661 ENDDO
1663 END IF
1665 IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1667 i = ide-2
1668 DO j = j_start, jtf
1669 DO k = 2, kte-1
1670 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1671 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1672 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1673 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1674 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1675 ENDDO
1676 ENDDO
1678 k = kte
1679 DO j = j_start, jtf
1680 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1681 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1682 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1683 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1684 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1685 ENDDO
1687 END IF
1689 !--------------------------
1691 ELSE IF (advective_order <= 6) THEN
1693 !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1694 !!! the branches covering the other advective_order cases have not. 20090923. JM
1696 ! y (v) advection
1698 i_start = its
1699 j_start = jts
1700 itf=MIN(ite,ide-1)
1701 jtf=MIN(jte,jde-1)
1703 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1704 IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
1706 DO j = j_start, jtf
1708 DO k = 2, kte-1
1709 DO i = i_start, itf
1710 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1711 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1712 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
1713 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1714 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1715 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1716 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1717 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1718 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1721 ENDDO
1722 ENDDO
1724 k = kte
1725 DO i = i_start, itf
1726 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1727 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1728 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
1729 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1730 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1731 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1732 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1733 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1734 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1736 ENDDO
1738 ENDDO
1740 ! 4th order stencil for open or specified conditions two in form the boundary
1742 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
1744 j = jds+2
1745 DO k = 2, kte-1
1746 DO i = i_start, itf
1747 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1748 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1749 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
1750 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1751 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1752 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1753 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1756 ENDDO
1757 ENDDO
1759 k = kte
1760 DO i = i_start, itf
1761 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1762 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1763 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1764 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1765 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1766 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1767 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1769 ENDDO
1771 END IF
1773 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
1774 j = jde-3
1775 DO k = 2, kte-1
1776 DO i = i_start, itf
1777 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1778 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1779 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
1780 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1781 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1782 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1783 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1786 ENDDO
1787 ENDDO
1789 k = kte
1790 DO i = i_start, itf
1791 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1792 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1793 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1794 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1795 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1796 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1797 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1799 ENDDO
1801 END IF
1803 ! 2nd order stencil for open and specified conditions one row in from boundary
1805 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
1807 j = jds+1
1808 DO k = 2, kte-1
1809 DO i = i_start, itf
1810 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1811 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1812 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1813 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1814 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1815 ENDDO
1816 ENDDO
1818 k = kte
1819 DO i = i_start, itf
1820 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1821 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1822 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1823 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1824 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1825 ENDDO
1827 END IF
1829 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
1831 j = jde-2
1832 DO k = 2, kte-1
1833 DO i = i_start, itf
1834 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1835 ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1836 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1837 +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1838 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1839 ENDDO
1840 ENDDO
1842 k = kte
1843 DO i = i_start, itf
1844 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1845 ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1846 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1847 +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1848 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1849 ENDDO
1851 END IF
1853 ! x (u) advection
1855 i_start = its
1856 j_start = jts
1857 itf=MIN(ite,ide-1)
1858 jtf=MIN(jte,jde-1)
1860 IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1861 IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
1863 DO j = j_start, jtf
1865 DO k = 2, kte-1
1866 DO i = i_start, itf
1867 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1868 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1869 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
1870 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1871 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1872 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1873 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1874 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1875 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1876 ENDDO
1877 ENDDO
1879 k = kte
1880 DO i = i_start, itf
1881 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1882 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1883 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
1884 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1885 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1886 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1887 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1888 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1889 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1890 ENDDO
1892 ENDDO
1894 ! 4th order stencil two in from the boundary for open and specified conditions
1896 IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1897 i = ids + 2
1898 DO j = j_start, jtf
1899 DO k = 2, kte-1
1900 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1901 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1902 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1903 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1904 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1905 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1906 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1907 ENDDO
1908 k = kte
1909 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1910 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1911 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1912 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1913 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1914 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1915 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1917 ENDDO
1918 END IF
1920 IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
1921 i = ide-3
1922 DO j = j_start, jtf
1923 DO k = 2, kte-1
1924 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1925 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1926 +muu(i,j )*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1927 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1928 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1929 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1930 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1931 ENDDO
1932 k = kte
1933 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1934 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1935 +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1936 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1937 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1938 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1939 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1941 ENDDO
1942 END IF
1944 ! 2nd order stencil for open and specified conditions one in from the boundary
1946 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
1948 i = ids + 1
1950 DO j = j_start, jtf
1951 DO k = 2, kte-1
1952 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1953 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1954 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1955 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1956 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1957 ENDDO
1958 ENDDO
1960 k = kte
1961 DO j = j_start, jtf
1962 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1963 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1964 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1965 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1966 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1967 ENDDO
1969 END IF
1971 IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
1973 i = ide-2
1974 DO j = j_start, jtf
1975 DO k = 2, kte-1
1976 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1977 ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1978 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1979 +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1980 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1981 ENDDO
1982 ENDDO
1984 k = kte
1985 DO j = j_start, jtf
1986 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1987 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1988 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1989 +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1990 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1991 ENDDO
1993 END IF
1995 END IF ! 6th order advection
1997 ! lateral open boundary conditions,
1998 ! start with north and south (y) boundaries
2000 i_start = its
2001 itf=MIN(ite,ide-1)
2003 ! south
2005 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2007 j=jts
2009 DO k=2,kde
2010 kz = min(k,kde-1)
2011 DO i = its,itf
2012 vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
2013 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
2014 vl=amin1(vb,0.)
2015 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2016 +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2017 ENDDO
2018 ENDDO
2020 END IF
2022 ! north
2024 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2026 j=jte-1
2028 DO k=2,kde
2029 kz = min(k,kde-1)
2030 DO i = its,itf
2031 vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
2032 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2033 vr=amax1(vb,0.)
2034 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( &
2035 +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2036 ENDDO
2037 ENDDO
2039 END IF
2041 ! now the east and west (y) boundaries
2043 j_start = its
2044 jtf=MIN(jte,jde-1)
2046 ! west
2048 IF ( (config_flags%open_xs) .and. its == ids ) THEN
2050 i=its
2052 DO j = jts,jtf
2053 DO k=2,kde-1
2054 kz = k
2055 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2056 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2057 ul=amin1(ub,0.)
2058 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2059 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2060 ENDDO
2062 k = kde
2063 kz = k
2064 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2065 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2066 ul=amin1(ub,0.)
2067 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2068 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2069 ENDDO
2071 END IF
2073 ! east
2075 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2077 i = ite-1
2079 DO j = jts,jtf
2080 DO k=2,kde-1
2081 kz = k
2082 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2083 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2084 ur=amax1(ub,0.)
2085 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2086 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2087 ENDDO
2089 k = kde
2090 kz = k-1
2091 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2092 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2093 ur=amax1(ub,0.)
2094 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2095 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2097 ENDDO
2099 END IF
2101 END SUBROUTINE rhs_ph
2104 !-------------------------------------------------------------------------------
2106 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
2107 ph,alt,p,pb,al,php,cqu,cqv, &
2108 muu,muv,mu,fnm,fnp,rdnw, &
2109 cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
2110 msfvx,msfvy,msftx,msfty, &
2111 config_flags, non_hydrostatic, &
2112 top_lid, &
2113 ids, ide, jds, jde, kds, kde, &
2114 ims, ime, jms, jme, kms, kme, &
2115 its, ite, jts, jte, kts, kte )
2117 IMPLICIT NONE
2119 ! Input data
2122 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2124 LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
2126 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2127 ims, ime, jms, jme, kms, kme, &
2128 its, ite, jts, jte, kts, kte
2130 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
2131 ph, &
2132 alt, &
2133 al, &
2134 p, &
2135 pb, &
2136 php, &
2137 cqu, &
2138 cqv
2141 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
2142 ru_tend, &
2143 rv_tend
2145 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
2146 msfux, msfuy, &
2147 msfvx, msfvy, &
2148 msftx, msfty
2150 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
2152 REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3
2154 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2155 REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2156 REAL :: dpx, dpy
2158 LOGICAL :: specified
2160 !<DESCRIPTION>
2161 !
2162 ! horizontal_pressure_gradient calculates the
2163 ! horizontal pressure gradient terms for the large-timestep tendency
2164 ! in the horizontal momentum equations (u,v).
2165 !
2166 !</DESCRIPTION>
2168 specified = .false.
2169 if(config_flags%specified .or. config_flags%nested) specified = .true.
2171 ! Notes on map scale factors:
2172 ! Calculates the pressure gradient terms in ADT eqns 44 and 45
2173 ! With upper rho -> 'mu', these are:
2174 ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
2175 ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
2176 !
2177 ! As we are on nu, rather than height, surfaces:
2178 !
2179 ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2180 ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2181 !
2182 ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2183 ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2185 ! start with the north-south (y) pressure gradient
2187 itf=MIN(ite,ide-1)
2188 jtf=jte
2189 ktf=MIN(kte,kde-1)
2190 i_start = its
2191 j_start = jts
2192 IF ( (config_flags%open_ys .or. specified .or. &
2193 config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2194 IF ( (config_flags%open_ye .or. specified .or. &
2195 config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2197 DO j = j_start, jtf
2199 IF ( non_hydrostatic ) THEN
2201 k=1
2203 DO i = i_start, itf
2204 dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
2205 +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2206 +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2207 dpn(i,kde) = 0.
2208 ENDDO
2209 IF (top_lid) THEN
2210 DO i = i_start, itf
2211 dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2212 +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j)) &
2213 +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j)) )
2214 ENDDO
2215 ENDIF
2217 DO k=2,ktf
2218 DO i = i_start, itf
2219 dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
2220 +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2221 END DO
2222 END DO
2224 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2225 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2226 DO K=1,ktf
2227 DO i = i_start, itf
2228 ! Here are mu dp/dy terms 1-3
2229 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2230 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2231 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2232 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2233 ! Here is mu dp/dy term 4
2234 dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2235 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2236 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2237 END DO
2238 END DO
2240 ELSE
2242 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2243 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2244 DO K=1,ktf
2245 DO i = i_start, itf
2246 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2247 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*( &
2248 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2249 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2250 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2251 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2252 END DO
2253 END DO
2255 END IF
2257 ENDDO
2259 ! now the east-west (x) pressure gradient
2261 itf=ite
2262 jtf=MIN(jte,jde-1)
2263 ktf=MIN(kte,kde-1)
2264 i_start = its
2265 j_start = jts
2266 IF ( (config_flags%open_xs .or. specified .or. &
2267 config_flags%nested ) .and. its == ids ) i_start = its+1
2268 IF ( (config_flags%open_xe .or. specified .or. &
2269 config_flags%nested ) .and. ite == ide ) itf = itf-1
2270 IF ( config_flags%periodic_x ) i_start = its
2271 IF ( config_flags%periodic_x ) itf=ite
2273 DO j = j_start, jtf
2275 IF ( non_hydrostatic ) THEN
2277 k=1
2279 DO i = i_start, itf
2280 dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
2281 +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2282 +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2283 dpn(i,kde) = 0.
2284 ENDDO
2285 IF (top_lid) THEN
2286 DO i = i_start, itf
2287 dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2288 +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j)) &
2289 +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j)) )
2290 ENDDO
2291 ENDIF
2293 DO k=2,ktf
2294 DO i = i_start, itf
2295 dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
2296 +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2297 END DO
2298 END DO
2300 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2301 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2302 DO K=1,ktf
2303 DO i = i_start, itf
2304 ! Here are mu dp/dy terms 1-3
2305 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2306 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2307 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2308 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2309 ! Here is mu dp/dy term 4
2310 dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2311 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2312 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2313 END DO
2314 END DO
2316 ELSE
2318 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2319 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2320 DO K=1,ktf
2321 DO i = i_start, itf
2322 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2323 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*( &
2324 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2325 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2326 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2327 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2328 END DO
2329 END DO
2331 END IF
2333 ENDDO
2335 END SUBROUTINE horizontal_pressure_gradient
2337 !-------------------------------------------------------------------------------
2339 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, &
2340 rdnw, rdn, g, msftx, msfty, &
2341 ids, ide, jds, jde, kds, kde, &
2342 ims, ime, jms, jme, kms, kme, &
2343 its, ite, jts, jte, kts, kte )
2345 IMPLICIT NONE
2347 ! Input data
2349 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2350 ims, ime, jms, jme, kms, kme, &
2351 its, ite, jts, jte, kts, kte
2353 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2354 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2357 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2359 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msftx, msfty
2361 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2363 REAL, INTENT(IN ) :: g
2365 INTEGER :: itf, jtf, i, j, k
2366 REAL :: cq1, cq2
2369 !<DESCRIPTION>
2370 !
2371 ! pg_buoy_w calculates the
2372 ! vertical pressure gradient and buoyancy terms for the large-timestep
2373 ! tendency in the vertical momentum equation.
2374 !
2375 !</DESCRIPTION>
2377 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2379 ! Map scale factor notes
2380 ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2381 ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2382 ! term 6: +(g/my) partial dp'/dnu
2383 ! term 7: -(g/my) mu'
2384 !
2385 ! For moisture-free atmosphere, cq1=1, cq2=0
2386 ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2388 itf=MIN(ite,ide-1)
2389 jtf=MIN(jte,jde-1)
2391 DO j = jts,jtf
2393 k=kde
2394 DO i=its,itf
2395 cq1 = 1./(1.+cqw(i,k-1,j))
2396 cq2 = cqw(i,k-1,j)*cq1
2397 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2398 cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
2399 -mu(i,j)-cq2*mub(i,j) )
2400 END DO
2402 DO k = 2, kde-1
2403 DO i = its,itf
2404 cq1 = 1./(1.+cqw(i,k,j))
2405 cq2 = cqw(i,k,j)*cq1
2406 cqw(i,k,j) = cq1
2407 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2408 cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
2409 -mu(i,j)-cq2*mub(i,j) )
2410 END DO
2411 ENDDO
2414 ENDDO
2416 END SUBROUTINE pg_buoy_w
2418 !-------------------------------------------------------------------------------
2420 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2421 u, v, ww, w, mut, rdnw, &
2422 rdx, rdy, msfux, msfuy, &
2423 msfvx, msfvy, dt, &
2424 config_flags, &
2425 ids, ide, jds, jde, kds, kde, &
2426 ims, ime, jms, jme, kms, kme, &
2427 its, ite, jts, jte, kts, kte )
2429 USE module_llxy
2430 IMPLICIT NONE
2432 ! Input data
2434 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2436 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2437 ims, ime, jms, jme, kms, kme, &
2438 its, ite, jts, jte, kts, kte
2440 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2442 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2444 REAL, INTENT(OUT) :: max_vert_cfl
2445 REAL, INTENT(OUT) :: max_horiz_cfl
2446 REAL :: horiz_cfl
2448 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2450 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2452 REAL, INTENT(IN) :: dt
2453 REAL, INTENT(IN) :: rdx, rdy
2454 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2455 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2457 REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2459 INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2460 INTEGER :: some
2461 CHARACTER*512 :: temp
2463 CHARACTER (LEN=256) :: time_str
2464 CHARACTER (LEN=256) :: grid_str
2466 integer :: total
2467 REAL :: msfuxt , msfxffl
2469 !<DESCRIPTION>
2470 !
2471 ! w_damp computes a damping term for the vertical velocity when the
2472 ! vertical Courant number is too large. This was found to be preferable to
2473 ! decreasing the timestep or increasing the diffusion in real-data applications
2474 ! that produced potentially-unstable large vertical velocities because of
2475 ! unphysically large heating rates coming from the cumulus parameterization
2476 ! schemes run at moderately high resolutions (dx ~ O(10) km).
2477 !
2478 ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2479 ! horizontal motion. These values are returned via the max_vert_cfl and
2480 ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2481 !
2482 !</DESCRIPTION>
2484 itf=MIN(ite,ide-1)
2485 jtf=MIN(jte,jde-1)
2487 some = 0
2488 max_vert_cfl = 0.
2489 max_horiz_cfl = 0.
2490 total = 0
2492 IF(config_flags%map_proj == PROJ_CASSINI ) then
2493 msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2494 END IF
2496 IF ( config_flags%w_damping == 1 ) THEN
2497 DO j = jts,jtf
2499 DO k = 2, kde-1
2500 DO i = its,itf
2501 #if 1
2502 IF(config_flags%map_proj == PROJ_CASSINI ) then
2503 msfuxt = MIN(msfux(i,j), msfxffl)
2504 ELSE
2505 msfuxt = msfux(i,j)
2506 END IF
2507 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2509 IF ( vert_cfl > max_vert_cfl ) THEN
2510 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2511 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2512 ENDIF
2514 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2515 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2516 if (horiz_cfl > max_horiz_cfl) then
2517 max_horiz_cfl = horiz_cfl
2518 endif
2520 if(vert_cfl .gt. w_beta)then
2521 #else
2522 ! restructure to get rid of divide
2523 !
2524 ! This had been used for efficiency, but with the addition of returning the cfl values,
2525 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2526 !
2527 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2528 cf_d = abs(mut(i,j))
2529 if(cf_n .gt. cf_d*w_beta )then
2530 #endif
2532 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2533 CALL wrf_debug ( 100 , TRIM(temp) )
2534 if ( vert_cfl > 2. ) some = some + 1
2535 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)
2536 endif
2537 END DO
2538 ENDDO
2539 ENDDO
2540 ELSE
2541 ! just print
2542 DO j = jts,jtf
2544 DO k = 2, kde-1
2545 DO i = its,itf
2547 #if 1
2548 IF(config_flags%map_proj == PROJ_CASSINI ) then
2549 msfuxt = MIN(msfux(i,j), msfxffl)
2550 ELSE
2551 msfuxt = msfux(i,j)
2552 END IF
2553 vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2555 IF ( vert_cfl > max_vert_cfl ) THEN
2556 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2557 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2558 ENDIF
2560 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2561 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2563 if (horiz_cfl > max_horiz_cfl) then
2564 max_horiz_cfl = horiz_cfl
2565 endif
2567 if(vert_cfl .gt. w_beta)then
2568 #else
2569 ! restructure to get rid of divide
2570 !
2571 ! This had been used for efficiency, but with the addition of returning the cfl values,
2572 ! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
2573 !
2574 cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2575 cf_d = abs(mut(i,j))
2576 if(cf_n .gt. cf_d*w_beta )then
2577 #endif
2578 WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2579 CALL wrf_debug ( 100 , TRIM(temp) )
2580 if ( vert_cfl > 2. ) some = some + 1
2581 endif
2582 END DO
2583 ENDDO
2584 ENDDO
2585 ENDIF
2586 IF ( some .GT. 0 ) THEN
2587 CALL get_current_time_string( time_str )
2588 CALL get_current_grid_name( grid_str )
2589 WRITE(wrf_err_message,*)some, &
2590 ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2591 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2592 WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2593 maxdub,maxdeta
2594 CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2595 ENDIF
2597 END SUBROUTINE w_damp
2599 !-------------------------------------------------------------------------------
2601 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, &
2602 config_flags, &
2603 msfux, msfuy, msfvx, msfvx_inv, &
2604 msfvy, msftx, msfty, &
2605 khdif, xkmhd, rdx, rdy, &
2606 ids, ide, jds, jde, kds, kde, &
2607 ims, ime, jms, jme, kms, kme, &
2608 its, ite, jts, jte, kts, kte )
2610 IMPLICIT NONE
2612 ! Input data
2614 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2616 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2617 ims, ime, jms, jme, kms, kme, &
2618 its, ite, jts, jte, kts, kte
2620 CHARACTER(LEN=1) , INTENT(IN ) :: name
2622 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2624 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2626 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2628 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2629 msfuy, &
2630 msfvx, &
2631 msfvx_inv, &
2632 msfvy, &
2633 msftx, &
2634 msfty
2636 REAL , INTENT(IN ) :: rdx, &
2637 rdy, &
2638 khdif
2640 ! Local data
2642 INTEGER :: i, j, k, itf, jtf, ktf
2644 INTEGER :: i_start, i_end, j_start, j_end
2646 REAL :: mrdx, mkrdxm, mkrdxp, &
2647 mrdy, mkrdym, mkrdyp
2649 LOGICAL :: specified
2651 !<DESCRIPTION>
2652 !
2653 ! horizontal_diffusion computes the horizontal diffusion tendency
2654 ! on model horizontal coordinate surfaces.
2655 !
2656 !</DESCRIPTION>
2658 specified = .false.
2659 if(config_flags%specified .or. config_flags%nested) specified = .true.
2661 ktf=MIN(kte,kde-1)
2663 IF (name .EQ. 'u') THEN
2665 i_start = its
2666 i_end = ite
2667 j_start = jts
2668 j_end = MIN(jte,jde-1)
2670 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2671 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2672 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2673 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2674 IF ( config_flags%periodic_x ) i_start = its
2675 IF ( config_flags%periodic_x ) i_end = ite
2678 DO j = j_start, j_end
2679 DO k=kts,ktf
2680 DO i = i_start, i_end
2682 ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2683 ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2685 mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2686 mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2687 mrdx=msfux(i,j)*msfuy(i,j)*rdx
2688 mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2689 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2690 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2691 mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2692 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2693 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2694 ! need to do four-corners (t) for diffusion coefficient as there are
2695 ! no values at u,v points
2696 ! msfuy - has to be y as part of d/dY
2697 ! has to be u as we're at a u point
2698 mrdy=msfux(i,j)*msfuy(i,j)*rdy
2700 ! correctly averaged version of rho~ * m^2 *
2701 ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2702 tendency(i,k,j)=tendency(i,k,j)+( &
2703 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2704 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2705 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2706 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2707 ENDDO
2708 ENDDO
2709 ENDDO
2711 ELSE IF (name .EQ. 'v')THEN
2713 i_start = its
2714 i_end = MIN(ite,ide-1)
2715 j_start = jts
2716 j_end = jte
2718 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2719 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2720 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2721 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2722 IF ( config_flags%periodic_x ) i_start = its
2723 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2724 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2725 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2727 DO j = j_start, j_end
2728 DO k=kts,ktf
2729 DO i = i_start, i_end
2731 mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2732 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2733 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2734 mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2735 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2736 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2737 mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2738 mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2739 mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2740 mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2742 tendency(i,k,j)=tendency(i,k,j)+( &
2743 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2744 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2745 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2746 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2747 ENDDO
2748 ENDDO
2749 ENDDO
2751 ELSE IF (name .EQ. 'w')THEN
2753 i_start = its
2754 i_end = MIN(ite,ide-1)
2755 j_start = jts
2756 j_end = MIN(jte,jde-1)
2758 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2759 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2760 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2761 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2762 IF ( config_flags%periodic_x ) i_start = its
2763 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2765 DO j = j_start, j_end
2766 DO k=kts+1,ktf
2767 DO i = i_start, i_end
2769 mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2770 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2771 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2772 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2773 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2774 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2775 mrdx=msftx(i,j)*msfty(i,j)*rdx
2776 ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2777 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
2778 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2779 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2780 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2781 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
2782 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2783 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2784 mrdy=msftx(i,j)*msfty(i,j)*rdy
2786 tendency(i,k,j)=tendency(i,k,j)+( &
2787 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2788 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2789 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2790 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2791 ENDDO
2792 ENDDO
2793 ENDDO
2795 ELSE
2798 i_start = its
2799 i_end = MIN(ite,ide-1)
2800 j_start = jts
2801 j_end = MIN(jte,jde-1)
2803 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2804 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2805 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2806 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2807 IF ( config_flags%periodic_x ) i_start = its
2808 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2810 DO j = j_start, j_end
2811 DO k=kts,ktf
2812 DO i = i_start, i_end
2814 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
2815 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
2816 mrdx=msftx(i,j)*msfty(i,j)*rdx
2817 ! 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
2818 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
2819 ! 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
2820 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
2821 mrdy=msftx(i,j)*msfty(i,j)*rdy
2823 tendency(i,k,j)=tendency(i,k,j)+( &
2824 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2825 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2826 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2827 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2828 ENDDO
2829 ENDDO
2830 ENDDO
2832 ENDIF
2834 END SUBROUTINE horizontal_diffusion
2836 !-----------------------------------------------------------------------------------------
2838 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, &
2839 config_flags, base_3d, &
2840 msfux, msfuy, msfvx, msfvx_inv, &
2841 msfvy, msftx, msfty, &
2842 khdif, xkmhd, rdx, rdy, &
2843 ids, ide, jds, jde, kds, kde, &
2844 ims, ime, jms, jme, kms, kme, &
2845 its, ite, jts, jte, kts, kte )
2847 IMPLICIT NONE
2849 ! Input data
2851 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2853 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2854 ims, ime, jms, jme, kms, kme, &
2855 its, ite, jts, jte, kts, kte
2857 CHARACTER(LEN=1) , INTENT(IN ) :: name
2859 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2860 xkmhd, &
2861 base_3d
2863 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2865 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
2867 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2868 msfuy, &
2869 msfvx, &
2870 msfvx_inv, &
2871 msfvy, &
2872 msftx, &
2873 msfty
2875 REAL , INTENT(IN ) :: rdx, &
2876 rdy, &
2877 khdif
2879 ! Local data
2881 INTEGER :: i, j, k, itf, jtf, ktf
2883 INTEGER :: i_start, i_end, j_start, j_end
2885 REAL :: mrdx, mkrdxm, mkrdxp, &
2886 mrdy, mkrdym, mkrdyp
2888 LOGICAL :: specified
2890 !<DESCRIPTION>
2891 !
2892 ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2893 ! on model horizontal coordinate surfaces. This routine computes diffusion
2894 ! a perturbation scalar (field-base_3d).
2895 !
2896 !</DESCRIPTION>
2898 specified = .false.
2899 if(config_flags%specified .or. config_flags%nested) specified = .true.
2901 ktf=MIN(kte,kde-1)
2903 i_start = its
2904 i_end = MIN(ite,ide-1)
2905 j_start = jts
2906 j_end = MIN(jte,jde-1)
2908 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2909 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2910 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2911 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2912 IF ( config_flags%periodic_x ) i_start = its
2913 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2915 DO j = j_start, j_end
2916 DO k=kts,ktf
2917 DO i = i_start, i_end
2919 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
2920 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
2921 mrdx=msftx(i,j)*msfty(i,j)*rdx
2922 ! 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
2923 ! 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
2924 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
2925 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
2926 mrdy=msftx(i,j)*msfty(i,j)*rdy
2928 tendency(i,k,j)=tendency(i,k,j)+( &
2929 mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
2930 -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
2931 -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
2932 -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
2933 +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
2934 -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
2935 -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
2936 -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
2937 )
2938 ENDDO
2939 ENDDO
2940 ENDDO
2942 END SUBROUTINE horizontal_diffusion_3dmp
2944 !-----------------------------------------------------------------------------------------
2946 SUBROUTINE vertical_diffusion ( name, field, tendency, &
2947 config_flags, &
2948 alt, mut, rdn, rdnw, kvdif, &
2949 ids, ide, jds, jde, kds, kde, &
2950 ims, ime, jms, jme, kms, kme, &
2951 its, ite, jts, jte, kts, kte )
2954 IMPLICIT NONE
2956 ! Input data
2958 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2960 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2961 ims, ime, jms, jme, kms, kme, &
2962 its, ite, jts, jte, kts, kte
2964 CHARACTER(LEN=1) , INTENT(IN ) :: name
2966 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2967 INTENT(IN ) :: field, &
2968 alt
2970 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2972 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2974 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
2976 REAL , INTENT(IN ) :: kvdif
2978 ! Local data
2980 INTEGER :: i, j, k, itf, jtf, ktf
2981 INTEGER :: i_start, i_end, j_start, j_end
2983 REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2984 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2986 REAL :: rdz
2988 LOGICAL :: specified
2990 !<DESCRIPTION>
2991 !
2992 ! vertical_diffusion
2993 ! computes vertical diffusion tendency.
2994 !
2995 !</DESCRIPTION>
2997 specified = .false.
2998 if(config_flags%specified .or. config_flags%nested) specified = .true.
3000 ktf=MIN(kte,kde-1)
3002 IF (name .EQ. 'w')THEN
3005 i_start = its
3006 i_end = MIN(ite,ide-1)
3007 j_start = jts
3008 j_end = MIN(jte,jde-1)
3010 j_loop_w : DO j = j_start, j_end
3012 DO k=kts,ktf-1
3013 DO i = i_start, i_end
3014 vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
3015 ENDDO
3016 ENDDO
3018 DO i = i_start, i_end
3019 vflux(i,ktf)=0.
3020 ENDDO
3022 DO k=kts+1,ktf
3023 DO i = i_start, i_end
3024 tendency(i,k,j)=tendency(i,k,j) &
3025 +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
3026 *(vflux(i,k)-vflux(i,k-1))
3027 ENDDO
3028 ENDDO
3030 ENDDO j_loop_w
3032 ELSE IF(name .EQ. 'm')THEN
3034 i_start = its
3035 i_end = MIN(ite,ide-1)
3036 j_start = jts
3037 j_end = MIN(jte,jde-1)
3039 j_loop_s : DO j = j_start, j_end
3041 DO k=kts,ktf-1
3042 DO i = i_start, i_end
3043 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3044 *(field(i,k+1,j)-field(i,k,j))
3045 ENDDO
3046 ENDDO
3048 DO i = i_start, i_end
3049 vflux(i,0)=vflux(i,1)
3050 ENDDO
3052 DO i = i_start, i_end
3053 vflux(i,ktf)=0.
3054 ENDDO
3056 DO k=kts,ktf
3057 DO i = i_start, i_end
3058 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3059 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3060 ENDDO
3061 ENDDO
3063 ENDDO j_loop_s
3065 ENDIF
3067 END SUBROUTINE vertical_diffusion
3070 !-------------------------------------------------------------------------------
3072 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3073 base, &
3074 alt, mut, rdn, rdnw, kvdif, &
3075 ids, ide, jds, jde, kds, kde, &
3076 ims, ime, jms, jme, kms, kme, &
3077 its, ite, jts, jte, kts, kte )
3080 IMPLICIT NONE
3082 ! Input data
3084 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3086 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3087 ims, ime, jms, jme, kms, kme, &
3088 its, ite, jts, jte, kts, kte
3090 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3091 INTENT(IN ) :: field, &
3092 alt
3094 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3096 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3098 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3099 rdnw, &
3100 base
3102 REAL , INTENT(IN ) :: kvdif
3104 ! Local data
3106 INTEGER :: i, j, k, itf, jtf, ktf
3107 INTEGER :: i_start, i_end, j_start, j_end
3109 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3111 REAL :: rdz
3113 LOGICAL :: specified
3115 !<DESCRIPTION>
3116 !
3117 ! vertical_diffusion_mp
3118 ! computes vertical diffusion tendency of a perturbation variable
3119 ! (field-base). Note that base as a 1D (k) field.
3120 !
3121 !</DESCRIPTION>
3123 specified = .false.
3124 if(config_flags%specified .or. config_flags%nested) specified = .true.
3126 ktf=MIN(kte,kde-1)
3128 i_start = its
3129 i_end = MIN(ite,ide-1)
3130 j_start = jts
3131 j_end = MIN(jte,jde-1)
3133 j_loop_s : DO j = j_start, j_end
3135 DO k=kts,ktf-1
3136 DO i = i_start, i_end
3137 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3138 *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3139 ENDDO
3140 ENDDO
3142 DO i = i_start, i_end
3143 vflux(i,0)=vflux(i,1)
3144 ENDDO
3146 DO i = i_start, i_end
3147 vflux(i,ktf)=0.
3148 ENDDO
3150 DO k=kts,ktf
3151 DO i = i_start, i_end
3152 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3153 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3154 ENDDO
3155 ENDDO
3157 ENDDO j_loop_s
3159 END SUBROUTINE vertical_diffusion_mp
3162 !-------------------------------------------------------------------------------
3164 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3165 base_3d, &
3166 alt, mut, rdn, rdnw, kvdif, &
3167 ids, ide, jds, jde, kds, kde, &
3168 ims, ime, jms, jme, kms, kme, &
3169 its, ite, jts, jte, kts, kte )
3172 IMPLICIT NONE
3174 ! Input data
3176 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3178 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3179 ims, ime, jms, jme, kms, kme, &
3180 its, ite, jts, jte, kts, kte
3182 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3183 INTENT(IN ) :: field, &
3184 alt, &
3185 base_3d
3187 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3189 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3191 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3192 rdnw
3194 REAL , INTENT(IN ) :: kvdif
3196 ! Local data
3198 INTEGER :: i, j, k, itf, jtf, ktf
3199 INTEGER :: i_start, i_end, j_start, j_end
3201 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3203 REAL :: rdz
3205 LOGICAL :: specified
3207 !<DESCRIPTION>
3208 !
3209 ! vertical_diffusion_3dmp
3210 ! computes vertical diffusion tendency of a perturbation variable
3211 ! (field-base_3d).
3212 !
3213 !</DESCRIPTION>
3215 specified = .false.
3216 if(config_flags%specified .or. config_flags%nested) specified = .true.
3218 ktf=MIN(kte,kde-1)
3220 i_start = its
3221 i_end = MIN(ite,ide-1)
3222 j_start = jts
3223 j_end = MIN(jte,jde-1)
3225 j_loop_s : DO j = j_start, j_end
3227 DO k=kts,ktf-1
3228 DO i = i_start, i_end
3229 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3230 *( field(i,k+1,j) -field(i,k,j) &
3231 -base_3d(i,k+1,j)+base_3d(i,k,j) )
3232 ENDDO
3233 ENDDO
3235 DO i = i_start, i_end
3236 vflux(i,0)=vflux(i,1)
3237 ENDDO
3239 DO i = i_start, i_end
3240 vflux(i,ktf)=0.
3241 ENDDO
3243 DO k=kts,ktf
3244 DO i = i_start, i_end
3245 tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) &
3246 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3247 ENDDO
3248 ENDDO
3250 ENDDO j_loop_s
3252 END SUBROUTINE vertical_diffusion_3dmp
3255 !-------------------------------------------------------------------------------
3258 SUBROUTINE vertical_diffusion_u ( field, tendency, &
3259 config_flags, u_base, &
3260 alt, muu, rdn, rdnw, kvdif, &
3261 ids, ide, jds, jde, kds, kde, &
3262 ims, ime, jms, jme, kms, kme, &
3263 its, ite, jts, jte, kts, kte )
3266 IMPLICIT NONE
3268 ! Input data
3270 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3272 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3273 ims, ime, jms, jme, kms, kme, &
3274 its, ite, jts, jte, kts, kte
3276 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3277 INTENT(IN ) :: field, &
3278 alt
3280 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3282 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3284 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3286 REAL , INTENT(IN ) :: kvdif
3288 ! Local data
3290 INTEGER :: i, j, k, itf, jtf, ktf
3291 INTEGER :: i_start, i_end, j_start, j_end
3293 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3295 REAL :: rdz, zz
3297 LOGICAL :: specified
3299 !<DESCRIPTION>
3300 !
3301 ! vertical_diffusion_u computes vertical diffusion tendency for
3302 ! the u momentum equation. This routine assumes a constant eddy
3303 ! viscosity kvdif.
3304 !
3305 !</DESCRIPTION>
3307 specified = .false.
3308 if(config_flags%specified .or. config_flags%nested) specified = .true.
3310 ktf=MIN(kte,kde-1)
3312 i_start = its
3313 i_end = ite
3314 j_start = jts
3315 j_end = MIN(jte,jde-1)
3317 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3318 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3319 IF ( config_flags%periodic_x ) i_start = its
3320 IF ( config_flags%periodic_x ) i_end = ite
3323 j_loop_u : DO j = j_start, j_end
3325 DO k=kts,ktf-1
3326 DO i = i_start, i_end
3327 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
3328 +alt(i-1,k ,j) &
3329 +alt(i ,k+1,j) &
3330 +alt(i-1,k+1,j) ) ) &
3331 *(field(i,k+1,j)-field(i,k,j) &
3332 -u_base(k+1) +u_base(k) )
3333 ENDDO
3334 ENDDO
3336 DO i = i_start, i_end
3337 vflux(i,0)=vflux(i,1)
3338 ENDDO
3340 DO i = i_start, i_end
3341 vflux(i,ktf)=0.
3342 ENDDO
3344 DO k=kts,ktf-1
3345 DO i = i_start, i_end
3346 tendency(i,k,j)=tendency(i,k,j)+ &
3347 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3348 (vflux(i,k)-vflux(i,k-1))
3349 ENDDO
3350 ENDDO
3352 ENDDO j_loop_u
3354 END SUBROUTINE vertical_diffusion_u
3356 !-------------------------------------------------------------------------------
3359 SUBROUTINE vertical_diffusion_v ( field, tendency, &
3360 config_flags, v_base, &
3361 alt, muv, rdn, rdnw, kvdif, &
3362 ids, ide, jds, jde, kds, kde, &
3363 ims, ime, jms, jme, kms, kme, &
3364 its, ite, jts, jte, kts, kte )
3367 IMPLICIT NONE
3369 ! Input data
3371 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3373 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3374 ims, ime, jms, jme, kms, kme, &
3375 its, ite, jts, jte, kts, kte
3377 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3378 INTENT(IN ) :: field, &
3379 alt
3380 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3382 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3384 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3386 REAL , INTENT(IN ) :: kvdif
3388 ! Local data
3390 INTEGER :: i, j, k, itf, jtf, ktf, jm1
3391 INTEGER :: i_start, i_end, j_start, j_end
3393 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3395 REAL :: rdz, zz
3397 LOGICAL :: specified
3399 !<DESCRIPTION>
3400 !
3401 ! vertical_diffusion_v computes vertical diffusion tendency for
3402 ! the v momentum equation. This routine assumes a constant eddy
3403 ! viscosity kvdif.
3404 !
3405 !</DESCRIPTION>
3407 specified = .false.
3408 if(config_flags%specified .or. config_flags%nested) specified = .true.
3410 ktf=MIN(kte,kde-1)
3412 i_start = its
3413 i_end = MIN(ite,ide-1)
3414 j_start = jts
3415 j_end = MIN(jte,jde-1)
3417 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3418 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3420 j_loop_v : DO j = j_start, j_end
3421 ! jm1 = max(j-1,1)
3422 jm1 = j-1
3424 DO k=kts,ktf-1
3425 DO i = i_start, i_end
3426 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
3427 +alt(i,k ,jm1) &
3428 +alt(i,k+1,j ) &
3429 +alt(i,k+1,jm1) ) ) &
3430 *(field(i,k+1,j)-field(i,k,j) &
3431 -v_base(k+1) +v_base(k) )
3432 ENDDO
3433 ENDDO
3435 DO i = i_start, i_end
3436 vflux(i,0)=vflux(i,1)
3437 ENDDO
3439 DO i = i_start, i_end
3440 vflux(i,ktf)=0.
3441 ENDDO
3443 DO k=kts,ktf-1
3444 DO i = i_start, i_end
3445 tendency(i,k,j)=tendency(i,k,j)+ &
3446 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
3447 (vflux(i,k)-vflux(i,k-1))
3448 ENDDO
3449 ENDDO
3451 ENDDO j_loop_v
3453 END SUBROUTINE vertical_diffusion_v
3455 !*************** end new mass coordinate routines
3457 !-------------------------------------------------------------------------------
3459 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3460 ids, ide, jds, jde, kds, kde, &
3461 ims, ime, jms, jme, kms, kme, &
3462 its, ite, jts, jte, kts, kte )
3464 IMPLICIT NONE
3466 ! Input data
3468 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3469 ims, ime, jms, jme, kms, kme, &
3470 its, ite, jts, jte, kts, kte
3472 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3473 rfieldp
3475 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3477 ! Local indices.
3479 INTEGER :: i, j, k, itf, jtf, ktf
3481 !<DESCRIPTION>
3482 !
3483 ! calculate_full
3484 ! calculates full 3D field from pertubation and base field.
3485 !
3486 !</DESCRIPTION>
3488 itf=MIN(ite,ide-1)
3489 jtf=MIN(jte,jde-1)
3490 ktf=MIN(kte,kde-1)
3492 DO j=jts,jtf
3493 DO k=kts,ktf
3494 DO i=its,itf
3495 rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3496 ENDDO
3497 ENDDO
3498 ENDDO
3500 END SUBROUTINE calculate_full
3502 !------------------------------------------------------------------------------
3504 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3505 config_flags, &
3506 msftx, msfty, msfux, msfuy, &
3507 msfvx, msfvy, &
3508 f, e, sina, cosa, fzm, fzp, &
3509 ids, ide, jds, jde, kds, kde, &
3510 ims, ime, jms, jme, kms, kme, &
3511 its, ite, jts, jte, kts, kte )
3513 IMPLICIT NONE
3515 ! Input data
3517 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3519 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3520 ims, ime, jms, jme, kms, kme, &
3521 its, ite, jts, jte, kts, kte
3523 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3524 rv_tend, &
3525 rw_tend
3526 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3527 rv, &
3528 rw
3530 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3531 msfuy, &
3532 msfvx, &
3533 msfvy, &
3534 msftx, &
3535 msfty
3537 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3538 e, &
3539 sina, &
3540 cosa
3542 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3543 fzp
3545 ! Local indices.
3547 INTEGER :: i, j , k, ktf
3548 INTEGER :: i_start, i_end, j_start, j_end
3550 LOGICAL :: specified
3552 !<DESCRIPTION>
3553 !
3554 ! coriolis calculates the large timestep tendency terms in the
3555 ! u, v, and w momentum equations arise from the coriolis force.
3556 !
3557 !</DESCRIPTION>
3559 specified = .false.
3560 if(config_flags%specified .or. config_flags%nested) specified = .true.
3562 ktf=MIN(kte,kde-1)
3564 ! coriolis for u-momentum equation
3566 ! Notes on map scale factor
3567 ! cosa, sina are related to rotating the coordinate frame if desired
3568 ! generally sina=0, cosa=1
3569 ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3570 ! + 2 mu v omega sin(lat)/my
3571 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3572 ! => terms are: -e mu w / my + f mu v / my
3573 ! rv = mu v / mx ; rw = mu w / my
3574 ! => terms are: -e rw + f rv *mx / my
3576 i_start = its
3577 i_end = ite
3578 IF ( config_flags%open_xs .or. specified .or. &
3579 config_flags%nested) i_start = MAX(ids+1,its)
3580 IF ( config_flags%open_xe .or. specified .or. &
3581 config_flags%nested) i_end = MIN(ide-1,ite)
3582 IF ( config_flags%periodic_x ) i_start = its
3583 IF ( config_flags%periodic_x ) i_end = ite
3585 DO j = jts, MIN(jte,jde-1)
3587 DO k=kts,ktf
3588 DO i = i_start, i_end
3590 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)) &
3591 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3592 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3593 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3595 ENDDO
3596 ENDDO
3598 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3600 DO k=kts,ktf
3602 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3603 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3604 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3605 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3607 ENDDO
3609 ENDIF
3611 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3613 DO k=kts,ktf
3615 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)) &
3616 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3617 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3618 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3620 ENDDO
3622 ENDIF
3624 ENDDO
3626 ! coriolis term for v-momentum equation
3628 ! Notes on map scale factors
3629 ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3630 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3631 ! => terms are: -f mu u / mx
3632 ! ru = mu u / my ; rw = mu w / my
3633 ! => terms are: -f ru *my / mx + ?
3635 j_start = jts
3636 j_end = jte
3638 IF ( config_flags%open_ys .or. specified .or. &
3639 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3640 IF ( config_flags%open_ye .or. specified .or. &
3641 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3643 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3645 DO k=kts,ktf
3646 DO i=its,MIN(ide-1,ite)
3648 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3649 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3650 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3651 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3653 ENDDO
3654 ENDDO
3656 ENDIF
3658 DO j=j_start, j_end
3659 DO k=kts,ktf
3660 DO i=its,MIN(ide-1,ite)
3662 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)) &
3663 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3664 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3665 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3667 ENDDO
3668 ENDDO
3669 ENDDO
3672 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3674 DO k=kts,ktf
3675 DO i=its,MIN(ide-1,ite)
3677 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)) &
3678 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3679 + (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)) &
3680 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3682 ENDDO
3683 ENDDO
3685 ENDIF
3687 ! coriolis term for w-mometum
3689 ! Notes on map scale factors
3690 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3691 ! Define e=2 omega cos(lat)
3692 ! => terms are: e mu u / my + ???
3693 ! ru = mu u / my ; ru = mu v / mx
3694 ! => terms are: e ru + ???
3696 DO j=jts,MIN(jte, jde-1)
3697 DO k=kts+1,ktf
3698 DO i=its,MIN(ite, ide-1)
3700 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3701 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3702 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3703 -(msftx(i,j)/msfty(i,j))* &
3704 sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3705 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3707 ENDDO
3708 ENDDO
3709 ENDDO
3711 END SUBROUTINE coriolis
3713 !------------------------------------------------------------------------------
3715 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3716 config_flags, &
3717 u_base, v_base, z_base, &
3718 muu, muv, phb, ph, &
3719 msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3720 f, e, sina, cosa, fzm, fzp, &
3721 ids, ide, jds, jde, kds, kde, &
3722 ims, ime, jms, jme, kms, kme, &
3723 its, ite, jts, jte, kts, kte )
3725 IMPLICIT NONE
3727 ! Input data
3729 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3731 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3732 ims, ime, jms, jme, kms, kme, &
3733 its, ite, jts, jte, kts, kte
3735 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3736 rv_tend, &
3737 rw_tend
3738 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3739 rv_in, &
3740 rw, &
3741 ph, &
3742 phb
3745 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3746 msfuy, &
3747 msfvx, &
3748 msfvy, &
3749 msftx, &
3750 msfty
3752 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3753 e, &
3754 sina, &
3755 cosa
3757 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3758 muv
3761 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3762 fzp
3764 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3765 v_base, &
3766 z_base
3768 ! Local storage
3770 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3771 rv
3773 REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3775 ! Local indices.
3777 INTEGER :: i, j , k, ktf
3778 INTEGER :: i_start, i_end, j_start, j_end
3780 LOGICAL :: specified
3782 !<DESCRIPTION>
3783 !
3784 ! perturbation_coriolis calculates the large timestep tendency terms in the
3785 ! u, v, and w momentum equations arise from the coriolis force. This version
3786 ! subtracts off the horizontal velocities from the initial sounding when
3787 ! computing the forcing terms, hence "perturbation" coriolis.
3788 !
3789 !</DESCRIPTION>
3791 specified = .false.
3792 if(config_flags%specified .or. config_flags%nested) specified = .true.
3794 ktf=MIN(kte,kde-1)
3796 ! coriolis for u-momentum equation
3798 i_start = its
3799 i_end = ite
3800 IF ( config_flags%open_xs .or. specified .or. &
3801 config_flags%nested) i_start = MAX(ids+1,its)
3802 IF ( config_flags%open_xe .or. specified .or. &
3803 config_flags%nested) i_end = MIN(ide-1,ite)
3804 IF ( config_flags%periodic_x ) i_start = its
3805 IF ( config_flags%periodic_x ) i_end = ite
3807 ! compute perturbation mu*v for use in u momentum equation
3809 DO j = jts, MIN(jte,jde-1)+1
3810 DO k=kts+1,ktf-1
3811 DO i = i_start-1, i_end
3812 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3813 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3814 +ph(i,k,j )+ph(i,k+1,j ) &
3815 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3816 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3817 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3818 wk = 1.-wkp1-wkm1
3819 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3820 wkm1*v_base(k-1) &
3821 +wk *v_base(k ) &
3822 +wkp1*v_base(k+1) )
3823 ENDDO
3824 ENDDO
3825 ENDDO
3828 ! pick up top and bottom v
3830 DO j = jts, MIN(jte,jde-1)+1
3831 DO i = i_start-1, i_end
3833 k = kts
3834 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3835 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3836 +ph(i,k,j )+ph(i,k+1,j ) &
3837 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3838 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3839 wk = 1.-wkp1
3840 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3841 +wk *v_base(k ) &
3842 +wkp1*v_base(k+1) )
3844 k = ktf
3845 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3846 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3847 +ph(i,k,j )+ph(i,k+1,j ) &
3848 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3849 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3850 wk = 1.-wkm1
3851 rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*( &
3852 wkm1*v_base(k-1) &
3853 +wk *v_base(k ) )
3855 ENDDO
3856 ENDDO
3858 ! compute coriolis forcing for u
3860 ! Map scale factors: see comments above for Coriolis
3862 DO j = jts, MIN(jte,jde-1)
3864 DO k=kts,ktf
3865 DO i = i_start, i_end
3866 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)) &
3867 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3868 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3869 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3870 ENDDO
3871 ENDDO
3873 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3875 DO k=kts,ktf
3877 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3878 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3879 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3880 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3882 ENDDO
3884 ENDIF
3886 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3888 DO k=kts,ktf
3890 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)) &
3891 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3892 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3893 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3895 ENDDO
3897 ENDIF
3899 ENDDO
3901 ! coriolis term for v-momentum equation
3902 ! Map scale factors: see comments above for Coriolis
3904 j_start = jts
3905 j_end = jte
3907 IF ( config_flags%open_ys .or. specified .or. &
3908 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3909 IF ( config_flags%open_ye .or. specified .or. &
3910 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3912 ! compute perturbation mu*u for use in v momentum equation
3914 DO j = j_start-1,j_end
3915 DO k=kts+1,ktf-1
3916 DO i = its, MIN(ite,ide-1)+1
3917 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3918 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3919 +ph(i ,k,j)+ph(i ,k+1,j) &
3920 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3921 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3922 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3923 wk = 1.-wkp1-wkm1
3924 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3925 wkm1*u_base(k-1) &
3926 +wk *u_base(k ) &
3927 +wkp1*u_base(k+1) )
3928 ENDDO
3929 ENDDO
3930 ENDDO
3932 ! pick up top and bottom u
3934 DO j = j_start-1,j_end
3935 DO i = its, MIN(ite,ide-1)+1
3937 k = kts
3938 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3939 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3940 +ph(i ,k,j)+ph(i ,k+1,j) &
3941 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3942 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3943 wk = 1.-wkp1
3944 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3945 +wk *u_base(k ) &
3946 +wkp1*u_base(k+1) )
3949 k = ktf
3950 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
3951 +phb(i-1,k,j)+phb(i-1,k+1,j) &
3952 +ph(i ,k,j)+ph(i ,k+1,j) &
3953 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3954 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3955 wk = 1.-wkm1
3956 ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*( &
3957 wkm1*u_base(k-1) &
3958 +wk *u_base(k ) )
3960 ENDDO
3961 ENDDO
3963 ! compute coriolis forcing for v momentum equation
3964 ! Map scale factors: see comments above for Coriolis
3966 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3968 DO k=kts,ktf
3969 DO i=its,MIN(ide-1,ite)
3971 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3972 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3973 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3974 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3976 ENDDO
3977 ENDDO
3979 ENDIF
3981 DO j=j_start, j_end
3982 DO k=kts,ktf
3983 DO i=its,MIN(ide-1,ite)
3985 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)) &
3986 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3987 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3988 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3990 ENDDO
3991 ENDDO
3992 ENDDO
3995 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3997 DO k=kts,ktf
3998 DO i=its,MIN(ide-1,ite)
4000 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)) &
4001 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
4002 + (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)) &
4003 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
4005 ENDDO
4006 ENDDO
4008 ENDIF
4010 ! coriolis term for w-mometum
4011 ! Map scale factors: see comments above for Coriolis
4013 DO j=jts,MIN(jte, jde-1)
4014 DO k=kts+1,ktf
4015 DO i=its,MIN(ite, ide-1)
4017 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
4018 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4019 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4020 -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
4021 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4023 ENDDO
4024 ENDDO
4025 ENDDO
4027 END SUBROUTINE perturbation_coriolis
4029 !------------------------------------------------------------------------------
4031 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4032 config_flags, &
4033 msfux, msfuy, msfvx, msfvy, msftx, msfty, &
4034 xlat, fzm, fzp, rdx, rdy, &
4035 ids, ide, jds, jde, kds, kde, &
4036 ims, ime, jms, jme, kms, kme, &
4037 its, ite, jts, jte, kts, kte )
4040 IMPLICIT NONE
4042 ! Input data
4044 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4046 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4047 ims, ime, jms, jme, kms, kme, &
4048 its, ite, jts, jte, kts, kte
4050 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4051 INTENT(INOUT) :: ru_tend, &
4052 rv_tend, &
4053 rw_tend
4055 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4056 INTENT(IN ) :: ru, &
4057 rv, &
4058 rw, &
4059 u, &
4060 v, &
4061 w
4063 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4064 msfuy, &
4065 msfvx, &
4066 msfvy, &
4067 msftx, &
4068 msfty, &
4069 xlat
4071 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4072 fzp
4074 REAL , INTENT(IN ) :: rdx, &
4075 rdy
4077 ! Local data
4079 ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4080 INTEGER :: i, j, k, itf, jtf, ktf
4081 INTEGER :: i_start, i_end, j_start, j_end
4082 ! INTEGER :: irmin, irmax, jrmin, jrmax
4084 REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4086 LOGICAL :: specified
4088 !<DESCRIPTION>
4089 !
4090 ! curvature calculates the large timestep tendency terms in the
4091 ! u, v, and w momentum equations arise from the curvature terms.
4092 !
4093 !</DESCRIPTION>
4095 specified = .false.
4096 if(config_flags%specified .or. config_flags%nested) specified = .true.
4098 itf=MIN(ite,ide-1)
4099 jtf=MIN(jte,jde-1)
4100 ktf=MIN(kte,kde-1)
4102 ! irmin = ims
4103 ! irmax = ime
4104 ! jrmin = jms
4105 ! jrmax = jme
4106 ! IF ( config_flags%open_xs ) irmin = ids
4107 ! IF ( config_flags%open_xe ) irmax = ide-1
4108 ! IF ( config_flags%open_ys ) jrmin = jds
4109 ! IF ( config_flags%open_ye ) jrmax = jde-1
4111 ! Define v cross grad m at scalar points - vxgm(i,j)
4113 i_start = its-1
4114 i_end = ite
4115 j_start = jts-1
4116 j_end = jte
4118 IF ( ( config_flags%open_xs .or. specified .or. &
4119 config_flags%nested) .and. (its == ids) ) i_start = its
4120 IF ( ( config_flags%open_xe .or. specified .or. &
4121 config_flags%nested) .and. (ite == ide) ) i_end = ite-1
4122 IF ( ( config_flags%open_ys .or. specified .or. &
4123 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4124 IF ( ( config_flags%open_ye .or. specified .or. &
4125 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
4126 IF ( config_flags%periodic_x ) i_start = its-1
4127 IF ( config_flags%periodic_x ) i_end = ite
4129 DO j=j_start, j_end
4130 DO k=kts,ktf
4131 DO i=i_start, i_end
4132 ! Map scale factor notes:
4133 ! msf...y is constant everywhere for cylindrical map projection
4134 ! msf...x varies with y only
4135 ! But we know that this is not = 0 for cylindrical,
4136 ! therefore use msfvX in 1st line
4137 ! which => by symmetry use msfuY in 2nd line - ???
4138 vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4139 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
4140 ENDDO
4141 ENDDO
4142 ENDDO
4144 ! Pick up the boundary rows for open (radiation) lateral b.c.
4145 ! Rather crude at present, we are assuming there is no
4146 ! variation in this term at the boundary.
4148 IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4149 config_flags%nested) .and. (its == ids) ) THEN
4151 DO j = jts, jte-1
4152 DO k = kts, ktf
4153 vxgm(its-1,k,j) = vxgm(its,k,j)
4154 ENDDO
4155 ENDDO
4157 ENDIF
4159 IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4160 config_flags%nested) .and. (ite == ide) ) THEN
4162 DO j = jts, jte-1
4163 DO k = kts, ktf
4164 vxgm(ite,k,j) = vxgm(ite-1,k,j)
4165 ENDDO
4166 ENDDO
4168 ENDIF
4170 ! Polar boundary condition:
4171 ! The following change is needed in case one tries using the vxgm route with
4172 ! polar B.C.'s in the future, but not needed if 'tan' used
4173 IF ( ( config_flags%open_ys .or. specified .or. &
4174 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4176 DO k = kts, ktf
4177 DO i = its-1, ite
4178 vxgm(i,k,jts-1) = vxgm(i,k,jts)
4179 ENDDO
4180 ENDDO
4182 ENDIF
4184 ! Polar boundary condition:
4185 ! The following change is needed in case one tries using the vxgm route with
4186 ! polar B.C.'s in the future, but not needed if 'tan' used
4187 IF ( ( config_flags%open_ye .or. specified .or. &
4188 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4190 DO k = kts, ktf
4191 DO i = its-1, ite
4192 vxgm(i,k,jte) = vxgm(i,k,jte-1)
4193 ENDDO
4194 ENDDO
4196 ENDIF
4198 ! curvature term for u momentum eqn.
4200 ! Map scale factor notes:
4201 ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4202 ! - mu u w /(a my)
4203 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4204 ! => terms are:
4205 ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4206 ! ru v tan(lat) / a - u rw / a
4207 ! xlat defined with end points half grid space from pole,
4208 ! hence are on u latitude points
4210 i_start = its
4211 IF ( config_flags%open_xs .or. specified .or. &
4212 config_flags%nested) i_start = MAX ( ids+1 , its )
4213 IF ( config_flags%open_xe .or. specified .or. &
4214 config_flags%nested) i_end = MIN ( ide-1 , ite )
4215 IF ( config_flags%periodic_x ) i_start = its
4216 IF ( config_flags%periodic_x ) i_end = ite
4218 ! Polar boundary condition
4219 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4221 DO j=jts,MIN(jde-1,jte)
4222 DO k=kts,ktf
4223 DO i=i_start,i_end
4225 ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
4226 (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4227 rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4228 - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4229 ENDDO
4230 ENDDO
4231 ENDDO
4233 ELSE ! normal code
4236 DO j=jts,MIN(jde-1,jte)
4237 DO k=kts,ktf
4238 DO i=i_start,i_end
4240 ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4241 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4242 - u(i,k,j)*reradius &
4243 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4245 ENDDO
4246 ENDDO
4247 ENDDO
4249 END IF
4251 ! curvature term for v momentum eqn.
4253 ! Map scale factor notes
4254 ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
4255 ! - mu v w /(a mx)
4256 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4257 ! terms are:
4258 ! - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
4259 ! = - [my/(mx*a)]*[u ru tan(lat) + v rw]
4260 ! - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
4261 ! xlat defined with end points half grid space from pole, hence are on
4262 ! u latitude points => av here
4263 !
4264 ! in original wrf, there was a sign error for the rw contribution
4266 j_start = jts
4267 IF ( config_flags%open_ys .or. specified .or. &
4268 config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4269 IF ( config_flags%open_ye .or. specified .or. &
4270 config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4272 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4274 DO j=j_start,j_end
4275 DO k=kts,ktf
4276 DO i=its,MIN(ite,ide-1)
4277 rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
4278 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
4279 tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
4280 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4281 + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4282 rw(i,k+1,j)+rw(i,k,j)) )
4283 ENDDO
4284 ENDDO
4285 ENDDO
4287 ELSE ! normal code
4289 DO j=j_start,j_end
4290 DO k=kts,ktf
4291 DO i=its,MIN(ite,ide-1)
4293 rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4294 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4295 - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
4296 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4298 ENDDO
4299 ENDDO
4300 ENDDO
4302 END IF
4304 ! curvature term for vertical momentum eqn.
4306 ! Notes on map scale factors:
4307 ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4308 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4309 ! terms are: u ru / a + (mx/my)v rv / a
4311 DO j=jts,MIN(jte,jde-1)
4312 DO k=MAX(2,kts),ktf
4313 DO i=its,MIN(ite,ide-1)
4315 rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4316 (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))) &
4317 *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))) &
4318 +(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))) &
4319 *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))))
4321 ENDDO
4322 ENDDO
4323 ENDDO
4325 END SUBROUTINE curvature
4327 !------------------------------------------------------------------------------
4329 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4330 fzm, fzp, &
4331 ids, ide, jds, jde, kds, kde, &
4332 ims, ime, jms, jme, kms, kme, &
4333 its, ite, jts, jte, kts, kte )
4335 IMPLICIT NONE
4337 ! Input data
4339 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4341 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4342 ims, ime, jms, jme, kms, kme, &
4343 its, ite, jts, jte, kts, kte
4345 CHARACTER(LEN=1) , INTENT(IN ) :: name
4347 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4349 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4351 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4353 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4355 ! Local data
4357 INTEGER :: i, j, k, itf, jtf, ktf
4359 !<DESCRIPTION>
4360 !
4361 ! decouple decouples a variable from the column dry-air mass.
4362 !
4363 !</DESCRIPTION>
4365 ktf=MIN(kte,kde-1)
4367 IF (name .EQ. 'u')THEN
4368 itf=ite
4369 jtf=MIN(jte,jde-1)
4371 DO j=jts,jtf
4372 DO k=kts,ktf
4373 DO i=its,itf
4374 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4375 ENDDO
4376 ENDDO
4377 ENDDO
4379 ELSE IF (name .EQ. 'v')THEN
4380 itf=MIN(ite,ide-1)
4381 jtf=jte
4383 DO j=jts,jtf
4384 DO k=kts,ktf
4385 DO i=its,itf
4386 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4387 ENDDO
4388 ENDDO
4389 ENDDO
4391 ELSE IF (name .EQ. 'w')THEN
4392 itf=MIN(ite,ide-1)
4393 jtf=MIN(jte,jde-1)
4394 DO j=jts,jtf
4395 DO k=kts+1,ktf
4396 DO i=its,itf
4397 field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4398 ENDDO
4399 ENDDO
4400 ENDDO
4402 DO j=jts,jtf
4403 DO i=its,itf
4404 field(i,kte,j) = 0.
4405 ENDDO
4406 ENDDO
4408 ELSE
4409 itf=MIN(ite,ide-1)
4410 jtf=MIN(jte,jde-1)
4411 ! For theta we will decouple tb and tp and add them to give t afterwards
4412 DO j=jts,jtf
4413 DO k=kts,ktf
4414 DO i=its,itf
4415 field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4416 ENDDO
4417 ENDDO
4418 ENDDO
4420 ENDIF
4422 END SUBROUTINE decouple
4424 !-------------------------------------------------------------------------------
4427 SUBROUTINE zero_tend ( tendency, &
4428 ids, ide, jds, jde, kds, kde, &
4429 ims, ime, jms, jme, kms, kme, &
4430 its, ite, jts, jte, kts, kte )
4433 IMPLICIT NONE
4435 ! Input data
4437 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4438 ims, ime, jms, jme, kms, kme, &
4439 its, ite, jts, jte, kts, kte
4441 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4443 ! Local data
4445 INTEGER :: i, j, k, itf, jtf, ktf
4447 !<DESCRIPTION>
4448 !
4449 ! zero_tend sets the input tendency array to zero.
4450 !
4451 !</DESCRIPTION>
4453 DO j = jts, jte
4454 DO k = kts, kte
4455 DO i = its, ite
4456 tendency(i,k,j) = 0.
4457 ENDDO
4458 ENDDO
4459 ENDDO
4461 END SUBROUTINE zero_tend
4463 !-------------------------------------------------------------------------------
4464 ! Sets the an array on the polar v point(s) to zero
4465 SUBROUTINE zero_pole ( field, &
4466 ids, ide, jds, jde, kds, kde, &
4467 ims, ime, jms, jme, kms, kme, &
4468 its, ite, jts, jte, kts, kte )
4471 IMPLICIT NONE
4473 ! Input data
4475 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4476 ims, ime, jms, jme, kms, kme, &
4477 its, ite, jts, jte, kts, kte
4479 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4481 ! Local data
4483 INTEGER :: i, k
4485 IF (jts == jds) THEN
4486 DO k = kts, kte
4487 DO i = its-1, ite+1
4488 field(i,k,jts) = 0.
4489 END DO
4490 END DO
4491 END IF
4492 IF (jte == jde) THEN
4493 DO k = kts, kte
4494 DO i = its-1, ite+1
4495 field(i,k,jte) = 0.
4496 END DO
4497 END DO
4498 END IF
4500 END SUBROUTINE zero_pole
4502 !-------------------------------------------------------------------------------
4503 ! Sets the an array on the polar v point(s)
4504 SUBROUTINE pole_point_bc ( field, &
4505 ids, ide, jds, jde, kds, kde, &
4506 ims, ime, jms, jme, kms, kme, &
4507 its, ite, jts, jte, kts, kte )
4510 IMPLICIT NONE
4512 ! Input data
4514 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4515 ims, ime, jms, jme, kms, kme, &
4516 its, ite, jts, jte, kts, kte
4518 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4520 ! Local data
4522 INTEGER :: i, k
4524 IF (jts == jds) THEN
4525 DO k = kts, kte
4526 DO i = its, ite
4527 ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4528 field(i,k,jts) = field(i,k,jts+1)
4529 END DO
4530 END DO
4531 END IF
4532 IF (jte == jde) THEN
4533 DO k = kts, kte
4534 DO i = its, ite
4535 ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4536 field(i,k,jte) = field(i,k,jte-1)
4537 END DO
4538 END DO
4539 END IF
4541 END SUBROUTINE pole_point_bc
4543 !======================================================================
4544 ! physics prep routines
4545 !======================================================================
4547 SUBROUTINE phy_prep ( config_flags, & ! input
4548 mu, muu, muv, u, v, p, pb, alt, ph, & ! input
4549 phb, t, tsk, moist, n_moist, & ! input
4550 rho, th_phy, p_phy , pi_phy , & ! output
4551 u_phy, v_phy, p8w, t_phy, t8w, & ! output
4552 z, z_at_w, dz8w, & ! output
4553 p_hyd, p_hyd_w, & ! output
4554 fzm, fzp, znw, p_top, & ! params
4555 RTHRATEN, &
4556 RTHBLTEN, RUBLTEN, RVBLTEN, &
4557 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4558 RTHCUTEN, RQVCUTEN, RQCCUTEN, &
4559 RQRCUTEN, RQICUTEN, RQSCUTEN, &
4560 RTHFTEN, RQVFTEN, &
4561 RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4562 RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
4563 ids, ide, jds, jde, kds, kde, &
4564 ims, ime, jms, jme, kms, kme, &
4565 its, ite, jts, jte, kts, kte )
4566 !----------------------------------------------------------------------
4567 IMPLICIT NONE
4568 !----------------------------------------------------------------------
4570 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4572 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4573 ims, ime, jms, jme, kms, kme, &
4574 its, ite, jts, jte, kts, kte
4575 INTEGER , INTENT(IN ) :: n_moist
4577 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4580 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu, muu, muv
4582 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4583 INTENT( OUT) :: u_phy, &
4584 v_phy, &
4585 pi_phy, &
4586 p_phy, &
4587 p8w, &
4588 t_phy, &
4589 th_phy, &
4590 t8w, &
4591 rho, &
4592 z, &
4593 dz8w, &
4594 z_at_w
4596 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4597 INTENT( OUT) :: p_hyd, &
4598 p_hyd_w
4600 REAL , INTENT(IN ) :: p_top
4602 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4603 INTENT(IN ) :: pb, &
4604 p, &
4605 u, &
4606 v, &
4607 alt, &
4608 ph, &
4609 phb, &
4610 t
4613 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4614 fzp
4616 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw
4618 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4619 INTENT(INOUT) :: RTHRATEN
4621 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
4622 INTENT(INOUT) :: RTHCUTEN, &
4623 RQVCUTEN, &
4624 RQCCUTEN, &
4625 RQRCUTEN, &
4626 RQICUTEN, &
4627 RQSCUTEN
4629 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4630 INTENT(INOUT) :: RUBLTEN, &
4631 RVBLTEN, &
4632 RTHBLTEN, &
4633 RQVBLTEN, &
4634 RQCBLTEN, &
4635 RQIBLTEN
4637 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4638 INTENT(INOUT) :: RTHFTEN, &
4639 RQVFTEN
4641 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
4642 INTENT(INOUT) :: RUNDGDTEN, &
4643 RVNDGDTEN, &
4644 RTHNDGDTEN, &
4645 RPHNDGDTEN, &
4646 RQVNDGDTEN
4648 REAL, DIMENSION( ims:ime, jms:jme ) , &
4649 INTENT(INOUT) :: RMUNDGDTEN
4651 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4652 INTEGER :: i, j, k
4653 REAL :: w1, w2, z0, z1, z2
4654 REAL :: e_vapor
4656 !-----------------------------------------------------------------------
4658 !<DESCRIPTION>
4659 !
4660 ! phys_prep calculates a number of diagnostic quantities needed by
4661 ! the physics routines. It also decouples the physics tendencies from
4662 ! the column dry-air mass (the physics routines expect to see/update the
4663 ! uncoupled tendencies).
4664 !
4665 !</DESCRIPTION>
4667 ! set up loop bounds for this grid's boundary conditions
4669 i_start = its
4670 i_end = min( ite,ide-1 )
4671 j_start = jts
4672 j_end = min( jte,jde-1 )
4674 k_start = kts
4675 k_end = min( kte, kde-1 )
4677 ! compute thermodynamics and velocities at pressure points (or half levels)
4679 do j = j_start,j_end
4680 do k = k_start, k_end
4681 do i = i_start, i_end
4683 th_phy(i,k,j) = t(i,k,j) + t0
4684 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4685 pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4686 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4687 rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4688 u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4689 v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4691 enddo
4692 enddo
4693 enddo
4695 ! compute z at w points
4697 do j = j_start,j_end
4698 do k = k_start, kte
4699 do i = i_start, i_end
4700 z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4701 enddo
4702 enddo
4703 enddo
4705 do j = j_start,j_end
4706 do k = k_start, kte-1
4707 do i = i_start, i_end
4708 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4709 enddo
4710 enddo
4711 enddo
4713 do j = j_start,j_end
4714 do i = i_start, i_end
4715 dz8w(i,kte,j) = 0.
4716 enddo
4717 enddo
4719 ! compute z at p points or half levels (average of z at full levels)
4721 do j = j_start,j_end
4722 do k = k_start, k_end
4723 do i = i_start, i_end
4724 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4725 enddo
4726 enddo
4727 enddo
4729 ! interp t and p to full levels
4731 do j = j_start,j_end
4732 do k = 2, k_end
4733 do i = i_start, i_end
4734 p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4735 t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4736 enddo
4737 enddo
4738 enddo
4740 ! extrapolate p and t to surface and top.
4741 ! we'll use an extrapolation in z for now
4743 do j = j_start,j_end
4744 do i = i_start, i_end
4746 ! bottom
4748 z0 = z_at_w(i,1,j)
4749 z1 = z(i,1,j)
4750 z2 = z(i,2,j)
4751 w1 = (z0 - z2)/(z1 - z2)
4752 w2 = 1. - w1
4753 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4754 t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4756 ! top
4758 z0 = z_at_w(i,kte,j)
4759 z1 = z(i,k_end,j)
4760 z2 = z(i,k_end-1,j)
4761 w1 = (z0 - z2)/(z1 - z2)
4762 w2 = 1. - w1
4764 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4765 !!! bug fix extrapolate ln(p) so p is positive definite
4766 p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4767 t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4769 enddo
4770 enddo
4772 ! calculate hydrostatic pressure at both full and half levels
4773 ! first, full level p: assuming dry over model top
4775 do j = j_start,j_end
4776 do i = i_start, i_end
4777 p_hyd_w(i,kte,j) = p_top
4778 enddo
4779 enddo
4781 e_vapor = 0.
4782 do j = j_start,j_end
4783 do k = kte-1, k_start, -1
4784 do i = i_start, i_end
4785 ! e_vapor = 1./alt(i,k,j)*moist(i,k,j,P_QV)*g*dz8w(i,k,j)
4786 ! p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+e_vapor
4787 p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))*g*dz8w(i,k,j)
4788 enddo
4789 enddo
4790 enddo
4792 ! now calculate hydrostatic pressure at half levels
4794 do j = j_start,j_end
4795 do k = k_start, k_end
4796 do i = i_start, i_end
4797 p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4798 enddo
4799 enddo
4800 enddo
4802 ! decouple all physics tendencies
4804 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4806 DO J=j_start,j_end
4807 DO K=k_start,k_end
4808 DO I=i_start,i_end
4809 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4810 ENDDO
4811 ENDDO
4812 ENDDO
4814 ENDIF
4816 IF (config_flags%cu_physics .gt. 0) THEN
4818 DO J=j_start,j_end
4819 DO I=i_start,i_end
4820 DO K=k_start,k_end
4821 RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4822 ENDDO
4823 ENDDO
4824 ENDDO
4826 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4827 DO J=j_start,j_end
4828 DO I=i_start,i_end
4829 DO K=k_start,k_end
4830 RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4831 ENDDO
4832 ENDDO
4833 ENDDO
4834 ENDIF
4836 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4837 DO J=j_start,j_end
4838 DO I=i_start,i_end
4839 DO K=k_start,k_end
4840 RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4841 ENDDO
4842 ENDDO
4843 ENDDO
4844 ENDIF
4846 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4847 DO J=j_start,j_end
4848 DO I=i_start,i_end
4849 DO K=k_start,k_end
4850 RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4851 ENDDO
4852 ENDDO
4853 ENDDO
4854 ENDIF
4856 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4857 DO J=j_start,j_end
4858 DO I=i_start,i_end
4859 DO K=k_start,k_end
4860 RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4861 ENDDO
4862 ENDDO
4863 ENDDO
4864 ENDIF
4866 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4867 DO J=j_start,j_end
4868 DO I=i_start,i_end
4869 DO K=k_start,k_end
4870 RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4871 ENDDO
4872 ENDDO
4873 ENDDO
4874 ENDIF
4876 ENDIF
4878 IF (config_flags%bl_pbl_physics .gt. 0) THEN
4880 DO J=j_start,j_end
4881 DO K=k_start,k_end
4882 DO I=i_start,i_end
4883 RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4884 RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4885 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4886 ENDDO
4887 ENDDO
4888 ENDDO
4890 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4891 DO J=j_start,j_end
4892 DO K=k_start,k_end
4893 DO I=i_start,i_end
4894 RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4895 ENDDO
4896 ENDDO
4897 ENDDO
4898 ENDIF
4900 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4901 DO J=j_start,j_end
4902 DO K=k_start,k_end
4903 DO I=i_start,i_end
4904 RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4905 ENDDO
4906 ENDDO
4907 ENDDO
4908 ENDIF
4910 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4911 DO J=j_start,j_end
4912 DO K=k_start,k_end
4913 DO I=i_start,i_end
4914 RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4915 ENDDO
4916 ENDDO
4917 ENDDO
4918 ENDIF
4920 ENDIF
4922 ! decouple advective forcing required by Grell-Devenyi scheme
4924 if(( config_flags%cu_physics == GDSCHEME ) .OR. &
4925 ( config_flags%cu_physics == G3SCHEME )) then
4927 DO J=j_start,j_end
4928 DO I=i_start,i_end
4929 DO K=k_start,k_end
4930 RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4931 ENDDO
4932 ENDDO
4933 ENDDO
4935 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4936 DO J=j_start,j_end
4937 DO I=i_start,i_end
4938 DO K=k_start,k_end
4939 RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4940 ENDDO
4941 ENDDO
4942 ENDDO
4943 ENDIF
4945 END IF
4947 ! fdda
4948 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4949 ! so only decouple those
4951 IF (config_flags%grid_fdda .gt. 0) THEN
4953 i_startu=MAX(its,ids+1)
4954 j_startv=MAX(jts,jds+1)
4956 DO J=j_start,j_end
4957 DO K=k_start,k_end
4958 DO I=i_startu,i_end
4959 RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4960 ENDDO
4961 ENDDO
4962 ENDDO
4963 DO J=j_startv,j_end
4964 DO K=k_start,k_end
4965 DO I=i_start,i_end
4966 RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4967 ENDDO
4968 ENDDO
4969 ENDDO
4970 DO J=j_start,j_end
4971 DO K=k_start,k_end
4972 DO I=i_start,i_end
4973 RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4974 ! RMUNDGDTEN(I,J) - no coupling
4975 ENDDO
4976 ENDDO
4977 ENDDO
4979 IF (config_flags%grid_fdda .EQ. 2) THEN
4980 DO J=j_start,j_end
4981 DO K=k_start,kte
4982 DO I=i_start,i_end
4983 RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
4984 ENDDO
4985 ENDDO
4986 ENDDO
4988 ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
4989 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4990 DO J=j_start,j_end
4991 DO K=k_start,k_end
4992 DO I=i_start,i_end
4993 RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
4994 ENDDO
4995 ENDDO
4996 ENDDO
4997 ENDIF
4998 ENDIF
5000 ENDIF
5002 END SUBROUTINE phy_prep
5004 !------------------------------------------------------------
5006 SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5007 p, p8w, p0, pb, ph, phb, &
5008 th_phy, pii, pf, &
5009 z, z_at_w, dz8w, &
5010 dt,h_diabatic, &
5011 config_flags,fzm, fzp, &
5012 ids,ide, jds,jde, kds,kde, &
5013 ims,ime, jms,jme, kms,kme, &
5014 its,ite, jts,jte, kts,kte )
5016 IMPLICIT NONE
5018 ! Here we construct full fields
5019 ! needed by the microphysics
5021 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5023 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5024 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5025 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5027 REAL, INTENT(IN ) :: dt
5029 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5030 INTENT(IN ) :: al, &
5031 alb, &
5032 p, &
5033 pb, &
5034 ph, &
5035 phb
5038 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5039 fzp
5041 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5042 INTENT( OUT) :: rho, &
5043 th_phy, &
5044 pii, &
5045 pf, &
5046 z, &
5047 z_at_w, &
5048 dz8w, &
5049 p8w
5051 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5052 INTENT(INOUT) :: h_diabatic
5054 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5055 INTENT(INOUT) :: t_new, &
5056 t_old
5058 REAL, INTENT(IN ) :: t0, p0
5059 REAL :: z0,z1,z2,w1,w2
5061 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5062 INTEGER :: i, j, k
5064 !--------------------------------------------------------------------
5066 !<DESCRIPTION>
5067 !
5068 ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
5069 ! the microphysics routines.
5070 !
5071 !</DESCRIPTION>
5073 ! set up loop bounds for this grid's boundary conditions
5075 i_start = its
5076 i_end = min( ite,ide-1 )
5077 j_start = jts
5078 j_end = min( jte,jde-1 )
5080 k_start = kts
5081 k_end = min( kte, kde-1 )
5083 DO j = j_start, j_end
5084 DO k = k_start, kte
5085 DO i = i_start, i_end
5086 z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5087 ENDDO
5088 ENDDO
5089 ENDDO
5091 do j = j_start,j_end
5092 do k = k_start, kte-1
5093 do i = i_start, i_end
5094 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5095 enddo
5096 enddo
5097 enddo
5099 do j = j_start,j_end
5100 do i = i_start, i_end
5101 dz8w(i,kte,j) = 0.
5102 enddo
5103 enddo
5106 ! compute full pii, rho, and z at the new time-level
5107 ! (needed for physics).
5108 ! convert perturbation theta to full theta (th_phy)
5109 ! use h_diabatic to temporarily save pre-microphysics full theta
5111 DO j = j_start, j_end
5112 DO k = k_start, k_end
5113 DO i = i_start, i_end
5115 #ifdef REVERT
5116 t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
5117 #endif
5118 th_phy(i,k,j) = t_new(i,k,j) + t0
5119 h_diabatic(i,k,j) = th_phy(i,k,j)
5120 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
5121 pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5122 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5123 pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5125 ENDDO
5126 ENDDO
5127 ENDDO
5129 ! interp t and p at w points
5131 do j = j_start,j_end
5132 do k = 2, k_end
5133 do i = i_start, i_end
5134 p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
5135 enddo
5136 enddo
5137 enddo
5139 ! extrapolate p and t to surface and top.
5140 ! we'll use an extrapolation in z for now
5142 do j = j_start,j_end
5143 do i = i_start, i_end
5145 ! bottom
5147 z0 = z_at_w(i,1,j)
5148 z1 = z(i,1,j)
5149 z2 = z(i,2,j)
5150 w1 = (z0 - z2)/(z1 - z2)
5151 w2 = 1. - w1
5152 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
5154 ! top
5156 z0 = z_at_w(i,kte,j)
5157 z1 = z(i,k_end,j)
5158 z2 = z(i,k_end-1,j)
5159 w1 = (z0 - z2)/(z1 - z2)
5160 w2 = 1. - w1
5161 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
5162 p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
5164 enddo
5165 enddo
5167 END SUBROUTINE moist_physics_prep_em
5169 !------------------------------------------------------------------------------
5171 SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
5172 th_phy, h_diabatic, dt, &
5173 config_flags, &
5174 #if ( WRF_DFI_RADAR == 1 )
5175 dfi_tten_rad,dfi_stage, &
5176 #endif
5177 ids,ide, jds,jde, kds,kde, &
5178 ims,ime, jms,jme, kms,kme, &
5179 its,ite, jts,jte, kts,kte )
5181 IMPLICIT NONE
5183 ! Here we construct full fields
5184 ! needed by the microphysics
5186 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5188 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5189 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5190 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5192 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5193 INTENT(INOUT) :: t_new, &
5194 t_old, &
5195 th_phy, &
5196 h_diabatic
5197 #if ( WRF_DFI_RADAR == 1 )
5198 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5199 INTENT(IN), OPTIONAL :: dfi_tten_rad
5200 INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
5201 ! REAL :: dfi_tten_max, old_max
5202 #endif
5204 REAL mpten, mptenmax, mptenmin
5206 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
5209 REAL, INTENT(IN ) :: t0, dt
5211 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5212 INTEGER :: i, j, k, imax, jmax, imin, jmin
5214 !--------------------------------------------------------------------
5216 !<DESCRIPTION>
5217 !
5218 ! moist_phys_finish_em resets theta to its perturbation value and
5219 ! computes and stores the microphysics diabatic heating term.
5220 !
5221 !</DESCRIPTION>
5223 ! set up loop bounds for this grid's boundary conditions
5225 i_start = its
5226 i_end = min( ite,ide-1 )
5227 j_start = jts
5228 j_end = min( jte,jde-1 )
5229 ! i_start=max(its,ids+4)
5230 ! i_end=min(ite,ide-5)
5231 ! j_start=max(jts,jds+4)
5232 ! j_end=min(jte,jde-5)
5234 k_start = kts
5235 k_end = min( kte, kde-1 )
5237 #if ( WRF_DFI_RADAR == 1 )
5238 if (config_flags%dfi_radar == 1 .and. PRESENT(dfi_stage) .and. dfi_stage ==DFI_FWD &
5239 .and. PRESENT(dfi_tten_rad) ) then
5240 WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5241 CALL wrf_debug ( 50 , TRIM(wrf_err_message) )
5242 endif
5243 dfi_tten_max=-999
5244 old_max=-999
5245 #endif
5247 ! add microphysics theta diff to perturbation theta, set h_diabatic
5249 IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5250 DO j = j_start, j_end
5251 DO k = k_start, k_end
5252 DO i = i_start, i_end
5253 #if ( WRF_DFI_RADAR == 1 )
5254 mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5256 mpten=min(config_flags%mp_tend_lim*dt, mpten)
5257 mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5259 if (config_flags%dfi_radar == 1 .and. PRESENT(dfi_stage) .and. &
5260 dfi_stage == DFI_FWD .and. PRESENT(dfi_tten_rad) .and. &
5261 dfi_tten_rad(i,k,j) >= 1.0e-7 .and. dfi_tten_rad(i,k,j) <= 10. &
5262 .and. k < k_end ) then
5263 ! add radar temp tendency
5264 !tgs if(dfi_tten_rad(i,k,j) > th_phy(i,k,j)-h_diabatic(i,k,j) ) then
5265 if(dfi_tten_rad(i,k,j) > mpten ) then
5266 t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))
5267 else
5268 !tgs t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
5269 t_new(i,k,j) = t_new(i,k,j) + mpten
5270 endif
5271 else
5272 !tgs t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
5273 t_new(i,k,j) = t_new(i,k,j) + mpten
5274 endif
5275 #else
5276 t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
5277 #endif
5278 h_diabatic(i,k,j) = (th_phy(i,k,j)-h_diabatic(i,k,j))/dt
5279 ENDDO
5280 ENDDO
5281 ENDDO
5283 ELSE
5285 DO j = j_start, j_end
5286 DO k = k_start, k_end
5287 DO i = i_start, i_end
5288 ! t_new(i,k,j) = t_new(i,k,j)
5289 h_diabatic(i,k,j) = 0.
5290 ENDDO
5291 ENDDO
5292 ENDDO
5293 ENDIF
5295 END SUBROUTINE moist_physics_finish_em
5297 !----------------------------------------------------------------
5300 SUBROUTINE init_module_big_step
5301 END SUBROUTINE init_module_big_step
5303 SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5304 ids, ide, jds, jde, kds, kde, &
5305 ims, ime, jms, jme, kms, kme, &
5306 its, ite, jts, jte, kts, kte )
5308 IMPLICIT NONE
5310 ! Input data
5312 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5313 ims, ime, jms, jme, kms, kme, &
5314 its, ite, jts, jte, kts, kte
5316 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5318 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5320 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5322 ! Local data
5324 INTEGER :: i, j, k, itf, jtf, ktf
5326 !<DESCRIPTION>
5327 !
5328 ! set_tend copies the advective tendency array into the tendency array.
5329 !
5330 !</DESCRIPTION>
5332 jtf = MIN(jte,jde-1)
5333 ktf = MIN(kte,kde-1)
5334 itf = MIN(ite,ide-1)
5335 DO j = jts, jtf
5336 DO k = kts, ktf
5337 DO i = its, itf
5338 field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5339 ENDDO
5340 ENDDO
5341 ENDDO
5343 END SUBROUTINE set_tend
5345 !------------------------------------------------------------------------------
5347 SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5348 rw_tendf, t_tendf, &
5349 u, v, w, t, t_init, &
5350 mut, muu, muv, ph, phb, &
5351 u_base, v_base, t_base, z_base, &
5352 dampcoef, zdamp, &
5353 ids, ide, jds, jde, kds, kde, &
5354 ims, ime, jms, jme, kms, kme, &
5355 its, ite, jts, jte, kts, kte )
5357 ! History: Apr 2005 Modifications by George Bryan, NCAR:
5358 ! - Generalized the code in a way that allows for
5359 ! simulations with steep terrain.
5360 !
5361 ! Jul 2004 Modifications by George Bryan, NCAR:
5362 ! - Modified the code to use u_base, v_base, and t_base
5363 ! arrays for the background state. Removed the hard-wired
5364 ! base-state values.
5365 ! - Modified the code to use dampcoef, zdamp, and damp_opt,
5366 ! i.e., the upper-level damper variables in namelist.input.
5367 ! Removed the hard-wired variables in the older version.
5368 ! This damper is used when damp_opt = 2.
5369 ! - Modified the code to account for the movement of the
5370 ! model surfaces with time. The code now obtains a base-
5371 ! state value by interpolation using the "_base" arrays.
5373 ! Nov 2003 Bug fix by Jason Knievel, NCAR
5375 ! Aug 2003 Meridional dimension, some comments, and
5376 ! changes in layout of the code added by
5377 ! Jason Knievel, NCAR
5379 ! Jul 2003 Original code by Bill Skamarock, NCAR
5381 ! Purpose: This routine applies Rayleigh damping to a layer at top
5382 ! of the model domain.
5384 !-----------------------------------------------------------------------
5385 ! Begin declarations.
5387 IMPLICIT NONE
5389 INTEGER, INTENT( IN ) &
5390 :: ids, ide, jds, jde, kds, kde, &
5391 ims, ime, jms, jme, kms, kme, &
5392 its, ite, jts, jte, kts, kte
5394 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5395 :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5397 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5398 :: u, v, w, t, t_init, ph, phb
5400 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5401 :: mut, muu, muv
5403 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5404 :: u_base, v_base, t_base, z_base
5406 REAL, INTENT(IN ) &
5407 :: dampcoef, zdamp
5409 ! Local variables.
5411 INTEGER &
5412 :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5414 REAL &
5415 :: pii, dcoef, z, ztop
5417 REAL :: wkp1, wk, wkm1
5419 REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5421 ! End declarations.
5422 !-----------------------------------------------------------------------
5424 pii = 2.0 * asin(1.0)
5426 ktf = MIN( kte, kde-1 )
5428 !-----------------------------------------------------------------------
5429 ! Adjust u to base state.
5431 DO j = jts, MIN( jte, jde-1 )
5432 DO i = its, MIN( ite, ide )
5434 ! Get height at top of model
5435 ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5436 +ph(i ,kde,j)+ph(i-1,kde,j) )/g
5438 ! Find bottom of damping layer
5439 k1 = ktf
5440 z = ztop
5441 DO WHILE( z >= (ztop-zdamp) )
5442 z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5443 +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5444 +ph(i ,k1,j)+ph(i ,k1+1,j) &
5445 +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5446 z00(k1) = z
5447 k1 = k1 - 1
5448 ENDDO
5449 k1 = k1 + 2
5451 ! Get reference state at model levels
5452 DO k = k1, ktf
5453 k2 = ktf
5454 DO WHILE( z_base(k2) .gt. z00(k) )
5455 k2 = k2 - 1
5456 ENDDO
5457 if(k2+1.gt.ktf)then
5458 u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5459 * ( z00(k) - z_base(k2) ) &
5460 / ( z_base(k2) - z_base(k2-1) )
5461 else
5462 u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5463 * ( z00(k) - z_base(k2) ) &
5464 / ( z_base(k2+1) - z_base(k2) )
5465 endif
5466 ENDDO
5468 ! Apply the Rayleigh damper
5469 DO k = k1, ktf
5470 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5471 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5472 ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5473 muu(i,j) * ( dcoef * dampcoef ) * &
5474 ( u(i,k,j) - u00(k) )
5475 END DO
5477 END DO
5478 END DO
5480 ! End adjustment of u.
5481 !-----------------------------------------------------------------------
5483 !-----------------------------------------------------------------------
5484 ! Adjust v to base state.
5486 DO j = jts, MIN( jte, jde )
5487 DO i = its, MIN( ite, ide-1 )
5489 ! Get height at top of model
5490 ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5491 +ph(i,kde,j )+ph(i,kde,j-1) )/g
5493 ! Find bottom of damping layer
5494 k1 = ktf
5495 z = ztop
5496 DO WHILE( z >= (ztop-zdamp) )
5497 z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5498 +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5499 +ph(i,k1,j )+ph(i,k1+1,j ) &
5500 +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5501 z00(k1) = z
5502 k1 = k1 - 1
5503 ENDDO
5504 k1 = k1 + 2
5506 ! Get reference state at model levels
5507 DO k = k1, ktf
5508 k2 = ktf
5509 DO WHILE( z_base(k2) .gt. z00(k) )
5510 k2 = k2 - 1
5511 ENDDO
5512 if(k2+1.gt.ktf)then
5513 v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
5514 * ( z00(k) - z_base(k2) ) &
5515 / ( z_base(k2) - z_base(k2-1) )
5516 else
5517 v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
5518 * ( z00(k) - z_base(k2) ) &
5519 / ( z_base(k2+1) - z_base(k2) )
5520 endif
5521 ENDDO
5523 ! Apply the Rayleigh damper
5524 DO k = k1, ktf
5525 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5526 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5527 rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
5528 muv(i,j) * ( dcoef * dampcoef ) * &
5529 ( v(i,k,j) - v00(k) )
5530 END DO
5532 END DO
5533 END DO
5535 ! End adjustment of v.
5536 !-----------------------------------------------------------------------
5538 !-----------------------------------------------------------------------
5539 ! Adjust w to base state.
5541 DO j = jts, MIN( jte, jde-1 )
5542 DO i = its, MIN( ite, ide-1 )
5543 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5544 DO k = kts, MIN( kte, kde )
5545 z = ( phb(i,k,j) + ph(i,k,j) ) / g
5546 IF ( z >= (ztop-zdamp) ) THEN
5547 dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5548 dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5549 rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
5550 mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5551 END IF
5552 END DO
5553 END DO
5554 END DO
5556 ! End adjustment of w.
5557 !-----------------------------------------------------------------------
5559 !-----------------------------------------------------------------------
5560 ! Adjust potential temperature to base state.
5562 DO j = jts, MIN( jte, jde-1 )
5563 DO i = its, MIN( ite, ide-1 )
5565 ! Get height at top of model
5566 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5568 ! Find bottom of damping layer
5569 k1 = ktf
5570 z = ztop
5571 DO WHILE( z >= (ztop-zdamp) )
5572 z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
5573 ph(i,k1,j) + ph(i,k1+1,j) ) / g
5574 z00(k1) = z
5575 k1 = k1 - 1
5576 ENDDO
5577 k1 = k1 + 2
5579 ! Get reference state at model levels
5580 DO k = k1, ktf
5581 k2 = ktf
5582 DO WHILE( z_base(k2) .gt. z00(k) )
5583 k2 = k2 - 1
5584 ENDDO
5585 if(k2+1.gt.ktf)then
5586 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
5587 * ( z00(k) - z_base(k2) ) &
5588 / ( z_base(k2) - z_base(k2-1) )
5589 else
5590 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
5591 * ( z00(k) - z_base(k2) ) &
5592 / ( z_base(k2+1) - z_base(k2) )
5593 endif
5594 ENDDO
5596 ! Apply the Rayleigh damper
5597 DO k = k1, ktf
5598 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5599 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5600 t_tendf(i,k,j) = t_tendf(i,k,j) - &
5601 mut(i,j) * ( dcoef * dampcoef ) * &
5602 ( t(i,k,j) - t00(k) )
5603 END DO
5605 END DO
5606 END DO
5608 ! End adjustment of potential temperature.
5609 !-----------------------------------------------------------------------
5611 END SUBROUTINE rk_rayleigh_damp
5613 !==============================================================================
5614 !==============================================================================
5616 SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt, &
5617 config_flags, &
5618 diff_6th_opt, diff_6th_factor, &
5619 ids, ide, jds, jde, kds, kde, &
5620 ims, ime, jms, jme, kms, kme, &
5621 its, ite, jts, jte, kts, kte )
5623 ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
5624 ! 07 Jun 2006 Revised and generalized by Jason Knievel
5625 ! 25 Apr 2005 Original code by Jason Knievel, NCAR
5627 ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
5628 ! diffusion to 3-d velocity and to scalars.
5630 ! References: Ming Xue (MWR Aug 2000)
5631 ! Durran ("Numerical Methods for Wave Equations..." 1999)
5632 ! George Bryan (personal communication)
5634 !------------------------------------------------------------------------------
5635 ! Begin: Declarations.
5637 IMPLICIT NONE
5639 INTEGER, INTENT(IN) &
5640 :: ids, ide, jds, jde, kds, kde, &
5641 ims, ime, jms, jme, kms, kme, &
5642 its, ite, jts, jte, kts, kte
5644 TYPE(grid_config_rec_type), INTENT(IN) &
5645 :: config_flags
5647 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
5648 :: tendency
5650 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
5651 :: field
5653 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
5654 :: mu
5656 REAL, INTENT(IN) &
5657 :: dt
5659 REAL, INTENT(IN) &
5660 :: diff_6th_factor
5662 INTEGER, INTENT(IN) &
5663 :: diff_6th_opt
5665 CHARACTER(LEN=1) , INTENT(IN) &
5666 :: name
5668 INTEGER &
5669 :: i, j, k, &
5670 i_start, i_end, &
5671 j_start, j_end, &
5672 k_start, k_end, &
5673 ktf
5675 REAL &
5676 :: dflux_x_p0, dflux_y_p0, &
5677 dflux_x_p1, dflux_y_p1, &
5678 tendency_x, tendency_y, &
5679 mu_avg_p0, mu_avg_p1, &
5680 diff_6th_coef
5682 LOGICAL &
5683 :: specified
5685 ! End: Declarations.
5686 !------------------------------------------------------------------------------
5688 !------------------------------------------------------------------------------
5689 ! Begin: Translate the diffusion factor into a diffusion coefficient. See
5690 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5691 ! fourth) and for diffusion in two dimensions (not one). For reference, a
5692 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5693 ! although application of the flux limiter reduces somewhat the effects of
5694 ! diffusion for a given coefficient.
5696 diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
5698 ! End: Translate diffusion factor.
5699 !------------------------------------------------------------------------------
5701 !------------------------------------------------------------------------------
5702 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5703 ! The halo regions are already filled with values by the time this subroutine
5704 ! is called, which allows the stencil to extend beyond the domains' edges.
5706 ktf = MIN( kte, kde-1 )
5708 IF ( name .EQ. 'u' ) THEN
5710 i_start = its
5711 i_end = ite
5712 j_start = jts
5713 j_end = MIN(jde-1,jte)
5714 k_start = kts
5715 k_end = ktf
5717 ELSE IF ( name .EQ. 'v' ) THEN
5719 i_start = its
5720 i_end = MIN(ide-1,ite)
5721 j_start = jts
5722 j_end = jte
5723 k_start = kts
5724 k_end = ktf
5726 ELSE IF ( name .EQ. 'w' ) THEN
5728 i_start = its
5729 i_end = MIN(ide-1,ite)
5730 j_start = jts
5731 j_end = MIN(jde-1,jte)
5732 k_start = kts+1
5733 k_end = ktf
5735 ELSE
5737 i_start = its
5738 i_end = MIN(ide-1,ite)
5739 j_start = jts
5740 j_end = MIN(jde-1,jte)
5741 k_start = kts
5742 k_end = ktf
5744 ENDIF
5746 ! End: Assignment of limits of spatial loops.
5747 !------------------------------------------------------------------------------
5749 !------------------------------------------------------------------------------
5750 ! Begin: Loop across spatial dimensions.
5752 DO j = j_start, j_end
5753 DO k = k_start, k_end
5754 DO i = i_start, i_end
5756 !------------------------------------------------------------------------------
5757 ! Begin: Diffusion in x (i index).
5759 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5761 dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
5762 - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
5763 + ( field(i+2,k,j) - field(i-3,k,j) ) )
5765 dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
5766 - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
5767 + ( field(i+3,k,j) - field(i-2,k,j) ) )
5769 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5770 ! (variation on Xue's eq. 10).
5772 IF ( diff_6th_opt .EQ. 2 ) THEN
5774 IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5775 dflux_x_p0 = 0.0
5776 END IF
5778 IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
5779 dflux_x_p1 = 0.0
5780 END IF
5782 END IF
5784 ! Apply 6th-order diffusion in x direction.
5786 IF ( name .EQ. 'u' ) THEN
5787 mu_avg_p0 = mu(i-1,j)
5788 mu_avg_p1 = mu(i ,j)
5789 ELSE IF ( name .EQ. 'v' ) THEN
5790 mu_avg_p0 = 0.25 * ( &
5791 mu(i-1,j-1) + &
5792 mu(i ,j-1) + &
5793 mu(i-1,j ) + &
5794 mu(i ,j ) )
5795 mu_avg_p1 = 0.25 * ( &
5796 mu(i ,j-1) + &
5797 mu(i+1,j-1) + &
5798 mu(i ,j ) + &
5799 mu(i+1,j ) )
5800 ELSE
5801 mu_avg_p0 = 0.5 * ( &
5802 mu(i-1,j) + &
5803 mu(i ,j) )
5804 mu_avg_p1 = 0.5 * ( &
5805 mu(i ,j) + &
5806 mu(i+1,j) )
5807 END IF
5809 tendency_x = diff_6th_coef * &
5810 ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5812 ! End: Diffusion in x.
5813 !------------------------------------------------------------------------------
5815 !------------------------------------------------------------------------------
5816 ! Begin: Diffusion in y (j index).
5818 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5820 dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
5821 - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
5822 + ( field(i,k,j+2) - field(i,k,j-3) ) )
5824 dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
5825 - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
5826 + ( field(i,k,j+3) - field(i,k,j-2) ) )
5828 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5829 ! (variation on Xue's eq. 10).
5831 IF ( diff_6th_opt .EQ. 2 ) THEN
5833 IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5834 dflux_y_p0 = 0.0
5835 END IF
5837 IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
5838 dflux_y_p1 = 0.0
5839 END IF
5841 END IF
5843 ! Apply 6th-order diffusion in y direction.
5845 IF ( name .EQ. 'u' ) THEN
5846 mu_avg_p0 = 0.25 * ( &
5847 mu(i-1,j-1) + &
5848 mu(i ,j-1) + &
5849 mu(i-1,j ) + &
5850 mu(i ,j ) )
5851 mu_avg_p1 = 0.25 * ( &
5852 mu(i-1,j ) + &
5853 mu(i ,j ) + &
5854 mu(i-1,j+1) + &
5855 mu(i ,j+1) )
5856 ELSE IF ( name .EQ. 'v' ) THEN
5857 mu_avg_p0 = mu(i,j-1)
5858 mu_avg_p1 = mu(i,j )
5859 ELSE
5860 mu_avg_p0 = 0.5 * ( &
5861 mu(i,j-1) + &
5862 mu(i,j ) )
5863 mu_avg_p1 = 0.5 * ( &
5864 mu(i,j ) + &
5865 mu(i,j+1) )
5866 END IF
5868 tendency_y = diff_6th_coef * &
5869 ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5871 ! End: Diffusion in y.
5872 !------------------------------------------------------------------------------
5874 !------------------------------------------------------------------------------
5875 ! Begin: Combine diffusion in x and y.
5877 tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5879 ! End: Combine diffusion in x and y.
5880 !------------------------------------------------------------------------------
5882 ENDDO
5883 ENDDO
5884 ENDDO
5886 ! End: Loop across spatial dimensions.
5887 !------------------------------------------------------------------------------
5889 END SUBROUTINE sixth_order_diffusion
5891 !==============================================================================
5892 !==============================================================================
5894 END MODULE module_big_step_utilities_em