1 !WRF:MODEL_LAYER:DYNAMICS
3 MODULE module_advect_em
6 USE module_model_constants
12 SUBROUTINE mass_flux_divergence ( field, field_old, tendency, &
15 msfux, msfuy, msfvx, msfvy, &
19 ids, ide, jds, jde, kds, kde, &
20 ims, ime, jms, jme, kms, kme, &
21 its, ite, jts, jte, kts, kte )
27 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
29 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
30 ims, ime, jms, jme, kms, kme, &
31 its, ite, jts, jte, kts, kte
33 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
39 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
40 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
42 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
49 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
53 REAL , INTENT(IN ) :: rdx, &
58 INTEGER :: i, j, k, itf, jtf, ktf
59 INTEGER :: i_start, i_end, j_start, j_end
60 INTEGER :: imin, imax, jmin, jmax
62 REAL :: mrdx, mrdy, ub, vb, uw, vw
63 REAL , DIMENSION(its:ite,kts:kte) :: vflux
67 !--------------- horizontal flux
70 if(config_flags%specified .or. config_flags%nested) specified = .true.
74 i_end = MIN(ite,ide-1)
76 j_end = MIN(jte,jde-1)
82 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
83 *(ru(i+1,k,j)*(field(i+1,k,j)+field(i ,k,j)) &
84 -ru(i ,k,j)*(field(i ,k,j)+field(i-1,k,j)))
93 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
94 *(rv(i,k,j+1)*(field(i,k,j+1)+field(i,k,j )) &
95 -rv(i,k,j )*(field(i,k,j )+field(i,k,j-1)))
100 !---------------- vertical flux divergence
103 DO i = i_start, i_end
108 DO j = j_start, j_end
111 DO i = i_start, i_end
112 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
117 DO i = i_start, i_end
118 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
124 END SUBROUTINE mass_flux_divergence
126 !-------------------------------------------------------------------------------
128 SUBROUTINE advect_u ( u, u_old, tendency, &
130 mut, time_step, config_flags, &
131 msfux, msfuy, msfvx, msfvy, &
135 ids, ide, jds, jde, kds, kde, &
136 ims, ime, jms, jme, kms, kme, &
137 its, ite, jts, jte, kts, kte )
143 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
145 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
146 ims, ime, jms, jme, kms, kme, &
147 its, ite, jts, jte, kts, kte
149 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, &
155 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
156 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
158 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
165 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
169 REAL , INTENT(IN ) :: rdx, &
171 INTEGER , INTENT(IN ) :: time_step
175 INTEGER :: i, j, k, itf, jtf, ktf
176 INTEGER :: i_start, i_end, j_start, j_end
177 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
178 INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
179 INTEGER :: jp1, jp0, jtmp
181 INTEGER :: horz_order, vert_order
183 REAL :: mrdx, mrdy, ub, vb, uw, vw, dvm, dvp
184 REAL , DIMENSION(its:ite, kts:kte) :: vflux
187 REAL, DIMENSION( its-1:ite+1, kts:kte ) :: fqx
188 REAL, DIMENSION( its:ite, kts:kte, 2) :: fqy
190 LOGICAL :: degrade_xs, degrade_ys
191 LOGICAL :: degrade_xe, degrade_ye
193 ! definition of flux operators, 3rd, 4th, 5th or 6th order
195 REAL :: flux3, flux4, flux5, flux6
196 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
198 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
199 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
201 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
202 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
203 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
205 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
206 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
207 +(q_ip2+q_im3) )/60.0
209 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
210 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
211 -sign(1,time_step)*sign(1.,ua)*( &
212 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
218 if(config_flags%specified .or. config_flags%nested) specified = .true.
220 ! set order for vertical and horzontal flux operators
222 horz_order = config_flags%h_mom_adv_order
223 vert_order = config_flags%v_mom_adv_order
227 ! begin with horizontal flux divergence
229 horizontal_order_test : IF( horz_order == 6 ) THEN
231 ! determine boundary mods for flux operators
232 ! We degrade the flux operators from 3rd/4th order
233 ! to second order one gridpoint in from the boundaries for
234 ! all boundary conditions except periodic and symmetry - these
235 ! conditions have boundary zone data fill for correct application
236 ! of the higher order flux stencils
243 IF( config_flags%periodic_x .or. &
244 config_flags%symmetric_xs .or. &
245 (its > ids+2) ) degrade_xs = .false.
246 IF( config_flags%periodic_x .or. &
247 config_flags%symmetric_xe .or. &
248 (ite < ide-2) ) degrade_xe = .false.
249 IF( config_flags%periodic_y .or. &
250 config_flags%symmetric_ys .or. &
251 (jts > jds+2) ) degrade_ys = .false.
252 IF( config_flags%periodic_y .or. &
253 config_flags%symmetric_ye .or. &
254 (jte < jde-3) ) degrade_ye = .false.
256 !--------------- y - advection first
260 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
261 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
262 IF ( config_flags%periodic_x ) i_start = its
263 IF ( config_flags%periodic_x ) i_end = ite
266 j_end = MIN(jte,jde-1)
268 ! higher order flux has a 5 or 7 point stencil, so compute
269 ! bounds so we can switch to second order flux close to the boundary
275 j_start = MAX(jts,jds+1)
280 j_end = MIN(jte,jde-2)
284 IF(config_flags%polar) j_end = MIN(jte,jde-1)
286 ! compute fluxes, 5th or 6th order
291 j_loop_y_flux_6 : DO j = j_start, j_end+1
293 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
296 DO i = i_start, i_end
297 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
298 fqy( i, k, jp1 ) = vel*flux6( &
299 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), &
300 u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel )
304 ! we must be close to some boundary where we need to reduce the order of the stencil
306 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
309 DO i = i_start, i_end
310 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
311 *(u(i,k,j)+u(i,k,j-1))
315 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
318 DO i = i_start, i_end
319 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
320 fqy( i, k, jp1 ) = vel*flux4( &
321 u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel )
325 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
328 DO i = i_start, i_end
329 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
330 *(u(i,k,j)+u(i,k,j-1))
334 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
337 DO i = i_start, i_end
338 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
339 fqy( i, k, jp1 ) = vel*flux4( &
340 u(i,k,j-2),u(i,k,j-1), &
341 u(i,k,j),u(i,k,j+1),vel )
349 ! y flux-divergence into tendency
351 ! Comments for polar boundary condition
352 ! Flow is only from one side for points next to poles
353 ! S. pole at j=jds, N. pole at j=jde for v-stagger points
354 ! Tendencies affected are held at j=jds and j=jde-1 (non-stagger)
355 ! jp0 will always hold the flux from the south, and
356 ! jp1 will hold the flux from the north.
358 ! When j=jds+1 we are 1 in from S. pole, and jp1 contains fqy(jds+1), jp0 has fqy(jds)
359 ! tendency(j-1) = - mx/dy * (u rho v (jds+1)/mx - u rho v (jds)/mx)
361 ! tendency(j-1) = - mx/dy * (u rho v (jds+1)/mx) = - mx/dy * fqy(jp1)
363 ! When j=jde-1 we are 1 in from N. pole, and jp1 contains fqy(jde-1), jp0 has fqy(jde-2)
364 ! tendency(j-1) = - mx/dy * (u rho v (jde)/mx - u rho v (jde-1)/mx)
366 ! tendency(j-1) = + mx/dy * (u rho v (jde-1)/mx) = + mx/dy * fqy(jp0)
368 ! (j > j_start) will miss the u(,,jds) tendency
369 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
371 DO i = i_start, i_end
372 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
373 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
376 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
377 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
379 DO i = i_start, i_end
380 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
381 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
389 DO i = i_start, i_end
390 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
391 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
404 ENDDO j_loop_y_flux_6
406 ! next, x - flux divergence
412 j_end = MIN(jte,jde-1)
414 ! higher order flux has a 5 or 7 point stencil, so compute
415 ! bounds so we can switch to second order flux close to the boundary
421 i_start = MAX(ids+1,its)
426 i_end = MIN(ide-1,ite)
432 DO j = j_start, j_end
434 ! 5th or 6th order flux
437 DO i = i_start_f, i_end_f
438 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
439 fqx( i,k ) = vel*flux6( u(i-3,k,j), u(i-2,k,j), &
440 u(i-1,k,j), u(i ,k,j), &
441 u(i+1,k,j), u(i+2,k,j), &
446 ! lower order fluxes close to boundaries (if not periodic or symmetric)
447 ! specified uses upstream normal wind at boundaries
449 IF( degrade_xs ) THEN
451 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
455 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
456 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
463 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
464 fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
465 u(i ,k,j), u(i+1,k,j), &
471 IF( degrade_xe ) THEN
473 IF( i_end == ide-1 ) THEN ! second order flux next to the boundary
477 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
478 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
485 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
486 fqx( i,k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
487 u(i ,k,j), u(i+1,k,j), &
493 ! x flux-divergence into tendency
496 DO i = i_start, i_end
497 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
498 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
504 ELSE IF( horz_order == 5 ) THEN
506 ! 5th order horizontal flux calculation
507 ! This code is EXACTLY the same as the 6th order code
508 ! EXCEPT the 5th order and 3rd operators are used in
509 ! place of the 6th and 4th order operators
511 ! determine boundary mods for flux operators
512 ! We degrade the flux operators from 3rd/4th order
513 ! to second order one gridpoint in from the boundaries for
514 ! all boundary conditions except periodic and symmetry - these
515 ! conditions have boundary zone data fill for correct application
516 ! of the higher order flux stencils
523 IF( config_flags%periodic_x .or. &
524 config_flags%symmetric_xs .or. &
525 (its > ids+2) ) degrade_xs = .false.
526 IF( config_flags%periodic_x .or. &
527 config_flags%symmetric_xe .or. &
528 (ite < ide-2) ) degrade_xe = .false.
529 IF( config_flags%periodic_y .or. &
530 config_flags%symmetric_ys .or. &
531 (jts > jds+2) ) degrade_ys = .false.
532 IF( config_flags%periodic_y .or. &
533 config_flags%symmetric_ye .or. &
534 (jte < jde-3) ) degrade_ye = .false.
536 !--------------- y - advection first
540 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
541 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
542 IF ( config_flags%periodic_x ) i_start = its
543 IF ( config_flags%periodic_x ) i_end = ite
546 j_end = MIN(jte,jde-1)
548 ! higher order flux has a 5 or 7 point stencil, so compute
549 ! bounds so we can switch to second order flux close to the boundary
555 j_start = MAX(jts,jds+1)
560 j_end = MIN(jte,jde-2)
564 IF(config_flags%polar) j_end = MIN(jte,jde-1)
566 ! compute fluxes, 5th or 6th order
571 j_loop_y_flux_5 : DO j = j_start, j_end+1
573 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
576 DO i = i_start, i_end
577 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
578 fqy( i, k, jp1 ) = vel*flux5( &
579 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), &
580 u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel )
584 ! we must be close to some boundary where we need to reduce the order of the stencil
586 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
589 DO i = i_start, i_end
590 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
591 *(u(i,k,j)+u(i,k,j-1))
595 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
598 DO i = i_start, i_end
599 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
600 fqy( i, k, jp1 ) = vel*flux3( &
601 u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel )
605 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
608 DO i = i_start, i_end
609 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
610 *(u(i,k,j)+u(i,k,j-1))
614 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
617 DO i = i_start, i_end
618 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
619 fqy( i, k, jp1 ) = vel*flux3( &
620 u(i,k,j-2),u(i,k,j-1), &
621 u(i,k,j),u(i,k,j+1),vel )
627 ! y flux-divergence into tendency
629 ! (j > j_start) will miss the u(,,jds) tendency
630 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
632 DO i = i_start, i_end
633 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
634 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
637 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
638 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
640 DO i = i_start, i_end
641 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
642 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
650 DO i = i_start, i_end
651 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
652 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
665 ENDDO j_loop_y_flux_5
667 ! next, x - flux divergence
673 j_end = MIN(jte,jde-1)
675 ! higher order flux has a 5 or 7 point stencil, so compute
676 ! bounds so we can switch to second order flux close to the boundary
682 i_start = MAX(ids+1,its)
687 i_end = MIN(ide-1,ite)
693 DO j = j_start, j_end
695 ! 5th or 6th order flux
698 DO i = i_start_f, i_end_f
699 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
700 fqx( i,k ) = vel*flux5( u(i-3,k,j), u(i-2,k,j), &
701 u(i-1,k,j), u(i ,k,j), &
702 u(i+1,k,j), u(i+2,k,j), &
707 ! lower order fluxes close to boundaries (if not periodic or symmetric)
708 ! specified uses upstream normal wind at boundaries
710 IF( degrade_xs ) THEN
712 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
716 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
717 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
724 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
725 fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
726 u(i ,k,j), u(i+1,k,j), &
732 IF( degrade_xe ) THEN
734 IF( i_end == ide-1 ) THEN ! second order flux next to the boundary
738 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
739 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
746 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
747 fqx( i,k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
748 u(i ,k,j), u(i+1,k,j), &
754 ! x flux-divergence into tendency
757 DO i = i_start, i_end
758 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
759 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
765 ELSE IF( horz_order == 4 ) THEN
767 ! determine boundary mods for flux operators
768 ! We degrade the flux operators from 3rd/4th order
769 ! to second order one gridpoint in from the boundaries for
770 ! all boundary conditions except periodic and symmetry - these
771 ! conditions have boundary zone data fill for correct application
772 ! of the higher order flux stencils
779 IF( config_flags%periodic_x .or. &
780 config_flags%symmetric_xs .or. &
781 (its > ids+1) ) degrade_xs = .false.
782 IF( config_flags%periodic_x .or. &
783 config_flags%symmetric_xe .or. &
784 (ite < ide-1) ) degrade_xe = .false.
785 IF( config_flags%periodic_y .or. &
786 config_flags%symmetric_ys .or. &
787 (jts > jds+1) ) degrade_ys = .false.
788 IF( config_flags%periodic_y .or. &
789 config_flags%symmetric_ye .or. &
790 (jte < jde-2) ) degrade_ye = .false.
792 !--------------- x - advection first
797 j_end = MIN(jte,jde-1)
799 ! 3rd or 4th order flux has a 5 point stencil, so compute
800 ! bounds so we can switch to second order flux close to the boundary
807 i_start_f = i_start+1
817 DO j = j_start, j_end
820 DO i = i_start_f, i_end_f
821 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
822 fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
823 u(i ,k,j), u(i+1,k,j), vel )
827 ! second order flux close to boundaries (if not periodic or symmetric)
828 ! specified uses upstream normal wind at boundaries
830 IF( degrade_xs ) THEN
834 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
835 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
840 IF( degrade_xe ) THEN
844 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
845 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
850 ! x flux-divergence into tendency
853 DO i = i_start, i_end
854 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
855 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
865 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
866 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
867 IF ( config_flags%periodic_x ) i_start = its
868 IF ( config_flags%periodic_x ) i_end = ite
871 j_end = MIN(jte,jde-1)
873 ! 3rd or 4th order flux has a 5 point stencil, so compute
874 ! bounds so we can switch to second order flux close to the boundary
879 !CJM these may not work with tiling because they define j_start and end in terms of domain dim
882 j_start_f = j_start+1
890 IF(config_flags%polar) j_end = MIN(jte,jde-1)
892 ! j flux loop for v flux of u momentum
897 DO j = j_start, j_end+1
899 IF ( (j < j_start_f) .and. degrade_ys) THEN
901 DO i = i_start, i_end
902 fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) &
903 *(u(i,k,j_start)+u(i,k,j_start-1))
906 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
908 DO i = i_start, i_end
909 ! Assumes j>j_end_f is ONLY j_end+1 ...
910 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) &
911 ! *(u(i,k,j_end+1)+u(i,k,j_end))
912 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
913 *(u(i,k,j)+u(i,k,j-1))
917 ! 3rd or 4th order flux
919 DO i = i_start, i_end
920 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
921 fqy( i, k, jp1 ) = vel*flux4( u(i,k,j-2), u(i,k,j-1), &
922 u(i,k,j ), u(i,k,j+1), &
929 ! y flux-divergence into tendency
931 ! (j > j_start) will miss the u(,,jds) tendency
932 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
934 DO i = i_start, i_end
935 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
936 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
939 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
940 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
942 DO i = i_start, i_end
943 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
944 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
949 IF (j > j_start) THEN
952 DO i = i_start, i_end
953 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
954 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
968 ELSE IF ( horz_order == 3 ) THEN
970 ! As with the 5th and 6th order flux chioces, the 3rd and 4th order
971 ! code is EXACTLY the same EXCEPT for the flux operator.
973 ! determine boundary mods for flux operators
974 ! We degrade the flux operators from 3rd/4th order
975 ! to second order one gridpoint in from the boundaries for
976 ! all boundary conditions except periodic and symmetry - these
977 ! conditions have boundary zone data fill for correct application
978 ! of the higher order flux stencils
985 IF( config_flags%periodic_x .or. &
986 config_flags%symmetric_xs .or. &
987 (its > ids+1) ) degrade_xs = .false.
988 IF( config_flags%periodic_x .or. &
989 config_flags%symmetric_xe .or. &
990 (ite < ide-1) ) degrade_xe = .false.
991 IF( config_flags%periodic_y .or. &
992 config_flags%symmetric_ys .or. &
993 (jts > jds+1) ) degrade_ys = .false.
994 IF( config_flags%periodic_y .or. &
995 config_flags%symmetric_ye .or. &
996 (jte < jde-2) ) degrade_ye = .false.
998 !--------------- x - advection first
1003 j_end = MIN(jte,jde-1)
1005 ! 3rd or 4th order flux has a 5 point stencil, so compute
1006 ! bounds so we can switch to second order flux close to the boundary
1013 i_start_f = i_start+1
1023 DO j = j_start, j_end
1026 DO i = i_start_f, i_end_f
1027 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
1028 fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
1029 u(i ,k,j), u(i+1,k,j), vel )
1033 ! second order flux close to boundaries (if not periodic or symmetric)
1034 ! specified uses upstream normal wind at boundaries
1036 IF( degrade_xs ) THEN
1040 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
1041 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
1046 IF( degrade_xe ) THEN
1050 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
1051 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
1056 ! x flux-divergence into tendency
1059 DO i = i_start, i_end
1060 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1061 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
1070 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
1071 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
1072 IF ( config_flags%periodic_x ) i_start = its
1073 IF ( config_flags%periodic_x ) i_end = ite
1076 j_end = MIN(jte,jde-1)
1078 ! 3rd or 4th order flux has a 5 point stencil, so compute
1079 ! bounds so we can switch to second order flux close to the boundary
1084 !CJM these may not work with tiling because they define j_start and end in terms of domain dim
1087 j_start_f = j_start+1
1095 IF(config_flags%polar) j_end = MIN(jte,jde-1)
1097 ! j flux loop for v flux of u momentum
1102 DO j = j_start, j_end+1
1104 IF ( (j < j_start_f) .and. degrade_ys) THEN
1106 DO i = i_start, i_end
1107 fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) &
1108 *(u(i,k,j_start)+u(i,k,j_start-1))
1111 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
1113 DO i = i_start, i_end
1114 ! Assumes j>j_end_f is ONLY j_end+1 ...
1115 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) &
1116 ! *(u(i,k,j_end+1)+u(i,k,j_end))
1117 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
1118 *(u(i,k,j)+u(i,k,j-1))
1122 ! 3rd or 4th order flux
1124 DO i = i_start, i_end
1125 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
1126 fqy( i, k, jp1 ) = vel*flux3( u(i,k,j-2), u(i,k,j-1), &
1127 u(i,k,j ), u(i,k,j+1), &
1134 ! y flux-divergence into tendency
1136 ! (j > j_start) will miss the u(,,jds) tendency
1137 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
1139 DO i = i_start, i_end
1140 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1141 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
1144 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
1145 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
1147 DO i = i_start, i_end
1148 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1149 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
1154 IF (j > j_start) THEN
1157 DO i = i_start, i_end
1158 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1159 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
1173 ELSE IF ( horz_order == 2 ) THEN
1178 j_end = MIN(jte,jde-1)
1180 IF ( config_flags%open_xs ) i_start = MAX(ids+1,its)
1181 IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite)
1182 IF ( specified ) i_start = MAX(ids+2,its)
1183 IF ( specified ) i_end = MIN(ide-2,ite)
1184 IF ( config_flags%periodic_x ) i_start = its
1185 IF ( config_flags%periodic_x ) i_end = ite
1187 DO j = j_start, j_end
1189 DO i = i_start, i_end
1190 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1191 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1192 *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) &
1193 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j)))
1198 IF ( specified .AND. its .LE. ids+1 .AND. .NOT. config_flags%periodic_x ) THEN
1199 DO j = j_start, j_end
1202 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1204 IF (u(i,k,j) .LT. 0.) ub = u(i,k,j)
1205 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1206 *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) &
1207 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+ub))
1211 IF ( specified .AND. ite .GE. ide-1 .AND. .NOT. config_flags%periodic_x ) THEN
1212 DO j = j_start, j_end
1215 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1217 IF (u(i,k,j) .GT. 0.) ub = u(i,k,j)
1218 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1219 *((ru(i+1,k,j)+ru(i,k,j))*(ub+u(i,k,j)) &
1220 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j)))
1225 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
1226 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
1228 DO j = j_start, j_end
1230 DO i = i_start, i_end
1231 mrdy=msfux(i,j)*rdy ! ADT eqn 44, 1st term on RHS
1232 ! Comments for polar boundary condition
1233 ! Flow is only from one side for points next to poles
1234 IF ( (config_flags%polar) .AND. (j == jds) ) THEN
1235 tendency(i,k,j)=tendency(i,k,j)-mrdy*0.25 &
1236 *(rv(i,k,j+1)+rv(i-1,k,j+1))*(u(i,k,j+1)+u(i,k,j))
1237 ELSE IF ( (config_flags%polar) .AND. (j == jde-1) ) THEN
1238 tendency(i,k,j)=tendency(i,k,j)+mrdy*0.25 &
1239 *(rv(i,k,j)+rv(i-1,k,j))*(u(i,k,j)+u(i,k,j-1))
1241 tendency(i,k,j)=tendency(i,k,j)-mrdy*0.25 &
1242 *((rv(i,k,j+1)+rv(i-1,k,j+1))*(u(i,k,j+1)+u(i,k,j)) &
1243 -(rv(i,k,j)+rv(i-1,k,j))*(u(i,k,j)+u(i,k,j-1)))
1249 ELSE IF ( horz_order == 0 ) THEN
1251 ! Just in case we want to turn horizontal advection off, we can do it
1255 WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: h_order not known ',horz_order
1256 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
1258 ENDIF horizontal_order_test
1260 ! radiative lateral boundary condition in x for normal velocity (u)
1262 IF ( (config_flags%open_xs) .and. its == ids ) THEN
1265 j_end = MIN(jte,jde-1)
1267 DO j = j_start, j_end
1269 ub = MIN(ru(its,k,j)-cb*mut(its,j), 0.)
1270 tendency(its,k,j) = tendency(its,k,j) &
1271 - rdx*ub*(u_old(its+1,k,j) - u_old(its,k,j))
1277 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1280 j_end = MIN(jte,jde-1)
1282 DO j = j_start, j_end
1284 ub = MAX(ru(ite,k,j)+cb*mut(ite-1,j), 0.)
1285 tendency(ite,k,j) = tendency(ite,k,j) &
1286 - rdx*ub*(u_old(ite,k,j) - u_old(ite-1,k,j))
1292 ! pick up the rest of the horizontal radiation boundary conditions.
1293 ! (these are the computations that don't require 'cb')
1294 ! first, set to index ranges
1297 i_end = MIN(ite,ide)
1301 IF (config_flags%open_xs) THEN
1302 i_start = MAX(ids+1, its)
1305 IF (config_flags%open_xe) THEN
1306 i_end = MIN(ite,ide-1)
1310 IF( (config_flags%open_ys) .and. (jts == jds)) THEN
1312 DO i = i_start, i_end
1314 mrdy=msfux(i,jts)*rdy ! ADT eqn 44, 2nd term on RHS
1316 im = MAX( imin, i-1 )
1320 vw = 0.5*(rv(ip,k,jts)+rv(im,k,jts))
1322 dvm = rv(ip,k,jts+1)-rv(ip,k,jts)
1323 dvp = rv(im,k,jts+1)-rv(im,k,jts)
1324 tendency(i,k,jts)=tendency(i,k,jts)-mrdy*( &
1325 vb*(u_old(i,k,jts+1)-u_old(i,k,jts)) &
1326 +0.5*u(i,k,jts)*(dvm+dvp))
1332 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
1334 DO i = i_start, i_end
1336 mrdy=msfux(i,jte-1)*rdy ! ADT eqn 44, 2nd term on RHS
1338 im = MAX( imin, i-1 )
1342 vw = 0.5*(rv(ip,k,jte)+rv(im,k,jte))
1344 dvm = rv(ip,k,jte)-rv(ip,k,jte-1)
1345 dvp = rv(im,k,jte)-rv(im,k,jte-1)
1346 tendency(i,k,jte-1)=tendency(i,k,jte-1)-mrdy*( &
1347 vb*(u_old(i,k,jte-1)-u_old(i,k,jte-2)) &
1348 +0.5*u(i,k,jte-1)*(dvm+dvp))
1354 !-------------------- vertical advection
1355 ! ADT eqn 44 has 3rd term on RHS = -(1/my) partial d/dz (rho u w)
1356 ! Here we have: - partial d/dz (u*rom) = - partial d/dz (u rho w / my)
1357 ! Since 'my' (map scale factor in y-direction) isn't a function of z,
1358 ! this is what we need, so leave unchanged in advect_u
1363 j_end = min(jte,jde-1)
1365 ! IF ( config_flags%open_xs ) i_start = MAX(ids+1,its)
1366 ! IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite)
1368 IF ( config_flags%open_ys .or. specified ) i_start = MAX(ids+1,its)
1369 IF ( config_flags%open_ye .or. specified ) i_end = MIN(ide-1,ite)
1370 IF ( config_flags%periodic_x ) i_start = its
1371 IF ( config_flags%periodic_x ) i_end = ite
1373 DO i = i_start, i_end
1378 vert_order_test : IF (vert_order == 6) THEN
1380 DO j = j_start, j_end
1383 DO i = i_start, i_end
1384 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1385 vflux(i,k) = vel*flux6( &
1386 u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), &
1387 u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel )
1391 DO i = i_start, i_end
1394 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1395 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1397 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1398 vflux(i,k) = vel*flux4( &
1399 u(i,k-2,j), u(i,k-1,j), &
1400 u(i,k ,j), u(i,k+1,j), -vel )
1402 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1403 vflux(i,k) = vel*flux4( &
1404 u(i,k-2,j), u(i,k-1,j), &
1405 u(i,k ,j), u(i,k+1,j), -vel )
1407 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1408 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1412 DO i = i_start, i_end
1413 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1418 ELSE IF (vert_order == 5) THEN
1420 DO j = j_start, j_end
1423 DO i = i_start, i_end
1424 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1425 vflux(i,k) = vel*flux5( &
1426 u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), &
1427 u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel )
1431 DO i = i_start, i_end
1434 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1435 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1437 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1438 vflux(i,k) = vel*flux3( &
1439 u(i,k-2,j), u(i,k-1,j), &
1440 u(i,k ,j), u(i,k+1,j), -vel )
1442 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1443 vflux(i,k) = vel*flux3( &
1444 u(i,k-2,j), u(i,k-1,j), &
1445 u(i,k ,j), u(i,k+1,j), -vel )
1447 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1448 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1452 DO i = i_start, i_end
1453 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1458 ELSE IF (vert_order == 4) THEN
1460 DO j = j_start, j_end
1463 DO i = i_start, i_end
1464 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1465 vflux(i,k) = vel*flux4( &
1466 u(i,k-2,j), u(i,k-1,j), &
1467 u(i,k ,j), u(i,k+1,j), -vel )
1471 DO i = i_start, i_end
1474 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1475 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1477 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1478 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1482 DO i = i_start, i_end
1483 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1488 ELSE IF (vert_order == 3) THEN
1490 DO j = j_start, j_end
1493 DO i = i_start, i_end
1494 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1495 vflux(i,k) = vel*flux3( &
1496 u(i,k-2,j), u(i,k-1,j), &
1497 u(i,k ,j), u(i,k+1,j), -vel )
1501 DO i = i_start, i_end
1504 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1505 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1507 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1508 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1512 DO i = i_start, i_end
1513 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1518 ELSE IF (vert_order == 2) THEN
1520 DO j = j_start, j_end
1522 DO i = i_start, i_end
1523 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1524 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1530 DO i = i_start, i_end
1531 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1539 WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: v_order not known ',vert_order
1540 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
1542 ENDIF vert_order_test
1544 END SUBROUTINE advect_u
1546 !-------------------------------------------------------------------------------
1548 SUBROUTINE advect_v ( v, v_old, tendency, &
1550 mut, time_step, config_flags, &
1551 msfux, msfuy, msfvx, msfvy, &
1555 ids, ide, jds, jde, kds, kde, &
1556 ims, ime, jms, jme, kms, kme, &
1557 its, ite, jts, jte, kts, kte )
1563 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1565 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1566 ims, ime, jms, jme, kms, kme, &
1567 its, ite, jts, jte, kts, kte
1569 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: v, &
1575 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
1576 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
1578 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
1585 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
1589 REAL , INTENT(IN ) :: rdx, &
1591 INTEGER , INTENT(IN ) :: time_step
1596 INTEGER :: i, j, k, itf, jtf, ktf
1597 INTEGER :: i_start, i_end, j_start, j_end
1598 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
1599 INTEGER :: jmin, jmax, jp, jm, imin, imax
1601 REAL :: mrdx, mrdy, ub, vb, uw, vw, dup, dum
1602 REAL , DIMENSION(its:ite, kts:kte) :: vflux
1605 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
1606 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
1608 INTEGER :: horz_order
1609 INTEGER :: vert_order
1611 LOGICAL :: degrade_xs, degrade_ys
1612 LOGICAL :: degrade_xe, degrade_ye
1614 INTEGER :: jp1, jp0, jtmp
1617 ! definition of flux operators, 3rd, 4th, 5th or 6th order
1619 REAL :: flux3, flux4, flux5, flux6
1620 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
1622 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
1623 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
1625 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
1626 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
1627 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
1629 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
1630 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
1631 +(q_ip2+q_im3) )/60.0
1633 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
1634 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
1635 -sign(1,time_step)*sign(1.,ua)*( &
1636 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
1640 LOGICAL :: specified
1643 if(config_flags%specified .or. config_flags%nested) specified = .true.
1645 ! set order for the advection schemes
1648 horz_order = config_flags%h_mom_adv_order
1649 vert_order = config_flags%v_mom_adv_order
1652 ! here is the choice of flux operators
1655 horizontal_order_test : IF( horz_order == 6 ) THEN
1657 ! determine boundary mods for flux operators
1658 ! We degrade the flux operators from 3rd/4th order
1659 ! to second order one gridpoint in from the boundaries for
1660 ! all boundary conditions except periodic and symmetry - these
1661 ! conditions have boundary zone data fill for correct application
1662 ! of the higher order flux stencils
1669 IF( config_flags%periodic_x .or. &
1670 config_flags%symmetric_xs .or. &
1671 (its > ids+2) ) degrade_xs = .false.
1672 IF( config_flags%periodic_x .or. &
1673 config_flags%symmetric_xe .or. &
1674 (ite < ide-3) ) degrade_xe = .false.
1675 IF( config_flags%periodic_y .or. &
1676 config_flags%symmetric_ys .or. &
1677 (jts > jds+2) ) degrade_ys = .false.
1678 IF( config_flags%periodic_y .or. &
1679 config_flags%symmetric_ye .or. &
1680 (jte < jde-2) ) degrade_ye = .false.
1682 !--------------- y - advection first
1687 i_end = MIN(ite,ide-1)
1691 ! higher order flux has a 5 or 7 point stencil, so compute
1692 ! bounds so we can switch to second order flux close to the boundary
1698 j_start = MAX(jts,jds+1)
1703 j_end = MIN(jte,jde-1)
1707 ! compute fluxes, 5th or 6th order
1712 j_loop_y_flux_6 : DO j = j_start, j_end+1
1714 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
1717 DO i = i_start, i_end
1718 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1719 fqy( i, k, jp1 ) = vel*flux6( &
1720 v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), &
1721 v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel )
1725 ! we must be close to some boundary where we need to reduce the order of the stencil
1726 ! specified uses upstream normal wind at boundaries
1728 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
1731 DO i = i_start, i_end
1733 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
1734 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1739 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
1742 DO i = i_start, i_end
1743 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1744 fqy( i, k, jp1 ) = vel*flux4( &
1745 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1750 ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary
1753 DO i = i_start, i_end
1755 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
1756 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1761 ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary
1764 DO i = i_start, i_end
1765 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1766 fqy( i, k, jp1 ) = vel*flux4( &
1767 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1773 ! y flux-divergence into tendency
1775 ! Comments on polar boundary conditions
1776 ! No advection over the poles means tendencies (held from jds [S. pole]
1777 ! to jde [N pole], i.e., on v grid) must be zero at poles
1778 ! [tendency(jds) and tendency(jde)=0]
1779 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
1781 DO i = i_start, i_end
1782 tendency(i,k,j-1) = 0.
1785 ! If j_end were set to jde in a special if statement apart from
1786 ! degrade_ye, then we would hit the next conditional. But since
1787 ! we want the tendency to be zero anyway, not looping to jde+1
1788 ! will produce the same effect.
1789 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
1791 DO i = i_start, i_end
1792 tendency(i,k,j-1) = 0.
1797 IF(j > j_start) THEN
1800 DO i = i_start, i_end
1801 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
1802 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
1814 ENDDO j_loop_y_flux_6
1816 ! next, x - flux divergence
1819 i_end = MIN(ite,ide-1)
1823 ! Polar boundary conditions are like open or specified
1824 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
1825 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
1827 ! higher order flux has a 5 or 7 point stencil, so compute
1828 ! bounds so we can switch to second order flux close to the boundary
1834 i_start = MAX(ids+1,its)
1835 i_start_f = i_start+2
1839 i_end = MIN(ide-2,ite)
1845 DO j = j_start, j_end
1847 ! 5th or 6th order flux
1850 DO i = i_start_f, i_end_f
1851 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1852 fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j), &
1853 v(i-1,k,j), v(i ,k,j), &
1854 v(i+1,k,j), v(i+2,k,j), &
1859 ! lower order fluxes close to boundaries (if not periodic or symmetric)
1861 IF( degrade_xs ) THEN
1863 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
1866 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
1867 *(v(i,k,j)+v(i-1,k,j))
1873 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1874 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
1875 v(i ,k,j), v(i+1,k,j), &
1881 IF( degrade_xe ) THEN
1883 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
1886 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
1887 *(v(i_end+1,k,j)+v(i_end,k,j))
1893 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1894 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
1895 v(i ,k,j), v(i+1,k,j), &
1901 ! x flux-divergence into tendency
1904 DO i = i_start, i_end
1905 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
1906 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
1912 ELSE IF( horz_order == 5 ) THEN
1914 ! 5th order horizontal flux calculation
1915 ! This code is EXACTLY the same as the 6th order code
1916 ! EXCEPT the 5th order and 3rd operators are used in
1917 ! place of the 6th and 4th order operators
1919 ! determine boundary mods for flux operators
1920 ! We degrade the flux operators from 3rd/4th order
1921 ! to second order one gridpoint in from the boundaries for
1922 ! all boundary conditions except periodic and symmetry - these
1923 ! conditions have boundary zone data fill for correct application
1924 ! of the higher order flux stencils
1931 IF( config_flags%periodic_x .or. &
1932 config_flags%symmetric_xs .or. &
1933 (its > ids+2) ) degrade_xs = .false.
1934 IF( config_flags%periodic_x .or. &
1935 config_flags%symmetric_xe .or. &
1936 (ite < ide-3) ) degrade_xe = .false.
1937 IF( config_flags%periodic_y .or. &
1938 config_flags%symmetric_ys .or. &
1939 (jts > jds+2) ) degrade_ys = .false.
1940 IF( config_flags%periodic_y .or. &
1941 config_flags%symmetric_ye .or. &
1942 (jte < jde-2) ) degrade_ye = .false.
1944 !--------------- y - advection first
1947 i_end = MIN(ite,ide-1)
1951 ! higher order flux has a 5 or 7 point stencil, so compute
1952 ! bounds so we can switch to second order flux close to the boundary
1958 j_start = MAX(jts,jds+1)
1963 j_end = MIN(jte,jde-1)
1967 ! compute fluxes, 5th or 6th order
1972 j_loop_y_flux_5 : DO j = j_start, j_end+1
1974 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
1977 DO i = i_start, i_end
1978 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1979 fqy( i, k, jp1 ) = vel*flux5( &
1980 v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), &
1981 v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel )
1985 ! we must be close to some boundary where we need to reduce the order of the stencil
1986 ! specified uses upstream normal wind at boundaries
1988 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
1991 DO i = i_start, i_end
1993 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
1994 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1999 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
2002 DO i = i_start, i_end
2003 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2004 fqy( i, k, jp1 ) = vel*flux3( &
2005 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
2010 ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary
2013 DO i = i_start, i_end
2015 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
2016 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2021 ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary
2024 DO i = i_start, i_end
2025 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2026 fqy( i, k, jp1 ) = vel*flux3( &
2027 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
2033 ! y flux-divergence into tendency
2035 ! Comments on polar boundary conditions
2036 ! No advection over the poles means tendencies (held from jds [S. pole]
2037 ! to jde [N pole], i.e., on v grid) must be zero at poles
2038 ! [tendency(jds) and tendency(jde)=0]
2039 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
2041 DO i = i_start, i_end
2042 tendency(i,k,j-1) = 0.
2045 ! If j_end were set to jde in a special if statement apart from
2046 ! degrade_ye, then we would hit the next conditional. But since
2047 ! we want the tendency to be zero anyway, not looping to jde+1
2048 ! will produce the same effect.
2049 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
2051 DO i = i_start, i_end
2052 tendency(i,k,j-1) = 0.
2057 IF(j > j_start) THEN
2060 DO i = i_start, i_end
2061 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
2062 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
2074 ENDDO j_loop_y_flux_5
2076 ! next, x - flux divergence
2079 i_end = MIN(ite,ide-1)
2083 ! Polar boundary conditions are like open or specified
2084 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2085 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2087 ! higher order flux has a 5 or 7 point stencil, so compute
2088 ! bounds so we can switch to second order flux close to the boundary
2094 i_start = MAX(ids+1,its)
2095 i_start_f = i_start+2
2099 i_end = MIN(ide-2,ite)
2105 DO j = j_start, j_end
2107 ! 5th or 6th order flux
2110 DO i = i_start_f, i_end_f
2111 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2112 fqx( i, k ) = vel*flux5( v(i-3,k,j), v(i-2,k,j), &
2113 v(i-1,k,j), v(i ,k,j), &
2114 v(i+1,k,j), v(i+2,k,j), &
2119 ! lower order fluxes close to boundaries (if not periodic or symmetric)
2121 IF( degrade_xs ) THEN
2123 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
2126 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
2127 *(v(i,k,j)+v(i-1,k,j))
2133 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2134 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2135 v(i ,k,j), v(i+1,k,j), &
2141 IF( degrade_xe ) THEN
2143 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
2146 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2147 *(v(i_end+1,k,j)+v(i_end,k,j))
2153 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2154 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2155 v(i ,k,j), v(i+1,k,j), &
2161 ! x flux-divergence into tendency
2164 DO i = i_start, i_end
2165 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2166 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2172 ELSE IF( horz_order == 4 ) THEN
2174 ! determine boundary mods for flux operators
2175 ! We degrade the flux operators from 3rd/4th order
2176 ! to second order one gridpoint in from the boundaries for
2177 ! all boundary conditions except periodic and symmetry - these
2178 ! conditions have boundary zone data fill for correct application
2179 ! of the higher order flux stencils
2186 IF( config_flags%periodic_x .or. &
2187 config_flags%symmetric_xs .or. &
2188 (its > ids+1) ) degrade_xs = .false.
2189 IF( config_flags%periodic_x .or. &
2190 config_flags%symmetric_xe .or. &
2191 (ite < ide-2) ) degrade_xe = .false.
2192 IF( config_flags%periodic_y .or. &
2193 config_flags%symmetric_ys .or. &
2194 (jts > jds+1) ) degrade_ys = .false.
2195 IF( config_flags%periodic_y .or. &
2196 config_flags%symmetric_ye .or. &
2197 (jte < jde-1) ) degrade_ye = .false.
2199 !--------------- y - advection first
2205 i_end = MIN(ite,ide-1)
2209 ! 3rd or 4th order flux has a 5 point stencil, so compute
2210 ! bounds so we can switch to second order flux close to the boundary
2215 !CJM May not work with tiling because defined in terms of domain dims
2218 j_start_f = j_start+1
2227 ! specified uses upstream normal wind at boundaries
2232 DO j = j_start, j_end+1
2234 IF ((j == j_start) .and. degrade_ys) THEN
2236 DO i = i_start, i_end
2238 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
2239 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2243 ELSE IF ((j == j_end+1) .and. degrade_ye) THEN
2245 DO i = i_start, i_end
2247 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
2248 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2254 DO i = i_start, i_end
2255 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2256 fqy( i,k,jp1 ) = vel*flux4( v(i,k,j-2), v(i,k,j-1), &
2257 v(i,k,j ), v(i,k,j+1), &
2263 ! Comments on polar boundary conditions
2264 ! No advection over the poles means tendencies (held from jds [S. pole]
2265 ! to jde [N pole], i.e., on v grid) must be zero at poles
2266 ! [tendency(jds) and tendency(jde)=0]
2267 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
2269 DO i = i_start, i_end
2270 tendency(i,k,j-1) = 0.
2273 ! If j_end were set to jde in a special if statement apart from
2274 ! degrade_ye, then we would hit the next conditional. But since
2275 ! we want the tendency to be zero anyway, not looping to jde+1
2276 ! will produce the same effect.
2277 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
2279 DO i = i_start, i_end
2280 tendency(i,k,j-1) = 0.
2285 IF( j > j_start) THEN
2287 DO i = i_start, i_end
2288 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
2289 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
2303 ! next, x - flux divergence
2307 i_end = MIN(ite,ide-1)
2311 ! Polar boundary conditions are like open or specified
2312 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2313 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2315 ! 3rd or 4th order flux has a 5 point stencil, so compute
2316 ! bounds so we can switch to second order flux close to the boundary
2323 i_start_f = i_start+1
2333 DO j = j_start, j_end
2335 ! 3rd or 4th order flux
2338 DO i = i_start_f, i_end_f
2339 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2340 fqx( i, k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
2341 v(i ,k,j), v(i+1,k,j), &
2346 ! second order flux close to boundaries (if not periodic or symmetric)
2348 IF( degrade_xs ) THEN
2350 fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) &
2351 *(v(i_start,k,j)+v(i_start-1,k,j))
2355 IF( degrade_xe ) THEN
2357 fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2358 *(v(i_end+1,k,j)+v(i_end,k,j))
2362 ! x flux-divergence into tendency
2365 DO i = i_start, i_end
2366 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2367 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2373 ELSE IF( horz_order == 3 ) THEN
2375 ! determine boundary mods for flux operators
2376 ! We degrade the flux operators from 3rd/4th order
2377 ! to second order one gridpoint in from the boundaries for
2378 ! all boundary conditions except periodic and symmetry - these
2379 ! conditions have boundary zone data fill for correct application
2380 ! of the higher order flux stencils
2387 IF( config_flags%periodic_x .or. &
2388 config_flags%symmetric_xs .or. &
2389 (its > ids+1) ) degrade_xs = .false.
2390 IF( config_flags%periodic_x .or. &
2391 config_flags%symmetric_xe .or. &
2392 (ite < ide-2) ) degrade_xe = .false.
2393 IF( config_flags%periodic_y .or. &
2394 config_flags%symmetric_ys .or. &
2395 (jts > jds+1) ) degrade_ys = .false.
2396 IF( config_flags%periodic_y .or. &
2397 config_flags%symmetric_ye .or. &
2398 (jte < jde-1) ) degrade_ye = .false.
2400 !--------------- y - advection first
2406 i_end = MIN(ite,ide-1)
2410 ! 3rd or 4th order flux has a 5 point stencil, so compute
2411 ! bounds so we can switch to second order flux close to the boundary
2416 !CJM May not work with tiling because defined in terms of domain dims
2419 j_start_f = j_start+1
2428 ! specified uses upstream normal wind at boundaries
2433 DO j = j_start, j_end+1
2435 IF ((j == j_start) .and. degrade_ys) THEN
2437 DO i = i_start, i_end
2439 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
2440 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2444 ELSE IF ((j == j_end+1) .and. degrade_ye) THEN
2446 DO i = i_start, i_end
2448 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
2449 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2455 DO i = i_start, i_end
2456 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2457 fqy( i,k,jp1 ) = vel*flux3( v(i,k,j-2), v(i,k,j-1), &
2458 v(i,k,j ), v(i,k,j+1), &
2464 ! Comments on polar boundary conditions
2465 ! No advection over the poles means tendencies (held from jds [S. pole]
2466 ! to jde [N pole], i.e., on v grid) must be zero at poles
2467 ! [tendency(jds) and tendency(jde)=0]
2468 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
2470 DO i = i_start, i_end
2471 tendency(i,k,j-1) = 0.
2474 ! If j_end were set to jde in a special if statement apart from
2475 ! degrade_ye, then we would hit the next conditional. But since
2476 ! we want the tendency to be zero anyway, not looping to jde+1
2477 ! will produce the same effect.
2478 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
2480 DO i = i_start, i_end
2481 tendency(i,k,j-1) = 0.
2486 IF( j > j_start) THEN
2488 DO i = i_start, i_end
2489 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
2490 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
2504 ! next, x - flux divergence
2508 i_end = MIN(ite,ide-1)
2512 ! Polar boundary conditions are like open or specified
2513 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2514 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2516 ! 3rd or 4th order flux has a 5 point stencil, so compute
2517 ! bounds so we can switch to second order flux close to the boundary
2524 i_start_f = i_start+1
2534 DO j = j_start, j_end
2536 ! 3rd or 4th order flux
2539 DO i = i_start_f, i_end_f
2540 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2541 fqx( i, k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2542 v(i ,k,j), v(i+1,k,j), &
2547 ! second order flux close to boundaries (if not periodic or symmetric)
2549 IF( degrade_xs ) THEN
2551 fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) &
2552 *(v(i_start,k,j)+v(i_start-1,k,j))
2556 IF( degrade_xe ) THEN
2558 fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2559 *(v(i_end+1,k,j)+v(i_end,k,j))
2563 ! x flux-divergence into tendency
2566 DO i = i_start, i_end
2567 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2568 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2574 ELSE IF( horz_order == 2 ) THEN
2578 i_end = MIN(ite,ide-1)
2582 IF ( config_flags%open_ys ) j_start = MAX(jds+1,jts)
2583 IF ( config_flags%open_ye ) j_end = MIN(jde-1,jte)
2584 IF ( specified ) j_start = MAX(jds+2,jts)
2585 IF ( specified ) j_end = MIN(jde-2,jte)
2586 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2587 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2589 DO j = j_start, j_end
2591 DO i = i_start, i_end
2593 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2595 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2596 *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) &
2597 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1)))
2603 ! Comments on polar boundary conditions
2604 ! tendencies = 0 at poles, and polar points do not contribute at points
2606 IF (config_flags%polar) THEN
2607 IF (jts == jds) THEN
2609 DO i = i_start, i_end
2610 tendency(i,k,jds) = 0.
2614 IF (jte == jde) THEN
2616 DO i = i_start, i_end
2617 tendency(i,k,jde) = 0.
2623 ! specified uses upstream normal wind at boundaries
2625 IF ( specified .AND. jts .LE. jds+1 ) THEN
2628 DO i = i_start, i_end
2629 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2631 IF (v(i,k,j) .LT. 0.) vb = v(i,k,j)
2633 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2634 *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) &
2635 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+vb))
2641 IF ( specified .AND. jte .GE. jde-1 ) THEN
2644 DO i = i_start, i_end
2646 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2648 IF (v(i,k,j) .GT. 0.) vb = v(i,k,j)
2650 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2651 *((rv(i,k,j+1)+rv(i,k,j ))*(vb+v(i,k,j )) &
2652 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1)))
2658 IF ( .NOT. config_flags%periodic_x ) THEN
2659 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2660 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2662 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2663 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2665 DO j = j_start, j_end
2667 DO i = i_start, i_end
2669 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2671 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
2672 *((ru(i+1,k,j)+ru(i+1,k,j-1))*(v(i+1,k,j)+v(i ,k,j)) &
2673 -(ru(i ,k,j)+ru(i ,k,j-1))*(v(i ,k,j)+v(i-1,k,j)))
2679 ELSE IF ( horz_order == 0 ) THEN
2681 ! Just in case we want to turn horizontal advection off, we can do it
2686 WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: h_order not known ',horz_order
2687 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
2689 ENDIF horizontal_order_test
2691 ! Comments on polar boundary condition
2692 ! Force tendency=0 at NP and SP
2693 ! We keep setting this everywhere, but it can't hurt...
2694 IF ( config_flags%polar .AND. (jts == jds) ) THEN
2697 tendency(i,k,jts)=0.
2701 IF ( config_flags%polar .AND. (jte == jde) ) THEN
2704 tendency(i,k,jte)=0.
2709 ! radiative lateral boundary condition in y for normal velocity (v)
2711 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2714 i_end = MIN(ite,ide-1)
2716 DO i = i_start, i_end
2718 vb = MIN(rv(i,k,jts)-cb*mut(i,jts), 0.)
2719 tendency(i,k,jts) = tendency(i,k,jts) &
2720 - rdy*vb*(v_old(i,k,jts+1) - v_old(i,k,jts))
2726 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2729 i_end = MIN(ite,ide-1)
2731 DO i = i_start, i_end
2733 vb = MAX(rv(i,k,jte)+cb*mut(i,jte-1), 0.)
2734 tendency(i,k,jte) = tendency(i,k,jte) &
2735 - rdy*vb*(v_old(i,k,jte) - v_old(i,k,jte-1))
2741 ! pick up the rest of the horizontal radiation boundary conditions.
2742 ! (these are the computations that don't require 'cb'.
2743 ! first, set to index ranges
2746 j_end = MIN(jte,jde)
2751 IF (config_flags%open_ys) THEN
2752 j_start = MAX(jds+1, jts)
2755 IF (config_flags%open_ye) THEN
2756 j_end = MIN(jte,jde-1)
2760 ! compute x (u) conditions for v, w, or scalar
2762 IF( (config_flags%open_xs) .and. (its == ids)) THEN
2764 DO j = j_start, j_end
2766 mrdx=msfvy(its,j)*rdx ! ADT eqn 45, 1st term on RHS
2768 jm = MAX( jmin, j-1 )
2772 uw = 0.5*(ru(its,k,jp)+ru(its,k,jm))
2774 dup = ru(its+1,k,jp)-ru(its,k,jp)
2775 dum = ru(its+1,k,jm)-ru(its,k,jm)
2776 tendency(its,k,j)=tendency(its,k,j)-mrdx*( &
2777 ub*(v_old(its+1,k,j)-v_old(its,k,j)) &
2778 +0.5*v(its,k,j)*(dup+dum))
2784 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
2785 DO j = j_start, j_end
2787 mrdx=msfvy(ite-1,j)*rdx ! ADT eqn 45, 1st term on RHS
2789 jm = MAX( jmin, j-1 )
2793 uw = 0.5*(ru(ite,k,jp)+ru(ite,k,jm))
2795 dup = ru(ite,k,jp)-ru(ite-1,k,jp)
2796 dum = ru(ite,k,jm)-ru(ite-1,k,jm)
2798 ! tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( &
2799 ! ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) &
2800 ! +0.5*v(ite-1,k,j)* &
2801 ! ( ru(ite,k,jp)-ru(ite-1,k,jp) &
2802 ! +ru(ite,k,jm)-ru(ite-1,k,jm)) )
2803 tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( &
2804 ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) &
2805 +0.5*v(ite-1,k,j)*(dup+dum))
2812 !-------------------- vertical advection
2813 ! ADT eqn 45 has 3rd term on RHS = -(1/mx) partial d/dz (rho v w)
2814 ! Here we have: - partial d/dz (v*rom) = - partial d/dz (v rho w / my)
2815 ! We therefore need to make a correction for advect_v
2816 ! since 'my' (map scale factor in y direction) isn't a function of z,
2817 ! we can do this using *(my/mx) (see eqn. 45 for example)
2821 i_end = MIN(ite,ide-1)
2825 DO i = i_start, i_end
2830 ! Polar boundary conditions are like open or specified
2831 ! We don't want to calculate vertical v tendencies at the N or S pole
2832 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2833 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2835 vert_order_test : IF (vert_order == 6) THEN
2837 DO j = j_start, j_end
2841 DO i = i_start, i_end
2842 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2843 vflux(i,k) = vel*flux6( &
2844 v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), &
2845 v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel )
2849 DO i = i_start, i_end
2851 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2852 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2854 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2855 vflux(i,k) = vel*flux4( &
2856 v(i,k-2,j), v(i,k-1,j), &
2857 v(i,k ,j), v(i,k+1,j), -vel )
2859 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2860 vflux(i,k) = vel*flux4( &
2861 v(i,k-2,j), v(i,k-1,j), &
2862 v(i,k ,j), v(i,k+1,j), -vel )
2864 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2865 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2871 DO i = i_start, i_end
2872 ! We are calculating vertical fluxes on v points,
2873 ! so we must mean msf_v_x/y variables
2874 tendency(i,k,j)=tendency(i,k,j)-(msfvy(i,j)/msfvx(i,j))*rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ! ADT eqn 45, 3rd term on RHS
2880 ELSE IF (vert_order == 5) THEN
2882 DO j = j_start, j_end
2886 DO i = i_start, i_end
2887 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2888 vflux(i,k) = vel*flux5( &
2889 v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), &
2890 v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel )
2894 DO i = i_start, i_end
2896 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2897 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2899 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2900 vflux(i,k) = vel*flux3( &
2901 v(i,k-2,j), v(i,k-1,j), &
2902 v(i,k ,j), v(i,k+1,j), -vel )
2904 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2905 vflux(i,k) = vel*flux3( &
2906 v(i,k-2,j), v(i,k-1,j), &
2907 v(i,k ,j), v(i,k+1,j), -vel )
2909 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2910 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2916 DO i = i_start, i_end
2917 ! We are calculating vertical fluxes on v points,
2918 ! so we must mean msf_v_x/y variables
2919 tendency(i,k,j)=tendency(i,k,j)-(msfvy(i,j)/msfvx(i,j))*rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ! ADT eqn 45, 3rd term on RHS
2925 ELSE IF (vert_order == 4) THEN
2927 DO j = j_start, j_end
2931 DO i = i_start, i_end
2932 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2933 vflux(i,k) = vel*flux4( &
2934 v(i,k-2,j), v(i,k-1,j), &
2935 v(i,k ,j), v(i,k+1,j), -vel )
2939 DO i = i_start, i_end
2941 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2942 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2944 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2945 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2951 DO i = i_start, i_end
2952 ! We are calculating vertical fluxes on v points,
2953 ! so we must mean msf_v_x/y variables
2954 tendency(i,k,j)=tendency(i,k,j)-(msfvy(i,j)/msfvx(i,j))*rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ! ADT eqn 45, 3rd term on RHS
2960 ELSE IF (vert_order == 3) THEN
2962 DO j = j_start, j_end
2966 DO i = i_start, i_end
2967 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2968 vflux(i,k) = vel*flux3( &
2969 v(i,k-2,j), v(i,k-1,j), &
2970 v(i,k ,j), v(i,k+1,j), -vel )
2974 DO i = i_start, i_end
2976 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2977 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2979 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2980 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2986 DO i = i_start, i_end
2987 ! We are calculating vertical fluxes on v points,
2988 ! so we must mean msf_v_x/y variables
2989 tendency(i,k,j)=tendency(i,k,j)-(msfvy(i,j)/msfvx(i,j))*rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ! ADT eqn 45, 3rd term on RHS
2996 ELSE IF (vert_order == 2) THEN
2998 DO j = j_start, j_end
3000 DO i = i_start, i_end
3002 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
3003 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
3008 DO i = i_start, i_end
3009 ! We are calculating vertical fluxes on v points,
3010 ! so we must mean msf_v_x/y variables
3011 tendency(i,k,j)=tendency(i,k,j)-(msfvy(i,j)/msfvx(i,j))*rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ! ADT eqn 45, 3rd term on RHS
3018 WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: v_order not known ',vert_order
3019 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
3021 ENDIF vert_order_test
3023 END SUBROUTINE advect_v
3025 !-------------------------------------------------------------------
3027 SUBROUTINE advect_scalar ( field, field_old, tendency, &
3029 mut, time_step, config_flags, &
3030 msfux, msfuy, msfvx, msfvy, &
3034 ids, ide, jds, jde, kds, kde, &
3035 ims, ime, jms, jme, kms, kme, &
3036 its, ite, jts, jte, kts, kte )
3042 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3044 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3045 ims, ime, jms, jme, kms, kme, &
3046 its, ite, jts, jte, kts, kte
3048 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
3054 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
3055 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3057 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3064 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3068 REAL , INTENT(IN ) :: rdx, &
3070 INTEGER , INTENT(IN ) :: time_step
3075 INTEGER :: i, j, k, itf, jtf, ktf
3076 INTEGER :: i_start, i_end, j_start, j_end
3077 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
3078 INTEGER :: jmin, jmax, jp, jm, imin, imax
3080 REAL :: mrdx, mrdy, ub, vb, uw, vw
3081 REAL , DIMENSION(its:ite, kts:kte) :: vflux
3084 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
3085 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
3087 INTEGER :: horz_order, vert_order
3089 LOGICAL :: degrade_xs, degrade_ys
3090 LOGICAL :: degrade_xe, degrade_ye
3092 INTEGER :: jp1, jp0, jtmp
3095 ! definition of flux operators, 3rd, 4th, 5th or 6th order
3097 REAL :: flux3, flux4, flux5, flux6
3098 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
3100 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
3101 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
3103 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
3104 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
3105 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
3107 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
3108 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
3109 +(q_ip2+q_im3) )/60.0
3111 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
3112 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
3113 -sign(1,time_step)*sign(1.,ua)*( &
3114 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
3117 LOGICAL :: specified
3120 if(config_flags%specified .or. config_flags%nested) specified = .true.
3122 ! set order for the advection schemes
3125 horz_order = config_flags%h_sca_adv_order
3126 vert_order = config_flags%v_sca_adv_order
3128 ! begin with horizontal flux divergence
3129 ! here is the choice of flux operators
3132 horizontal_order_test : IF( horz_order == 6 ) THEN
3134 ! determine boundary mods for flux operators
3135 ! We degrade the flux operators from 3rd/4th order
3136 ! to second order one gridpoint in from the boundaries for
3137 ! all boundary conditions except periodic and symmetry - these
3138 ! conditions have boundary zone data fill for correct application
3139 ! of the higher order flux stencils
3146 IF( config_flags%periodic_x .or. &
3147 config_flags%symmetric_xs .or. &
3148 (its > ids+2) ) degrade_xs = .false.
3149 IF( config_flags%periodic_x .or. &
3150 config_flags%symmetric_xe .or. &
3151 (ite < ide-3) ) degrade_xe = .false.
3152 IF( config_flags%periodic_y .or. &
3153 config_flags%symmetric_ys .or. &
3154 (jts > jds+2) ) degrade_ys = .false.
3155 IF( config_flags%periodic_y .or. &
3156 config_flags%symmetric_ye .or. &
3157 (jte < jde-3) ) degrade_ye = .false.
3159 !--------------- y - advection first
3163 i_end = MIN(ite,ide-1)
3165 j_end = MIN(jte,jde-1)
3167 ! higher order flux has a 5 or 7 point stencil, so compute
3168 ! bounds so we can switch to second order flux close to the boundary
3174 j_start = MAX(jts,jds+1)
3179 j_end = MIN(jte,jde-2)
3183 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3185 ! compute fluxes, 5th or 6th order
3190 j_loop_y_flux_6 : DO j = j_start, j_end+1
3192 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
3195 DO i = i_start, i_end
3197 fqy( i, k, jp1 ) = vel*flux6( &
3198 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
3199 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
3203 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
3206 DO i = i_start, i_end
3207 fqy(i,k, jp1) = 0.5*rv(i,k,j)* &
3208 (field(i,k,j)+field(i,k,j-1))
3213 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
3216 DO i = i_start, i_end
3218 fqy( i, k, jp1 ) = vel*flux4( &
3219 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
3223 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
3226 DO i = i_start, i_end
3227 fqy(i, k, jp1) = 0.5*rv(i,k,j)* &
3228 (field(i,k,j)+field(i,k,j-1))
3232 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
3235 DO i = i_start, i_end
3237 fqy( i, k, jp1) = vel*flux4( &
3238 field(i,k,j-2),field(i,k,j-1), &
3239 field(i,k,j),field(i,k,j+1),vel )
3245 ! y flux-divergence into tendency
3247 ! Comments on polar boundary conditions
3248 ! Same process as for advect_u - tendencies run from jds to jde-1
3249 ! (latitudes are as for u grid, longitudes are displaced)
3250 ! Therefore: flow is only from one side for points next to poles
3251 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3253 DO i = i_start, i_end
3254 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3255 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3258 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3260 DO i = i_start, i_end
3261 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3262 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3267 IF(j > j_start) THEN
3270 DO i = i_start, i_end
3271 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3272 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3284 ENDDO j_loop_y_flux_6
3286 ! next, x - flux divergence
3289 i_end = MIN(ite,ide-1)
3292 j_end = MIN(jte,jde-1)
3294 ! higher order flux has a 5 or 7 point stencil, so compute
3295 ! bounds so we can switch to second order flux close to the boundary
3301 i_start = MAX(ids+1,its)
3302 i_start_f = i_start+2
3306 i_end = MIN(ide-2,ite)
3312 DO j = j_start, j_end
3314 ! 5th or 6th order flux
3317 DO i = i_start_f, i_end_f
3319 fqx( i,k ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), &
3320 field(i-1,k,j), field(i ,k,j), &
3321 field(i+1,k,j), field(i+2,k,j), &
3326 ! lower order fluxes close to boundaries (if not periodic or symmetric)
3328 IF( degrade_xs ) THEN
3330 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
3333 fqx(i,k) = 0.5*(ru(i,k,j)) &
3334 *(field(i,k,j)+field(i-1,k,j))
3342 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
3343 field(i ,k,j), field(i+1,k,j), &
3349 IF( degrade_xe ) THEN
3351 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
3354 fqx(i,k) = 0.5*(ru(i,k,j)) &
3355 *(field(i,k,j)+field(i-1,k,j))
3362 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
3363 field(i ,k,j), field(i+1,k,j), &
3369 ! x flux-divergence into tendency
3372 DO i = i_start, i_end
3373 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3374 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3380 ELSE IF( horz_order == 5 ) THEN
3382 ! determine boundary mods for flux operators
3383 ! We degrade the flux operators from 3rd/4th order
3384 ! to second order one gridpoint in from the boundaries for
3385 ! all boundary conditions except periodic and symmetry - these
3386 ! conditions have boundary zone data fill for correct application
3387 ! of the higher order flux stencils
3394 IF( config_flags%periodic_x .or. &
3395 config_flags%symmetric_xs .or. &
3396 (its > ids+2) ) degrade_xs = .false.
3397 IF( config_flags%periodic_x .or. &
3398 config_flags%symmetric_xe .or. &
3399 (ite < ide-3) ) degrade_xe = .false.
3400 IF( config_flags%periodic_y .or. &
3401 config_flags%symmetric_ys .or. &
3402 (jts > jds+2) ) degrade_ys = .false.
3403 IF( config_flags%periodic_y .or. &
3404 config_flags%symmetric_ye .or. &
3405 (jte < jde-3) ) degrade_ye = .false.
3407 !--------------- y - advection first
3411 i_end = MIN(ite,ide-1)
3413 j_end = MIN(jte,jde-1)
3415 ! higher order flux has a 5 or 7 point stencil, so compute
3416 ! bounds so we can switch to second order flux close to the boundary
3422 j_start = MAX(jts,jds+1)
3427 j_end = MIN(jte,jde-2)
3431 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3433 ! compute fluxes, 5th or 6th order
3438 j_loop_y_flux_5 : DO j = j_start, j_end+1
3440 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
3443 DO i = i_start, i_end
3445 fqy( i, k, jp1 ) = vel*flux5( &
3446 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
3447 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
3451 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
3454 DO i = i_start, i_end
3455 fqy(i,k, jp1) = 0.5*rv(i,k,j)* &
3456 (field(i,k,j)+field(i,k,j-1))
3461 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
3464 DO i = i_start, i_end
3466 fqy( i, k, jp1 ) = vel*flux3( &
3467 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
3471 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
3474 DO i = i_start, i_end
3475 fqy(i, k, jp1) = 0.5*rv(i,k,j)* &
3476 (field(i,k,j)+field(i,k,j-1))
3480 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
3483 DO i = i_start, i_end
3485 fqy( i, k, jp1) = vel*flux3( &
3486 field(i,k,j-2),field(i,k,j-1), &
3487 field(i,k,j),field(i,k,j+1),vel )
3493 ! y flux-divergence into tendency
3495 ! Comments on polar boundary conditions
3496 ! Same process as for advect_u - tendencies run from jds to jde-1
3497 ! (latitudes are as for u grid, longitudes are displaced)
3498 ! Therefore: flow is only from one side for points next to poles
3499 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3501 DO i = i_start, i_end
3502 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3503 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3506 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3508 DO i = i_start, i_end
3509 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3510 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3515 IF(j > j_start) THEN
3518 DO i = i_start, i_end
3519 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3520 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3532 ENDDO j_loop_y_flux_5
3534 ! next, x - flux divergence
3537 i_end = MIN(ite,ide-1)
3540 j_end = MIN(jte,jde-1)
3542 ! higher order flux has a 5 or 7 point stencil, so compute
3543 ! bounds so we can switch to second order flux close to the boundary
3549 i_start = MAX(ids+1,its)
3550 i_start_f = i_start+2
3554 i_end = MIN(ide-2,ite)
3560 DO j = j_start, j_end
3562 ! 5th or 6th order flux
3565 DO i = i_start_f, i_end_f
3567 fqx( i,k ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), &
3568 field(i-1,k,j), field(i ,k,j), &
3569 field(i+1,k,j), field(i+2,k,j), &
3574 ! lower order fluxes close to boundaries (if not periodic or symmetric)
3576 IF( degrade_xs ) THEN
3578 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
3581 fqx(i,k) = 0.5*(ru(i,k,j)) &
3582 *(field(i,k,j)+field(i-1,k,j))
3590 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
3591 field(i ,k,j), field(i+1,k,j), &
3597 IF( degrade_xe ) THEN
3599 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
3602 fqx(i,k) = 0.5*(ru(i,k,j)) &
3603 *(field(i,k,j)+field(i-1,k,j))
3610 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
3611 field(i ,k,j), field(i+1,k,j), &
3617 ! x flux-divergence into tendency
3620 DO i = i_start, i_end
3621 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3622 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3629 ELSE IF( horz_order == 4 ) THEN
3636 IF( config_flags%periodic_x .or. &
3637 config_flags%symmetric_xs .or. &
3638 (its > ids+1) ) degrade_xs = .false.
3639 IF( config_flags%periodic_x .or. &
3640 config_flags%symmetric_xe .or. &
3641 (ite < ide-2) ) degrade_xe = .false.
3642 IF( config_flags%periodic_y .or. &
3643 config_flags%symmetric_ys .or. &
3644 (jts > jds+1) ) degrade_ys = .false.
3645 IF( config_flags%periodic_y .or. &
3646 config_flags%symmetric_ye .or. &
3647 (jte < jde-2) ) degrade_ye = .false.
3649 ! begin flux computations
3650 ! start with x flux divergence
3655 i_end = MIN(ite,ide-1)
3657 j_end = MIN(jte,jde-1)
3659 ! 3rd or 4th order flux has a 5 point stencil, so compute
3660 ! bounds so we can switch to second order flux close to the boundary
3667 i_start_f = i_start+1
3677 DO j = j_start, j_end
3679 ! 3rd or 4th order flux
3682 DO i = i_start_f, i_end_f
3684 fqx( i, k) = ru(i,k,j)*flux4( field(i-2,k,j), field(i-1,k,j), &
3685 field(i ,k,j), field(i+1,k,j), &
3690 ! second order flux close to boundaries (if not periodic or symmetric)
3692 IF( degrade_xs ) THEN
3694 fqx(i_start, k) = 0.5*ru(i_start,k,j) &
3695 *(field(i_start,k,j)+field(i_start-1,k,j))
3699 IF( degrade_xe ) THEN
3701 fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) &
3702 *(field(i_end+1,k,j)+field(i_end,k,j))
3706 ! x flux-divergence into tendency
3709 DO i = i_start, i_end
3710 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3711 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3718 ! next -> y flux divergence calculation
3721 i_end = MIN(ite,ide-1)
3723 j_end = MIN(jte,jde-1)
3725 ! 3rd or 4th order flux has a 5 point stencil, so compute
3726 ! bounds so we can switch to second order flux close to the boundary
3733 j_start_f = j_start+1
3741 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3746 DO j = j_start, j_end+1
3748 IF ((j < j_start_f) .and. degrade_ys) THEN
3750 DO i = i_start, i_end
3751 fqy(i,k,jp1) = 0.5*rv(i,k,j_start) &
3752 *(field(i,k,j_start)+field(i,k,j_start-1))
3755 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
3757 DO i = i_start, i_end
3758 ! Assumes j>j_end_f is ONLY j_end+1 ...
3759 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) &
3760 ! *(field(i,k,j_end+1)+field(i,k,j_end))
3761 fqy(i,k,jp1) = 0.5*rv(i,k,j) &
3762 *(field(i,k,j)+field(i,k,j-1))
3766 ! 3rd or 4th order flux
3768 DO i = i_start, i_end
3769 fqy( i, k, jp1 ) = rv(i,k,j)*flux4( field(i,k,j-2), field(i,k,j-1), &
3770 field(i,k,j ), field(i,k,j+1), &
3776 ! y flux-divergence into tendency
3778 ! Comments on polar boundary conditions
3779 ! Same process as for advect_u - tendencies run from jds to jde-1
3780 ! (latitudes are as for u grid, longitudes are displaced)
3781 ! Therefore: flow is only from one side for points next to poles
3782 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3784 DO i = i_start, i_end
3785 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3786 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3789 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3791 DO i = i_start, i_end
3792 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3793 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3798 IF ( j > j_start ) THEN
3801 DO i = i_start, i_end
3802 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3803 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3818 ELSE IF( horz_order == 3 ) THEN
3825 IF( config_flags%periodic_x .or. &
3826 config_flags%symmetric_xs .or. &
3827 (its > ids+1) ) degrade_xs = .false.
3828 IF( config_flags%periodic_x .or. &
3829 config_flags%symmetric_xe .or. &
3830 (ite < ide-2) ) degrade_xe = .false.
3831 IF( config_flags%periodic_y .or. &
3832 config_flags%symmetric_ys .or. &
3833 (jts > jds+1) ) degrade_ys = .false.
3834 IF( config_flags%periodic_y .or. &
3835 config_flags%symmetric_ye .or. &
3836 (jte < jde-2) ) degrade_ye = .false.
3838 ! begin flux computations
3839 ! start with x flux divergence
3844 i_end = MIN(ite,ide-1)
3846 j_end = MIN(jte,jde-1)
3848 ! 3rd or 4th order flux has a 5 point stencil, so compute
3849 ! bounds so we can switch to second order flux close to the boundary
3856 i_start_f = i_start+1
3866 DO j = j_start, j_end
3868 ! 3rd or 4th order flux
3871 DO i = i_start_f, i_end_f
3873 fqx( i, k) = ru(i,k,j)*flux3( field(i-2,k,j), field(i-1,k,j), &
3874 field(i ,k,j), field(i+1,k,j), &
3879 ! second order flux close to boundaries (if not periodic or symmetric)
3881 IF( degrade_xs ) THEN
3883 fqx(i_start, k) = 0.5*ru(i_start,k,j) &
3884 *(field(i_start,k,j)+field(i_start-1,k,j))
3888 IF( degrade_xe ) THEN
3890 fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) &
3891 *(field(i_end+1,k,j)+field(i_end,k,j))
3895 ! x flux-divergence into tendency
3898 DO i = i_start, i_end
3899 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3900 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3907 ! next -> y flux divergence calculation
3910 i_end = MIN(ite,ide-1)
3912 j_end = MIN(jte,jde-1)
3914 ! 3rd or 4th order flux has a 5 point stencil, so compute
3915 ! bounds so we can switch to second order flux close to the boundary
3922 j_start_f = j_start+1
3930 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3935 DO j = j_start, j_end+1
3937 IF ((j < j_start_f) .and. degrade_ys) THEN
3939 DO i = i_start, i_end
3940 fqy(i,k,jp1) = 0.5*rv(i,k,j_start) &
3941 *(field(i,k,j_start)+field(i,k,j_start-1))
3944 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
3946 DO i = i_start, i_end
3947 ! Assumes j>j_end_f is ONLY j_end+1 ...
3948 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) &
3949 ! *(field(i,k,j_end+1)+field(i,k,j_end))
3950 fqy(i,k,jp1) = 0.5*rv(i,k,j) &
3951 *(field(i,k,j)+field(i,k,j-1))
3955 ! 3rd or 4th order flux
3957 DO i = i_start, i_end
3958 fqy( i, k, jp1 ) = rv(i,k,j)*flux3( field(i,k,j-2), field(i,k,j-1), &
3959 field(i,k,j ), field(i,k,j+1), &
3965 ! y flux-divergence into tendency
3967 ! Comments on polar boundary conditions
3968 ! Same process as for advect_u - tendencies run from jds to jde-1
3969 ! (latitudes are as for u grid, longitudes are displaced)
3970 ! Therefore: flow is only from one side for points next to poles
3971 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3973 DO i = i_start, i_end
3974 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3975 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3978 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3980 DO i = i_start, i_end
3981 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3982 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3987 IF ( j > j_start ) THEN
3990 DO i = i_start, i_end
3991 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3992 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
4006 ELSE IF( horz_order == 2 ) THEN
4009 i_end = MIN(ite,ide-1)
4011 j_end = MIN(jte,jde-1)
4013 IF ( .NOT. config_flags%periodic_x ) THEN
4014 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
4015 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
4018 DO j = j_start, j_end
4020 DO i = i_start, i_end
4021 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
4022 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
4023 *(ru(i+1,k,j)*(field(i+1,k,j)+field(i ,k,j)) &
4024 -ru(i ,k,j)*(field(i ,k,j)+field(i-1,k,j)))
4030 i_end = MIN(ite,ide-1)
4032 ! Polar boundary conditions are like open or specified
4033 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
4034 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-2,jte)
4036 DO j = j_start, j_end
4038 DO i = i_start, i_end
4039 mrdy=msftx(i,j)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
4040 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
4041 *(rv(i,k,j+1)*(field(i,k,j+1)+field(i,k,j )) &
4042 -rv(i,k,j )*(field(i,k,j )+field(i,k,j-1)))
4047 ! Polar boundary condtions
4048 ! These won't be covered in the loop above...
4049 IF (config_flags%polar) THEN
4050 IF (jts == jds) THEN
4052 DO i = i_start, i_end
4053 mrdy=msftx(i,jds)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
4054 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
4055 *rv(i,k,jds+1)*(field(i,k,jds+1)+field(i,k,jds))
4059 IF (jte == jde) THEN
4061 DO i = i_start, i_end
4062 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
4063 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
4064 *rv(i,k,jde-1)*(field(i,k,jde-1)+field(i,k,jde-2))
4070 ELSE IF ( horz_order == 0 ) THEN
4072 ! Just in case we want to turn horizontal advection off, we can do it
4076 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_6a, h_order not known ',horz_order
4077 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
4079 ENDIF horizontal_order_test
4081 ! pick up the rest of the horizontal radiation boundary conditions.
4082 ! (these are the computations that don't require 'cb'.
4083 ! first, set to index ranges
4086 i_end = MIN(ite,ide-1)
4088 j_end = MIN(jte,jde-1)
4090 ! compute x (u) conditions for v, w, or scalar
4092 IF( (config_flags%open_xs) .and. (its == ids) ) THEN
4094 DO j = j_start, j_end
4096 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
4097 tendency(its,k,j) = tendency(its,k,j) &
4099 ub*( field_old(its+1,k,j) &
4100 - field_old(its ,k,j) ) + &
4101 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) &
4108 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
4110 DO j = j_start, j_end
4112 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
4113 tendency(i_end,k,j) = tendency(i_end,k,j) &
4115 ub*( field_old(i_end ,k,j) &
4116 - field_old(i_end-1,k,j) ) + &
4117 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) &
4124 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
4126 DO i = i_start, i_end
4128 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
4129 tendency(i,k,jts) = tendency(i,k,jts) &
4131 vb*( field_old(i,k,jts+1) &
4132 - field_old(i,k,jts ) ) + &
4133 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) &
4140 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
4142 DO i = i_start, i_end
4144 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
4145 tendency(i,k,j_end) = tendency(i,k,j_end) &
4147 vb*( field_old(i,k,j_end ) &
4148 - field_old(i,k,j_end-1) ) + &
4149 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) &
4157 !-------------------- vertical advection
4158 ! Scalar equation has 3rd term on RHS = - partial d/dz (q rho w /my)
4159 ! Here we have: - partial d/dz (q*rom) = - partial d/dz (q rho w / my)
4160 ! So we don't need to make a correction for advect_scalar
4163 i_end = MIN(ite,ide-1)
4165 j_end = MIN(jte,jde-1)
4167 DO i = i_start, i_end
4172 vert_order_test : IF (vert_order == 6) THEN
4174 DO j = j_start, j_end
4177 DO i = i_start, i_end
4179 vflux(i,k) = vel*flux6( &
4180 field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
4181 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
4185 DO i = i_start, i_end
4188 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4192 vflux(i,k) = vel*flux4( &
4193 field(i,k-2,j), field(i,k-1,j), &
4194 field(i,k ,j), field(i,k+1,j), -vel )
4197 vflux(i,k) = vel*flux4( &
4198 field(i,k-2,j), field(i,k-1,j), &
4199 field(i,k ,j), field(i,k+1,j), -vel )
4202 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4206 DO i = i_start, i_end
4207 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4213 ELSE IF (vert_order == 5) THEN
4215 DO j = j_start, j_end
4218 DO i = i_start, i_end
4220 vflux(i,k) = vel*flux5( &
4221 field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
4222 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
4226 DO i = i_start, i_end
4229 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4233 vflux(i,k) = vel*flux3( &
4234 field(i,k-2,j), field(i,k-1,j), &
4235 field(i,k ,j), field(i,k+1,j), -vel )
4238 vflux(i,k) = vel*flux3( &
4239 field(i,k-2,j), field(i,k-1,j), &
4240 field(i,k ,j), field(i,k+1,j), -vel )
4243 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4247 DO i = i_start, i_end
4248 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4254 ELSE IF (vert_order == 4) THEN
4256 DO j = j_start, j_end
4259 DO i = i_start, i_end
4261 vflux(i,k) = vel*flux4( &
4262 field(i,k-2,j), field(i,k-1,j), &
4263 field(i,k ,j), field(i,k+1,j), -vel )
4267 DO i = i_start, i_end
4270 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4272 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4276 DO i = i_start, i_end
4277 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4283 ELSE IF (vert_order == 3) THEN
4285 DO j = j_start, j_end
4288 DO i = i_start, i_end
4290 vflux(i,k) = vel*flux3( &
4291 field(i,k-2,j), field(i,k-1,j), &
4292 field(i,k ,j), field(i,k+1,j), -vel )
4296 DO i = i_start, i_end
4299 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4301 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4305 DO i = i_start, i_end
4306 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4313 ELSE IF (vert_order == 2) THEN
4315 DO j = j_start, j_end
4317 DO i = i_start, i_end
4318 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4323 DO i = i_start, i_end
4324 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4332 WRITE (wrf_err_message,*) ' advect_scalar_6a, v_order not known ',vert_order
4333 CALL wrf_error_fatal ( wrf_err_message )
4335 ENDIF vert_order_test
4337 END SUBROUTINE advect_scalar
4339 !---------------------------------------------------------------------------------
4341 SUBROUTINE advect_w ( w, w_old, tendency, &
4343 mut, time_step, config_flags, &
4344 msfux, msfuy, msfvx, msfvy, &
4348 ids, ide, jds, jde, kds, kde, &
4349 ims, ime, jms, jme, kms, kme, &
4350 its, ite, jts, jte, kts, kte )
4356 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4358 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4359 ims, ime, jms, jme, kms, kme, &
4360 its, ite, jts, jte, kts, kte
4362 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: w, &
4368 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
4369 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4371 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4378 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4382 REAL , INTENT(IN ) :: rdx, &
4384 INTEGER , INTENT(IN ) :: time_step
4389 INTEGER :: i, j, k, itf, jtf, ktf
4390 INTEGER :: i_start, i_end, j_start, j_end
4391 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
4392 INTEGER :: jmin, jmax, jp, jm, imin, imax
4394 REAL :: mrdx, mrdy, ub, vb, uw, vw
4395 REAL , DIMENSION(its:ite, kts:kte) :: vflux
4397 INTEGER :: horz_order, vert_order
4399 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
4400 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
4402 LOGICAL :: degrade_xs, degrade_ys
4403 LOGICAL :: degrade_xe, degrade_ye
4405 INTEGER :: jp1, jp0, jtmp
4407 ! definition of flux operators, 3rd, 4th, 5th or 6th order
4409 REAL :: flux3, flux4, flux5, flux6
4410 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
4412 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
4413 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
4415 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
4416 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
4417 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
4419 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
4420 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
4421 +(q_ip2+q_im3) )/60.0
4423 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
4424 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
4425 -sign(1,time_step)*sign(1.,ua)*( &
4426 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
4429 LOGICAL :: specified
4432 if(config_flags%specified .or. config_flags%nested) specified = .true.
4434 ! set order for the advection scheme
4437 horz_order = config_flags%h_sca_adv_order
4438 vert_order = config_flags%v_sca_adv_order
4440 ! here is the choice of flux operators
4442 ! begin with horizontal flux divergence
4444 horizontal_order_test : IF( horz_order == 6 ) THEN
4446 ! determine boundary mods for flux operators
4447 ! We degrade the flux operators from 3rd/4th order
4448 ! to second order one gridpoint in from the boundaries for
4449 ! all boundary conditions except periodic and symmetry - these
4450 ! conditions have boundary zone data fill for correct application
4451 ! of the higher order flux stencils
4458 IF( config_flags%periodic_x .or. &
4459 config_flags%symmetric_xs .or. &
4460 (its > ids+2) ) degrade_xs = .false.
4461 IF( config_flags%periodic_x .or. &
4462 config_flags%symmetric_xe .or. &
4463 (ite < ide-3) ) degrade_xe = .false.
4464 IF( config_flags%periodic_y .or. &
4465 config_flags%symmetric_ys .or. &
4466 (jts > jds+2) ) degrade_ys = .false.
4467 IF( config_flags%periodic_y .or. &
4468 config_flags%symmetric_ye .or. &
4469 (jte < jde-3) ) degrade_ye = .false.
4471 !--------------- y - advection first
4474 i_end = MIN(ite,ide-1)
4476 j_end = MIN(jte,jde-1)
4478 ! higher order flux has a 5 or 7 point stencil, so compute
4479 ! bounds so we can switch to second order flux close to the boundary
4485 j_start = MAX(jts,jds+1)
4490 j_end = MIN(jte,jde-2)
4494 IF(config_flags%polar) j_end = MIN(jte,jde-1)
4496 ! compute fluxes, 5th or 6th order
4501 j_loop_y_flux_6 : DO j = j_start, j_end+1
4503 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
4506 DO i = i_start, i_end
4507 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4508 fqy( i, k, jp1 ) = vel*flux6( &
4509 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4510 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4515 DO i = i_start, i_end
4516 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4517 fqy( i, k, jp1 ) = vel*flux6( &
4518 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4519 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4522 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
4525 DO i = i_start, i_end
4526 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4527 (w(i,k,j)+w(i,k,j-1))
4532 DO i = i_start, i_end
4533 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4534 (w(i,k,j)+w(i,k,j-1))
4537 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
4540 DO i = i_start, i_end
4541 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4542 fqy( i, k, jp1 ) = vel*flux4( &
4543 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4548 DO i = i_start, i_end
4549 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4550 fqy( i, k, jp1 ) = vel*flux4( &
4551 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4554 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
4557 DO i = i_start, i_end
4558 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4559 (w(i,k,j)+w(i,k,j-1))
4564 DO i = i_start, i_end
4565 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4566 (w(i,k,j)+w(i,k,j-1))
4569 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
4572 DO i = i_start, i_end
4573 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4574 fqy( i, k, jp1 ) = vel*flux4( &
4575 w(i,k,j-2),w(i,k,j-1), &
4576 w(i,k,j),w(i,k,j+1),vel )
4581 DO i = i_start, i_end
4582 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4583 fqy( i, k, jp1 ) = vel*flux4( &
4584 w(i,k,j-2),w(i,k,j-1), &
4585 w(i,k,j),w(i,k,j+1),vel )
4590 ! y flux-divergence into tendency
4592 ! Comments for polar boundary conditions
4593 ! Same process as for advect_u - tendencies run from jds to jde-1
4594 ! (latitudes are as for u grid, longitudes are displaced)
4595 ! Therefore: flow is only from one side for points next to poles
4596 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
4598 DO i = i_start, i_end
4599 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4600 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
4603 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
4605 DO i = i_start, i_end
4606 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4607 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
4612 IF(j > j_start) THEN
4615 DO i = i_start, i_end
4616 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4617 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
4629 ENDDO j_loop_y_flux_6
4631 ! next, x - flux divergence
4634 i_end = MIN(ite,ide-1)
4637 j_end = MIN(jte,jde-1)
4639 ! higher order flux has a 5 or 7 point stencil, so compute
4640 ! bounds so we can switch to second order flux close to the boundary
4646 i_start = MAX(ids+1,its)
4647 i_start_f = i_start+2
4651 i_end = MIN(ide-2,ite)
4657 DO j = j_start, j_end
4659 ! 5th or 6th order flux
4662 DO i = i_start_f, i_end_f
4663 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4664 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &
4665 w(i-1,k,j), w(i ,k,j), &
4666 w(i+1,k,j), w(i+2,k,j), &
4672 DO i = i_start_f, i_end_f
4673 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4674 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &
4675 w(i-1,k,j), w(i ,k,j), &
4676 w(i+1,k,j), w(i+2,k,j), &
4680 ! lower order fluxes close to boundaries (if not periodic or symmetric)
4682 IF( degrade_xs ) THEN
4684 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
4687 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4688 *(w(i,k,j)+w(i-1,k,j))
4691 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4692 *(w(i,k,j)+w(i-1,k,j))
4697 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4698 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4699 w(i ,k,j), w(i+1,k,j), &
4704 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4705 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4706 w(i ,k,j), w(i+1,k,j), &
4710 IF( degrade_xe ) THEN
4712 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
4715 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4716 *(w(i,k,j)+w(i-1,k,j))
4719 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4720 *(w(i,k,j)+w(i-1,k,j))
4725 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4726 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4727 w(i ,k,j), w(i+1,k,j), &
4732 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4733 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4734 w(i ,k,j), w(i+1,k,j), &
4738 ! x flux-divergence into tendency
4741 DO i = i_start, i_end
4742 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
4743 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
4750 ELSE IF (horz_order == 5 ) THEN
4752 ! determine boundary mods for flux operators
4753 ! We degrade the flux operators from 3rd/4th order
4754 ! to second order one gridpoint in from the boundaries for
4755 ! all boundary conditions except periodic and symmetry - these
4756 ! conditions have boundary zone data fill for correct application
4757 ! of the higher order flux stencils
4764 IF( config_flags%periodic_x .or. &
4765 config_flags%symmetric_xs .or. &
4766 (its > ids+2) ) degrade_xs = .false.
4767 IF( config_flags%periodic_x .or. &
4768 config_flags%symmetric_xe .or. &
4769 (ite < ide-3) ) degrade_xe = .false.
4770 IF( config_flags%periodic_y .or. &
4771 config_flags%symmetric_ys .or. &
4772 (jts > jds+2) ) degrade_ys = .false.
4773 IF( config_flags%periodic_y .or. &
4774 config_flags%symmetric_ye .or. &
4775 (jte < jde-3) ) degrade_ye = .false.
4777 !--------------- y - advection first
4780 i_end = MIN(ite,ide-1)
4782 j_end = MIN(jte,jde-1)
4784 ! higher order flux has a 5 or 7 point stencil, so compute
4785 ! bounds so we can switch to second order flux close to the boundary
4791 j_start = MAX(jts,jds+1)
4796 j_end = MIN(jte,jde-2)
4800 IF(config_flags%polar) j_end = MIN(jte,jde-1)
4802 ! compute fluxes, 5th or 6th order
4807 j_loop_y_flux_5 : DO j = j_start, j_end+1
4809 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
4812 DO i = i_start, i_end
4813 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4814 fqy( i, k, jp1 ) = vel*flux5( &
4815 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4816 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4821 DO i = i_start, i_end
4822 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4823 fqy( i, k, jp1 ) = vel*flux5( &
4824 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4825 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4828 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
4831 DO i = i_start, i_end
4832 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4833 (w(i,k,j)+w(i,k,j-1))
4838 DO i = i_start, i_end
4839 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4840 (w(i,k,j)+w(i,k,j-1))
4843 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
4846 DO i = i_start, i_end
4847 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4848 fqy( i, k, jp1 ) = vel*flux3( &
4849 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4854 DO i = i_start, i_end
4855 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4856 fqy( i, k, jp1 ) = vel*flux3( &
4857 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4860 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
4863 DO i = i_start, i_end
4864 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4865 (w(i,k,j)+w(i,k,j-1))
4870 DO i = i_start, i_end
4871 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4872 (w(i,k,j)+w(i,k,j-1))
4875 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
4878 DO i = i_start, i_end
4879 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4880 fqy( i, k, jp1 ) = vel*flux3( &
4881 w(i,k,j-2),w(i,k,j-1), &
4882 w(i,k,j),w(i,k,j+1),vel )
4887 DO i = i_start, i_end
4888 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4889 fqy( i, k, jp1 ) = vel*flux3( &
4890 w(i,k,j-2),w(i,k,j-1), &
4891 w(i,k,j),w(i,k,j+1),vel )
4896 ! y flux-divergence into tendency
4898 ! Comments for polar boundary conditions
4899 ! Same process as for advect_u - tendencies run from jds to jde-1
4900 ! (latitudes are as for u grid, longitudes are displaced)
4901 ! Therefore: flow is only from one side for points next to poles
4902 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
4904 DO i = i_start, i_end
4905 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4906 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
4909 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
4911 DO i = i_start, i_end
4912 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4913 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
4918 IF(j > j_start) THEN
4921 DO i = i_start, i_end
4922 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4923 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
4935 ENDDO j_loop_y_flux_5
4937 ! next, x - flux divergence
4940 i_end = MIN(ite,ide-1)
4943 j_end = MIN(jte,jde-1)
4945 ! higher order flux has a 5 or 7 point stencil, so compute
4946 ! bounds so we can switch to second order flux close to the boundary
4952 i_start = MAX(ids+1,its)
4953 i_start_f = i_start+2
4957 i_end = MIN(ide-2,ite)
4963 DO j = j_start, j_end
4965 ! 5th or 6th order flux
4968 DO i = i_start_f, i_end_f
4969 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4970 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &
4971 w(i-1,k,j), w(i ,k,j), &
4972 w(i+1,k,j), w(i+2,k,j), &
4978 DO i = i_start_f, i_end_f
4979 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4980 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &
4981 w(i-1,k,j), w(i ,k,j), &
4982 w(i+1,k,j), w(i+2,k,j), &
4986 ! lower order fluxes close to boundaries (if not periodic or symmetric)
4988 IF( degrade_xs ) THEN
4990 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
4993 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4994 *(w(i,k,j)+w(i-1,k,j))
4997 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4998 *(w(i,k,j)+w(i-1,k,j))
5003 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5004 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5005 w(i ,k,j), w(i+1,k,j), &
5009 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5010 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5011 w(i ,k,j), w(i+1,k,j), &
5016 IF( degrade_xe ) THEN
5018 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
5021 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
5022 *(w(i,k,j)+w(i-1,k,j))
5025 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
5026 *(w(i,k,j)+w(i-1,k,j))
5031 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5032 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5033 w(i ,k,j), w(i+1,k,j), &
5037 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5038 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5039 w(i ,k,j), w(i+1,k,j), &
5043 ! x flux-divergence into tendency
5046 DO i = i_start, i_end
5047 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5048 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
5054 ELSE IF ( horz_order == 4 ) THEN
5061 IF( config_flags%periodic_x .or. &
5062 config_flags%symmetric_xs .or. &
5063 (its > ids+1) ) degrade_xs = .false.
5064 IF( config_flags%periodic_x .or. &
5065 config_flags%symmetric_xe .or. &
5066 (ite < ide-2) ) degrade_xe = .false.
5067 IF( config_flags%periodic_y .or. &
5068 config_flags%symmetric_ys .or. &
5069 (jts > jds+1) ) degrade_ys = .false.
5070 IF( config_flags%periodic_y .or. &
5071 config_flags%symmetric_ye .or. &
5072 (jte < jde-2) ) degrade_ye = .false.
5074 ! begin flux computations
5075 ! start with x flux divergence
5082 i_end = MIN(ite,ide-1)
5084 j_end = MIN(jte,jde-1)
5086 ! 3rd or 4th order flux has a 5 point stencil, so compute
5087 ! bounds so we can switch to second order flux close to the boundary
5094 i_start_f = i_start+1
5104 DO j = j_start, j_end
5107 DO i = i_start_f, i_end_f
5108 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5109 fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
5110 w(i ,k,j), w(i+1,k,j), &
5116 DO i = i_start_f, i_end_f
5117 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5118 fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
5119 w(i ,k,j), w(i+1,k,j), &
5122 ! second order flux close to boundaries (if not periodic or symmetric)
5124 IF( degrade_xs ) THEN
5127 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) &
5128 *(w(i_start,k,j)+w(i_start-1,k,j))
5132 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) &
5133 *(w(i_start,k,j)+w(i_start-1,k,j))
5136 IF( degrade_xe ) THEN
5139 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) &
5140 *(w(i_end+1,k,j)+w(i_end,k,j))
5144 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) &
5145 *(w(i_end+1,k,j)+w(i_end,k,j))
5148 ! x flux-divergence into tendency
5151 DO i = i_start, i_end
5152 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5153 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
5159 ! next -> y flux divergence calculation
5162 i_end = MIN(ite,ide-1)
5164 j_end = MIN(jte,jde-1)
5167 ! 3rd or 4th order flux has a 5 point stencil, so compute
5168 ! bounds so we can switch to second order flux close to the boundary
5175 j_start_f = j_start+1
5183 IF(config_flags%polar) j_end = MIN(jte,jde-1)
5188 DO j = j_start, j_end+1
5190 IF ((j < j_start_f) .and. degrade_ys) THEN
5192 DO i = i_start, i_end
5194 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) &
5195 *(w(i,k,j_start)+w(i,k,j_start-1))
5199 DO i = i_start, i_end
5201 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) &
5202 *(w(i,k,j_start)+w(i,k,j_start-1))
5204 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
5206 DO i = i_start, i_end
5207 ! Assumes j>j_end_f is ONLY j_end+1 ...
5208 ! fqy(i, k, jp1) = &
5209 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) &
5210 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5212 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5213 *(w(i,k,j)+w(i,k,j-1))
5217 DO i = i_start, i_end
5218 ! Assumes j>j_end_f is ONLY j_end+1 ...
5219 ! fqy(i, k, jp1) = &
5220 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &
5221 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5223 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5224 *(w(i,k,j)+w(i,k,j-1))
5227 ! 3rd or 4th order flux
5229 DO i = i_start, i_end
5230 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
5231 fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), &
5232 w(i,k,j ), w(i,k,j+1), &
5237 DO i = i_start, i_end
5238 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
5239 fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), &
5240 w(i,k,j ), w(i,k,j+1), &
5245 ! y flux-divergence into tendency
5247 ! Comments for polar boundary conditions
5248 ! Same process as for advect_u - tendencies run from jds to jde-1
5249 ! (latitudes are as for u grid, longitudes are displaced)
5250 ! Therefore: flow is only from one side for points next to poles
5251 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
5253 DO i = i_start, i_end
5254 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5255 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
5258 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
5260 DO i = i_start, i_end
5261 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5262 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
5267 IF( j > j_start ) THEN
5270 DO i = i_start, i_end
5271 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5272 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
5286 ELSE IF ( horz_order == 3 ) THEN
5293 IF( config_flags%periodic_x .or. &
5294 config_flags%symmetric_xs .or. &
5295 (its > ids+1) ) degrade_xs = .false.
5296 IF( config_flags%periodic_x .or. &
5297 config_flags%symmetric_xe .or. &
5298 (ite < ide-2) ) degrade_xe = .false.
5299 IF( config_flags%periodic_y .or. &
5300 config_flags%symmetric_ys .or. &
5301 (jts > jds+1) ) degrade_ys = .false.
5302 IF( config_flags%periodic_y .or. &
5303 config_flags%symmetric_ye .or. &
5304 (jte < jde-2) ) degrade_ye = .false.
5306 ! begin flux computations
5307 ! start with x flux divergence
5314 i_end = MIN(ite,ide-1)
5316 j_end = MIN(jte,jde-1)
5318 ! 3rd or 4th order flux has a 5 point stencil, so compute
5319 ! bounds so we can switch to second order flux close to the boundary
5326 i_start_f = i_start+1
5336 DO j = j_start, j_end
5339 DO i = i_start_f, i_end_f
5340 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5341 fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5342 w(i ,k,j), w(i+1,k,j), &
5347 DO i = i_start_f, i_end_f
5348 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5349 fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5350 w(i ,k,j), w(i+1,k,j), &
5354 ! second order flux close to boundaries (if not periodic or symmetric)
5356 IF( degrade_xs ) THEN
5359 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) &
5360 *(w(i_start,k,j)+w(i_start-1,k,j))
5364 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) &
5365 *(w(i_start,k,j)+w(i_start-1,k,j))
5368 IF( degrade_xe ) THEN
5371 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) &
5372 *(w(i_end+1,k,j)+w(i_end,k,j))
5376 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) &
5377 *(w(i_end+1,k,j)+w(i_end,k,j))
5380 ! x flux-divergence into tendency
5383 DO i = i_start, i_end
5384 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5385 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
5391 ! next -> y flux divergence calculation
5394 i_end = MIN(ite,ide-1)
5396 j_end = MIN(jte,jde-1)
5399 ! 3rd or 4th order flux has a 5 point stencil, so compute
5400 ! bounds so we can switch to second order flux close to the boundary
5407 j_start_f = j_start+1
5415 IF(config_flags%polar) j_end = MIN(jte,jde-1)
5420 DO j = j_start, j_end+1
5422 IF ((j < j_start_f) .and. degrade_ys) THEN
5424 DO i = i_start, i_end
5426 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) &
5427 *(w(i,k,j_start)+w(i,k,j_start-1))
5431 DO i = i_start, i_end
5433 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) &
5434 *(w(i,k,j_start)+w(i,k,j_start-1))
5436 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
5438 DO i = i_start, i_end
5439 ! Assumes j>j_end_f is ONLY j_end+1 ...
5440 ! fqy(i, k, jp1) = &
5441 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) &
5442 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5444 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5445 *(w(i,k,j)+w(i,k,j-1))
5449 DO i = i_start, i_end
5450 ! Assumes j>j_end_f is ONLY j_end+1 ...
5451 ! fqy(i, k, jp1) = &
5452 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &
5453 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5455 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5456 *(w(i,k,j)+w(i,k,j-1))
5459 ! 3rd or 4th order flux
5461 DO i = i_start, i_end
5462 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
5463 fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), &
5464 w(i,k,j ), w(i,k,j+1), &
5469 DO i = i_start, i_end
5470 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
5471 fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), &
5472 w(i,k,j ), w(i,k,j+1), &
5477 ! y flux-divergence into tendency
5479 ! Comments for polar boundary conditions
5480 ! Same process as for advect_u - tendencies run from jds to jde-1
5481 ! (latitudes are as for u grid, longitudes are displaced)
5482 ! Therefore: flow is only from one side for points next to poles
5483 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
5485 DO i = i_start, i_end
5486 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5487 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
5490 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
5492 DO i = i_start, i_end
5493 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5494 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
5499 IF( j > j_start ) THEN
5502 DO i = i_start, i_end
5503 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5504 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
5518 ELSE IF (horz_order == 2 ) THEN
5521 i_end = MIN(ite,ide-1)
5523 j_end = MIN(jte,jde-1)
5525 IF ( .NOT. config_flags%periodic_x ) THEN
5526 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
5527 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
5530 DO j = j_start, j_end
5532 DO i = i_start, i_end
5534 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5536 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
5537 *((fzm(k)*ru(i+1,k,j)+fzp(k)*ru(i+1,k-1,j)) &
5538 *(w(i+1,k,j)+w(i,k,j)) &
5539 -(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
5540 *(w(i,k,j)+w(i-1,k,j)))
5546 DO i = i_start, i_end
5548 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5550 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
5551 *(((2.-fzm(k-1))*ru(i+1,k-1,j)-fzp(k-1)*ru(i+1,k-2,j)) &
5552 *(w(i+1,k,j)+w(i,k,j)) &
5553 -((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
5554 *(w(i,k,j)+w(i-1,k,j)))
5561 i_end = MIN(ite,ide-1)
5562 ! Polar boundary conditions are like open or specified
5563 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
5564 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-2,jte)
5566 DO j = j_start, j_end
5568 DO i = i_start, i_end
5570 mrdy=msftx(i,j)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5572 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
5573 *((fzm(k)*rv(i,k,j+1)+fzp(k)*rv(i,k-1,j+1))* &
5574 (w(i,k,j+1)+w(i,k,j)) &
5575 -(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5576 *(w(i,k,j)+w(i,k,j-1)))
5582 DO i = i_start, i_end
5584 mrdy=msftx(i,j)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5586 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
5587 *(((2.-fzm(k-1))*rv(i,k-1,j+1)-fzp(k-1)*rv(i,k-2,j+1))* &
5588 (w(i,k,j+1)+w(i,k,j)) &
5589 -((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5590 *(w(i,k,j)+w(i,k,j-1)))
5596 ! Polar boundary condition ... not covered in above j-loop
5597 IF (config_flags%polar) THEN
5598 IF (jts == jds) THEN
5600 DO i = i_start, i_end
5601 mrdy=msftx(i,jds)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5602 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
5603 *((fzm(k)*rv(i,k,jds+1)+fzp(k)*rv(i,k-1,jds+1))* &
5604 (w(i,k,jds+1)+w(i,k,jds)))
5608 DO i = i_start, i_end
5609 mrdy=msftx(i,jds)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5610 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
5611 *((2.-fzm(k-1))*rv(i,k-1,jds+1)-fzp(k-1)*rv(i,k-2,jds+1))* &
5612 (w(i,k,jds+1)+w(i,k,jds))
5615 IF (jte == jde) THEN
5617 DO i = i_start, i_end
5618 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5619 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
5620 *((fzm(k)*rv(i,k,jde-1)+fzp(k)*rv(i,k-1,jde-1))* &
5621 (w(i,k,jde-1)+w(i,k,jde-2)))
5625 DO i = i_start, i_end
5626 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5627 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
5628 *((2.-fzm(k-1))*rv(i,k-1,jde-1)-fzp(k-1)*rv(i,k-2,jde-1)) &
5629 *(w(i,k,jde-1)+w(i,k,jde-2))
5634 ELSE IF ( horz_order == 0 ) THEN
5636 ! Just in case we want to turn horizontal advection off, we can do it
5640 WRITE ( wrf_err_message ,*) ' advect_w_6a, h_order not known ',horz_order
5641 CALL wrf_error_fatal ( wrf_err_message )
5643 ENDIF horizontal_order_test
5646 ! pick up the the horizontal radiation boundary conditions.
5647 ! (these are the computations that don't require 'cb'.
5648 ! first, set to index ranges
5652 i_end = MIN(ite,ide-1)
5654 j_end = MIN(jte,jde-1)
5656 IF( (config_flags%open_xs) .and. (its == ids)) THEN
5658 DO j = j_start, j_end
5661 uw = 0.5*(fzm(k)*(ru(its,k ,j)+ru(its+1,k ,j)) + &
5662 fzp(k)*(ru(its,k-1,j)+ru(its+1,k-1,j)) )
5665 tendency(its,k,j) = tendency(its,k,j) &
5667 ub*(w_old(its+1,k,j) - w_old(its,k,j)) + &
5669 fzm(k)*(ru(its+1,k ,j)-ru(its,k ,j))+ &
5670 fzp(k)*(ru(its+1,k-1,j)-ru(its,k-1,j))) &
5676 DO j = j_start, j_end
5678 uw = 0.5*( (2.-fzm(k-1))*(ru(its,k-1,j)+ru(its+1,k-1,j)) &
5679 -fzp(k-1)*(ru(its,k-2,j)+ru(its+1,k-2,j)) )
5682 tendency(its,k,j) = tendency(its,k,j) &
5684 ub*(w_old(its+1,k,j) - w_old(its,k,j)) + &
5686 (2.-fzm(k-1))*(ru(its+1,k-1,j)-ru(its,k-1,j))- &
5687 fzp(k-1)*(ru(its+1,k-2,j)-ru(its,k-2,j))) &
5693 IF( (config_flags%open_xe) .and. (ite == ide)) THEN
5695 DO j = j_start, j_end
5698 uw = 0.5*(fzm(k)*(ru(ite-1,k ,j)+ru(ite,k ,j)) + &
5699 fzp(k)*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) )
5702 tendency(i_end,k,j) = tendency(i_end,k,j) &
5704 ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + &
5706 fzm(k)*(ru(ite,k ,j)-ru(ite-1,k ,j)) + &
5707 fzp(k)*(ru(ite,k-1,j)-ru(ite-1,k-1,j))) &
5713 DO j = j_start, j_end
5715 uw = 0.5*( (2.-fzm(k-1))*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) &
5716 -fzp(k-1)*(ru(ite-1,k-2,j)+ru(ite,k-2,j)) )
5719 tendency(i_end,k,j) = tendency(i_end,k,j) &
5721 ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + &
5723 (2.-fzm(k-1))*(ru(ite,k-1,j)-ru(ite-1,k-1,j)) - &
5724 fzp(k-1)*(ru(ite,k-2,j)-ru(ite-1,k-2,j))) &
5731 IF( (config_flags%open_ys) .and. (jts == jds)) THEN
5733 DO i = i_start, i_end
5736 vw = 0.5*( fzm(k)*(rv(i,k ,jts)+rv(i,k ,jts+1)) + &
5737 fzp(k)*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) )
5740 tendency(i,k,jts) = tendency(i,k,jts) &
5742 vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + &
5744 fzm(k)*(rv(i,k ,jts+1)-rv(i,k ,jts))+ &
5745 fzp(k)*(rv(i,k-1,jts+1)-rv(i,k-1,jts))) &
5751 DO i = i_start, i_end
5752 vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) &
5753 -fzp(k-1)*(rv(i,k-2,jts)+rv(i,k-2,jts+1)) )
5756 tendency(i,k,jts) = tendency(i,k,jts) &
5758 vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + &
5760 (2.-fzm(k-1))*(rv(i,k-1,jts+1)-rv(i,k-1,jts))- &
5761 fzp(k-1)*(rv(i,k-2,jts+1)-rv(i,k-2,jts))) &
5767 IF( (config_flags%open_ye) .and. (jte == jde) ) THEN
5769 DO i = i_start, i_end
5772 vw = 0.5*( fzm(k)*(rv(i,k ,jte-1)+rv(i,k ,jte)) + &
5773 fzp(k)*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) )
5776 tendency(i,k,j_end) = tendency(i,k,j_end) &
5778 vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + &
5780 fzm(k)*(rv(i,k ,jte)-rv(i,k ,jte-1))+ &
5781 fzp(k)*(rv(i,k-1,jte)-rv(i,k-1,jte-1))) &
5787 DO i = i_start, i_end
5789 vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) &
5790 -fzp(k-1)*(rv(i,k-2,jte-1)+rv(i,k-2,jte)) )
5793 tendency(i,k,j_end) = tendency(i,k,j_end) &
5795 vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + &
5797 (2.-fzm(k-1))*(rv(i,k-1,jte)-rv(i,k-1,jte-1))- &
5798 fzp(k-1)*(rv(i,k-2,jte)-rv(i,k-2,jte-1))) &
5804 !-------------------- vertical advection
5805 ! ADT eqn 46 has 3rd term on RHS (dividing through by my) = - partial d/dz (w rho w /my)
5806 ! Here we have: - partial d/dz (w*rom) = - partial d/dz (w rho w / my)
5807 ! Therefore we don't need to make a correction for advect_w
5810 i_end = MIN(ite,ide-1)
5812 j_end = MIN(jte,jde-1)
5814 DO i = i_start, i_end
5819 vert_order_test : IF (vert_order == 6) THEN
5821 DO j = j_start, j_end
5824 DO i = i_start, i_end
5825 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5826 vflux(i,k) = vel*flux6( &
5827 w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), &
5828 w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel )
5832 DO i = i_start, i_end
5835 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5838 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5839 vflux(i,k) = vel*flux4( &
5840 w(i,k-2,j), w(i,k-1,j), &
5841 w(i,k ,j), w(i,k+1,j), -vel )
5844 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5845 vflux(i,k) = vel*flux4( &
5846 w(i,k-2,j), w(i,k-1,j), &
5847 w(i,k ,j), w(i,k+1,j), -vel )
5850 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5855 DO i = i_start, i_end
5856 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5860 ! pick up flux contribution for w at the lid. wcs, 13 march 2004
5862 DO i = i_start, i_end
5863 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5868 ELSE IF (vert_order == 5) THEN
5870 DO j = j_start, j_end
5873 DO i = i_start, i_end
5874 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5875 vflux(i,k) = vel*flux5( &
5876 w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), &
5877 w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel )
5881 DO i = i_start, i_end
5884 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5887 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5888 vflux(i,k) = vel*flux3( &
5889 w(i,k-2,j), w(i,k-1,j), &
5890 w(i,k ,j), w(i,k+1,j), -vel )
5892 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5893 vflux(i,k) = vel*flux3( &
5894 w(i,k-2,j), w(i,k-1,j), &
5895 w(i,k ,j), w(i,k+1,j), -vel )
5898 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5903 DO i = i_start, i_end
5904 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5908 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5910 DO i = i_start, i_end
5911 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5916 ELSE IF (vert_order == 4) THEN
5918 DO j = j_start, j_end
5921 DO i = i_start, i_end
5922 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5923 vflux(i,k) = vel*flux4( &
5924 w(i,k-2,j), w(i,k-1,j), &
5925 w(i,k ,j), w(i,k+1,j), -vel )
5929 DO i = i_start, i_end
5932 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5934 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5939 DO i = i_start, i_end
5940 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5944 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5946 DO i = i_start, i_end
5947 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5952 ELSE IF (vert_order == 3) THEN
5954 DO j = j_start, j_end
5957 DO i = i_start, i_end
5958 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5959 vflux(i,k) = vel*flux3( &
5960 w(i,k-2,j), w(i,k-1,j), &
5961 w(i,k ,j), w(i,k+1,j), -vel )
5965 DO i = i_start, i_end
5968 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5970 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5975 DO i = i_start, i_end
5976 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5980 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5982 DO i = i_start, i_end
5983 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5988 ELSE IF (vert_order == 2) THEN
5990 DO j = j_start, j_end
5992 DO i = i_start, i_end
5994 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5998 DO i = i_start, i_end
5999 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
6004 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
6006 DO i = i_start, i_end
6007 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
6014 WRITE (wrf_err_message ,*) ' advect_w, v_order not known ',vert_order
6015 CALL wrf_error_fatal ( wrf_err_message )
6017 ENDIF vert_order_test
6019 END SUBROUTINE advect_w
6021 !----------------------------------------------------------------
6023 SUBROUTINE advect_scalar_pd ( field, field_old, tendency, &
6027 msfux, msfuy, msfvx, msfvy, &
6030 rdx, rdy, rdzw, dt, &
6031 ids, ide, jds, jde, kds, kde, &
6032 ims, ime, jms, jme, kms, kme, &
6033 its, ite, jts, jte, kts, kte )
6035 ! this is a first cut at a positive definite advection option
6036 ! for scalars in WRF. This version is memory intensive ->
6037 ! we save 3d arrays of x, y and z both high and low order fluxes
6038 ! (six in all). Alternatively, we could sweep in a direction
6039 ! and lower the cost considerably.
6041 ! uses the Smolarkiewicz MWR 1989 approach, with addition of first-order
6044 ! WCS, 3 December 2002, 24 February 2003
6050 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
6052 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
6053 ims, ime, jms, jme, kms, kme, &
6054 its, ite, jts, jte, kts, kte
6056 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
6062 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, mub, mu_old
6063 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
6065 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
6072 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
6076 REAL , INTENT(IN ) :: rdx, &
6082 INTEGER :: i, j, k, itf, jtf, ktf
6083 INTEGER :: i_start, i_end, j_start, j_end
6084 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
6085 INTEGER :: jmin, jmax, jp, jm, imin, imax
6087 REAL :: mrdx, mrdy, ub, vb, uw, vw, mu
6089 ! storage for high and low order fluxes
6091 REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqx, fqy, fqz
6092 REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqxl, fqyl, fqzl
6094 INTEGER :: horz_order, vert_order
6096 LOGICAL :: degrade_xs, degrade_ys
6097 LOGICAL :: degrade_xe, degrade_ye
6099 INTEGER :: jp1, jp0, jtmp
6101 REAL :: flux_out, ph_low, scale
6102 REAL, PARAMETER :: eps=1.e-20
6105 ! definition of flux operators, 3rd, 4th, 5th or 6th order
6107 REAL :: flux3, flux4, flux5, flux6, flux_upwind
6108 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr
6110 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
6111 (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2)
6113 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
6114 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
6115 sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1))
6117 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
6118 (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2) &
6119 +(1./60.)*(q_ip2+q_im3)
6121 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
6122 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
6123 -sign(1.,ua)*(1./60.)*( &
6124 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )
6126 flux_upwind(q_im1, q_i, cr ) = 0.5*min( 1.0,(cr+abs(cr)))*q_im1 &
6127 +0.5*max(-1.0,(cr-abs(cr)))*q_i
6128 ! flux_upwind(q_im1, q_i, cr ) = 0.
6132 LOGICAL, PARAMETER :: pd_limit = .true.
6134 ! set order for the advection schemes
6136 ! write(6,*) ' in pd advection routine '
6138 ! Empty arrays just in case:
6139 IF (config_flags%polar) THEN
6149 horz_order = config_flags%h_sca_adv_order
6150 vert_order = config_flags%v_sca_adv_order
6152 ! determine boundary mods for flux operators
6153 ! We degrade the flux operators from 3rd/4th order
6154 ! to second order one gridpoint in from the boundaries for
6155 ! all boundary conditions except periodic and symmetry - these
6156 ! conditions have boundary zone data fill for correct application
6157 ! of the higher order flux stencils
6164 ! begin with horizontal flux divergence
6165 ! here is the choice of flux operators
6168 horizontal_order_test : IF( horz_order == 6 ) THEN
6170 IF( config_flags%periodic_x .or. &
6171 config_flags%symmetric_xs .or. &
6172 (its > ids+2) ) degrade_xs = .false.
6173 IF( config_flags%periodic_x .or. &
6174 config_flags%symmetric_xe .or. &
6175 (ite < ide-3) ) degrade_xe = .false.
6176 IF( config_flags%periodic_y .or. &
6177 config_flags%symmetric_ys .or. &
6178 (jts > jds+2) ) degrade_ys = .false.
6179 IF( config_flags%periodic_y .or. &
6180 config_flags%symmetric_ye .or. &
6181 (jte < jde-3) ) degrade_ye = .false.
6183 !--------------- y - advection first
6185 !-- y flux compute; these bounds are for periodic and sym b.c.
6189 i_end = MIN(ite,ide-1)+1
6191 j_end = MIN(jte,jde-1)+1
6195 !-- modify loop bounds if open or specified
6197 IF(degrade_xs) i_start = its
6198 IF(degrade_xe) i_end = MIN(ite,ide-1)
6201 j_start = MAX(jts,jds+1)
6206 j_end = MIN(jte,jde-2)
6210 ! compute fluxes, 6th order
6212 j_loop_y_flux_6 : DO j = j_start, j_end+1
6214 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6217 DO i = i_start, i_end
6219 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6220 mu = 0.5*(mut(i,j)+mut(i,j-1))
6223 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6225 fqy( i, k, j ) = vel*flux6( &
6226 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
6227 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
6229 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6234 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6237 DO i = i_start, i_end
6239 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6240 mu = 0.5*(mut(i,j)+mut(i,j-1))
6243 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6245 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6246 (field(i,k,j)+field(i,k,j-1))
6248 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6253 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
6256 DO i = i_start, i_end
6258 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6259 mu = 0.5*(mut(i,j)+mut(i,j-1))
6262 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6264 fqy( i, k, j ) = vel*flux4( &
6265 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
6266 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6271 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6274 DO i = i_start, i_end
6276 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6277 mu = 0.5*(mut(i,j)+mut(i,j-1))
6280 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6282 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6283 (field(i,k,j)+field(i,k,j-1))
6284 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6289 ELSE IF ( j == jde-2 ) THEN ! 4th order flux 2 in from north boundary
6292 DO i = i_start, i_end
6294 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6295 mu = 0.5*(mut(i,j)+mut(i,j-1))
6298 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6300 fqy( i, k, j) = vel*flux4( &
6301 field(i,k,j-2),field(i,k,j-1), &
6302 field(i,k,j),field(i,k,j+1),vel )
6303 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6310 ENDDO j_loop_y_flux_6
6314 !-- these bounds are for periodic and sym conditions
6317 i_end = MIN(ite,ide-1)+1
6322 j_end = MIN(jte,jde-1)+1
6324 !-- modify loop bounds for open and specified b.c
6326 IF(degrade_ys) j_start = jts
6327 IF(degrade_ye) j_end = MIN(jte,jde-1)
6330 i_start = MAX(ids+1,its)
6331 i_start_f = i_start+2
6335 i_end = MIN(ide-2,ite)
6341 DO j = j_start, j_end
6346 DO i = i_start_f, i_end_f
6348 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6349 mu = 0.5*(mut(i,j)+mut(i-1,j))
6352 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6354 fqx( i,k,j ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), &
6355 field(i-1,k,j), field(i ,k,j), &
6356 field(i+1,k,j), field(i+2,k,j), &
6358 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6363 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6365 IF( degrade_xs ) THEN
6367 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
6371 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6372 mu = 0.5*(mut(i,j)+mut(i-1,j))
6375 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6377 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6378 *(field(i,k,j)+field(i-1,k,j))
6380 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6387 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6388 mu = 0.5*(mut(i,j)+mut(i-1,j))
6391 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6392 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6393 field(i ,k,j), field(i+1,k,j), &
6395 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6401 IF( degrade_xe ) THEN
6403 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
6406 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6407 mu = 0.5*(mut(i,j)+mut(i-1,j))
6410 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6411 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6412 *(field(i,k,j)+field(i-1,k,j))
6413 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6421 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6422 mu = 0.5*(mut(i,j)+mut(i-1,j))
6425 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6426 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6427 field(i ,k,j), field(i+1,k,j), &
6429 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6435 ENDDO ! enddo for outer J loop
6437 !--- end of 6th order horizontal flux calculation
6439 ELSE IF( horz_order == 5 ) THEN
6441 IF( config_flags%periodic_x .or. &
6442 config_flags%symmetric_xs .or. &
6443 (its > ids+2) ) degrade_xs = .false.
6444 IF( config_flags%periodic_x .or. &
6445 config_flags%symmetric_xe .or. &
6446 (ite < ide-3) ) degrade_xe = .false.
6447 IF( config_flags%periodic_y .or. &
6448 config_flags%symmetric_ys .or. &
6449 (jts > jds+2) ) degrade_ys = .false.
6450 IF( config_flags%periodic_y .or. &
6451 config_flags%symmetric_ye .or. &
6452 (jte < jde-3) ) degrade_ye = .false.
6454 !--------------- y - advection first
6456 !-- y flux compute; these bounds are for periodic and sym b.c.
6460 i_end = MIN(ite,ide-1)+1
6462 j_end = MIN(jte,jde-1)+1
6466 !-- modify loop bounds if open or specified
6468 IF(degrade_xs) i_start = its
6469 IF(degrade_xe) i_end = MIN(ite,ide-1)
6472 j_start = MAX(jts,jds+1)
6477 j_end = MIN(jte,jde-2)
6481 ! compute fluxes, 5th order
6483 j_loop_y_flux_5 : DO j = j_start, j_end+1
6485 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6488 DO i = i_start, i_end
6490 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6491 mu = 0.5*(mut(i,j)+mut(i,j-1))
6494 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6496 fqy( i, k, j ) = vel*flux5( &
6497 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
6498 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
6500 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6505 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6508 DO i = i_start, i_end
6510 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6511 mu = 0.5*(mut(i,j)+mut(i,j-1))
6514 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6516 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6517 (field(i,k,j)+field(i,k,j-1))
6519 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6524 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
6527 DO i = i_start, i_end
6529 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6530 mu = 0.5*(mut(i,j)+mut(i,j-1))
6533 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6535 fqy( i, k, j ) = vel*flux3( &
6536 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
6537 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6542 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6545 DO i = i_start, i_end
6547 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6548 mu = 0.5*(mut(i,j)+mut(i,j-1))
6551 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6553 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6554 (field(i,k,j)+field(i,k,j-1))
6555 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6560 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
6563 DO i = i_start, i_end
6565 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6566 mu = 0.5*(mut(i,j)+mut(i,j-1))
6569 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6571 fqy( i, k, j) = vel*flux3( &
6572 field(i,k,j-2),field(i,k,j-1), &
6573 field(i,k,j),field(i,k,j+1),vel )
6574 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6581 ENDDO j_loop_y_flux_5
6585 !-- these bounds are for periodic and sym conditions
6588 i_end = MIN(ite,ide-1)+1
6593 j_end = MIN(jte,jde-1)+1
6595 !-- modify loop bounds for open and specified b.c
6597 IF(degrade_ys) j_start = jts
6598 IF(degrade_ye) j_end = MIN(jte,jde-1)
6601 i_start = MAX(ids+1,its)
6602 i_start_f = i_start+2
6606 i_end = MIN(ide-2,ite)
6612 DO j = j_start, j_end
6617 DO i = i_start_f, i_end_f
6619 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6620 mu = 0.5*(mut(i,j)+mut(i-1,j))
6623 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6625 fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), &
6626 field(i-1,k,j), field(i ,k,j), &
6627 field(i+1,k,j), field(i+2,k,j), &
6629 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6634 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6636 IF( degrade_xs ) THEN
6638 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
6642 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6643 mu = 0.5*(mut(i,j)+mut(i-1,j))
6646 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6648 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6649 *(field(i,k,j)+field(i-1,k,j))
6651 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6658 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6659 mu = 0.5*(mut(i,j)+mut(i-1,j))
6662 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6663 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
6664 field(i ,k,j), field(i+1,k,j), &
6666 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6672 IF( degrade_xe ) THEN
6674 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
6677 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6678 mu = 0.5*(mut(i,j)+mut(i-1,j))
6681 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6682 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6683 *(field(i,k,j)+field(i-1,k,j))
6684 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6692 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6693 mu = 0.5*(mut(i,j)+mut(i-1,j))
6696 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6697 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
6698 field(i ,k,j), field(i+1,k,j), &
6700 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6706 ENDDO ! enddo for outer J loop
6708 !--- end of 5th order horizontal flux calculation
6710 ELSE IF( horz_order == 4 ) THEN
6712 IF( config_flags%periodic_x .or. &
6713 config_flags%symmetric_xs .or. &
6714 (its > ids+1) ) degrade_xs = .false.
6715 IF( config_flags%periodic_x .or. &
6716 config_flags%symmetric_xe .or. &
6717 (ite < ide-2) ) degrade_xe = .false.
6718 IF( config_flags%periodic_y .or. &
6719 config_flags%symmetric_ys .or. &
6720 (jts > jds+1) ) degrade_ys = .false.
6721 IF( config_flags%periodic_y .or. &
6722 config_flags%symmetric_ye .or. &
6723 (jte < jde-2) ) degrade_ye = .false.
6725 !--------------- y - advection first
6727 !-- y flux compute; these bounds are for periodic and sym b.c.
6731 i_end = MIN(ite,ide-1)+1
6733 j_end = MIN(jte,jde-1)+1
6737 !-- modify loop bounds if open or specified
6739 IF(degrade_xs) i_start = its
6740 IF(degrade_xe) i_end = MIN(ite,ide-1)
6743 j_start = MAX(jts,jds+1)
6748 j_end = MIN(jte,jde-2)
6752 ! compute fluxes, 4th order
6754 j_loop_y_flux_4 : DO j = j_start, j_end+1
6756 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6759 DO i = i_start, i_end
6761 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6762 mu = 0.5*(mut(i,j)+mut(i,j-1))
6765 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6767 fqy( i, k, j ) = vel*flux4( field(i,k,j-2), field(i,k,j-1), &
6768 field(i,k,j ), field(i,k,j+1), vel )
6770 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6775 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6778 DO i = i_start, i_end
6780 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6781 mu = 0.5*(mut(i,j)+mut(i,j-1))
6784 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6786 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6787 (field(i,k,j)+field(i,k,j-1))
6789 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6794 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6797 DO i = i_start, i_end
6799 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6800 mu = 0.5*(mut(i,j)+mut(i,j-1))
6803 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6805 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6806 (field(i,k,j)+field(i,k,j-1))
6807 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6814 ENDDO j_loop_y_flux_4
6818 !-- these bounds are for periodic and sym conditions
6821 i_end = MIN(ite,ide-1)+1
6826 j_end = MIN(jte,jde-1)+1
6828 !-- modify loop bounds for open and specified b.c
6830 IF(degrade_ys) j_start = jts
6831 IF(degrade_ye) j_end = MIN(jte,jde-1)
6834 i_start = MAX(ids+1,its)
6835 i_start_f = i_start+1
6839 i_end = MIN(ide-2,ite)
6845 DO j = j_start, j_end
6850 DO i = i_start_f, i_end_f
6852 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6853 mu = 0.5*(mut(i,j)+mut(i-1,j))
6856 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6858 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6859 field(i ,k,j), field(i+1,k,j), vel )
6860 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6865 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6867 IF( degrade_xs ) THEN
6868 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
6872 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6873 mu = 0.5*(mut(i,j)+mut(i-1,j))
6876 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6878 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6879 *(field(i,k,j)+field(i-1,k,j))
6881 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6887 IF( degrade_xe ) THEN
6888 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
6891 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6892 mu = 0.5*(mut(i,j)+mut(i-1,j))
6895 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6896 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6897 *(field(i,k,j)+field(i-1,k,j))
6898 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6904 ENDDO ! enddo for outer J loop
6906 !--- end of 4th order horizontal flux calculation
6908 ELSE IF( horz_order == 3 ) THEN
6910 IF( config_flags%periodic_x .or. &
6911 config_flags%symmetric_xs .or. &
6912 (its > ids+1) ) degrade_xs = .false.
6913 IF( config_flags%periodic_x .or. &
6914 config_flags%symmetric_xe .or. &
6915 (ite < ide-2) ) degrade_xe = .false.
6916 IF( config_flags%periodic_y .or. &
6917 config_flags%symmetric_ys .or. &
6918 (jts > jds+1) ) degrade_ys = .false.
6919 IF( config_flags%periodic_y .or. &
6920 config_flags%symmetric_ye .or. &
6921 (jte < jde-2) ) degrade_ye = .false.
6923 !--------------- y - advection first
6925 !-- y flux compute; these bounds are for periodic and sym b.c.
6929 i_end = MIN(ite,ide-1)+1
6931 j_end = MIN(jte,jde-1)+1
6935 !-- modify loop bounds if open or specified
6937 IF(degrade_xs) i_start = its
6938 IF(degrade_xe) i_end = MIN(ite,ide-1)
6941 j_start = MAX(jts,jds+1)
6946 j_end = MIN(jte,jde-2)
6950 ! compute fluxes, 3rd order
6952 j_loop_y_flux_3 : DO j = j_start, j_end+1
6954 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6957 DO i = i_start, i_end
6959 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6960 mu = 0.5*(mut(i,j)+mut(i,j-1))
6963 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6965 fqy( i, k, j ) = vel*flux3( field(i,k,j-2), field(i,k,j-1), &
6966 field(i,k,j ), field(i,k,j+1), vel )
6968 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6973 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6976 DO i = i_start, i_end
6978 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6979 mu = 0.5*(mut(i,j)+mut(i,j-1))
6982 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6984 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6985 (field(i,k,j)+field(i,k,j-1))
6987 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6992 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6995 DO i = i_start, i_end
6997 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6998 mu = 0.5*(mut(i,j)+mut(i,j-1))
7001 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
7003 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
7004 (field(i,k,j)+field(i,k,j-1))
7005 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7012 ENDDO j_loop_y_flux_3
7016 !-- these bounds are for periodic and sym conditions
7019 i_end = MIN(ite,ide-1)+1
7024 j_end = MIN(jte,jde-1)+1
7026 !-- modify loop bounds for open and specified b.c
7028 IF(degrade_ys) j_start = jts
7029 IF(degrade_ye) j_end = MIN(jte,jde-1)
7032 i_start = MAX(ids+1,its)
7033 i_start_f = i_start+1
7037 i_end = MIN(ide-2,ite)
7043 DO j = j_start, j_end
7048 DO i = i_start_f, i_end_f
7050 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7051 mu = 0.5*(mut(i,j)+mut(i-1,j))
7054 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7056 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
7057 field(i ,k,j), field(i+1,k,j), vel )
7058 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7063 ! lower order fluxes close to boundaries (if not periodic or symmetric)
7065 IF( degrade_xs ) THEN
7067 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
7071 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7072 mu = 0.5*(mut(i,j)+mut(i-1,j))
7075 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7077 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
7078 *(field(i,k,j)+field(i-1,k,j))
7080 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7086 IF( degrade_xe ) THEN
7087 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
7090 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7091 mu = 0.5*(mut(i,j)+mut(i-1,j))
7094 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7095 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
7096 *(field(i,k,j)+field(i-1,k,j))
7097 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7103 ENDDO ! enddo for outer J loop
7105 !--- end of 3rd order horizontal flux calculation
7108 ELSE IF( horz_order == 2 ) THEN
7110 IF( config_flags%periodic_x .or. &
7111 config_flags%symmetric_xs .or. &
7112 (its > ids) ) degrade_xs = .false.
7113 IF( config_flags%periodic_x .or. &
7114 config_flags%symmetric_xe .or. &
7115 (ite < ide-1) ) degrade_xe = .false.
7116 IF( config_flags%periodic_y .or. &
7117 config_flags%symmetric_ys .or. &
7118 (jts > jds) ) degrade_ys = .false.
7119 IF( config_flags%periodic_y .or. &
7120 config_flags%symmetric_ye .or. &
7121 (jte < jde-1) ) degrade_ye = .false.
7123 !-- y flux compute; these bounds are for periodic and sym b.c.
7127 i_end = MIN(ite,ide-1)+1
7129 j_end = MIN(jte,jde-1)+1
7131 !-- modify loop bounds if open or specified
7133 IF(degrade_xs) i_start = its
7134 IF(degrade_xe) i_end = MIN(ite,ide-1)
7135 IF(degrade_ys) j_start = MAX(jts,jds+1)
7136 IF(degrade_ye) j_end = MIN(jte,jde-2)
7138 ! compute fluxes, 2nd order, y flux
7140 DO j = j_start, j_end+1
7142 DO i = i_start, i_end
7143 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
7144 mu = 0.5*(mut(i,j)+mut(i,j-1))
7147 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
7149 fqy(i,k, j) = 0.5*rv(i,k,j)* &
7150 (field(i,k,j)+field(i,k,j-1))
7152 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7159 DO j = j_start, j_end
7161 DO i = i_start, i_end+1
7162 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7163 mu = 0.5*(mut(i,j)+mut(i-1,j))
7166 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7167 fqx( i,k,j ) = 0.5*ru(i,k,j)* &
7168 (field(i,k,j)+field(i-1,k,j))
7170 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7175 !--- end of 2nd order horizontal flux calculation
7179 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_pd, h_order not known ',horz_order
7180 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
7182 ENDIF horizontal_order_test
7184 ! pick up the rest of the horizontal radiation boundary conditions.
7185 ! (these are the computations that don't require 'cb'.
7186 ! first, set to index ranges
7189 i_end = MIN(ite,ide-1)
7191 j_end = MIN(jte,jde-1)
7193 ! compute x (u) conditions for v, w, or scalar
7195 IF( (config_flags%open_xs) .and. (its == ids) ) THEN
7197 DO j = j_start, j_end
7199 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
7200 tendency(its,k,j) = tendency(its,k,j) &
7202 ub*( field_old(its+1,k,j) &
7203 - field_old(its ,k,j) ) + &
7204 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) &
7211 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
7213 DO j = j_start, j_end
7215 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
7216 tendency(i_end,k,j) = tendency(i_end,k,j) &
7218 ub*( field_old(i_end ,k,j) &
7219 - field_old(i_end-1,k,j) ) + &
7220 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) &
7227 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
7229 DO i = i_start, i_end
7231 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
7232 tendency(i,k,jts) = tendency(i,k,jts) &
7234 vb*( field_old(i,k,jts+1) &
7235 - field_old(i,k,jts ) ) + &
7236 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) &
7243 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
7245 DO i = i_start, i_end
7247 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
7248 tendency(i,k,j_end) = tendency(i,k,j_end) &
7250 vb*( field_old(i,k,j_end ) &
7251 - field_old(i,k,j_end-1) ) + &
7252 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) &
7259 IF( (config_flags%polar) .and. (jts == jds) ) THEN
7261 ! Assuming rv(i,k,jds) = 0.
7262 DO i = i_start, i_end
7264 vb = MIN( 0.5*rv(i,k,jts+1), 0. )
7265 tendency(i,k,jts) = tendency(i,k,jts) &
7267 vb*( field_old(i,k,jts+1) &
7268 - field_old(i,k,jts ) ) + &
7269 field(i,k,jts)*rv(i,k,jts+1) &
7276 IF( (config_flags%polar) .and. (jte == jde)) THEN
7278 ! Assuming rv(i,k,jde) = 0.
7279 DO i = i_start, i_end
7281 vb = MAX( 0.5*rv(i,k,jte-1), 0. )
7282 tendency(i,k,j_end) = tendency(i,k,j_end) &
7284 vb*( field_old(i,k,j_end ) &
7285 - field_old(i,k,j_end-1) ) + &
7286 field(i,k,j_end)*(-rv(i,k,jte-1)) &
7293 !-------------------- vertical advection
7295 !-- loop bounds for periodic or sym conditions
7298 i_end = MIN(ite,ide-1)+1
7300 j_end = MIN(jte,jde-1)+1
7302 !-- loop bounds for open or specified conditions
7304 IF(degrade_xs) i_start = its
7305 IF(degrade_xe) i_end = MIN(ite,ide-1)
7306 IF(degrade_ys) j_start = jts
7307 IF(degrade_ye) j_end = MIN(jte,jde-1)
7309 vert_order_test : IF (vert_order == 6) THEN
7311 DO j = j_start, j_end
7313 DO i = i_start, i_end
7321 DO i = i_start, i_end
7322 dz = 2./(rdzw(k)+rdzw(k-1))
7323 mu = 0.5*(mut(i,j)+mut(i,j))
7326 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7328 fqz(i,k,j) = vel*flux6( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
7329 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
7330 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7334 DO i = i_start, i_end
7337 dz = 2./(rdzw(k)+rdzw(k-1))
7338 mu = 0.5*(mut(i,j)+mut(i,j))
7341 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7342 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7343 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7346 dz = 2./(rdzw(k)+rdzw(k-1))
7347 mu = 0.5*(mut(i,j)+mut(i,j))
7350 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7352 fqz(i,k,j) = vel*flux4( &
7353 field(i,k-2,j), field(i,k-1,j), &
7354 field(i,k ,j), field(i,k+1,j), -vel )
7355 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7358 dz = 2./(rdzw(k)+rdzw(k-1))
7359 mu = 0.5*(mut(i,j)+mut(i,j))
7362 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7364 fqz(i,k,j) = vel*flux4( &
7365 field(i,k-2,j), field(i,k-1,j), &
7366 field(i,k ,j), field(i,k+1,j), -vel )
7367 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7370 dz = 2./(rdzw(k)+rdzw(k-1))
7371 mu = 0.5*(mut(i,j)+mut(i,j))
7374 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7375 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7376 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7382 ELSE IF (vert_order == 5) THEN
7384 DO j = j_start, j_end
7386 DO i = i_start, i_end
7394 DO i = i_start, i_end
7395 dz = 2./(rdzw(k)+rdzw(k-1))
7396 mu = 0.5*(mut(i,j)+mut(i,j))
7399 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7401 fqz(i,k,j) = vel*flux5( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
7402 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
7403 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7407 DO i = i_start, i_end
7410 dz = 2./(rdzw(k)+rdzw(k-1))
7411 mu = 0.5*(mut(i,j)+mut(i,j))
7414 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7415 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7416 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7419 dz = 2./(rdzw(k)+rdzw(k-1))
7420 mu = 0.5*(mut(i,j)+mut(i,j))
7423 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7425 fqz(i,k,j) = vel*flux3( &
7426 field(i,k-2,j), field(i,k-1,j), &
7427 field(i,k ,j), field(i,k+1,j), -vel )
7428 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7431 dz = 2./(rdzw(k)+rdzw(k-1))
7432 mu = 0.5*(mut(i,j)+mut(i,j))
7435 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7437 fqz(i,k,j) = vel*flux3( &
7438 field(i,k-2,j), field(i,k-1,j), &
7439 field(i,k ,j), field(i,k+1,j), -vel )
7440 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7443 dz = 2./(rdzw(k)+rdzw(k-1))
7444 mu = 0.5*(mut(i,j)+mut(i,j))
7447 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7448 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7449 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7455 ELSE IF (vert_order == 4) THEN
7457 DO j = j_start, j_end
7459 DO i = i_start, i_end
7467 DO i = i_start, i_end
7469 dz = 2./(rdzw(k)+rdzw(k-1))
7470 mu = 0.5*(mut(i,j)+mut(i,j))
7473 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7475 fqz(i,k,j) = vel*flux4( &
7476 field(i,k-2,j), field(i,k-1,j), &
7477 field(i,k ,j), field(i,k+1,j), -vel )
7478 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7482 DO i = i_start, i_end
7485 dz = 2./(rdzw(k)+rdzw(k-1))
7486 mu = 0.5*(mut(i,j)+mut(i,j))
7489 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7490 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7491 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7494 dz = 2./(rdzw(k)+rdzw(k-1))
7495 mu = 0.5*(mut(i,j)+mut(i,j))
7498 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7499 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7500 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7506 ELSE IF (vert_order == 3) THEN
7508 DO j = j_start, j_end
7510 DO i = i_start, i_end
7518 DO i = i_start, i_end
7520 dz = 2./(rdzw(k)+rdzw(k-1))
7521 mu = 0.5*(mut(i,j)+mut(i,j))
7524 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7526 fqz(i,k,j) = vel*flux3( &
7527 field(i,k-2,j), field(i,k-1,j), &
7528 field(i,k ,j), field(i,k+1,j), -vel )
7529 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7533 DO i = i_start, i_end
7536 dz = 2./(rdzw(k)+rdzw(k-1))
7537 mu = 0.5*(mut(i,j)+mut(i,j))
7540 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7541 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7542 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7545 dz = 2./(rdzw(k)+rdzw(k-1))
7546 mu = 0.5*(mut(i,j)+mut(i,j))
7549 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7550 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7551 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7557 ELSE IF (vert_order == 2) THEN
7559 DO j = j_start, j_end
7561 DO i = i_start, i_end
7569 DO i = i_start, i_end
7571 dz = 2./(rdzw(k)+rdzw(k-1))
7572 mu = 0.5*(mut(i,j)+mut(i,j))
7575 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7576 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7577 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7586 WRITE (wrf_err_message,*) ' advect_scalar_pd, v_order not known ',vert_order
7587 CALL wrf_error_fatal ( wrf_err_message )
7589 ENDIF vert_order_test
7593 ! positive definite filter
7596 i_end = MIN(ite,ide-1)+1
7598 j_end = MIN(jte,jde-1)+1
7600 !-- loop bounds for open or specified conditions
7602 IF(degrade_xs) i_start = its
7603 IF(degrade_xe) i_end = MIN(ite,ide-1)
7604 IF(degrade_ys) j_start = jts
7605 IF(degrade_ye) j_end = MIN(jte,jde-1)
7607 IF(config_flags%specified .or. config_flags%nested) THEN
7608 IF (degrade_xs) i_start = MAX(its,ids+1)
7609 IF (degrade_xe) i_end = MIN(ite,ide-2)
7610 IF (degrade_ys) j_start = MAX(jts,jds+1)
7611 IF (degrade_ye) j_end = MIN(jte,jde-2)
7614 IF(config_flags%open_xs) THEN
7615 IF (degrade_xs) i_start = MAX(its,ids+1)
7617 IF(config_flags%open_xe) THEN
7618 IF (degrade_xe) i_end = MIN(ite,ide-2)
7620 IF(config_flags%open_ys) THEN
7621 IF (degrade_ys) j_start = MAX(jts,jds+1)
7623 IF(config_flags%open_ye) THEN
7624 IF (degrade_ye) j_end = MIN(jte,jde-2)
7627 ! We don't want to change j_start and j_end
7628 ! for polar BC's since we want to calculate
7629 ! fluxes for directions other than y at the
7632 !-- here is the limiter...
7638 ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) &
7639 - dt*( msftx(i,j)*msfty(i,j)*( &
7640 rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) + &
7641 rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) &
7642 +msfty(i,j)*rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) )
7644 flux_out = dt*( (msftx(i,j)*msfty(i,j))*( &
7645 rdx*( max(0.,fqx (i+1,k,j)) &
7646 -min(0.,fqx (i ,k,j)) ) &
7647 +rdy*( max(0.,fqy (i,k,j+1)) &
7648 -min(0.,fqy (i,k,j )) ) ) &
7649 +msfty(i,j)*rdzw(k)*( min(0.,fqz (i,k+1,j)) &
7650 -max(0.,fqz (i,k ,j)) ) )
7652 IF( flux_out .gt. ph_low ) THEN
7654 scale = max(0.,ph_low/(flux_out+eps))
7655 IF( fqx (i+1,k,j) .gt. 0.) fqx(i+1,k,j) = scale*fqx(i+1,k,j)
7656 IF( fqx (i ,k,j) .lt. 0.) fqx(i ,k,j) = scale*fqx(i ,k,j)
7657 IF( fqy (i,k,j+1) .gt. 0.) fqy(i,k,j+1) = scale*fqy(i,k,j+1)
7658 IF( fqy (i,k,j ) .lt. 0.) fqy(i,k,j ) = scale*fqy(i,k,j )
7659 ! note: z flux is opposite sign in mass coordinate because
7660 ! vertical coordinate decreases with increasing k
7661 IF( fqz (i,k+1,j) .lt. 0.) fqz(i,k+1,j) = scale*fqz(i,k+1,j)
7662 IF( fqz (i,k ,j) .gt. 0.) fqz(i,k ,j) = scale*fqz(i,k ,j)
7672 ! add in the pd-limited flux divergence
7675 i_end = MIN(ite,ide-1)
7677 j_end = MIN(jte,jde-1)
7679 DO j = j_start, j_end
7681 DO i = i_start, i_end
7683 tendency (i,k,j) = tendency(i,k,j) &
7684 -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) &
7685 +fqzl(i,k+1,j)-fqzl(i,k,j))
7693 IF(degrade_xs) i_start = i_start + 1
7694 IF(degrade_xe) i_end = i_end - 1
7696 DO j = j_start, j_end
7698 DO i = i_start, i_end
7700 ! Un-"canceled" map scale factor, ADT Eq. 48
7701 tendency (i,k,j) = tendency(i,k,j) &
7702 - msftx(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) &
7703 +fqxl(i+1,k,j)-fqxl(i,k,j)) )
7712 i_end = MIN(ite,ide-1)
7713 IF(degrade_ys) j_start = j_start + 1
7714 IF(degrade_ye) j_end = j_end - 1
7716 DO j = j_start, j_end
7718 DO i = i_start, i_end
7720 ! Un-"canceled" map scale factor, ADT Eq. 48
7721 ! It is correct to use msftx (and not msfty), per W. Skamarock, 20080606
7722 tendency (i,k,j) = tendency(i,k,j) &
7723 - msftx(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) &
7724 +fqyl(i,k,j+1)-fqyl(i,k,j)) )
7730 END SUBROUTINE advect_scalar_pd
7732 !----------------------------------------------------------------
7734 END MODULE module_advect_em