1 !WRF:MODEL_LAYER:DYNAMICS
3 MODULE module_advect_em
6 USE module_model_constants
11 !-------------------------------------------------------------------------------
13 SUBROUTINE advect_u ( u, u_old, tendency, &
15 mut, time_step, config_flags, &
16 msfux, msfuy, msfvx, msfvy, &
20 ids, ide, jds, jde, kds, kde, &
21 ims, ime, jms, jme, kms, kme, &
22 its, ite, jts, jte, kts, kte )
28 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
30 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
31 ims, ime, jms, jme, kms, kme, &
32 its, ite, jts, jte, kts, kte
34 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, &
40 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
41 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
43 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
50 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
54 REAL , INTENT(IN ) :: rdx, &
56 INTEGER , INTENT(IN ) :: time_step
60 INTEGER :: i, j, k, itf, jtf, ktf
61 INTEGER :: i_start, i_end, j_start, j_end
62 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
63 INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
64 INTEGER :: jp1, jp0, jtmp
66 INTEGER :: horz_order, vert_order
68 REAL :: mrdx, mrdy, ub, vb, uw, vw, dvm, dvp
69 REAL , DIMENSION(its:ite, kts:kte) :: vflux
72 REAL, DIMENSION( its-1:ite+1, kts:kte ) :: fqx
73 REAL, DIMENSION( its:ite, kts:kte, 2) :: fqy
75 LOGICAL :: degrade_xs, degrade_ys
76 LOGICAL :: degrade_xe, degrade_ye
78 ! definition of flux operators, 3rd, 4th, 5th or 6th order
80 REAL :: flux3, flux4, flux5, flux6
81 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
83 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
84 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
86 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
87 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
88 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
90 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
91 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
94 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
95 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
96 -sign(1,time_step)*sign(1.,ua)*( &
97 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
103 if(config_flags%specified .or. config_flags%nested) specified = .true.
105 ! set order for vertical and horzontal flux operators
107 horz_order = config_flags%h_mom_adv_order
108 vert_order = config_flags%v_mom_adv_order
112 ! begin with horizontal flux divergence
114 horizontal_order_test : IF( horz_order == 6 ) THEN
116 ! determine boundary mods for flux operators
117 ! We degrade the flux operators from 3rd/4th order
118 ! to second order one gridpoint in from the boundaries for
119 ! all boundary conditions except periodic and symmetry - these
120 ! conditions have boundary zone data fill for correct application
121 ! of the higher order flux stencils
128 IF( config_flags%periodic_x .or. &
129 config_flags%symmetric_xs .or. &
130 (its > ids+3) ) degrade_xs = .false.
131 IF( config_flags%periodic_x .or. &
132 config_flags%symmetric_xe .or. &
133 (ite < ide-2) ) degrade_xe = .false.
134 IF( config_flags%periodic_y .or. &
135 config_flags%symmetric_ys .or. &
136 (jts > jds+3) ) degrade_ys = .false.
137 IF( config_flags%periodic_y .or. &
138 config_flags%symmetric_ye .or. &
139 (jte < jde-4) ) degrade_ye = .false.
141 !--------------- y - advection first
145 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
146 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
147 IF ( config_flags%periodic_x ) i_start = its
148 IF ( config_flags%periodic_x ) i_end = ite
151 j_end = MIN(jte,jde-1)
153 ! higher order flux has a 5 or 7 point stencil, so compute
154 ! bounds so we can switch to second order flux close to the boundary
160 j_start = MAX(jts,jds+1)
165 j_end = MIN(jte,jde-2)
169 IF(config_flags%polar) j_end = MIN(jte,jde-1)
171 ! compute fluxes, 5th or 6th order
176 j_loop_y_flux_6 : DO j = j_start, j_end+1
178 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
181 DO i = i_start, i_end
182 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
183 fqy( i, k, jp1 ) = vel*flux6( &
184 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), &
185 u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel )
189 ! we must be close to some boundary where we need to reduce the order of the stencil
191 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
194 DO i = i_start, i_end
195 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
196 *(u(i,k,j)+u(i,k,j-1))
200 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
203 DO i = i_start, i_end
204 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
205 fqy( i, k, jp1 ) = vel*flux4( &
206 u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel )
210 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
213 DO i = i_start, i_end
214 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
215 *(u(i,k,j)+u(i,k,j-1))
219 ELSE IF ( j == jde-2 ) THEN ! 3rd order flux 2 in from north boundary
222 DO i = i_start, i_end
223 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
224 fqy( i, k, jp1 ) = vel*flux4( &
225 u(i,k,j-2),u(i,k,j-1), &
226 u(i,k,j),u(i,k,j+1),vel )
232 ! y flux-divergence into tendency
234 ! (j > j_start) will miss the u(,,jds) tendency
235 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
237 DO i = i_start, i_end
238 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
239 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
242 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
243 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
245 DO i = i_start, i_end
246 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
247 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
255 DO i = i_start, i_end
256 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
257 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
270 ENDDO j_loop_y_flux_6
272 ! next, x - flux divergence
278 j_end = MIN(jte,jde-1)
280 ! higher order flux has a 5 or 7 point stencil, so compute
281 ! bounds so we can switch to second order flux close to the boundary
287 i_start = MAX(ids+1,its)
292 i_end = MIN(ide-1,ite)
298 DO j = j_start, j_end
300 ! 5th or 6th order flux
303 DO i = i_start_f, i_end_f
304 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
305 fqx( i,k ) = vel*flux6( u(i-3,k,j), u(i-2,k,j), &
306 u(i-1,k,j), u(i ,k,j), &
307 u(i+1,k,j), u(i+2,k,j), &
312 ! lower order fluxes close to boundaries (if not periodic or symmetric)
313 ! specified uses upstream normal wind at boundaries
315 IF( degrade_xs ) THEN
317 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
321 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
322 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
329 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
330 fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
331 u(i ,k,j), u(i+1,k,j), &
337 IF( degrade_xe ) THEN
339 IF( i_end == ide-1 ) THEN ! second order flux next to the boundary
343 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
344 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
351 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
352 fqx( i,k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
353 u(i ,k,j), u(i+1,k,j), &
359 ! x flux-divergence into tendency
362 DO i = i_start, i_end
363 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
364 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
370 ELSE IF( horz_order == 5 ) THEN
372 ! 5th order horizontal flux calculation
373 ! This code is EXACTLY the same as the 6th order code
374 ! EXCEPT the 5th order and 3rd operators are used in
375 ! place of the 6th and 4th order operators
377 ! determine boundary mods for flux operators
378 ! We degrade the flux operators from 3rd/4th order
379 ! to second order one gridpoint in from the boundaries for
380 ! all boundary conditions except periodic and symmetry - these
381 ! conditions have boundary zone data fill for correct application
382 ! of the higher order flux stencils
389 IF( config_flags%periodic_x .or. &
390 config_flags%symmetric_xs .or. &
391 (its > ids+3) ) degrade_xs = .false.
392 IF( config_flags%periodic_x .or. &
393 config_flags%symmetric_xe .or. &
394 (ite < ide-2) ) degrade_xe = .false.
395 IF( config_flags%periodic_y .or. &
396 config_flags%symmetric_ys .or. &
397 (jts > jds+3) ) degrade_ys = .false.
398 IF( config_flags%periodic_y .or. &
399 config_flags%symmetric_ye .or. &
400 (jte < jde-4) ) degrade_ye = .false.
402 !--------------- y - advection first
406 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
407 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
408 IF ( config_flags%periodic_x ) i_start = its
409 IF ( config_flags%periodic_x ) i_end = ite
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 j_start = MAX(jts,jds+1)
426 j_end = MIN(jte,jde-2)
430 IF(config_flags%polar) j_end = MIN(jte,jde-1)
432 ! compute fluxes, 5th or 6th order
437 j_loop_y_flux_5 : DO j = j_start, j_end+1
439 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
442 DO i = i_start, i_end
443 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
444 fqy( i, k, jp1 ) = vel*flux5( &
445 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), &
446 u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel )
450 ! we must be close to some boundary where we need to reduce the order of the stencil
452 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
455 DO i = i_start, i_end
456 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
457 *(u(i,k,j)+u(i,k,j-1))
461 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
464 DO i = i_start, i_end
465 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
466 fqy( i, k, jp1 ) = vel*flux3( &
467 u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel )
471 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
474 DO i = i_start, i_end
475 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
476 *(u(i,k,j)+u(i,k,j-1))
480 ELSE IF ( j == jde-2 ) THEN ! 3rd order flux 2 in from north boundary
483 DO i = i_start, i_end
484 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
485 fqy( i, k, jp1 ) = vel*flux3( &
486 u(i,k,j-2),u(i,k,j-1), &
487 u(i,k,j),u(i,k,j+1),vel )
493 ! y flux-divergence into tendency
495 ! (j > j_start) will miss the u(,,jds) tendency
496 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
498 DO i = i_start, i_end
499 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
500 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
503 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
504 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
506 DO i = i_start, i_end
507 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
508 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
516 DO i = i_start, i_end
517 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
518 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
531 ENDDO j_loop_y_flux_5
533 ! next, x - flux divergence
539 j_end = MIN(jte,jde-1)
541 ! higher order flux has a 5 or 7 point stencil, so compute
542 ! bounds so we can switch to second order flux close to the boundary
548 i_start = MAX(ids+1,its)
553 i_end = MIN(ide-1,ite)
559 DO j = j_start, j_end
561 ! 5th or 6th order flux
564 DO i = i_start_f, i_end_f
565 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
566 fqx( i,k ) = vel*flux5( u(i-3,k,j), u(i-2,k,j), &
567 u(i-1,k,j), u(i ,k,j), &
568 u(i+1,k,j), u(i+2,k,j), &
573 ! lower order fluxes close to boundaries (if not periodic or symmetric)
574 ! specified uses upstream normal wind at boundaries
576 IF( degrade_xs ) THEN
578 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
582 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
583 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
590 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
591 fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
592 u(i ,k,j), u(i+1,k,j), &
598 IF( degrade_xe ) THEN
600 IF( i_end == ide-1 ) THEN ! second order flux next to the boundary
604 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
605 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
612 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
613 fqx( i,k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
614 u(i ,k,j), u(i+1,k,j), &
620 ! x flux-divergence into tendency
623 DO i = i_start, i_end
624 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
625 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
631 ELSE IF( horz_order == 4 ) THEN
633 ! determine boundary mods for flux operators
634 ! We degrade the flux operators from 3rd/4th order
635 ! to second order one gridpoint in from the boundaries for
636 ! all boundary conditions except periodic and symmetry - these
637 ! conditions have boundary zone data fill for correct application
638 ! of the higher order flux stencils
645 IF( config_flags%periodic_x .or. &
646 config_flags%symmetric_xs .or. &
647 (its > ids+2) ) degrade_xs = .false.
648 IF( config_flags%periodic_x .or. &
649 config_flags%symmetric_xe .or. &
650 (ite < ide-1) ) degrade_xe = .false.
651 IF( config_flags%periodic_y .or. &
652 config_flags%symmetric_ys .or. &
653 (jts > jds+2) ) degrade_ys = .false.
654 IF( config_flags%periodic_y .or. &
655 config_flags%symmetric_ye .or. &
656 (jte < jde-3) ) degrade_ye = .false.
658 !--------------- x - advection first
663 j_end = MIN(jte,jde-1)
665 ! 3rd or 4th order flux has a 5 point stencil, so compute
666 ! bounds so we can switch to second order flux close to the boundary
673 i_start_f = i_start+1
683 DO j = j_start, j_end
686 DO i = i_start_f, i_end_f
687 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
688 fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), &
689 u(i ,k,j), u(i+1,k,j), vel )
693 ! second order flux close to boundaries (if not periodic or symmetric)
694 ! specified uses upstream normal wind at boundaries
696 IF( degrade_xs ) THEN
700 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
701 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
706 IF( degrade_xe ) THEN
710 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
711 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
716 ! x flux-divergence into tendency
719 DO i = i_start, i_end
720 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
721 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
731 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
732 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
733 IF ( config_flags%periodic_x ) i_start = its
734 IF ( config_flags%periodic_x ) i_end = ite
737 j_end = MIN(jte,jde-1)
739 ! 3rd or 4th order flux has a 5 point stencil, so compute
740 ! bounds so we can switch to second order flux close to the boundary
745 !CJM these may not work with tiling because they define j_start and end in terms of domain dim
748 j_start_f = j_start+1
756 IF(config_flags%polar) j_end = MIN(jte,jde-1)
758 ! j flux loop for v flux of u momentum
763 DO j = j_start, j_end+1
765 IF ( (j < j_start_f) .and. degrade_ys) THEN
767 DO i = i_start, i_end
768 fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) &
769 *(u(i,k,j_start)+u(i,k,j_start-1))
772 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
774 DO i = i_start, i_end
775 ! Assumes j>j_end_f is ONLY j_end+1 ...
776 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) &
777 ! *(u(i,k,j_end+1)+u(i,k,j_end))
778 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
779 *(u(i,k,j)+u(i,k,j-1))
783 ! 3rd or 4th order flux
785 DO i = i_start, i_end
786 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
787 fqy( i, k, jp1 ) = vel*flux4( u(i,k,j-2), u(i,k,j-1), &
788 u(i,k,j ), u(i,k,j+1), &
795 ! y flux-divergence into tendency
797 ! (j > j_start) will miss the u(,,jds) tendency
798 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
800 DO i = i_start, i_end
801 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
802 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
805 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
806 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
808 DO i = i_start, i_end
809 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
810 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
815 IF (j > j_start) THEN
818 DO i = i_start, i_end
819 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
820 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
834 ELSE IF ( horz_order == 3 ) THEN
836 ! As with the 5th and 6th order flux chioces, the 3rd and 4th order
837 ! code is EXACTLY the same EXCEPT for the flux operator.
839 ! determine boundary mods for flux operators
840 ! We degrade the flux operators from 3rd/4th order
841 ! to second order one gridpoint in from the boundaries for
842 ! all boundary conditions except periodic and symmetry - these
843 ! conditions have boundary zone data fill for correct application
844 ! of the higher order flux stencils
851 IF( config_flags%periodic_x .or. &
852 config_flags%symmetric_xs .or. &
853 (its > ids+2) ) degrade_xs = .false.
854 IF( config_flags%periodic_x .or. &
855 config_flags%symmetric_xe .or. &
856 (ite < ide-1) ) degrade_xe = .false.
857 IF( config_flags%periodic_y .or. &
858 config_flags%symmetric_ys .or. &
859 (jts > jds+2) ) degrade_ys = .false.
860 IF( config_flags%periodic_y .or. &
861 config_flags%symmetric_ye .or. &
862 (jte < jde-3) ) degrade_ye = .false.
864 !--------------- x - advection first
869 j_end = MIN(jte,jde-1)
871 ! 3rd or 4th order flux has a 5 point stencil, so compute
872 ! bounds so we can switch to second order flux close to the boundary
879 i_start_f = i_start+1
889 DO j = j_start, j_end
892 DO i = i_start_f, i_end_f
893 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j))
894 fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), &
895 u(i ,k,j), u(i+1,k,j), vel )
899 ! second order flux close to boundaries (if not periodic or symmetric)
900 ! specified uses upstream normal wind at boundaries
902 IF( degrade_xs ) THEN
906 IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j)
907 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
912 IF( degrade_xe ) THEN
916 IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j)
917 fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) &
922 ! x flux-divergence into tendency
925 DO i = i_start, i_end
926 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
927 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
936 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
937 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
938 IF ( config_flags%periodic_x ) i_start = its
939 IF ( config_flags%periodic_x ) i_end = ite
942 j_end = MIN(jte,jde-1)
944 ! 3rd or 4th order flux has a 5 point stencil, so compute
945 ! bounds so we can switch to second order flux close to the boundary
950 !CJM these may not work with tiling because they define j_start and end in terms of domain dim
953 j_start_f = j_start+1
961 IF(config_flags%polar) j_end = MIN(jte,jde-1)
963 ! j flux loop for v flux of u momentum
968 DO j = j_start, j_end+1
970 IF ( (j < j_start_f) .and. degrade_ys) THEN
972 DO i = i_start, i_end
973 fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) &
974 *(u(i,k,j_start)+u(i,k,j_start-1))
977 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
979 DO i = i_start, i_end
980 ! Assumes j>j_end_f is ONLY j_end+1 ...
981 ! fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) &
982 ! *(u(i,k,j_end+1)+u(i,k,j_end))
983 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) &
984 *(u(i,k,j)+u(i,k,j-1))
988 ! 3rd or 4th order flux
990 DO i = i_start, i_end
991 vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
992 fqy( i, k, jp1 ) = vel*flux3( u(i,k,j-2), u(i,k,j-1), &
993 u(i,k,j ), u(i,k,j+1), &
1000 ! y flux-divergence into tendency
1002 ! (j > j_start) will miss the u(,,jds) tendency
1003 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
1005 DO i = i_start, i_end
1006 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1007 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
1010 ! This would be seen by (j > j_start) but we need to zero out the NP tendency
1011 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
1013 DO i = i_start, i_end
1014 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1015 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
1020 IF (j > j_start) THEN
1023 DO i = i_start, i_end
1024 mrdy=msfux(i,j-1)*rdy ! ADT eqn 44, 2nd term on RHS
1025 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
1039 ELSE IF ( horz_order == 2 ) THEN
1044 j_end = MIN(jte,jde-1)
1046 IF ( config_flags%open_xs ) i_start = MAX(ids+1,its)
1047 IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite)
1048 IF ( specified ) i_start = MAX(ids+2,its)
1049 IF ( specified ) i_end = MIN(ide-2,ite)
1050 IF ( config_flags%periodic_x ) i_start = its
1051 IF ( config_flags%periodic_x ) i_end = ite
1053 DO j = j_start, j_end
1055 DO i = i_start, i_end
1056 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1057 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1058 *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) &
1059 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j)))
1064 IF ( specified .AND. its .LE. ids+1 .AND. .NOT. config_flags%periodic_x ) THEN
1065 DO j = j_start, j_end
1068 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1070 IF (u(i,k,j) .LT. 0.) ub = u(i,k,j)
1071 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1072 *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) &
1073 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+ub))
1077 IF ( specified .AND. ite .GE. ide-1 .AND. .NOT. config_flags%periodic_x ) THEN
1078 DO j = j_start, j_end
1081 mrdx=msfux(i,j)*rdx ! ADT eqn 44, 1st term on RHS
1083 IF (u(i,k,j) .GT. 0.) ub = u(i,k,j)
1084 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
1085 *((ru(i+1,k,j)+ru(i,k,j))*(ub+u(i,k,j)) &
1086 -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j)))
1091 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
1092 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
1094 DO j = j_start, j_end
1096 DO i = i_start, i_end
1097 mrdy=msfux(i,j)*rdy ! ADT eqn 44, 1st term on RHS
1098 ! Comments for polar boundary condition
1099 ! Flow is only from one side for points next to poles
1100 IF ( (config_flags%polar) .AND. (j == jds) ) THEN
1101 tendency(i,k,j)=tendency(i,k,j)-mrdy*0.25 &
1102 *(rv(i,k,j+1)+rv(i-1,k,j+1))*(u(i,k,j+1)+u(i,k,j))
1103 ELSE IF ( (config_flags%polar) .AND. (j == jde-1) ) THEN
1104 tendency(i,k,j)=tendency(i,k,j)+mrdy*0.25 &
1105 *(rv(i,k,j)+rv(i-1,k,j))*(u(i,k,j)+u(i,k,j-1))
1107 tendency(i,k,j)=tendency(i,k,j)-mrdy*0.25 &
1108 *((rv(i,k,j+1)+rv(i-1,k,j+1))*(u(i,k,j+1)+u(i,k,j)) &
1109 -(rv(i,k,j)+rv(i-1,k,j))*(u(i,k,j)+u(i,k,j-1)))
1115 ELSE IF ( horz_order == 0 ) THEN
1117 ! Just in case we want to turn horizontal advection off, we can do it
1121 WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: h_order not known ',horz_order
1122 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
1124 ENDIF horizontal_order_test
1126 ! radiative lateral boundary condition in x for normal velocity (u)
1128 IF ( (config_flags%open_xs) .and. its == ids ) THEN
1131 j_end = MIN(jte,jde-1)
1133 DO j = j_start, j_end
1135 ub = MIN(ru(its,k,j)-cb*mut(its,j), 0.)
1136 tendency(its,k,j) = tendency(its,k,j) &
1137 - rdx*ub*(u_old(its+1,k,j) - u_old(its,k,j))
1143 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1146 j_end = MIN(jte,jde-1)
1148 DO j = j_start, j_end
1150 ub = MAX(ru(ite,k,j)+cb*mut(ite-1,j), 0.)
1151 tendency(ite,k,j) = tendency(ite,k,j) &
1152 - rdx*ub*(u_old(ite,k,j) - u_old(ite-1,k,j))
1158 ! pick up the rest of the horizontal radiation boundary conditions.
1159 ! (these are the computations that don't require 'cb')
1160 ! first, set to index ranges
1163 i_end = MIN(ite,ide)
1167 IF (config_flags%open_xs) THEN
1168 i_start = MAX(ids+1, its)
1171 IF (config_flags%open_xe) THEN
1172 i_end = MIN(ite,ide-1)
1176 IF( (config_flags%open_ys) .and. (jts == jds)) THEN
1178 DO i = i_start, i_end
1180 mrdy=msfux(i,jts)*rdy ! ADT eqn 44, 2nd term on RHS
1182 im = MAX( imin, i-1 )
1186 vw = 0.5*(rv(ip,k,jts)+rv(im,k,jts))
1188 dvm = rv(ip,k,jts+1)-rv(ip,k,jts)
1189 dvp = rv(im,k,jts+1)-rv(im,k,jts)
1190 tendency(i,k,jts)=tendency(i,k,jts)-mrdy*( &
1191 vb*(u_old(i,k,jts+1)-u_old(i,k,jts)) &
1192 +0.5*u(i,k,jts)*(dvm+dvp))
1198 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
1200 DO i = i_start, i_end
1202 mrdy=msfux(i,jte-1)*rdy ! ADT eqn 44, 2nd term on RHS
1204 im = MAX( imin, i-1 )
1208 vw = 0.5*(rv(ip,k,jte)+rv(im,k,jte))
1210 dvm = rv(ip,k,jte)-rv(ip,k,jte-1)
1211 dvp = rv(im,k,jte)-rv(im,k,jte-1)
1212 tendency(i,k,jte-1)=tendency(i,k,jte-1)-mrdy*( &
1213 vb*(u_old(i,k,jte-1)-u_old(i,k,jte-2)) &
1214 +0.5*u(i,k,jte-1)*(dvm+dvp))
1220 !-------------------- vertical advection
1221 ! ADT eqn 44 has 3rd term on RHS = -(1/my) partial d/dz (rho u w)
1222 ! Here we have: - partial d/dz (u*rom) = - partial d/dz (u rho w / my)
1223 ! Since 'my' (map scale factor in y-direction) isn't a function of z,
1224 ! this is what we need, so leave unchanged in advect_u
1229 j_end = min(jte,jde-1)
1231 ! IF ( config_flags%open_xs ) i_start = MAX(ids+1,its)
1232 ! IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite)
1234 IF ( config_flags%open_ys .or. specified ) i_start = MAX(ids+1,its)
1235 IF ( config_flags%open_ye .or. specified ) i_end = MIN(ide-1,ite)
1236 IF ( config_flags%periodic_x ) i_start = its
1237 IF ( config_flags%periodic_x ) i_end = ite
1239 DO i = i_start, i_end
1244 vert_order_test : IF (vert_order == 6) THEN
1246 DO j = j_start, j_end
1249 DO i = i_start, i_end
1250 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1251 vflux(i,k) = vel*flux6( &
1252 u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), &
1253 u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel )
1257 DO i = i_start, i_end
1260 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1261 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1263 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1264 vflux(i,k) = vel*flux4( &
1265 u(i,k-2,j), u(i,k-1,j), &
1266 u(i,k ,j), u(i,k+1,j), -vel )
1268 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1269 vflux(i,k) = vel*flux4( &
1270 u(i,k-2,j), u(i,k-1,j), &
1271 u(i,k ,j), u(i,k+1,j), -vel )
1273 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1274 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1278 DO i = i_start, i_end
1279 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1284 ELSE IF (vert_order == 5) THEN
1286 DO j = j_start, j_end
1289 DO i = i_start, i_end
1290 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1291 vflux(i,k) = vel*flux5( &
1292 u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), &
1293 u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel )
1297 DO i = i_start, i_end
1300 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1301 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1303 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1304 vflux(i,k) = vel*flux3( &
1305 u(i,k-2,j), u(i,k-1,j), &
1306 u(i,k ,j), u(i,k+1,j), -vel )
1308 vel=0.5*(rom(i,k,j)+rom(i-1,k,j))
1309 vflux(i,k) = vel*flux3( &
1310 u(i,k-2,j), u(i,k-1,j), &
1311 u(i,k ,j), u(i,k+1,j), -vel )
1313 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1314 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1318 DO i = i_start, i_end
1319 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1324 ELSE IF (vert_order == 4) THEN
1326 DO j = j_start, j_end
1329 DO i = i_start, i_end
1330 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1331 vflux(i,k) = vel*flux4( &
1332 u(i,k-2,j), u(i,k-1,j), &
1333 u(i,k ,j), u(i,k+1,j), -vel )
1337 DO i = i_start, i_end
1340 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1341 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1343 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1344 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1348 DO i = i_start, i_end
1349 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1354 ELSE IF (vert_order == 3) THEN
1356 DO j = j_start, j_end
1359 DO i = i_start, i_end
1360 vel=0.5*(rom(i-1,k,j)+rom(i,k,j))
1361 vflux(i,k) = vel*flux3( &
1362 u(i,k-2,j), u(i,k-1,j), &
1363 u(i,k ,j), u(i,k+1,j), -vel )
1367 DO i = i_start, i_end
1370 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1371 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1373 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1374 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1378 DO i = i_start, i_end
1379 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1384 ELSE IF (vert_order == 2) THEN
1386 DO j = j_start, j_end
1388 DO i = i_start, i_end
1389 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) &
1390 *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j))
1396 DO i = i_start, i_end
1397 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
1405 WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: v_order not known ',vert_order
1406 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
1408 ENDIF vert_order_test
1410 END SUBROUTINE advect_u
1412 !-------------------------------------------------------------------------------
1414 SUBROUTINE advect_v ( v, v_old, tendency, &
1416 mut, time_step, config_flags, &
1417 msfux, msfuy, msfvx, msfvy, &
1421 ids, ide, jds, jde, kds, kde, &
1422 ims, ime, jms, jme, kms, kme, &
1423 its, ite, jts, jte, kts, kte )
1429 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1431 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1432 ims, ime, jms, jme, kms, kme, &
1433 its, ite, jts, jte, kts, kte
1435 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: v, &
1441 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
1442 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
1444 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
1451 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
1455 REAL , INTENT(IN ) :: rdx, &
1457 INTEGER , INTENT(IN ) :: time_step
1462 INTEGER :: i, j, k, itf, jtf, ktf
1463 INTEGER :: i_start, i_end, j_start, j_end
1464 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
1465 INTEGER :: jmin, jmax, jp, jm, imin, imax
1467 REAL :: mrdx, mrdy, ub, vb, uw, vw, dup, dum
1468 REAL , DIMENSION(its:ite, kts:kte) :: vflux
1471 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
1472 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
1474 INTEGER :: horz_order
1475 INTEGER :: vert_order
1477 LOGICAL :: degrade_xs, degrade_ys
1478 LOGICAL :: degrade_xe, degrade_ye
1480 INTEGER :: jp1, jp0, jtmp
1483 ! definition of flux operators, 3rd, 4th, 5th or 6th order
1485 REAL :: flux3, flux4, flux5, flux6
1486 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
1488 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
1489 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
1491 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
1492 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
1493 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
1495 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
1496 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
1497 +(q_ip2+q_im3) )/60.0
1499 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
1500 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
1501 -sign(1,time_step)*sign(1.,ua)*( &
1502 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
1506 LOGICAL :: specified
1509 if(config_flags%specified .or. config_flags%nested) specified = .true.
1511 ! set order for the advection schemes
1514 horz_order = config_flags%h_mom_adv_order
1515 vert_order = config_flags%v_mom_adv_order
1518 ! here is the choice of flux operators
1521 horizontal_order_test : IF( horz_order == 6 ) THEN
1523 ! determine boundary mods for flux operators
1524 ! We degrade the flux operators from 3rd/4th order
1525 ! to second order one gridpoint in from the boundaries for
1526 ! all boundary conditions except periodic and symmetry - these
1527 ! conditions have boundary zone data fill for correct application
1528 ! of the higher order flux stencils
1535 IF( config_flags%periodic_x .or. &
1536 config_flags%symmetric_xs .or. &
1537 (its > ids+3) ) degrade_xs = .false.
1538 IF( config_flags%periodic_x .or. &
1539 config_flags%symmetric_xe .or. &
1540 (ite < ide-3) ) degrade_xe = .false.
1541 IF( config_flags%periodic_y .or. &
1542 config_flags%symmetric_ys .or. &
1543 (jts > jds+3) ) degrade_ys = .false.
1544 IF( config_flags%periodic_y .or. &
1545 config_flags%symmetric_ye .or. &
1546 (jte < jde-3) ) degrade_ye = .false.
1548 !--------------- y - advection first
1551 i_end = MIN(ite,ide-1)
1555 ! higher order flux has a 5 or 7 point stencil, so compute
1556 ! bounds so we can switch to second order flux close to the boundary
1562 j_start = MAX(jts,jds+1)
1567 j_end = MIN(jte,jde-1)
1571 ! compute fluxes, 5th or 6th order
1576 j_loop_y_flux_6 : DO j = j_start, j_end+1
1578 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
1581 DO i = i_start, i_end
1582 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1583 fqy( i, k, jp1 ) = vel*flux6( &
1584 v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), &
1585 v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel )
1589 ! we must be close to some boundary where we need to reduce the order of the stencil
1590 ! specified uses upstream normal wind at boundaries
1592 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
1595 DO i = i_start, i_end
1597 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
1598 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1603 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
1606 DO i = i_start, i_end
1607 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1608 fqy( i, k, jp1 ) = vel*flux4( &
1609 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1614 ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary
1617 DO i = i_start, i_end
1619 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
1620 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1625 ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary
1628 DO i = i_start, i_end
1629 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1630 fqy( i, k, jp1 ) = vel*flux4( &
1631 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1637 ! y flux-divergence into tendency
1639 ! Comments on polar boundary conditions
1640 ! No advection over the poles means tendencies (held from jds [S. pole]
1641 ! to jde [N pole], i.e., on v grid) must be zero at poles
1642 ! [tendency(jds) and tendency(jde)=0]
1643 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
1645 DO i = i_start, i_end
1646 tendency(i,k,j-1) = 0.
1649 ! If j_end were set to jde in a special if statement apart from
1650 ! degrade_ye, then we would hit the next conditional. But since
1651 ! we want the tendency to be zero anyway, not looping to jde+1
1652 ! will produce the same effect.
1653 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
1655 DO i = i_start, i_end
1656 tendency(i,k,j-1) = 0.
1661 IF(j > j_start) THEN
1664 DO i = i_start, i_end
1665 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
1666 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
1678 ENDDO j_loop_y_flux_6
1680 ! next, x - flux divergence
1683 i_end = MIN(ite,ide-1)
1687 ! Polar boundary conditions are like open or specified
1688 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
1689 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
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 i_start = MAX(ids+1,its)
1699 ! i_start_f = i_start+2
1700 i_start_f = MIN(i_start+2,ids+3)
1704 i_end = MIN(ide-2,ite)
1710 DO j = j_start, j_end
1712 ! 5th or 6th order flux
1715 DO i = i_start_f, i_end_f
1716 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1717 fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j), &
1718 v(i-1,k,j), v(i ,k,j), &
1719 v(i+1,k,j), v(i+2,k,j), &
1724 ! lower order fluxes close to boundaries (if not periodic or symmetric)
1726 IF( degrade_xs ) THEN
1728 DO i=i_start,i_start_f-1
1730 IF(i == ids+1) THEN ! second order
1732 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
1733 *(v(i,k,j)+v(i-1,k,j))
1737 IF(i == ids+2) THEN ! third order
1739 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1740 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
1741 v(i ,k,j), v(i+1,k,j), &
1750 IF( degrade_xe ) THEN
1752 DO i = i_end_f+1, i_end+1
1754 IF( i == ide-1 ) THEN ! second order flux next to the boundary
1756 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
1757 *(v(i_end+1,k,j)+v(i_end,k,j))
1761 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
1763 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1764 fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
1765 v(i ,k,j), v(i+1,k,j), &
1774 ! x flux-divergence into tendency
1777 DO i = i_start, i_end
1778 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
1779 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
1785 ELSE IF( horz_order == 5 ) THEN
1787 ! 5th order horizontal flux calculation
1788 ! This code is EXACTLY the same as the 6th order code
1789 ! EXCEPT the 5th order and 3rd operators are used in
1790 ! place of the 6th and 4th order operators
1792 ! determine boundary mods for flux operators
1793 ! We degrade the flux operators from 3rd/4th order
1794 ! to second order one gridpoint in from the boundaries for
1795 ! all boundary conditions except periodic and symmetry - these
1796 ! conditions have boundary zone data fill for correct application
1797 ! of the higher order flux stencils
1804 IF( config_flags%periodic_x .or. &
1805 config_flags%symmetric_xs .or. &
1806 (its > ids+3) ) degrade_xs = .false.
1807 IF( config_flags%periodic_x .or. &
1808 config_flags%symmetric_xe .or. &
1809 (ite < ide-3) ) degrade_xe = .false.
1810 IF( config_flags%periodic_y .or. &
1811 config_flags%symmetric_ys .or. &
1812 (jts > jds+3) ) degrade_ys = .false.
1813 IF( config_flags%periodic_y .or. &
1814 config_flags%symmetric_ye .or. &
1815 (jte < jde-3) ) degrade_ye = .false.
1817 !--------------- y - advection first
1820 i_end = MIN(ite,ide-1)
1824 ! higher order flux has a 5 or 7 point stencil, so compute
1825 ! bounds so we can switch to second order flux close to the boundary
1831 j_start = MAX(jts,jds+1)
1836 j_end = MIN(jte,jde-1)
1840 ! compute fluxes, 5th or 6th order
1845 j_loop_y_flux_5 : DO j = j_start, j_end+1
1847 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
1850 DO i = i_start, i_end
1851 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1852 fqy( i, k, jp1 ) = vel*flux5( &
1853 v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), &
1854 v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel )
1858 ! we must be close to some boundary where we need to reduce the order of the stencil
1859 ! specified uses upstream normal wind at boundaries
1861 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
1864 DO i = i_start, i_end
1866 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
1867 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1872 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
1875 DO i = i_start, i_end
1876 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1877 fqy( i, k, jp1 ) = vel*flux3( &
1878 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1883 ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary
1886 DO i = i_start, i_end
1888 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
1889 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
1894 ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary
1897 DO i = i_start, i_end
1898 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
1899 fqy( i, k, jp1 ) = vel*flux3( &
1900 v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
1906 ! y flux-divergence into tendency
1908 ! Comments on polar boundary conditions
1909 ! No advection over the poles means tendencies (held from jds [S. pole]
1910 ! to jde [N pole], i.e., on v grid) must be zero at poles
1911 ! [tendency(jds) and tendency(jde)=0]
1912 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
1914 DO i = i_start, i_end
1915 tendency(i,k,j-1) = 0.
1918 ! If j_end were set to jde in a special if statement apart from
1919 ! degrade_ye, then we would hit the next conditional. But since
1920 ! we want the tendency to be zero anyway, not looping to jde+1
1921 ! will produce the same effect.
1922 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
1924 DO i = i_start, i_end
1925 tendency(i,k,j-1) = 0.
1930 IF(j > j_start) THEN
1933 DO i = i_start, i_end
1934 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
1935 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
1947 ENDDO j_loop_y_flux_5
1949 ! next, x - flux divergence
1952 i_end = MIN(ite,ide-1)
1956 ! Polar boundary conditions are like open or specified
1957 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
1958 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
1960 ! higher order flux has a 5 or 7 point stencil, so compute
1961 ! bounds so we can switch to second order flux close to the boundary
1967 i_start = MAX(ids+1,its)
1968 ! i_start_f = i_start+2
1969 i_start_f = MIN(i_start+2,ids+3)
1973 i_end = MIN(ide-2,ite)
1979 DO j = j_start, j_end
1981 ! 5th or 6th order flux
1984 DO i = i_start_f, i_end_f
1985 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
1986 fqx( i, k ) = vel*flux5( v(i-3,k,j), v(i-2,k,j), &
1987 v(i-1,k,j), v(i ,k,j), &
1988 v(i+1,k,j), v(i+2,k,j), &
1993 ! lower order fluxes close to boundaries (if not periodic or symmetric)
1995 IF( degrade_xs ) THEN
1997 DO i=i_start,i_start_f-1
1999 IF(i == ids+1) THEN ! second order
2001 fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
2002 *(v(i,k,j)+v(i-1,k,j))
2006 IF(i == ids+2) THEN ! third order
2008 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2009 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2010 v(i ,k,j), v(i+1,k,j), &
2019 IF( degrade_xe ) THEN
2021 DO i = i_end_f+1, i_end+1
2023 IF( i == ide-1 ) THEN ! second order flux next to the boundary
2025 fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2026 *(v(i_end+1,k,j)+v(i_end,k,j))
2030 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
2032 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2033 fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2034 v(i ,k,j), v(i+1,k,j), &
2043 ! x flux-divergence into tendency
2046 DO i = i_start, i_end
2047 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2048 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2054 ELSE IF( horz_order == 4 ) THEN
2056 ! determine boundary mods for flux operators
2057 ! We degrade the flux operators from 3rd/4th order
2058 ! to second order one gridpoint in from the boundaries for
2059 ! all boundary conditions except periodic and symmetry - these
2060 ! conditions have boundary zone data fill for correct application
2061 ! of the higher order flux stencils
2068 IF( config_flags%periodic_x .or. &
2069 config_flags%symmetric_xs .or. &
2070 (its > ids+2) ) degrade_xs = .false.
2071 IF( config_flags%periodic_x .or. &
2072 config_flags%symmetric_xe .or. &
2073 (ite < ide-2) ) degrade_xe = .false.
2074 IF( config_flags%periodic_y .or. &
2075 config_flags%symmetric_ys .or. &
2076 (jts > jds+2) ) degrade_ys = .false.
2077 IF( config_flags%periodic_y .or. &
2078 config_flags%symmetric_ye .or. &
2079 (jte < jde-2) ) degrade_ye = .false.
2081 !--------------- y - advection first
2087 i_end = MIN(ite,ide-1)
2091 ! 3rd or 4th order flux has a 5 point stencil, so compute
2092 ! bounds so we can switch to second order flux close to the boundary
2097 !CJM May not work with tiling because defined in terms of domain dims
2100 j_start_f = j_start+1
2109 ! specified uses upstream normal wind at boundaries
2114 DO j = j_start, j_end+1
2116 IF ((j == j_start) .and. degrade_ys) THEN
2118 DO i = i_start, i_end
2120 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
2121 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2125 ELSE IF ((j == j_end+1) .and. degrade_ye) THEN
2127 DO i = i_start, i_end
2129 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
2130 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2136 DO i = i_start, i_end
2137 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2138 fqy( i,k,jp1 ) = vel*flux4( v(i,k,j-2), v(i,k,j-1), &
2139 v(i,k,j ), v(i,k,j+1), &
2145 ! Comments on polar boundary conditions
2146 ! No advection over the poles means tendencies (held from jds [S. pole]
2147 ! to jde [N pole], i.e., on v grid) must be zero at poles
2148 ! [tendency(jds) and tendency(jde)=0]
2149 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
2151 DO i = i_start, i_end
2152 tendency(i,k,j-1) = 0.
2155 ! If j_end were set to jde in a special if statement apart from
2156 ! degrade_ye, then we would hit the next conditional. But since
2157 ! we want the tendency to be zero anyway, not looping to jde+1
2158 ! will produce the same effect.
2159 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
2161 DO i = i_start, i_end
2162 tendency(i,k,j-1) = 0.
2167 IF( j > j_start) THEN
2169 DO i = i_start, i_end
2170 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
2171 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
2185 ! next, x - flux divergence
2189 i_end = MIN(ite,ide-1)
2193 ! Polar boundary conditions are like open or specified
2194 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2195 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2197 ! 3rd or 4th order flux has a 5 point stencil, so compute
2198 ! bounds so we can switch to second order flux close to the boundary
2205 i_start_f = i_start+1
2215 DO j = j_start, j_end
2217 ! 3rd or 4th order flux
2220 DO i = i_start_f, i_end_f
2221 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2222 fqx( i, k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), &
2223 v(i ,k,j), v(i+1,k,j), &
2228 ! second order flux close to boundaries (if not periodic or symmetric)
2230 IF( degrade_xs ) THEN
2232 fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) &
2233 *(v(i_start,k,j)+v(i_start-1,k,j))
2237 IF( degrade_xe ) THEN
2239 fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2240 *(v(i_end+1,k,j)+v(i_end,k,j))
2244 ! x flux-divergence into tendency
2247 DO i = i_start, i_end
2248 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2249 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2255 ELSE IF( horz_order == 3 ) THEN
2257 ! determine boundary mods for flux operators
2258 ! We degrade the flux operators from 3rd/4th order
2259 ! to second order one gridpoint in from the boundaries for
2260 ! all boundary conditions except periodic and symmetry - these
2261 ! conditions have boundary zone data fill for correct application
2262 ! of the higher order flux stencils
2269 IF( config_flags%periodic_x .or. &
2270 config_flags%symmetric_xs .or. &
2271 (its > ids+2) ) degrade_xs = .false.
2272 IF( config_flags%periodic_x .or. &
2273 config_flags%symmetric_xe .or. &
2274 (ite < ide-2) ) degrade_xe = .false.
2275 IF( config_flags%periodic_y .or. &
2276 config_flags%symmetric_ys .or. &
2277 (jts > jds+2) ) degrade_ys = .false.
2278 IF( config_flags%periodic_y .or. &
2279 config_flags%symmetric_ye .or. &
2280 (jte < jde-2) ) degrade_ye = .false.
2282 !--------------- y - advection first
2288 i_end = MIN(ite,ide-1)
2292 ! 3rd or 4th order flux has a 5 point stencil, so compute
2293 ! bounds so we can switch to second order flux close to the boundary
2298 !CJM May not work with tiling because defined in terms of domain dims
2301 j_start_f = j_start+1
2310 ! specified uses upstream normal wind at boundaries
2315 DO j = j_start, j_end+1
2317 IF ((j == j_start) .and. degrade_ys) THEN
2319 DO i = i_start, i_end
2321 IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j)
2322 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2326 ELSE IF ((j == j_end+1) .and. degrade_ye) THEN
2328 DO i = i_start, i_end
2330 IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1)
2331 fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) &
2337 DO i = i_start, i_end
2338 vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
2339 fqy( i,k,jp1 ) = vel*flux3( v(i,k,j-2), v(i,k,j-1), &
2340 v(i,k,j ), v(i,k,j+1), &
2346 ! Comments on polar boundary conditions
2347 ! No advection over the poles means tendencies (held from jds [S. pole]
2348 ! to jde [N pole], i.e., on v grid) must be zero at poles
2349 ! [tendency(jds) and tendency(jde)=0]
2350 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
2352 DO i = i_start, i_end
2353 tendency(i,k,j-1) = 0.
2356 ! If j_end were set to jde in a special if statement apart from
2357 ! degrade_ye, then we would hit the next conditional. But since
2358 ! we want the tendency to be zero anyway, not looping to jde+1
2359 ! will produce the same effect.
2360 ELSE IF( config_flags%polar .AND. (j == jde+1) ) THEN
2362 DO i = i_start, i_end
2363 tendency(i,k,j-1) = 0.
2368 IF( j > j_start) THEN
2370 DO i = i_start, i_end
2371 mrdy=msfvy(i,j-1)*rdy ! ADT eqn 45, 2nd term on RHS
2372 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
2386 ! next, x - flux divergence
2390 i_end = MIN(ite,ide-1)
2394 ! Polar boundary conditions are like open or specified
2395 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2396 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2398 ! 3rd or 4th order flux has a 5 point stencil, so compute
2399 ! bounds so we can switch to second order flux close to the boundary
2406 i_start_f = i_start+1
2416 DO j = j_start, j_end
2418 ! 3rd or 4th order flux
2421 DO i = i_start_f, i_end_f
2422 vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
2423 fqx( i, k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), &
2424 v(i ,k,j), v(i+1,k,j), &
2429 ! second order flux close to boundaries (if not periodic or symmetric)
2431 IF( degrade_xs ) THEN
2433 fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) &
2434 *(v(i_start,k,j)+v(i_start-1,k,j))
2438 IF( degrade_xe ) THEN
2440 fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) &
2441 *(v(i_end+1,k,j)+v(i_end,k,j))
2445 ! x flux-divergence into tendency
2448 DO i = i_start, i_end
2449 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2450 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
2456 ELSE IF( horz_order == 2 ) THEN
2460 i_end = MIN(ite,ide-1)
2464 IF ( config_flags%open_ys ) j_start = MAX(jds+1,jts)
2465 IF ( config_flags%open_ye ) j_end = MIN(jde-1,jte)
2466 IF ( specified ) j_start = MAX(jds+2,jts)
2467 IF ( specified ) j_end = MIN(jde-2,jte)
2468 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2469 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2471 DO j = j_start, j_end
2473 DO i = i_start, i_end
2475 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2477 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2478 *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) &
2479 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1)))
2485 ! Comments on polar boundary conditions
2486 ! tendencies = 0 at poles, and polar points do not contribute at points
2488 IF (config_flags%polar) THEN
2489 IF (jts == jds) THEN
2491 DO i = i_start, i_end
2492 tendency(i,k,jds) = 0.
2496 IF (jte == jde) THEN
2498 DO i = i_start, i_end
2499 tendency(i,k,jde) = 0.
2505 ! specified uses upstream normal wind at boundaries
2507 IF ( specified .AND. jts .LE. jds+1 ) THEN
2510 DO i = i_start, i_end
2511 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2513 IF (v(i,k,j) .LT. 0.) vb = v(i,k,j)
2515 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2516 *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) &
2517 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+vb))
2523 IF ( specified .AND. jte .GE. jde-1 ) THEN
2526 DO i = i_start, i_end
2528 mrdy=msfvy(i,j)*rdy ! ADT eqn 45, 2nd term on RHS
2530 IF (v(i,k,j) .GT. 0.) vb = v(i,k,j)
2532 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 &
2533 *((rv(i,k,j+1)+rv(i,k,j ))*(vb+v(i,k,j )) &
2534 -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1)))
2540 IF ( .NOT. config_flags%periodic_x ) THEN
2541 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2542 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2544 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2545 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2547 DO j = j_start, j_end
2549 DO i = i_start, i_end
2551 mrdx=msfvy(i,j)*rdx ! ADT eqn 45, 1st term on RHS
2553 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 &
2554 *((ru(i+1,k,j)+ru(i+1,k,j-1))*(v(i+1,k,j)+v(i ,k,j)) &
2555 -(ru(i ,k,j)+ru(i ,k,j-1))*(v(i ,k,j)+v(i-1,k,j)))
2561 ELSE IF ( horz_order == 0 ) THEN
2563 ! Just in case we want to turn horizontal advection off, we can do it
2568 WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: h_order not known ',horz_order
2569 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
2571 ENDIF horizontal_order_test
2573 ! Comments on polar boundary condition
2574 ! Force tendency=0 at NP and SP
2575 ! We keep setting this everywhere, but it can't hurt...
2576 IF ( config_flags%polar .AND. (jts == jds) ) THEN
2579 tendency(i,k,jts)=0.
2583 IF ( config_flags%polar .AND. (jte == jde) ) THEN
2586 tendency(i,k,jte)=0.
2591 ! radiative lateral boundary condition in y for normal velocity (v)
2593 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2596 i_end = MIN(ite,ide-1)
2598 DO i = i_start, i_end
2600 vb = MIN(rv(i,k,jts)-cb*mut(i,jts), 0.)
2601 tendency(i,k,jts) = tendency(i,k,jts) &
2602 - rdy*vb*(v_old(i,k,jts+1) - v_old(i,k,jts))
2608 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2611 i_end = MIN(ite,ide-1)
2613 DO i = i_start, i_end
2615 vb = MAX(rv(i,k,jte)+cb*mut(i,jte-1), 0.)
2616 tendency(i,k,jte) = tendency(i,k,jte) &
2617 - rdy*vb*(v_old(i,k,jte) - v_old(i,k,jte-1))
2623 ! pick up the rest of the horizontal radiation boundary conditions.
2624 ! (these are the computations that don't require 'cb'.
2625 ! first, set to index ranges
2628 j_end = MIN(jte,jde)
2633 IF (config_flags%open_ys) THEN
2634 j_start = MAX(jds+1, jts)
2637 IF (config_flags%open_ye) THEN
2638 j_end = MIN(jte,jde-1)
2642 ! compute x (u) conditions for v, w, or scalar
2644 IF( (config_flags%open_xs) .and. (its == ids)) THEN
2646 DO j = j_start, j_end
2648 mrdx=msfvy(its,j)*rdx ! ADT eqn 45, 1st term on RHS
2650 jm = MAX( jmin, j-1 )
2654 uw = 0.5*(ru(its,k,jp)+ru(its,k,jm))
2656 dup = ru(its+1,k,jp)-ru(its,k,jp)
2657 dum = ru(its+1,k,jm)-ru(its,k,jm)
2658 tendency(its,k,j)=tendency(its,k,j)-mrdx*( &
2659 ub*(v_old(its+1,k,j)-v_old(its,k,j)) &
2660 +0.5*v(its,k,j)*(dup+dum))
2666 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
2667 DO j = j_start, j_end
2669 mrdx=msfvy(ite-1,j)*rdx ! ADT eqn 45, 1st term on RHS
2671 jm = MAX( jmin, j-1 )
2675 uw = 0.5*(ru(ite,k,jp)+ru(ite,k,jm))
2677 dup = ru(ite,k,jp)-ru(ite-1,k,jp)
2678 dum = ru(ite,k,jm)-ru(ite-1,k,jm)
2680 ! tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( &
2681 ! ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) &
2682 ! +0.5*v(ite-1,k,j)* &
2683 ! ( ru(ite,k,jp)-ru(ite-1,k,jp) &
2684 ! +ru(ite,k,jm)-ru(ite-1,k,jm)) )
2685 tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( &
2686 ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) &
2687 +0.5*v(ite-1,k,j)*(dup+dum))
2694 !-------------------- vertical advection
2695 ! ADT eqn 45 has 3rd term on RHS = -(1/mx) partial d/dz (rho v w)
2696 ! Here we have: - partial d/dz (v*rom) = - partial d/dz (v rho w / my)
2697 ! We therefore need to make a correction for advect_v
2698 ! since 'my' (map scale factor in y direction) isn't a function of z,
2699 ! we can do this using *(my/mx) (see eqn. 45 for example)
2703 i_end = MIN(ite,ide-1)
2707 DO i = i_start, i_end
2712 ! Polar boundary conditions are like open or specified
2713 ! We don't want to calculate vertical v tendencies at the N or S pole
2714 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
2715 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-1,jte)
2717 vert_order_test : IF (vert_order == 6) THEN
2719 DO j = j_start, j_end
2723 DO i = i_start, i_end
2724 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2725 vflux(i,k) = vel*flux6( &
2726 v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), &
2727 v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel )
2731 DO i = i_start, i_end
2733 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2734 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2736 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2737 vflux(i,k) = vel*flux4( &
2738 v(i,k-2,j), v(i,k-1,j), &
2739 v(i,k ,j), v(i,k+1,j), -vel )
2741 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2742 vflux(i,k) = vel*flux4( &
2743 v(i,k-2,j), v(i,k-1,j), &
2744 v(i,k ,j), v(i,k+1,j), -vel )
2746 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2747 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2753 DO i = i_start, i_end
2754 ! We are calculating vertical fluxes on v points,
2755 ! so we must mean msf_v_x/y variables
2756 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
2762 ELSE IF (vert_order == 5) THEN
2764 DO j = j_start, j_end
2768 DO i = i_start, i_end
2769 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2770 vflux(i,k) = vel*flux5( &
2771 v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), &
2772 v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel )
2776 DO i = i_start, i_end
2778 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2779 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2781 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2782 vflux(i,k) = vel*flux3( &
2783 v(i,k-2,j), v(i,k-1,j), &
2784 v(i,k ,j), v(i,k+1,j), -vel )
2786 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2787 vflux(i,k) = vel*flux3( &
2788 v(i,k-2,j), v(i,k-1,j), &
2789 v(i,k ,j), v(i,k+1,j), -vel )
2791 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2792 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2798 DO i = i_start, i_end
2799 ! We are calculating vertical fluxes on v points,
2800 ! so we must mean msf_v_x/y variables
2801 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
2807 ELSE IF (vert_order == 4) THEN
2809 DO j = j_start, j_end
2813 DO i = i_start, i_end
2814 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2815 vflux(i,k) = vel*flux4( &
2816 v(i,k-2,j), v(i,k-1,j), &
2817 v(i,k ,j), v(i,k+1,j), -vel )
2821 DO i = i_start, i_end
2823 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2824 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2826 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2827 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2833 DO i = i_start, i_end
2834 ! We are calculating vertical fluxes on v points,
2835 ! so we must mean msf_v_x/y variables
2836 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
2842 ELSE IF (vert_order == 3) THEN
2844 DO j = j_start, j_end
2848 DO i = i_start, i_end
2849 vel=0.5*(rom(i,k,j)+rom(i,k,j-1))
2850 vflux(i,k) = vel*flux3( &
2851 v(i,k-2,j), v(i,k-1,j), &
2852 v(i,k ,j), v(i,k+1,j), -vel )
2856 DO i = i_start, i_end
2858 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2859 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2861 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2862 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2868 DO i = i_start, i_end
2869 ! We are calculating vertical fluxes on v points,
2870 ! so we must mean msf_v_x/y variables
2871 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
2878 ELSE IF (vert_order == 2) THEN
2880 DO j = j_start, j_end
2882 DO i = i_start, i_end
2884 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) &
2885 *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j))
2890 DO i = i_start, i_end
2891 ! We are calculating vertical fluxes on v points,
2892 ! so we must mean msf_v_x/y variables
2893 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
2900 WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: v_order not known ',vert_order
2901 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
2903 ENDIF vert_order_test
2905 END SUBROUTINE advect_v
2907 !-------------------------------------------------------------------
2909 SUBROUTINE advect_scalar ( field, field_old, tendency, &
2911 mut, time_step, config_flags, &
2912 msfux, msfuy, msfvx, msfvy, &
2916 ids, ide, jds, jde, kds, kde, &
2917 ims, ime, jms, jme, kms, kme, &
2918 its, ite, jts, jte, kts, kte )
2924 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2926 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2927 ims, ime, jms, jme, kms, kme, &
2928 its, ite, jts, jte, kts, kte
2930 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2936 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
2937 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2939 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2946 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
2950 REAL , INTENT(IN ) :: rdx, &
2952 INTEGER , INTENT(IN ) :: time_step
2957 INTEGER :: i, j, k, itf, jtf, ktf
2958 INTEGER :: i_start, i_end, j_start, j_end
2959 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
2960 INTEGER :: jmin, jmax, jp, jm, imin, imax
2962 REAL :: mrdx, mrdy, ub, vb, uw, vw
2963 REAL , DIMENSION(its:ite, kts:kte) :: vflux
2966 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
2967 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
2969 INTEGER :: horz_order, vert_order
2971 LOGICAL :: degrade_xs, degrade_ys
2972 LOGICAL :: degrade_xe, degrade_ye
2974 INTEGER :: jp1, jp0, jtmp
2977 ! definition of flux operators, 3rd, 4th, 5th or 6th order
2979 REAL :: flux3, flux4, flux5, flux6
2980 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
2982 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
2983 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
2985 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
2986 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
2987 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
2989 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
2990 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
2991 +(q_ip2+q_im3) )/60.0
2993 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
2994 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
2995 -sign(1,time_step)*sign(1.,ua)*( &
2996 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
2999 LOGICAL :: specified
3002 if(config_flags%specified .or. config_flags%nested) specified = .true.
3004 ! set order for the advection schemes
3007 horz_order = config_flags%h_sca_adv_order
3008 vert_order = config_flags%v_sca_adv_order
3010 ! begin with horizontal flux divergence
3011 ! here is the choice of flux operators
3014 horizontal_order_test : IF( horz_order == 6 ) THEN
3016 ! determine boundary mods for flux operators
3017 ! We degrade the flux operators from 3rd/4th order
3018 ! to second order one gridpoint in from the boundaries for
3019 ! all boundary conditions except periodic and symmetry - these
3020 ! conditions have boundary zone data fill for correct application
3021 ! of the higher order flux stencils
3028 IF( config_flags%periodic_x .or. &
3029 config_flags%symmetric_xs .or. &
3030 (its > ids+3) ) degrade_xs = .false.
3031 IF( config_flags%periodic_x .or. &
3032 config_flags%symmetric_xe .or. &
3033 (ite < ide-3) ) degrade_xe = .false.
3034 IF( config_flags%periodic_y .or. &
3035 config_flags%symmetric_ys .or. &
3036 (jts > jds+3) ) degrade_ys = .false.
3037 IF( config_flags%periodic_y .or. &
3038 config_flags%symmetric_ye .or. &
3039 (jte < jde-4) ) degrade_ye = .false.
3041 !--------------- y - advection first
3045 i_end = MIN(ite,ide-1)
3047 j_end = MIN(jte,jde-1)
3049 ! higher order flux has a 5 or 7 point stencil, so compute
3050 ! bounds so we can switch to second order flux close to the boundary
3056 j_start = MAX(jts,jds+1)
3061 j_end = MIN(jte,jde-2)
3065 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3067 ! compute fluxes, 5th or 6th order
3072 j_loop_y_flux_6 : DO j = j_start, j_end+1
3074 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
3077 DO i = i_start, i_end
3079 fqy( i, k, jp1 ) = vel*flux6( &
3080 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
3081 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
3086 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
3089 DO i = i_start, i_end
3090 fqy(i,k, jp1) = 0.5*rv(i,k,j)* &
3091 (field(i,k,j)+field(i,k,j-1))
3096 ELSE IF ( j == jds+2 ) THEN ! 4th order flux 2 in from south boundary
3099 DO i = i_start, i_end
3101 fqy( i, k, jp1 ) = vel*flux4( &
3102 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
3106 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
3109 DO i = i_start, i_end
3110 fqy(i, k, jp1) = 0.5*rv(i,k,j)* &
3111 (field(i,k,j)+field(i,k,j-1))
3115 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
3118 DO i = i_start, i_end
3120 fqy( i, k, jp1) = vel*flux4( &
3121 field(i,k,j-2),field(i,k,j-1), &
3122 field(i,k,j),field(i,k,j+1),vel )
3128 ! y flux-divergence into tendency
3130 ! Comments on polar boundary conditions
3131 ! Same process as for advect_u - tendencies run from jds to jde-1
3132 ! (latitudes are as for u grid, longitudes are displaced)
3133 ! Therefore: flow is only from one side for points next to poles
3134 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3136 DO i = i_start, i_end
3137 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3138 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3141 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3143 DO i = i_start, i_end
3144 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3145 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3150 IF(j > j_start) THEN
3153 DO i = i_start, i_end
3154 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3155 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3167 ENDDO j_loop_y_flux_6
3169 ! next, x - flux divergence
3172 i_end = MIN(ite,ide-1)
3175 j_end = MIN(jte,jde-1)
3177 ! higher order flux has a 5 or 7 point stencil, so compute
3178 ! bounds so we can switch to second order flux close to the boundary
3184 i_start = MAX(ids+1,its)
3185 ! i_start_f = i_start+2
3186 i_start_f = MIN(i_start+2,ids+3)
3190 i_end = MIN(ide-2,ite)
3196 DO j = j_start, j_end
3198 ! 5th or 6th order flux
3201 DO i = i_start_f, i_end_f
3203 fqx( i,k ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), &
3204 field(i-1,k,j), field(i ,k,j), &
3205 field(i+1,k,j), field(i+2,k,j), &
3210 ! lower order fluxes close to boundaries (if not periodic or symmetric)
3212 IF( degrade_xs ) THEN
3214 DO i=i_start,i_start_f-1
3216 IF(i == ids+1) THEN ! second order
3218 fqx(i,k) = 0.5*(ru(i,k,j)) &
3219 *(field(i,k,j)+field(i-1,k,j))
3223 IF(i == ids+2) THEN ! third order
3226 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
3227 field(i ,k,j), field(i+1,k,j), &
3236 IF( degrade_xe ) THEN
3238 DO i = i_end_f+1, i_end+1
3240 IF( i == ide-1 ) THEN ! second order flux next to the boundary
3242 fqx(i,k) = 0.5*(ru(i,k,j)) &
3243 *(field(i,k,j)+field(i-1,k,j))
3247 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
3250 fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
3251 field(i ,k,j), field(i+1,k,j), &
3260 ! x flux-divergence into tendency
3263 DO i = i_start, i_end
3264 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3265 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3271 ELSE IF( horz_order == 5 ) THEN
3273 ! determine boundary mods for flux operators
3274 ! We degrade the flux operators from 3rd/4th order
3275 ! to second order one gridpoint in from the boundaries for
3276 ! all boundary conditions except periodic and symmetry - these
3277 ! conditions have boundary zone data fill for correct application
3278 ! of the higher order flux stencils
3285 IF( config_flags%periodic_x .or. &
3286 config_flags%symmetric_xs .or. &
3287 (its > ids+3) ) degrade_xs = .false.
3288 IF( config_flags%periodic_x .or. &
3289 config_flags%symmetric_xe .or. &
3290 (ite < ide-3) ) degrade_xe = .false.
3291 IF( config_flags%periodic_y .or. &
3292 config_flags%symmetric_ys .or. &
3293 (jts > jds+3) ) degrade_ys = .false.
3294 IF( config_flags%periodic_y .or. &
3295 config_flags%symmetric_ye .or. &
3296 (jte < jde-4) ) degrade_ye = .false.
3298 !--------------- y - advection first
3302 i_end = MIN(ite,ide-1)
3304 j_end = MIN(jte,jde-1)
3306 ! higher order flux has a 5 or 7 point stencil, so compute
3307 ! bounds so we can switch to second order flux close to the boundary
3313 j_start = MAX(jts,jds+1)
3318 j_end = MIN(jte,jde-2)
3322 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3324 ! compute fluxes, 5th or 6th order
3329 j_loop_y_flux_5 : DO j = j_start, j_end+1
3331 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
3334 DO i = i_start, i_end
3336 fqy( i, k, jp1 ) = vel*flux5( &
3337 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
3338 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
3343 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
3346 DO i = i_start, i_end
3347 fqy(i,k, jp1) = 0.5*rv(i,k,j)* &
3348 (field(i,k,j)+field(i,k,j-1))
3353 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
3356 DO i = i_start, i_end
3358 fqy( i, k, jp1 ) = vel*flux3( &
3359 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
3363 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
3366 DO i = i_start, i_end
3367 fqy(i, k, jp1) = 0.5*rv(i,k,j)* &
3368 (field(i,k,j)+field(i,k,j-1))
3372 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
3375 DO i = i_start, i_end
3377 fqy( i, k, jp1) = vel*flux3( &
3378 field(i,k,j-2),field(i,k,j-1), &
3379 field(i,k,j),field(i,k,j+1),vel )
3385 ! y flux-divergence into tendency
3387 ! Comments on polar boundary conditions
3388 ! Same process as for advect_u - tendencies run from jds to jde-1
3389 ! (latitudes are as for u grid, longitudes are displaced)
3390 ! Therefore: flow is only from one side for points next to poles
3391 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3393 DO i = i_start, i_end
3394 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3395 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3398 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3400 DO i = i_start, i_end
3401 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3402 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3407 IF(j > j_start) THEN
3410 DO i = i_start, i_end
3411 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3412 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3424 ENDDO j_loop_y_flux_5
3426 ! next, x - flux divergence
3429 i_end = MIN(ite,ide-1)
3432 j_end = MIN(jte,jde-1)
3434 ! higher order flux has a 5 or 7 point stencil, so compute
3435 ! bounds so we can switch to second order flux close to the boundary
3441 i_start = MAX(ids+1,its)
3442 ! i_start_f = i_start+2
3443 i_start_f = MIN(i_start+2,ids+3)
3447 i_end = MIN(ide-2,ite)
3453 DO j = j_start, j_end
3455 ! 5th or 6th order flux
3458 DO i = i_start_f, i_end_f
3460 fqx( i,k ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), &
3461 field(i-1,k,j), field(i ,k,j), &
3462 field(i+1,k,j), field(i+2,k,j), &
3467 ! lower order fluxes close to boundaries (if not periodic or symmetric)
3469 IF( degrade_xs ) THEN
3471 DO i=i_start,i_start_f-1
3473 IF(i == ids+1) THEN ! second order
3475 fqx(i,k) = 0.5*(ru(i,k,j)) &
3476 *(field(i,k,j)+field(i-1,k,j))
3480 IF(i == ids+2) THEN ! third order
3483 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
3484 field(i ,k,j), field(i+1,k,j), &
3493 IF( degrade_xe ) THEN
3495 DO i = i_end_f+1, i_end+1
3497 IF( i == ide-1 ) THEN ! second order flux next to the boundary
3499 fqx(i,k) = 0.5*(ru(i,k,j)) &
3500 *(field(i,k,j)+field(i-1,k,j))
3504 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
3507 fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
3508 field(i ,k,j), field(i+1,k,j), &
3517 ! x flux-divergence into tendency
3520 DO i = i_start, i_end
3521 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3522 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3529 ELSE IF( horz_order == 4 ) THEN
3536 IF( config_flags%periodic_x .or. &
3537 config_flags%symmetric_xs .or. &
3538 (its > ids+2) ) degrade_xs = .false.
3539 IF( config_flags%periodic_x .or. &
3540 config_flags%symmetric_xe .or. &
3541 (ite < ide-2) ) degrade_xe = .false.
3542 IF( config_flags%periodic_y .or. &
3543 config_flags%symmetric_ys .or. &
3544 (jts > jds+2) ) degrade_ys = .false.
3545 IF( config_flags%periodic_y .or. &
3546 config_flags%symmetric_ye .or. &
3547 (jte < jde-3) ) degrade_ye = .false.
3549 ! begin flux computations
3550 ! start with x flux divergence
3555 i_end = MIN(ite,ide-1)
3557 j_end = MIN(jte,jde-1)
3559 ! 3rd or 4th order flux has a 5 point stencil, so compute
3560 ! bounds so we can switch to second order flux close to the boundary
3567 i_start_f = i_start+1
3577 DO j = j_start, j_end
3579 ! 3rd or 4th order flux
3582 DO i = i_start_f, i_end_f
3584 fqx( i, k) = ru(i,k,j)*flux4( field(i-2,k,j), field(i-1,k,j), &
3585 field(i ,k,j), field(i+1,k,j), &
3590 ! second order flux close to boundaries (if not periodic or symmetric)
3592 IF( degrade_xs ) THEN
3594 fqx(i_start, k) = 0.5*ru(i_start,k,j) &
3595 *(field(i_start,k,j)+field(i_start-1,k,j))
3599 IF( degrade_xe ) THEN
3601 fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) &
3602 *(field(i_end+1,k,j)+field(i_end,k,j))
3606 ! x flux-divergence into tendency
3609 DO i = i_start, i_end
3610 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3611 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3618 ! next -> y flux divergence calculation
3621 i_end = MIN(ite,ide-1)
3623 j_end = MIN(jte,jde-1)
3625 ! 3rd or 4th order flux has a 5 point stencil, so compute
3626 ! bounds so we can switch to second order flux close to the boundary
3633 j_start_f = j_start+1
3641 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3646 DO j = j_start, j_end+1
3648 IF ((j < j_start_f) .and. degrade_ys) THEN
3650 DO i = i_start, i_end
3651 fqy(i,k,jp1) = 0.5*rv(i,k,j_start) &
3652 *(field(i,k,j_start)+field(i,k,j_start-1))
3655 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
3657 DO i = i_start, i_end
3658 ! Assumes j>j_end_f is ONLY j_end+1 ...
3659 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) &
3660 ! *(field(i,k,j_end+1)+field(i,k,j_end))
3661 fqy(i,k,jp1) = 0.5*rv(i,k,j) &
3662 *(field(i,k,j)+field(i,k,j-1))
3666 ! 3rd or 4th order flux
3668 DO i = i_start, i_end
3669 fqy( i, k, jp1 ) = rv(i,k,j)*flux4( field(i,k,j-2), field(i,k,j-1), &
3670 field(i,k,j ), field(i,k,j+1), &
3676 ! y flux-divergence into tendency
3678 ! Comments on polar boundary conditions
3679 ! Same process as for advect_u - tendencies run from jds to jde-1
3680 ! (latitudes are as for u grid, longitudes are displaced)
3681 ! Therefore: flow is only from one side for points next to poles
3682 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3684 DO i = i_start, i_end
3685 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3686 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3689 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3691 DO i = i_start, i_end
3692 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3693 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3698 IF ( j > j_start ) THEN
3701 DO i = i_start, i_end
3702 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3703 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3718 ELSE IF( horz_order == 3 ) THEN
3725 IF( config_flags%periodic_x .or. &
3726 config_flags%symmetric_xs .or. &
3727 (its > ids+2) ) degrade_xs = .false.
3728 IF( config_flags%periodic_x .or. &
3729 config_flags%symmetric_xe .or. &
3730 (ite < ide-2) ) degrade_xe = .false.
3731 IF( config_flags%periodic_y .or. &
3732 config_flags%symmetric_ys .or. &
3733 (jts > jds+2) ) degrade_ys = .false.
3734 IF( config_flags%periodic_y .or. &
3735 config_flags%symmetric_ye .or. &
3736 (jte < jde-3) ) degrade_ye = .false.
3738 ! begin flux computations
3739 ! start with x flux divergence
3744 i_end = MIN(ite,ide-1)
3746 j_end = MIN(jte,jde-1)
3748 ! 3rd or 4th order flux has a 5 point stencil, so compute
3749 ! bounds so we can switch to second order flux close to the boundary
3756 i_start_f = i_start+1
3766 DO j = j_start, j_end
3768 ! 3rd or 4th order flux
3771 DO i = i_start_f, i_end_f
3773 fqx( i, k) = ru(i,k,j)*flux3( field(i-2,k,j), field(i-1,k,j), &
3774 field(i ,k,j), field(i+1,k,j), &
3779 ! second order flux close to boundaries (if not periodic or symmetric)
3781 IF( degrade_xs ) THEN
3783 fqx(i_start, k) = 0.5*ru(i_start,k,j) &
3784 *(field(i_start,k,j)+field(i_start-1,k,j))
3788 IF( degrade_xe ) THEN
3790 fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) &
3791 *(field(i_end+1,k,j)+field(i_end,k,j))
3795 ! x flux-divergence into tendency
3798 DO i = i_start, i_end
3799 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3800 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
3807 ! next -> y flux divergence calculation
3810 i_end = MIN(ite,ide-1)
3812 j_end = MIN(jte,jde-1)
3814 ! 3rd or 4th order flux has a 5 point stencil, so compute
3815 ! bounds so we can switch to second order flux close to the boundary
3822 j_start_f = j_start+1
3830 IF(config_flags%polar) j_end = MIN(jte,jde-1)
3835 DO j = j_start, j_end+1
3837 IF ((j < j_start_f) .and. degrade_ys) THEN
3839 DO i = i_start, i_end
3840 fqy(i,k,jp1) = 0.5*rv(i,k,j_start) &
3841 *(field(i,k,j_start)+field(i,k,j_start-1))
3844 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
3846 DO i = i_start, i_end
3847 ! Assumes j>j_end_f is ONLY j_end+1 ...
3848 ! fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) &
3849 ! *(field(i,k,j_end+1)+field(i,k,j_end))
3850 fqy(i,k,jp1) = 0.5*rv(i,k,j) &
3851 *(field(i,k,j)+field(i,k,j-1))
3855 ! 3rd or 4th order flux
3857 DO i = i_start, i_end
3858 fqy( i, k, jp1 ) = rv(i,k,j)*flux3( field(i,k,j-2), field(i,k,j-1), &
3859 field(i,k,j ), field(i,k,j+1), &
3865 ! y flux-divergence into tendency
3867 ! Comments on polar boundary conditions
3868 ! Same process as for advect_u - tendencies run from jds to jde-1
3869 ! (latitudes are as for u grid, longitudes are displaced)
3870 ! Therefore: flow is only from one side for points next to poles
3871 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
3873 DO i = i_start, i_end
3874 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3875 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
3878 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
3880 DO i = i_start, i_end
3881 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3882 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
3887 IF ( j > j_start ) THEN
3890 DO i = i_start, i_end
3891 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3892 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
3906 ELSE IF( horz_order == 2 ) THEN
3909 i_end = MIN(ite,ide-1)
3911 j_end = MIN(jte,jde-1)
3913 IF ( .NOT. config_flags%periodic_x ) THEN
3914 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3915 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
3918 DO j = j_start, j_end
3920 DO i = i_start, i_end
3921 mrdx=msftx(i,j)*rdx ! see ADT eqn 48 [rho->rho*q] dividing by my, 1st term RHS
3922 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
3923 *(ru(i+1,k,j)*(field(i+1,k,j)+field(i ,k,j)) &
3924 -ru(i ,k,j)*(field(i ,k,j)+field(i-1,k,j)))
3930 i_end = MIN(ite,ide-1)
3932 ! Polar boundary conditions are like open or specified
3933 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
3934 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-2,jte)
3936 DO j = j_start, j_end
3938 DO i = i_start, i_end
3939 mrdy=msftx(i,j)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3940 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
3941 *(rv(i,k,j+1)*(field(i,k,j+1)+field(i,k,j )) &
3942 -rv(i,k,j )*(field(i,k,j )+field(i,k,j-1)))
3947 ! Polar boundary condtions
3948 ! These won't be covered in the loop above...
3949 IF (config_flags%polar) THEN
3950 IF (jts == jds) THEN
3952 DO i = i_start, i_end
3953 mrdy=msftx(i,jds)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3954 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
3955 *rv(i,k,jds+1)*(field(i,k,jds+1)+field(i,k,jds))
3959 IF (jte == jde) THEN
3961 DO i = i_start, i_end
3962 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 48 [rho->rho*q] dividing by my, 2nd term RHS
3963 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
3964 *rv(i,k,jde-1)*(field(i,k,jde-1)+field(i,k,jde-2))
3970 ELSE IF ( horz_order == 0 ) THEN
3972 ! Just in case we want to turn horizontal advection off, we can do it
3976 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_6a, h_order not known ',horz_order
3977 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
3979 ENDIF horizontal_order_test
3981 ! pick up the rest of the horizontal radiation boundary conditions.
3982 ! (these are the computations that don't require 'cb'.
3983 ! first, set to index ranges
3986 i_end = MIN(ite,ide-1)
3988 j_end = MIN(jte,jde-1)
3990 ! compute x (u) conditions for v, w, or scalar
3992 IF( (config_flags%open_xs) .and. (its == ids) ) THEN
3994 DO j = j_start, j_end
3996 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
3997 tendency(its,k,j) = tendency(its,k,j) &
3999 ub*( field_old(its+1,k,j) &
4000 - field_old(its ,k,j) ) + &
4001 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) &
4008 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
4010 DO j = j_start, j_end
4012 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
4013 tendency(i_end,k,j) = tendency(i_end,k,j) &
4015 ub*( field_old(i_end ,k,j) &
4016 - field_old(i_end-1,k,j) ) + &
4017 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) &
4024 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
4026 DO i = i_start, i_end
4028 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
4029 tendency(i,k,jts) = tendency(i,k,jts) &
4031 vb*( field_old(i,k,jts+1) &
4032 - field_old(i,k,jts ) ) + &
4033 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) &
4040 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
4042 DO i = i_start, i_end
4044 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
4045 tendency(i,k,j_end) = tendency(i,k,j_end) &
4047 vb*( field_old(i,k,j_end ) &
4048 - field_old(i,k,j_end-1) ) + &
4049 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) &
4057 !-------------------- vertical advection
4058 ! Scalar equation has 3rd term on RHS = - partial d/dz (q rho w /my)
4059 ! Here we have: - partial d/dz (q*rom) = - partial d/dz (q rho w / my)
4060 ! So we don't need to make a correction for advect_scalar
4063 i_end = MIN(ite,ide-1)
4065 j_end = MIN(jte,jde-1)
4067 DO i = i_start, i_end
4072 vert_order_test : IF (vert_order == 6) THEN
4074 DO j = j_start, j_end
4077 DO i = i_start, i_end
4079 vflux(i,k) = vel*flux6( &
4080 field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
4081 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
4085 DO i = i_start, i_end
4088 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4092 vflux(i,k) = vel*flux4( &
4093 field(i,k-2,j), field(i,k-1,j), &
4094 field(i,k ,j), field(i,k+1,j), -vel )
4097 vflux(i,k) = vel*flux4( &
4098 field(i,k-2,j), field(i,k-1,j), &
4099 field(i,k ,j), field(i,k+1,j), -vel )
4102 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4106 DO i = i_start, i_end
4107 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4113 ELSE IF (vert_order == 5) THEN
4115 DO j = j_start, j_end
4118 DO i = i_start, i_end
4120 vflux(i,k) = vel*flux5( &
4121 field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
4122 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
4126 DO i = i_start, i_end
4129 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4133 vflux(i,k) = vel*flux3( &
4134 field(i,k-2,j), field(i,k-1,j), &
4135 field(i,k ,j), field(i,k+1,j), -vel )
4138 vflux(i,k) = vel*flux3( &
4139 field(i,k-2,j), field(i,k-1,j), &
4140 field(i,k ,j), field(i,k+1,j), -vel )
4143 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4147 DO i = i_start, i_end
4148 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4154 ELSE IF (vert_order == 4) THEN
4156 DO j = j_start, j_end
4159 DO i = i_start, i_end
4161 vflux(i,k) = vel*flux4( &
4162 field(i,k-2,j), field(i,k-1,j), &
4163 field(i,k ,j), field(i,k+1,j), -vel )
4167 DO i = i_start, i_end
4170 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4172 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4176 DO i = i_start, i_end
4177 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4183 ELSE IF (vert_order == 3) THEN
4185 DO j = j_start, j_end
4188 DO i = i_start, i_end
4190 vflux(i,k) = vel*flux3( &
4191 field(i,k-2,j), field(i,k-1,j), &
4192 field(i,k ,j), field(i,k+1,j), -vel )
4196 DO i = i_start, i_end
4199 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4201 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4205 DO i = i_start, i_end
4206 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4212 ELSE IF (vert_order == 2) THEN
4214 DO j = j_start, j_end
4216 DO i = i_start, i_end
4217 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
4222 DO i = i_start, i_end
4223 tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
4231 WRITE (wrf_err_message,*) ' advect_scalar_6a, v_order not known ',vert_order
4232 CALL wrf_error_fatal ( wrf_err_message )
4234 ENDIF vert_order_test
4236 END SUBROUTINE advect_scalar
4238 !---------------------------------------------------------------------------------
4240 SUBROUTINE advect_w ( w, w_old, tendency, &
4242 mut, time_step, config_flags, &
4243 msfux, msfuy, msfvx, msfvy, &
4247 ids, ide, jds, jde, kds, kde, &
4248 ims, ime, jms, jme, kms, kme, &
4249 its, ite, jts, jte, kts, kte )
4255 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4257 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4258 ims, ime, jms, jme, kms, kme, &
4259 its, ite, jts, jte, kts, kte
4261 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: w, &
4267 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut
4268 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4270 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4277 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4281 REAL , INTENT(IN ) :: rdx, &
4283 INTEGER , INTENT(IN ) :: time_step
4288 INTEGER :: i, j, k, itf, jtf, ktf
4289 INTEGER :: i_start, i_end, j_start, j_end
4290 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
4291 INTEGER :: jmin, jmax, jp, jm, imin, imax
4293 REAL :: mrdx, mrdy, ub, vb, uw, vw
4294 REAL , DIMENSION(its:ite, kts:kte) :: vflux
4296 INTEGER :: horz_order, vert_order
4298 REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx
4299 REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy
4301 LOGICAL :: degrade_xs, degrade_ys
4302 LOGICAL :: degrade_xe, degrade_ye
4304 INTEGER :: jp1, jp0, jtmp
4306 ! definition of flux operators, 3rd, 4th, 5th or 6th order
4308 REAL :: flux3, flux4, flux5, flux6
4309 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel
4311 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
4312 ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
4314 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
4315 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
4316 sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
4318 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
4319 ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) &
4320 +(q_ip2+q_im3) )/60.0
4322 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
4323 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
4324 -sign(1,time_step)*sign(1.,ua)*( &
4325 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
4328 LOGICAL :: specified
4331 if(config_flags%specified .or. config_flags%nested) specified = .true.
4333 ! set order for the advection scheme
4336 horz_order = config_flags%h_sca_adv_order
4337 vert_order = config_flags%v_sca_adv_order
4339 ! here is the choice of flux operators
4341 ! begin with horizontal flux divergence
4343 horizontal_order_test : IF( horz_order == 6 ) THEN
4345 ! determine boundary mods for flux operators
4346 ! We degrade the flux operators from 3rd/4th order
4347 ! to second order one gridpoint in from the boundaries for
4348 ! all boundary conditions except periodic and symmetry - these
4349 ! conditions have boundary zone data fill for correct application
4350 ! of the higher order flux stencils
4357 IF( config_flags%periodic_x .or. &
4358 config_flags%symmetric_xs .or. &
4359 (its > ids+3) ) degrade_xs = .false.
4360 IF( config_flags%periodic_x .or. &
4361 config_flags%symmetric_xe .or. &
4362 (ite < ide-3) ) degrade_xe = .false.
4363 IF( config_flags%periodic_y .or. &
4364 config_flags%symmetric_ys .or. &
4365 (jts > jds+3) ) degrade_ys = .false.
4366 IF( config_flags%periodic_y .or. &
4367 config_flags%symmetric_ye .or. &
4368 (jte < jde-4) ) degrade_ye = .false.
4370 !--------------- y - advection first
4373 i_end = MIN(ite,ide-1)
4375 j_end = MIN(jte,jde-1)
4377 ! higher order flux has a 5 or 7 point stencil, so compute
4378 ! bounds so we can switch to second order flux close to the boundary
4384 j_start = MAX(jts,jds+1)
4389 j_end = MIN(jte,jde-2)
4393 IF(config_flags%polar) j_end = MIN(jte,jde-1)
4395 ! compute fluxes, 5th or 6th order
4400 j_loop_y_flux_6 : DO j = j_start, j_end+1
4402 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
4405 DO i = i_start, i_end
4406 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4407 fqy( i, k, jp1 ) = vel*flux6( &
4408 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4409 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4414 DO i = i_start, i_end
4415 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4416 fqy( i, k, jp1 ) = vel*flux6( &
4417 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4418 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4421 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
4424 DO i = i_start, i_end
4425 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4426 (w(i,k,j)+w(i,k,j-1))
4431 DO i = i_start, i_end
4432 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4433 (w(i,k,j)+w(i,k,j-1))
4436 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
4439 DO i = i_start, i_end
4440 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4441 fqy( i, k, jp1 ) = vel*flux4( &
4442 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4447 DO i = i_start, i_end
4448 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4449 fqy( i, k, jp1 ) = vel*flux4( &
4450 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4453 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
4456 DO i = i_start, i_end
4457 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4458 (w(i,k,j)+w(i,k,j-1))
4463 DO i = i_start, i_end
4464 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4465 (w(i,k,j)+w(i,k,j-1))
4468 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
4471 DO i = i_start, i_end
4472 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4473 fqy( i, k, jp1 ) = vel*flux4( &
4474 w(i,k,j-2),w(i,k,j-1), &
4475 w(i,k,j),w(i,k,j+1),vel )
4480 DO i = i_start, i_end
4481 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4482 fqy( i, k, jp1 ) = vel*flux4( &
4483 w(i,k,j-2),w(i,k,j-1), &
4484 w(i,k,j),w(i,k,j+1),vel )
4489 ! y flux-divergence into tendency
4491 ! Comments for polar boundary conditions
4492 ! Same process as for advect_u - tendencies run from jds to jde-1
4493 ! (latitudes are as for u grid, longitudes are displaced)
4494 ! Therefore: flow is only from one side for points next to poles
4495 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
4497 DO i = i_start, i_end
4498 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4499 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
4502 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
4504 DO i = i_start, i_end
4505 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4506 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
4511 IF(j > j_start) THEN
4514 DO i = i_start, i_end
4515 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4516 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
4528 ENDDO j_loop_y_flux_6
4530 ! next, x - flux divergence
4533 i_end = MIN(ite,ide-1)
4536 j_end = MIN(jte,jde-1)
4538 ! higher order flux has a 5 or 7 point stencil, so compute
4539 ! bounds so we can switch to second order flux close to the boundary
4545 i_start = MAX(ids+1,its)
4546 ! i_start_f = i_start+2
4547 i_start_f = MIN(i_start+2,ids+3)
4551 i_end = MIN(ide-2,ite)
4557 DO j = j_start, j_end
4559 ! 5th or 6th order flux
4562 DO i = i_start_f, i_end_f
4563 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4564 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &
4565 w(i-1,k,j), w(i ,k,j), &
4566 w(i+1,k,j), w(i+2,k,j), &
4572 DO i = i_start_f, i_end_f
4573 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4574 fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), &
4575 w(i-1,k,j), w(i ,k,j), &
4576 w(i+1,k,j), w(i+2,k,j), &
4580 ! lower order fluxes close to boundaries (if not periodic or symmetric)
4582 IF( degrade_xs ) THEN
4584 DO i=i_start,i_start_f-1
4586 IF(i == ids+1) THEN ! second order
4588 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4589 *(w(i,k,j)+w(i-1,k,j))
4592 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4593 *(w(i,k,j)+w(i-1,k,j))
4596 IF(i == ids+2) THEN ! third order
4598 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4599 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4600 w(i ,k,j), w(i+1,k,j), &
4604 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4605 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4606 w(i ,k,j), w(i+1,k,j), &
4614 IF( degrade_xe ) THEN
4616 DO i = i_end_f+1, i_end+1
4618 IF( i == ide-1 ) THEN ! second order flux next to the boundary
4620 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4621 *(w(i,k,j)+w(i-1,k,j))
4624 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4625 *(w(i,k,j)+w(i-1,k,j))
4628 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
4630 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4631 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4632 w(i ,k,j), w(i+1,k,j), &
4636 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4637 fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
4638 w(i ,k,j), w(i+1,k,j), &
4646 ! x flux-divergence into tendency
4649 DO i = i_start, i_end
4650 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
4651 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
4657 ELSE IF (horz_order == 5 ) THEN
4659 ! determine boundary mods for flux operators
4660 ! We degrade the flux operators from 3rd/4th order
4661 ! to second order one gridpoint in from the boundaries for
4662 ! all boundary conditions except periodic and symmetry - these
4663 ! conditions have boundary zone data fill for correct application
4664 ! of the higher order flux stencils
4671 IF( config_flags%periodic_x .or. &
4672 config_flags%symmetric_xs .or. &
4673 (its > ids+3) ) degrade_xs = .false.
4674 IF( config_flags%periodic_x .or. &
4675 config_flags%symmetric_xe .or. &
4676 (ite < ide-3) ) degrade_xe = .false.
4677 IF( config_flags%periodic_y .or. &
4678 config_flags%symmetric_ys .or. &
4679 (jts > jds+3) ) degrade_ys = .false.
4680 IF( config_flags%periodic_y .or. &
4681 config_flags%symmetric_ye .or. &
4682 (jte < jde-4) ) degrade_ye = .false.
4684 !--------------- y - advection first
4687 i_end = MIN(ite,ide-1)
4689 j_end = MIN(jte,jde-1)
4691 ! higher order flux has a 5 or 7 point stencil, so compute
4692 ! bounds so we can switch to second order flux close to the boundary
4698 j_start = MAX(jts,jds+1)
4703 j_end = MIN(jte,jde-2)
4707 IF(config_flags%polar) j_end = MIN(jte,jde-1)
4709 ! compute fluxes, 5th or 6th order
4714 j_loop_y_flux_5 : DO j = j_start, j_end+1
4716 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
4719 DO i = i_start, i_end
4720 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4721 fqy( i, k, jp1 ) = vel*flux5( &
4722 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4723 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4728 DO i = i_start, i_end
4729 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4730 fqy( i, k, jp1 ) = vel*flux5( &
4731 w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), &
4732 w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel )
4735 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
4738 DO i = i_start, i_end
4739 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4740 (w(i,k,j)+w(i,k,j-1))
4745 DO i = i_start, i_end
4746 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4747 (w(i,k,j)+w(i,k,j-1))
4750 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
4753 DO i = i_start, i_end
4754 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4755 fqy( i, k, jp1 ) = vel*flux3( &
4756 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4761 DO i = i_start, i_end
4762 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4763 fqy( i, k, jp1 ) = vel*flux3( &
4764 w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
4767 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
4770 DO i = i_start, i_end
4771 fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* &
4772 (w(i,k,j)+w(i,k,j-1))
4777 DO i = i_start, i_end
4778 fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
4779 (w(i,k,j)+w(i,k,j-1))
4782 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
4785 DO i = i_start, i_end
4786 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
4787 fqy( i, k, jp1 ) = vel*flux3( &
4788 w(i,k,j-2),w(i,k,j-1), &
4789 w(i,k,j),w(i,k,j+1),vel )
4794 DO i = i_start, i_end
4795 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
4796 fqy( i, k, jp1 ) = vel*flux3( &
4797 w(i,k,j-2),w(i,k,j-1), &
4798 w(i,k,j),w(i,k,j+1),vel )
4803 ! y flux-divergence into tendency
4805 ! Comments for polar boundary conditions
4806 ! Same process as for advect_u - tendencies run from jds to jde-1
4807 ! (latitudes are as for u grid, longitudes are displaced)
4808 ! Therefore: flow is only from one side for points next to poles
4809 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
4811 DO i = i_start, i_end
4812 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4813 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
4816 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
4818 DO i = i_start, i_end
4819 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4820 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
4825 IF(j > j_start) THEN
4828 DO i = i_start, i_end
4829 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
4830 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
4842 ENDDO j_loop_y_flux_5
4844 ! next, x - flux divergence
4847 i_end = MIN(ite,ide-1)
4850 j_end = MIN(jte,jde-1)
4852 ! higher order flux has a 5 or 7 point stencil, so compute
4853 ! bounds so we can switch to second order flux close to the boundary
4859 i_start = MAX(ids+1,its)
4860 ! i_start_f = i_start+2
4861 i_start_f = MIN(i_start+2,ids+3)
4865 i_end = MIN(ide-2,ite)
4871 DO j = j_start, j_end
4873 ! 5th or 6th order flux
4876 DO i = i_start_f, i_end_f
4877 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4878 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &
4879 w(i-1,k,j), w(i ,k,j), &
4880 w(i+1,k,j), w(i+2,k,j), &
4886 DO i = i_start_f, i_end_f
4887 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4888 fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), &
4889 w(i-1,k,j), w(i ,k,j), &
4890 w(i+1,k,j), w(i+2,k,j), &
4894 ! lower order fluxes close to boundaries (if not periodic or symmetric)
4896 IF( degrade_xs ) THEN
4898 DO i=i_start,i_start_f-1
4900 IF(i == ids+1) THEN ! second order
4902 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4903 *(w(i,k,j)+w(i-1,k,j))
4906 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4907 *(w(i,k,j)+w(i-1,k,j))
4910 IF(i == ids+2) THEN ! third order
4912 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4913 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
4914 w(i ,k,j), w(i+1,k,j), &
4918 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4919 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
4920 w(i ,k,j), w(i+1,k,j), &
4928 IF( degrade_xe ) THEN
4930 DO i = i_end_f+1, i_end+1
4932 IF( i == ide-1 ) THEN ! second order flux next to the boundary
4934 fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
4935 *(w(i,k,j)+w(i-1,k,j))
4938 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
4939 *(w(i,k,j)+w(i-1,k,j))
4942 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
4944 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
4945 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
4946 w(i ,k,j), w(i+1,k,j), &
4950 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
4951 fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
4952 w(i ,k,j), w(i+1,k,j), &
4960 ! x flux-divergence into tendency
4963 DO i = i_start, i_end
4964 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
4965 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
4971 ELSE IF ( horz_order == 4 ) THEN
4978 IF( config_flags%periodic_x .or. &
4979 config_flags%symmetric_xs .or. &
4980 (its > ids+2) ) degrade_xs = .false.
4981 IF( config_flags%periodic_x .or. &
4982 config_flags%symmetric_xe .or. &
4983 (ite < ide-2) ) degrade_xe = .false.
4984 IF( config_flags%periodic_y .or. &
4985 config_flags%symmetric_ys .or. &
4986 (jts > jds+2) ) degrade_ys = .false.
4987 IF( config_flags%periodic_y .or. &
4988 config_flags%symmetric_ye .or. &
4989 (jte < jde-3) ) degrade_ye = .false.
4991 ! begin flux computations
4992 ! start with x flux divergence
4999 i_end = MIN(ite,ide-1)
5001 j_end = MIN(jte,jde-1)
5003 ! 3rd or 4th order flux has a 5 point stencil, so compute
5004 ! bounds so we can switch to second order flux close to the boundary
5011 i_start_f = i_start+1
5021 DO j = j_start, j_end
5024 DO i = i_start_f, i_end_f
5025 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5026 fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
5027 w(i ,k,j), w(i+1,k,j), &
5033 DO i = i_start_f, i_end_f
5034 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5035 fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), &
5036 w(i ,k,j), w(i+1,k,j), &
5039 ! second order flux close to boundaries (if not periodic or symmetric)
5041 IF( degrade_xs ) THEN
5044 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) &
5045 *(w(i_start,k,j)+w(i_start-1,k,j))
5049 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) &
5050 *(w(i_start,k,j)+w(i_start-1,k,j))
5053 IF( degrade_xe ) THEN
5056 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) &
5057 *(w(i_end+1,k,j)+w(i_end,k,j))
5061 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) &
5062 *(w(i_end+1,k,j)+w(i_end,k,j))
5065 ! x flux-divergence into tendency
5068 DO i = i_start, i_end
5069 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5070 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
5076 ! next -> y flux divergence calculation
5079 i_end = MIN(ite,ide-1)
5081 j_end = MIN(jte,jde-1)
5084 ! 3rd or 4th order flux has a 5 point stencil, so compute
5085 ! bounds so we can switch to second order flux close to the boundary
5092 j_start_f = j_start+1
5100 IF(config_flags%polar) j_end = MIN(jte,jde-1)
5105 DO j = j_start, j_end+1
5107 IF ((j < j_start_f) .and. degrade_ys) THEN
5109 DO i = i_start, i_end
5111 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) &
5112 *(w(i,k,j_start)+w(i,k,j_start-1))
5116 DO i = i_start, i_end
5118 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) &
5119 *(w(i,k,j_start)+w(i,k,j_start-1))
5121 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
5123 DO i = i_start, i_end
5124 ! Assumes j>j_end_f is ONLY j_end+1 ...
5125 ! fqy(i, k, jp1) = &
5126 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) &
5127 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5129 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5130 *(w(i,k,j)+w(i,k,j-1))
5134 DO i = i_start, i_end
5135 ! Assumes j>j_end_f is ONLY j_end+1 ...
5136 ! fqy(i, k, jp1) = &
5137 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &
5138 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5140 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5141 *(w(i,k,j)+w(i,k,j-1))
5144 ! 3rd or 4th order flux
5146 DO i = i_start, i_end
5147 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
5148 fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), &
5149 w(i,k,j ), w(i,k,j+1), &
5154 DO i = i_start, i_end
5155 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
5156 fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), &
5157 w(i,k,j ), w(i,k,j+1), &
5162 ! y flux-divergence into tendency
5164 ! Comments for polar boundary conditions
5165 ! Same process as for advect_u - tendencies run from jds to jde-1
5166 ! (latitudes are as for u grid, longitudes are displaced)
5167 ! Therefore: flow is only from one side for points next to poles
5168 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
5170 DO i = i_start, i_end
5171 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5172 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
5175 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
5177 DO i = i_start, i_end
5178 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5179 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
5184 IF( j > j_start ) THEN
5187 DO i = i_start, i_end
5188 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5189 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
5203 ELSE IF ( horz_order == 3 ) THEN
5210 IF( config_flags%periodic_x .or. &
5211 config_flags%symmetric_xs .or. &
5212 (its > ids+2) ) degrade_xs = .false.
5213 IF( config_flags%periodic_x .or. &
5214 config_flags%symmetric_xe .or. &
5215 (ite < ide-2) ) degrade_xe = .false.
5216 IF( config_flags%periodic_y .or. &
5217 config_flags%symmetric_ys .or. &
5218 (jts > jds+2) ) degrade_ys = .false.
5219 IF( config_flags%periodic_y .or. &
5220 config_flags%symmetric_ye .or. &
5221 (jte < jde-3) ) degrade_ye = .false.
5223 ! begin flux computations
5224 ! start with x flux divergence
5231 i_end = MIN(ite,ide-1)
5233 j_end = MIN(jte,jde-1)
5235 ! 3rd or 4th order flux has a 5 point stencil, so compute
5236 ! bounds so we can switch to second order flux close to the boundary
5243 i_start_f = i_start+1
5253 DO j = j_start, j_end
5256 DO i = i_start_f, i_end_f
5257 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
5258 fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5259 w(i ,k,j), w(i+1,k,j), &
5264 DO i = i_start_f, i_end_f
5265 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
5266 fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), &
5267 w(i ,k,j), w(i+1,k,j), &
5271 ! second order flux close to boundaries (if not periodic or symmetric)
5273 IF( degrade_xs ) THEN
5276 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) &
5277 *(w(i_start,k,j)+w(i_start-1,k,j))
5281 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) &
5282 *(w(i_start,k,j)+w(i_start-1,k,j))
5285 IF( degrade_xe ) THEN
5288 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) &
5289 *(w(i_end+1,k,j)+w(i_end,k,j))
5293 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) &
5294 *(w(i_end+1,k,j)+w(i_end,k,j))
5297 ! x flux-divergence into tendency
5300 DO i = i_start, i_end
5301 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5302 tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
5308 ! next -> y flux divergence calculation
5311 i_end = MIN(ite,ide-1)
5313 j_end = MIN(jte,jde-1)
5316 ! 3rd or 4th order flux has a 5 point stencil, so compute
5317 ! bounds so we can switch to second order flux close to the boundary
5324 j_start_f = j_start+1
5332 IF(config_flags%polar) j_end = MIN(jte,jde-1)
5337 DO j = j_start, j_end+1
5339 IF ((j < j_start_f) .and. degrade_ys) THEN
5341 DO i = i_start, i_end
5343 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) &
5344 *(w(i,k,j_start)+w(i,k,j_start-1))
5348 DO i = i_start, i_end
5350 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) &
5351 *(w(i,k,j_start)+w(i,k,j_start-1))
5353 ELSE IF ((j > j_end_f) .and. degrade_ye) THEN
5355 DO i = i_start, i_end
5356 ! Assumes j>j_end_f is ONLY j_end+1 ...
5357 ! fqy(i, k, jp1) = &
5358 ! 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) &
5359 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5361 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5362 *(w(i,k,j)+w(i,k,j-1))
5366 DO i = i_start, i_end
5367 ! Assumes j>j_end_f is ONLY j_end+1 ...
5368 ! fqy(i, k, jp1) = &
5369 ! 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) &
5370 ! *(w(i,k,j_end+1)+w(i,k,j_end))
5372 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5373 *(w(i,k,j)+w(i,k,j-1))
5376 ! 3rd or 4th order flux
5378 DO i = i_start, i_end
5379 vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
5380 fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), &
5381 w(i,k,j ), w(i,k,j+1), &
5386 DO i = i_start, i_end
5387 vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
5388 fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), &
5389 w(i,k,j ), w(i,k,j+1), &
5394 ! y flux-divergence into tendency
5396 ! Comments for polar boundary conditions
5397 ! Same process as for advect_u - tendencies run from jds to jde-1
5398 ! (latitudes are as for u grid, longitudes are displaced)
5399 ! Therefore: flow is only from one side for points next to poles
5400 IF ( config_flags%polar .AND. (j == jds+1) ) THEN
5402 DO i = i_start, i_end
5403 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5404 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*fqy(i,k,jp1)
5407 ELSE IF( config_flags%polar .AND. (j == jde) ) THEN
5409 DO i = i_start, i_end
5410 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5411 tendency(i,k,j-1) = tendency(i,k,j-1) + mrdy*fqy(i,k,jp0)
5416 IF( j > j_start ) THEN
5419 DO i = i_start, i_end
5420 mrdy=msftx(i,j-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5421 tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
5435 ELSE IF (horz_order == 2 ) THEN
5438 i_end = MIN(ite,ide-1)
5440 j_end = MIN(jte,jde-1)
5442 IF ( .NOT. config_flags%periodic_x ) THEN
5443 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
5444 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
5447 DO j = j_start, j_end
5449 DO i = i_start, i_end
5451 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5453 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
5454 *((fzm(k)*ru(i+1,k,j)+fzp(k)*ru(i+1,k-1,j)) &
5455 *(w(i+1,k,j)+w(i,k,j)) &
5456 -(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
5457 *(w(i,k,j)+w(i-1,k,j)))
5463 DO i = i_start, i_end
5465 mrdx=msftx(i,j)*rdx ! see ADT eqn 46 dividing by my, 1st term RHS
5467 tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 &
5468 *(((2.-fzm(k-1))*ru(i+1,k-1,j)-fzp(k-1)*ru(i+1,k-2,j)) &
5469 *(w(i+1,k,j)+w(i,k,j)) &
5470 -((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
5471 *(w(i,k,j)+w(i-1,k,j)))
5478 i_end = MIN(ite,ide-1)
5479 ! Polar boundary conditions are like open or specified
5480 IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
5481 IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end = MIN(jde-2,jte)
5483 DO j = j_start, j_end
5485 DO i = i_start, i_end
5487 mrdy=msftx(i,j)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5489 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
5490 *((fzm(k)*rv(i,k,j+1)+fzp(k)*rv(i,k-1,j+1))* &
5491 (w(i,k,j+1)+w(i,k,j)) &
5492 -(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) &
5493 *(w(i,k,j)+w(i,k,j-1)))
5499 DO i = i_start, i_end
5501 mrdy=msftx(i,j)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5503 tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 &
5504 *(((2.-fzm(k-1))*rv(i,k-1,j+1)-fzp(k-1)*rv(i,k-2,j+1))* &
5505 (w(i,k,j+1)+w(i,k,j)) &
5506 -((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) &
5507 *(w(i,k,j)+w(i,k,j-1)))
5513 ! Polar boundary condition ... not covered in above j-loop
5514 IF (config_flags%polar) THEN
5515 IF (jts == jds) THEN
5517 DO i = i_start, i_end
5518 mrdy=msftx(i,jds)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5519 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
5520 *((fzm(k)*rv(i,k,jds+1)+fzp(k)*rv(i,k-1,jds+1))* &
5521 (w(i,k,jds+1)+w(i,k,jds)))
5525 DO i = i_start, i_end
5526 mrdy=msftx(i,jds)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5527 tendency(i,k,jds)=tendency(i,k,jds) -mrdy*0.5 &
5528 *((2.-fzm(k-1))*rv(i,k-1,jds+1)-fzp(k-1)*rv(i,k-2,jds+1))* &
5529 (w(i,k,jds+1)+w(i,k,jds))
5532 IF (jte == jde) THEN
5534 DO i = i_start, i_end
5535 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5536 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
5537 *((fzm(k)*rv(i,k,jde-1)+fzp(k)*rv(i,k-1,jde-1))* &
5538 (w(i,k,jde-1)+w(i,k,jde-2)))
5542 DO i = i_start, i_end
5543 mrdy=msftx(i,jde-1)*rdy ! see ADT eqn 46 dividing by my, 2nd term RHS
5544 tendency(i,k,jde-1)=tendency(i,k,jde-1) +mrdy*0.5 &
5545 *((2.-fzm(k-1))*rv(i,k-1,jde-1)-fzp(k-1)*rv(i,k-2,jde-1)) &
5546 *(w(i,k,jde-1)+w(i,k,jde-2))
5551 ELSE IF ( horz_order == 0 ) THEN
5553 ! Just in case we want to turn horizontal advection off, we can do it
5557 WRITE ( wrf_err_message ,*) ' advect_w_6a, h_order not known ',horz_order
5558 CALL wrf_error_fatal ( wrf_err_message )
5560 ENDIF horizontal_order_test
5563 ! pick up the the horizontal radiation boundary conditions.
5564 ! (these are the computations that don't require 'cb'.
5565 ! first, set to index ranges
5569 i_end = MIN(ite,ide-1)
5571 j_end = MIN(jte,jde-1)
5573 IF( (config_flags%open_xs) .and. (its == ids)) THEN
5575 DO j = j_start, j_end
5578 uw = 0.5*(fzm(k)*(ru(its,k ,j)+ru(its+1,k ,j)) + &
5579 fzp(k)*(ru(its,k-1,j)+ru(its+1,k-1,j)) )
5582 tendency(its,k,j) = tendency(its,k,j) &
5584 ub*(w_old(its+1,k,j) - w_old(its,k,j)) + &
5586 fzm(k)*(ru(its+1,k ,j)-ru(its,k ,j))+ &
5587 fzp(k)*(ru(its+1,k-1,j)-ru(its,k-1,j))) &
5593 DO j = j_start, j_end
5595 uw = 0.5*( (2.-fzm(k-1))*(ru(its,k-1,j)+ru(its+1,k-1,j)) &
5596 -fzp(k-1)*(ru(its,k-2,j)+ru(its+1,k-2,j)) )
5599 tendency(its,k,j) = tendency(its,k,j) &
5601 ub*(w_old(its+1,k,j) - w_old(its,k,j)) + &
5603 (2.-fzm(k-1))*(ru(its+1,k-1,j)-ru(its,k-1,j))- &
5604 fzp(k-1)*(ru(its+1,k-2,j)-ru(its,k-2,j))) &
5610 IF( (config_flags%open_xe) .and. (ite == ide)) THEN
5612 DO j = j_start, j_end
5615 uw = 0.5*(fzm(k)*(ru(ite-1,k ,j)+ru(ite,k ,j)) + &
5616 fzp(k)*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) )
5619 tendency(i_end,k,j) = tendency(i_end,k,j) &
5621 ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + &
5623 fzm(k)*(ru(ite,k ,j)-ru(ite-1,k ,j)) + &
5624 fzp(k)*(ru(ite,k-1,j)-ru(ite-1,k-1,j))) &
5630 DO j = j_start, j_end
5632 uw = 0.5*( (2.-fzm(k-1))*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) &
5633 -fzp(k-1)*(ru(ite-1,k-2,j)+ru(ite,k-2,j)) )
5636 tendency(i_end,k,j) = tendency(i_end,k,j) &
5638 ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + &
5640 (2.-fzm(k-1))*(ru(ite,k-1,j)-ru(ite-1,k-1,j)) - &
5641 fzp(k-1)*(ru(ite,k-2,j)-ru(ite-1,k-2,j))) &
5648 IF( (config_flags%open_ys) .and. (jts == jds)) THEN
5650 DO i = i_start, i_end
5653 vw = 0.5*( fzm(k)*(rv(i,k ,jts)+rv(i,k ,jts+1)) + &
5654 fzp(k)*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) )
5657 tendency(i,k,jts) = tendency(i,k,jts) &
5659 vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + &
5661 fzm(k)*(rv(i,k ,jts+1)-rv(i,k ,jts))+ &
5662 fzp(k)*(rv(i,k-1,jts+1)-rv(i,k-1,jts))) &
5668 DO i = i_start, i_end
5669 vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) &
5670 -fzp(k-1)*(rv(i,k-2,jts)+rv(i,k-2,jts+1)) )
5673 tendency(i,k,jts) = tendency(i,k,jts) &
5675 vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + &
5677 (2.-fzm(k-1))*(rv(i,k-1,jts+1)-rv(i,k-1,jts))- &
5678 fzp(k-1)*(rv(i,k-2,jts+1)-rv(i,k-2,jts))) &
5684 IF( (config_flags%open_ye) .and. (jte == jde) ) THEN
5686 DO i = i_start, i_end
5689 vw = 0.5*( fzm(k)*(rv(i,k ,jte-1)+rv(i,k ,jte)) + &
5690 fzp(k)*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) )
5693 tendency(i,k,j_end) = tendency(i,k,j_end) &
5695 vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + &
5697 fzm(k)*(rv(i,k ,jte)-rv(i,k ,jte-1))+ &
5698 fzp(k)*(rv(i,k-1,jte)-rv(i,k-1,jte-1))) &
5704 DO i = i_start, i_end
5706 vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) &
5707 -fzp(k-1)*(rv(i,k-2,jte-1)+rv(i,k-2,jte)) )
5710 tendency(i,k,j_end) = tendency(i,k,j_end) &
5712 vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + &
5714 (2.-fzm(k-1))*(rv(i,k-1,jte)-rv(i,k-1,jte-1))- &
5715 fzp(k-1)*(rv(i,k-2,jte)-rv(i,k-2,jte-1))) &
5721 !-------------------- vertical advection
5722 ! ADT eqn 46 has 3rd term on RHS (dividing through by my) = - partial d/dz (w rho w /my)
5723 ! Here we have: - partial d/dz (w*rom) = - partial d/dz (w rho w / my)
5724 ! Therefore we don't need to make a correction for advect_w
5727 i_end = MIN(ite,ide-1)
5729 j_end = MIN(jte,jde-1)
5731 DO i = i_start, i_end
5736 vert_order_test : IF (vert_order == 6) THEN
5738 DO j = j_start, j_end
5741 DO i = i_start, i_end
5742 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5743 vflux(i,k) = vel*flux6( &
5744 w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), &
5745 w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel )
5749 DO i = i_start, i_end
5752 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5755 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5756 vflux(i,k) = vel*flux4( &
5757 w(i,k-2,j), w(i,k-1,j), &
5758 w(i,k ,j), w(i,k+1,j), -vel )
5761 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5762 vflux(i,k) = vel*flux4( &
5763 w(i,k-2,j), w(i,k-1,j), &
5764 w(i,k ,j), w(i,k+1,j), -vel )
5767 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5772 DO i = i_start, i_end
5773 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5777 ! pick up flux contribution for w at the lid. wcs, 13 march 2004
5779 DO i = i_start, i_end
5780 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5785 ELSE IF (vert_order == 5) THEN
5787 DO j = j_start, j_end
5790 DO i = i_start, i_end
5791 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5792 vflux(i,k) = vel*flux5( &
5793 w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), &
5794 w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel )
5798 DO i = i_start, i_end
5801 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5804 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5805 vflux(i,k) = vel*flux3( &
5806 w(i,k-2,j), w(i,k-1,j), &
5807 w(i,k ,j), w(i,k+1,j), -vel )
5809 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5810 vflux(i,k) = vel*flux3( &
5811 w(i,k-2,j), w(i,k-1,j), &
5812 w(i,k ,j), w(i,k+1,j), -vel )
5815 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5820 DO i = i_start, i_end
5821 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5825 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5827 DO i = i_start, i_end
5828 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5833 ELSE IF (vert_order == 4) THEN
5835 DO j = j_start, j_end
5838 DO i = i_start, i_end
5839 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5840 vflux(i,k) = vel*flux4( &
5841 w(i,k-2,j), w(i,k-1,j), &
5842 w(i,k ,j), w(i,k+1,j), -vel )
5846 DO i = i_start, i_end
5849 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5851 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5856 DO i = i_start, i_end
5857 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5861 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5863 DO i = i_start, i_end
5864 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5869 ELSE IF (vert_order == 3) THEN
5871 DO j = j_start, j_end
5874 DO i = i_start, i_end
5875 vel=0.5*(rom(i,k,j)+rom(i,k-1,j))
5876 vflux(i,k) = vel*flux3( &
5877 w(i,k-2,j), w(i,k-1,j), &
5878 w(i,k ,j), w(i,k+1,j), -vel )
5882 DO i = i_start, i_end
5885 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5887 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5892 DO i = i_start, i_end
5893 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5897 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5899 DO i = i_start, i_end
5900 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5905 ELSE IF (vert_order == 2) THEN
5907 DO j = j_start, j_end
5909 DO i = i_start, i_end
5911 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j))
5915 DO i = i_start, i_end
5916 tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k))
5921 ! pick up flux contribution for w at the lid, wcs. 13 march 2004
5923 DO i = i_start, i_end
5924 tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k))
5931 WRITE (wrf_err_message ,*) ' advect_w, v_order not known ',vert_order
5932 CALL wrf_error_fatal ( wrf_err_message )
5934 ENDIF vert_order_test
5936 END SUBROUTINE advect_w
5938 !----------------------------------------------------------------
5940 SUBROUTINE advect_scalar_pd ( field, field_old, tendency, &
5941 h_tendency, z_tendency, &
5944 time_step, config_flags, &
5946 msfux, msfuy, msfvx, msfvy, &
5949 rdx, rdy, rdzw, dt, &
5950 ids, ide, jds, jde, kds, kde, &
5951 ims, ime, jms, jme, kms, kme, &
5952 its, ite, jts, jte, kts, kte )
5954 ! this is a first cut at a positive definite advection option
5955 ! for scalars in WRF. This version is memory intensive ->
5956 ! we save 3d arrays of x, y and z both high and low order fluxes
5957 ! (six in all). Alternatively, we could sweep in a direction
5958 ! and lower the cost considerably.
5960 ! uses the Smolarkiewicz MWR 1989 approach, with addition of first-order
5963 ! WCS, 3 December 2002, 24 February 2003
5969 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5971 LOGICAL , INTENT(IN ) :: tenddec ! tendency flag
5973 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5974 ims, ime, jms, jme, kms, kme, &
5975 its, ite, jts, jte, kts, kte
5977 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
5983 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, mub, mu_old
5984 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
5985 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: h_tendency, z_tendency
5987 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
5994 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5998 REAL , INTENT(IN ) :: rdx, &
6001 INTEGER , INTENT(IN ) :: time_step
6005 INTEGER :: i, j, k, itf, jtf, ktf
6006 INTEGER :: i_start, i_end, j_start, j_end
6007 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
6008 INTEGER :: jmin, jmax, jp, jm, imin, imax
6010 REAL :: mrdx, mrdy, ub, vb, uw, vw, mu
6012 ! storage for high and low order fluxes
6014 REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqx, fqy, fqz
6015 REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqxl, fqyl, fqzl
6017 INTEGER :: horz_order, vert_order
6019 LOGICAL :: degrade_xs, degrade_ys
6020 LOGICAL :: degrade_xe, degrade_ye
6022 INTEGER :: jp1, jp0, jtmp
6024 REAL :: flux_out, ph_low, scale
6025 REAL, PARAMETER :: eps=1.e-20
6028 ! definition of flux operators, 3rd, 4th, 5th or 6th order
6030 REAL :: flux3, flux4, flux5, flux6, flux_upwind
6031 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr
6033 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
6034 (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2)
6036 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
6037 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
6038 sign(1,time_step)*sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1))
6040 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
6041 (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2) &
6042 +(1./60.)*(q_ip2+q_im3)
6044 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
6045 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
6046 -sign(1,time_step)*sign(1.,ua)*(1./60.)*( &
6047 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )
6049 flux_upwind(q_im1, q_i, cr ) = 0.5*min( 1.0,(cr+abs(cr)))*q_im1 &
6050 +0.5*max(-1.0,(cr-abs(cr)))*q_i
6052 ! flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 &
6053 ! +0.5*(1.-sign(1.,cr))*q_i
6054 ! flux_upwind(q_im1, q_i, cr ) = 0.
6058 LOGICAL, PARAMETER :: pd_limit = .true.
6060 ! set order for the advection schemes
6062 ! write(6,*) ' in pd advection routine '
6064 ! Empty arrays just in case:
6065 IF (config_flags%polar) THEN
6075 horz_order = config_flags%h_sca_adv_order
6076 vert_order = config_flags%v_sca_adv_order
6078 ! determine boundary mods for flux operators
6079 ! We degrade the flux operators from 3rd/4th order
6080 ! to second order one gridpoint in from the boundaries for
6081 ! all boundary conditions except periodic and symmetry - these
6082 ! conditions have boundary zone data fill for correct application
6083 ! of the higher order flux stencils
6090 ! begin with horizontal flux divergence
6091 ! here is the choice of flux operators
6094 horizontal_order_test : IF( horz_order == 6 ) THEN
6096 IF( config_flags%periodic_x .or. &
6097 config_flags%symmetric_xs .or. &
6098 (its > ids+3) ) degrade_xs = .false.
6099 IF( config_flags%periodic_x .or. &
6100 config_flags%symmetric_xe .or. &
6101 (ite < ide-4) ) degrade_xe = .false.
6102 IF( config_flags%periodic_y .or. &
6103 config_flags%symmetric_ys .or. &
6104 (jts > jds+3) ) degrade_ys = .false.
6105 IF( config_flags%periodic_y .or. &
6106 config_flags%symmetric_ye .or. &
6107 (jte < jde-4) ) degrade_ye = .false.
6109 !--------------- y - advection first
6111 !-- y flux compute; these bounds are for periodic and sym b.c.
6115 i_end = MIN(ite,ide-1)+1
6117 j_end = MIN(jte,jde-1)+1
6121 !-- modify loop bounds if open or specified
6123 ! IF(degrade_xs) i_start = MAX(its-1,ids-1)
6124 ! IF(degrade_xe) i_end = MIN(ite+1,ide-2)
6125 IF(degrade_xs) i_start = MAX(its-1,ids)
6126 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
6129 j_start = MAX(jts-1,jds+1)
6134 j_end = MIN(jte+1,jde-2)
6138 ! compute fluxes, 6th order
6140 j_loop_y_flux_6 : DO j = j_start, j_end+1
6142 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6145 DO i = i_start, i_end
6147 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6148 mu = 0.5*(mut(i,j)+mut(i,j-1))
6151 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6153 fqy( i, k, j ) = vel*flux6( &
6154 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
6155 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
6157 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6162 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6165 DO i = i_start, i_end
6167 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6168 mu = 0.5*(mut(i,j)+mut(i,j-1))
6171 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6173 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6174 (field(i,k,j)+field(i,k,j-1))
6176 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6181 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
6184 DO i = i_start, i_end
6186 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6187 mu = 0.5*(mut(i,j)+mut(i,j-1))
6190 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6192 fqy( i, k, j ) = vel*flux4( &
6193 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
6194 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6199 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6202 DO i = i_start, i_end
6204 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6205 mu = 0.5*(mut(i,j)+mut(i,j-1))
6208 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6210 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6211 (field(i,k,j)+field(i,k,j-1))
6212 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6217 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
6220 DO i = i_start, i_end
6222 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6223 mu = 0.5*(mut(i,j)+mut(i,j-1))
6226 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6228 fqy( i, k, j) = vel*flux4( &
6229 field(i,k,j-2),field(i,k,j-1), &
6230 field(i,k,j),field(i,k,j+1),vel )
6231 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6238 ENDDO j_loop_y_flux_6
6242 !-- these bounds are for periodic and sym conditions
6245 i_end = MIN(ite,ide-1)+1
6250 j_end = MIN(jte,jde-1)+1
6252 !-- modify loop bounds for open and specified b.c
6254 ! IF(degrade_ys) j_start = MAX(jts-1,jds+1)
6255 ! IF(degrade_ye) j_end = MIN(jte+1,jde-2)
6256 IF(degrade_ys) j_start = MAX(jts-1,jds)
6257 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
6260 i_start = MAX(ids+1,its-1)
6265 i_end = MIN(ide-2,ite+1)
6271 DO j = j_start, j_end
6276 DO i = i_start_f, i_end_f
6278 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6279 mu = 0.5*(mut(i,j)+mut(i-1,j))
6282 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6284 fqx( i,k,j ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), &
6285 field(i-1,k,j), field(i ,k,j), &
6286 field(i+1,k,j), field(i+2,k,j), &
6288 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6293 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6295 IF( degrade_xs ) THEN
6297 DO i=i_start,i_start_f-1
6299 IF(i == ids+1) THEN ! second order
6301 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6302 mu = 0.5*(mut(i,j)+mut(i-1,j))
6305 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6306 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6307 *(field(i,k,j)+field(i-1,k,j))
6308 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6312 IF(i == ids+2) THEN ! fourth order
6314 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6315 mu = 0.5*(mut(i,j)+mut(i-1,j))
6318 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6319 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6320 field(i ,k,j), field(i+1,k,j), &
6322 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6330 IF( degrade_xe ) THEN
6332 DO i = i_end_f+1, i_end+1
6334 IF( i == ide-1 ) THEN ! second order flux next to the boundary
6336 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6337 mu = 0.5*(mut(i,j)+mut(i-1,j))
6340 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6341 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6342 *(field(i,k,j)+field(i-1,k,j))
6343 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6348 IF( i == ide-2 ) THEN ! fourth order flux one in from the boundary
6350 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6351 mu = 0.5*(mut(i,j)+mut(i-1,j))
6354 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6355 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6356 field(i ,k,j), field(i+1,k,j), &
6358 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6366 ENDDO ! enddo for outer J loop
6368 !--- end of 6th order horizontal flux calculation
6370 ELSE IF( horz_order == 5 ) THEN
6372 IF( config_flags%periodic_x .or. &
6373 config_flags%symmetric_xs .or. &
6374 (its > ids+3) ) degrade_xs = .false.
6375 IF( config_flags%periodic_x .or. &
6376 config_flags%symmetric_xe .or. &
6377 (ite < ide-4) ) degrade_xe = .false.
6378 IF( config_flags%periodic_y .or. &
6379 config_flags%symmetric_ys .or. &
6380 (jts > jds+3) ) degrade_ys = .false.
6381 IF( config_flags%periodic_y .or. &
6382 config_flags%symmetric_ye .or. &
6383 (jte < jde-4) ) degrade_ye = .false.
6385 !--------------- y - advection first
6387 !-- y flux compute; these bounds are for periodic and sym b.c.
6391 i_end = MIN(ite,ide-1)+1
6393 j_end = MIN(jte,jde-1)+1
6397 !-- modify loop bounds if open or specified
6399 ! IF(degrade_xs) i_start = MAX(its-1,ids-1)
6400 ! IF(degrade_xe) i_end = MIN(ite+1,ide-2)
6401 IF(degrade_xs) i_start = MAX(its-1,ids)
6402 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
6405 j_start = MAX(jts-1,jds+1)
6410 j_end = MIN(jte+1,jde-2)
6414 ! compute fluxes, 5th order
6416 j_loop_y_flux_5 : DO j = j_start, j_end+1
6418 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6421 DO i = i_start, i_end
6423 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6424 mu = 0.5*(mut(i,j)+mut(i,j-1))
6427 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6429 fqy( i, k, j ) = vel*flux5( &
6430 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
6431 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
6433 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6438 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6441 DO i = i_start, i_end
6443 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6444 mu = 0.5*(mut(i,j)+mut(i,j-1))
6447 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6449 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6450 (field(i,k,j)+field(i,k,j-1))
6452 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6457 ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary
6460 DO i = i_start, i_end
6462 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6463 mu = 0.5*(mut(i,j)+mut(i,j-1))
6466 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6468 fqy( i, k, j ) = vel*flux3( &
6469 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
6470 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6475 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6478 DO i = i_start, i_end
6480 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6481 mu = 0.5*(mut(i,j)+mut(i,j-1))
6484 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6486 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6487 (field(i,k,j)+field(i,k,j-1))
6488 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6493 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary
6496 DO i = i_start, i_end
6498 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6499 mu = 0.5*(mut(i,j)+mut(i,j-1))
6502 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6504 fqy( i, k, j) = vel*flux3( &
6505 field(i,k,j-2),field(i,k,j-1), &
6506 field(i,k,j),field(i,k,j+1),vel )
6507 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6514 ENDDO j_loop_y_flux_5
6518 !-- these bounds are for periodic and sym conditions
6521 i_end = MIN(ite,ide-1)+1
6526 j_end = MIN(jte,jde-1)+1
6528 !-- modify loop bounds for open and specified b.c
6530 ! IF(degrade_ys) j_start = MAX(jts-1,jds+1)
6531 ! IF(degrade_ye) j_end = MIN(jte+1,jde-2)
6532 IF(degrade_ys) j_start = MAX(jts-1,jds)
6533 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
6536 i_start = MAX(ids+1,its-1)
6541 i_end = MIN(ide-2,ite+1)
6547 DO j = j_start, j_end
6552 DO i = i_start_f, i_end_f
6554 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6555 mu = 0.5*(mut(i,j)+mut(i-1,j))
6558 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6560 fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), &
6561 field(i-1,k,j), field(i ,k,j), &
6562 field(i+1,k,j), field(i+2,k,j), &
6564 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6569 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6571 IF( degrade_xs ) THEN
6573 DO i=i_start,i_start_f-1
6575 IF(i == ids+1) THEN ! second order
6577 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6578 mu = 0.5*(mut(i,j)+mut(i-1,j))
6581 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6582 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6583 *(field(i,k,j)+field(i-1,k,j))
6584 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6588 IF(i == ids+2) THEN ! third order
6590 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6591 mu = 0.5*(mut(i,j)+mut(i-1,j))
6594 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6595 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
6596 field(i ,k,j), field(i+1,k,j), &
6598 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6606 IF( degrade_xe ) THEN
6608 DO i = i_end_f+1, i_end+1
6610 IF( i == ide-1 ) THEN ! second order flux next to the boundary
6612 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6613 mu = 0.5*(mut(i,j)+mut(i-1,j))
6616 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6617 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6618 *(field(i,k,j)+field(i-1,k,j))
6619 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6624 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
6626 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6627 mu = 0.5*(mut(i,j)+mut(i-1,j))
6630 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6631 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
6632 field(i ,k,j), field(i+1,k,j), &
6634 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6642 ENDDO ! enddo for outer J loop
6644 !--- end of 5th order horizontal flux calculation
6646 ELSE IF( horz_order == 4 ) THEN
6648 IF( config_flags%periodic_x .or. &
6649 config_flags%symmetric_xs .or. &
6650 (its > ids+1) ) degrade_xs = .false.
6651 IF( config_flags%periodic_x .or. &
6652 config_flags%symmetric_xe .or. &
6653 (ite < ide-2) ) degrade_xe = .false.
6654 IF( config_flags%periodic_y .or. &
6655 config_flags%symmetric_ys .or. &
6656 (jts > jds+1) ) degrade_ys = .false.
6657 IF( config_flags%periodic_y .or. &
6658 config_flags%symmetric_ye .or. &
6659 (jte < jde-2) ) degrade_ye = .false.
6661 !--------------- y - advection first
6663 !-- y flux compute; these bounds are for periodic and sym b.c.
6667 i_end = MIN(ite,ide-1)+1
6669 j_end = MIN(jte,jde-1)+1
6673 !-- modify loop bounds if open or specified
6675 IF(degrade_xs) i_start = its
6676 IF(degrade_xe) i_end = MIN(ite,ide-1)
6679 j_start = MAX(jts,jds+1)
6684 j_end = MIN(jte,jde-2)
6688 ! compute fluxes, 4th order
6690 j_loop_y_flux_4 : DO j = j_start, j_end+1
6692 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6695 DO i = i_start, i_end
6697 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6698 mu = 0.5*(mut(i,j)+mut(i,j-1))
6701 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6703 fqy( i, k, j ) = vel*flux4( field(i,k,j-2), field(i,k,j-1), &
6704 field(i,k,j ), field(i,k,j+1), vel )
6706 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6711 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6714 DO i = i_start, i_end
6716 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6717 mu = 0.5*(mut(i,j)+mut(i,j-1))
6720 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6722 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6723 (field(i,k,j)+field(i,k,j-1))
6725 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6730 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6733 DO i = i_start, i_end
6735 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6736 mu = 0.5*(mut(i,j)+mut(i,j-1))
6739 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6741 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6742 (field(i,k,j)+field(i,k,j-1))
6743 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6750 ENDDO j_loop_y_flux_4
6754 !-- these bounds are for periodic and sym conditions
6757 i_end = MIN(ite,ide-1)+1
6762 j_end = MIN(jte,jde-1)+1
6764 !-- modify loop bounds for open and specified b.c
6766 IF(degrade_ys) j_start = jts
6767 IF(degrade_ye) j_end = MIN(jte,jde-1)
6770 i_start = MAX(ids+1,its)
6771 i_start_f = i_start+1
6775 i_end = MIN(ide-2,ite)
6781 DO j = j_start, j_end
6786 DO i = i_start_f, i_end_f
6788 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6789 mu = 0.5*(mut(i,j)+mut(i-1,j))
6792 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6794 fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), &
6795 field(i ,k,j), field(i+1,k,j), vel )
6796 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6801 ! lower order fluxes close to boundaries (if not periodic or symmetric)
6803 IF( degrade_xs ) THEN
6804 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
6808 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6809 mu = 0.5*(mut(i,j)+mut(i-1,j))
6812 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6814 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6815 *(field(i,k,j)+field(i-1,k,j))
6817 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6823 IF( degrade_xe ) THEN
6824 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
6827 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6828 mu = 0.5*(mut(i,j)+mut(i-1,j))
6831 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6832 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
6833 *(field(i,k,j)+field(i-1,k,j))
6834 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6840 ENDDO ! enddo for outer J loop
6842 !--- end of 4th order horizontal flux calculation
6844 ELSE IF( horz_order == 3 ) THEN
6846 IF( config_flags%periodic_x .or. &
6847 config_flags%symmetric_xs .or. &
6848 (its > ids+2) ) degrade_xs = .false.
6849 IF( config_flags%periodic_x .or. &
6850 config_flags%symmetric_xe .or. &
6851 (ite < ide-1) ) degrade_xe = .false.
6852 IF( config_flags%periodic_y .or. &
6853 config_flags%symmetric_ys .or. &
6854 (jts > jds+2) ) degrade_ys = .false.
6855 IF( config_flags%periodic_y .or. &
6856 config_flags%symmetric_ye .or. &
6857 (jte < jde-1) ) degrade_ye = .false.
6859 !--------------- y - advection first
6861 !-- y flux compute; these bounds are for periodic and sym b.c.
6865 i_end = MIN(ite,ide-1)+1
6867 j_end = MIN(jte,jde-1)+1
6871 !-- modify loop bounds if open or specified
6873 IF(degrade_xs) i_start = its
6874 IF(degrade_xe) i_end = MIN(ite,ide-1)
6877 j_start = MAX(jts,jds+1)
6882 j_end = MIN(jte,jde-2)
6886 ! compute fluxes, 3rd order
6888 j_loop_y_flux_3 : DO j = j_start, j_end+1
6890 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
6893 DO i = i_start, i_end
6895 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6896 mu = 0.5*(mut(i,j)+mut(i,j-1))
6899 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6901 fqy( i, k, j ) = vel*flux3( field(i,k,j-2), field(i,k,j-1), &
6902 field(i,k,j ), field(i,k,j+1), vel )
6904 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6909 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
6912 DO i = i_start, i_end
6914 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6915 mu = 0.5*(mut(i,j)+mut(i,j-1))
6918 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6920 fqy(i,k, j) = 0.5*rv(i,k,j)* &
6921 (field(i,k,j)+field(i,k,j-1))
6923 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6928 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
6931 DO i = i_start, i_end
6933 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
6934 mu = 0.5*(mut(i,j)+mut(i,j-1))
6937 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
6939 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
6940 (field(i,k,j)+field(i,k,j-1))
6941 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
6948 ENDDO j_loop_y_flux_3
6952 !-- these bounds are for periodic and sym conditions
6955 i_end = MIN(ite,ide-1)+1
6960 j_end = MIN(jte,jde-1)+1
6962 !-- modify loop bounds for open and specified b.c
6964 IF(degrade_ys) j_start = jts
6965 IF(degrade_ye) j_end = MIN(jte,jde-1)
6968 i_start = MAX(ids+1,its)
6969 i_start_f = i_start+1
6973 i_end = MIN(ide-2,ite)
6979 DO j = j_start, j_end
6984 DO i = i_start_f, i_end_f
6986 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
6987 mu = 0.5*(mut(i,j)+mut(i-1,j))
6990 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
6992 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
6993 field(i ,k,j), field(i+1,k,j), vel )
6994 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
6999 ! lower order fluxes close to boundaries (if not periodic or symmetric)
7001 IF( degrade_xs ) THEN
7003 IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
7007 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7008 mu = 0.5*(mut(i,j)+mut(i-1,j))
7011 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7013 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
7014 *(field(i,k,j)+field(i-1,k,j))
7016 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7022 IF( degrade_xe ) THEN
7023 IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
7026 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7027 mu = 0.5*(mut(i,j)+mut(i-1,j))
7030 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7031 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
7032 *(field(i,k,j)+field(i-1,k,j))
7033 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7039 ENDDO ! enddo for outer J loop
7041 !--- end of 3rd order horizontal flux calculation
7044 ELSE IF( horz_order == 2 ) THEN
7046 IF( config_flags%periodic_x .or. &
7047 config_flags%symmetric_xs .or. &
7048 (its > ids+1) ) degrade_xs = .false.
7049 IF( config_flags%periodic_x .or. &
7050 config_flags%symmetric_xe .or. &
7051 (ite < ide-2) ) degrade_xe = .false.
7052 IF( config_flags%periodic_y .or. &
7053 config_flags%symmetric_ys .or. &
7054 (jts > jds+1) ) degrade_ys = .false.
7055 IF( config_flags%periodic_y .or. &
7056 config_flags%symmetric_ye .or. &
7057 (jte < jde-2) ) degrade_ye = .false.
7059 !-- y flux compute; these bounds are for periodic and sym b.c.
7063 i_end = MIN(ite,ide-1)+1
7065 j_end = MIN(jte,jde-1)+1
7067 !-- modify loop bounds if open or specified
7069 IF(degrade_xs) i_start = its
7070 IF(degrade_xe) i_end = MIN(ite,ide-1)
7071 IF(degrade_ys) j_start = MAX(jts,jds+1)
7072 IF(degrade_ye) j_end = MIN(jte,jde-2)
7074 ! compute fluxes, 2nd order, y flux
7076 DO j = j_start, j_end+1
7078 DO i = i_start, i_end
7079 dy = 2./(msftx(i,j)+msftx(i,j-1))/rdy ! ADT eqn 48 d/dy
7080 mu = 0.5*(mut(i,j)+mut(i,j-1))
7083 fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
7085 fqy(i,k, j) = 0.5*rv(i,k,j)* &
7086 (field(i,k,j)+field(i,k,j-1))
7088 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7095 DO j = j_start, j_end
7097 DO i = i_start, i_end+1
7098 dx = 2./(msfty(i,j)+msfty(i-1,j))/rdx ! ADT eqn 48 d/dx
7099 mu = 0.5*(mut(i,j)+mut(i-1,j))
7102 fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
7103 fqx( i,k,j ) = 0.5*ru(i,k,j)* &
7104 (field(i,k,j)+field(i-1,k,j))
7106 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
7111 !--- end of 2nd order horizontal flux calculation
7115 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_pd, h_order not known ',horz_order
7116 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
7118 ENDIF horizontal_order_test
7120 ! pick up the rest of the horizontal radiation boundary conditions.
7121 ! (these are the computations that don't require 'cb'.
7122 ! first, set to index ranges
7125 i_end = MIN(ite,ide-1)
7127 j_end = MIN(jte,jde-1)
7129 ! compute x (u) conditions for v, w, or scalar
7131 IF( (config_flags%open_xs) .and. (its == ids) ) THEN
7133 DO j = j_start, j_end
7135 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
7136 tendency(its,k,j) = tendency(its,k,j) &
7138 ub*( field_old(its+1,k,j) &
7139 - field_old(its ,k,j) ) + &
7140 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) &
7147 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
7149 DO j = j_start, j_end
7151 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
7152 tendency(i_end,k,j) = tendency(i_end,k,j) &
7154 ub*( field_old(i_end ,k,j) &
7155 - field_old(i_end-1,k,j) ) + &
7156 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) &
7163 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
7165 DO i = i_start, i_end
7167 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
7168 tendency(i,k,jts) = tendency(i,k,jts) &
7170 vb*( field_old(i,k,jts+1) &
7171 - field_old(i,k,jts ) ) + &
7172 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) &
7179 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
7181 DO i = i_start, i_end
7183 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
7184 tendency(i,k,j_end) = tendency(i,k,j_end) &
7186 vb*( field_old(i,k,j_end ) &
7187 - field_old(i,k,j_end-1) ) + &
7188 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) &
7195 IF( (config_flags%polar) .and. (jts == jds) ) THEN
7197 ! Assuming rv(i,k,jds) = 0.
7198 DO i = i_start, i_end
7200 vb = MIN( 0.5*rv(i,k,jts+1), 0. )
7201 tendency(i,k,jts) = tendency(i,k,jts) &
7203 vb*( field_old(i,k,jts+1) &
7204 - field_old(i,k,jts ) ) + &
7205 field(i,k,jts)*rv(i,k,jts+1) &
7212 IF( (config_flags%polar) .and. (jte == jde)) THEN
7214 ! Assuming rv(i,k,jde) = 0.
7215 DO i = i_start, i_end
7217 vb = MAX( 0.5*rv(i,k,jte-1), 0. )
7218 tendency(i,k,j_end) = tendency(i,k,j_end) &
7220 vb*( field_old(i,k,j_end ) &
7221 - field_old(i,k,j_end-1) ) + &
7222 field(i,k,j_end)*(-rv(i,k,jte-1)) &
7229 !-------------------- vertical advection
7231 !-- loop bounds for periodic or sym conditions
7234 i_end = MIN(ite,ide-1)+1
7236 j_end = MIN(jte,jde-1)+1
7238 !-- loop bounds for open or specified conditions
7240 IF(degrade_xs) i_start = MAX(its-1,ids)
7241 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
7242 IF(degrade_ys) j_start = MAX(jts-1,jds)
7243 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
7245 vert_order_test : IF (vert_order == 6) THEN
7247 DO j = j_start, j_end
7249 DO i = i_start, i_end
7257 DO i = i_start, i_end
7258 dz = 2./(rdzw(k)+rdzw(k-1))
7259 mu = 0.5*(mut(i,j)+mut(i,j))
7262 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7264 fqz(i,k,j) = vel*flux6( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
7265 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
7266 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7270 DO i = i_start, i_end
7273 dz = 2./(rdzw(k)+rdzw(k-1))
7274 mu = 0.5*(mut(i,j)+mut(i,j))
7277 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7278 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7279 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7282 dz = 2./(rdzw(k)+rdzw(k-1))
7283 mu = 0.5*(mut(i,j)+mut(i,j))
7286 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7288 fqz(i,k,j) = vel*flux4( &
7289 field(i,k-2,j), field(i,k-1,j), &
7290 field(i,k ,j), field(i,k+1,j), -vel )
7291 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7294 dz = 2./(rdzw(k)+rdzw(k-1))
7295 mu = 0.5*(mut(i,j)+mut(i,j))
7298 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7300 fqz(i,k,j) = vel*flux4( &
7301 field(i,k-2,j), field(i,k-1,j), &
7302 field(i,k ,j), field(i,k+1,j), -vel )
7303 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7306 dz = 2./(rdzw(k)+rdzw(k-1))
7307 mu = 0.5*(mut(i,j)+mut(i,j))
7310 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7311 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7312 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7318 ELSE IF (vert_order == 5) THEN
7320 DO j = j_start, j_end
7322 DO i = i_start, i_end
7330 DO i = i_start, i_end
7331 dz = 2./(rdzw(k)+rdzw(k-1))
7332 mu = 0.5*(mut(i,j)+mut(i,j))
7335 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7337 fqz(i,k,j) = vel*flux5( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), &
7338 field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel )
7339 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7343 DO i = i_start, i_end
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)
7351 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7352 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7355 dz = 2./(rdzw(k)+rdzw(k-1))
7356 mu = 0.5*(mut(i,j)+mut(i,j))
7359 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7361 fqz(i,k,j) = vel*flux3( &
7362 field(i,k-2,j), field(i,k-1,j), &
7363 field(i,k ,j), field(i,k+1,j), -vel )
7364 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7367 dz = 2./(rdzw(k)+rdzw(k-1))
7368 mu = 0.5*(mut(i,j)+mut(i,j))
7371 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7373 fqz(i,k,j) = vel*flux3( &
7374 field(i,k-2,j), field(i,k-1,j), &
7375 field(i,k ,j), field(i,k+1,j), -vel )
7376 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7379 dz = 2./(rdzw(k)+rdzw(k-1))
7380 mu = 0.5*(mut(i,j)+mut(i,j))
7383 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7384 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7385 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7391 ELSE IF (vert_order == 4) THEN
7393 DO j = j_start, j_end
7395 DO i = i_start, i_end
7403 DO i = i_start, i_end
7405 dz = 2./(rdzw(k)+rdzw(k-1))
7406 mu = 0.5*(mut(i,j)+mut(i,j))
7409 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7411 fqz(i,k,j) = vel*flux4( &
7412 field(i,k-2,j), field(i,k-1,j), &
7413 field(i,k ,j), field(i,k+1,j), -vel )
7414 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7418 DO i = i_start, i_end
7421 dz = 2./(rdzw(k)+rdzw(k-1))
7422 mu = 0.5*(mut(i,j)+mut(i,j))
7425 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7426 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7427 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7430 dz = 2./(rdzw(k)+rdzw(k-1))
7431 mu = 0.5*(mut(i,j)+mut(i,j))
7434 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7435 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7436 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7442 ELSE IF (vert_order == 3) THEN
7444 DO j = j_start, j_end
7446 DO i = i_start, i_end
7454 DO i = i_start, i_end
7456 dz = 2./(rdzw(k)+rdzw(k-1))
7457 mu = 0.5*(mut(i,j)+mut(i,j))
7460 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7462 fqz(i,k,j) = vel*flux3( &
7463 field(i,k-2,j), field(i,k-1,j), &
7464 field(i,k ,j), field(i,k+1,j), -vel )
7465 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7469 DO i = i_start, i_end
7472 dz = 2./(rdzw(k)+rdzw(k-1))
7473 mu = 0.5*(mut(i,j)+mut(i,j))
7476 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7477 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7478 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7481 dz = 2./(rdzw(k)+rdzw(k-1))
7482 mu = 0.5*(mut(i,j)+mut(i,j))
7485 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7486 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7487 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7493 ELSE IF (vert_order == 2) THEN
7495 DO j = j_start, j_end
7497 DO i = i_start, i_end
7505 DO i = i_start, i_end
7507 dz = 2./(rdzw(k)+rdzw(k-1))
7508 mu = 0.5*(mut(i,j)+mut(i,j))
7511 fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
7512 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
7513 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
7522 WRITE (wrf_err_message,*) ' advect_scalar_pd, v_order not known ',vert_order
7523 CALL wrf_error_fatal ( wrf_err_message )
7525 ENDIF vert_order_test
7529 ! positive definite filter
7532 i_end = MIN(ite,ide-1)+1
7534 j_end = MIN(jte,jde-1)+1
7536 !-- loop bounds for open or specified conditions
7538 IF(degrade_xs) i_start = MAX(its-1,ids)
7539 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
7540 IF(degrade_ys) j_start = MAX(jts-1,jds)
7541 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
7543 IF(config_flags%specified .or. config_flags%nested) THEN
7544 IF (degrade_xs) i_start = MAX(its-1,ids+1)
7545 IF (degrade_xe) i_end = MIN(ite+1,ide-2)
7546 IF (degrade_ys) j_start = MAX(jts-1,jds+1)
7547 IF (degrade_ye) j_end = MIN(jte+1,jde-2)
7550 IF(config_flags%open_xs) THEN
7551 IF (degrade_xs) i_start = MAX(its-1,ids+1)
7553 IF(config_flags%open_xe) THEN
7554 IF (degrade_xe) i_end = MIN(ite+1,ide-2)
7556 IF(config_flags%open_ys) THEN
7557 IF (degrade_ys) j_start = MAX(jts-1,jds+1)
7559 IF(config_flags%open_ye) THEN
7560 IF (degrade_ye) j_end = MIN(jte+1,jde-2)
7563 ! We don't want to change j_start and j_end
7564 ! for polar BC's since we want to calculate
7565 ! fluxes for directions other than y at the
7568 !-- here is the limiter...
7574 ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) &
7575 - dt*( msftx(i,j)*msfty(i,j)*( &
7576 rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) + &
7577 rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) &
7578 +msfty(i,j)*rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) )
7580 flux_out = dt*( (msftx(i,j)*msfty(i,j))*( &
7581 rdx*( max(0.,fqx (i+1,k,j)) &
7582 -min(0.,fqx (i ,k,j)) ) &
7583 +rdy*( max(0.,fqy (i,k,j+1)) &
7584 -min(0.,fqy (i,k,j )) ) ) &
7585 +msfty(i,j)*rdzw(k)*( min(0.,fqz (i,k+1,j)) &
7586 -max(0.,fqz (i,k ,j)) ) )
7588 IF( flux_out .gt. ph_low ) THEN
7590 scale = max(0.,ph_low/(flux_out+eps))
7591 IF( fqx (i+1,k,j) .gt. 0.) fqx(i+1,k,j) = scale*fqx(i+1,k,j)
7592 IF( fqx (i ,k,j) .lt. 0.) fqx(i ,k,j) = scale*fqx(i ,k,j)
7593 IF( fqy (i,k,j+1) .gt. 0.) fqy(i,k,j+1) = scale*fqy(i,k,j+1)
7594 IF( fqy (i,k,j ) .lt. 0.) fqy(i,k,j ) = scale*fqy(i,k,j )
7595 ! note: z flux is opposite sign in mass coordinate because
7596 ! vertical coordinate decreases with increasing k
7597 IF( fqz (i,k+1,j) .lt. 0.) fqz(i,k+1,j) = scale*fqz(i,k+1,j)
7598 IF( fqz (i,k ,j) .gt. 0.) fqz(i,k ,j) = scale*fqz(i,k ,j)
7608 ! add in the pd-limited flux divergence
7611 i_end = MIN(ite,ide-1)
7613 j_end = MIN(jte,jde-1)
7615 DO j = j_start, j_end
7617 DO i = i_start, i_end
7619 tendency (i,k,j) = tendency(i,k,j) &
7620 -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) &
7621 +fqzl(i,k+1,j)-fqzl(i,k,j))
7628 DO j = j_start, j_end
7630 DO i = i_start, i_end
7632 z_tendency (i,k,j) = 0. -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) &
7633 +fqzl(i,k+1,j)-fqzl(i,k,j))
7642 IF(degrade_xs) i_start = MAX(its,ids+1)
7643 IF(degrade_xe) i_end = MIN(ite,ide-2)
7645 DO j = j_start, j_end
7647 DO i = i_start, i_end
7649 ! Un-"canceled" map scale factor, ADT Eq. 48
7650 tendency (i,k,j) = tendency(i,k,j) &
7651 - msftx(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) &
7652 +fqxl(i+1,k,j)-fqxl(i,k,j)) )
7659 DO j = j_start, j_end
7661 DO i = i_start, i_end
7663 h_tendency (i,k,j) = 0. &
7664 - msftx(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) &
7665 +fqxl(i+1,k,j)-fqxl(i,k,j)) )
7675 i_end = MIN(ite,ide-1)
7676 IF(degrade_ys) j_start = MAX(jts,jds+1)
7677 IF(degrade_ye) j_end = MIN(jte,jde-2)
7679 DO j = j_start, j_end
7681 DO i = i_start, i_end
7683 ! Un-"canceled" map scale factor, ADT Eq. 48
7684 ! It is correct to use msftx (and not msfty), per W. Skamarock, 20080606
7685 tendency (i,k,j) = tendency(i,k,j) &
7686 - msftx(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) &
7687 +fqyl(i,k,j+1)-fqyl(i,k,j)) )
7694 DO j = j_start, j_end
7696 DO i = i_start, i_end
7698 h_tendency (i,k,j) = h_tendency (i,k,j) &
7699 - msftx(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) &
7700 +fqyl(i,k,j+1)-fqyl(i,k,j)) )
7707 END SUBROUTINE advect_scalar_pd
7709 !----------------------------------------------------------------
7711 SUBROUTINE advect_scalar_mono ( field, field_old, tendency, &
7712 h_tendency, z_tendency, &
7717 msfux, msfuy, msfvx, msfvy, &
7720 rdx, rdy, rdzw, dt, &
7721 ids, ide, jds, jde, kds, kde, &
7722 ims, ime, jms, jme, kms, kme, &
7723 its, ite, jts, jte, kts, kte )
7725 ! monotonic advection option
7726 ! for scalars in WRF RK3 advection. This version is memory intensive ->
7727 ! we save 3d arrays of x, y and z both high and low order fluxes
7728 ! (six in all). Alternatively, we could sweep in a direction
7729 ! and lower the cost considerably.
7731 ! uses the Smolarkiewicz MWR 1989 approach, with addition of first-order
7738 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
7740 LOGICAL , INTENT(IN ) :: tenddec ! tendency flag
7742 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
7743 ims, ime, jms, jme, kms, kme, &
7744 its, ite, jts, jte, kts, kte
7746 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
7752 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, mub, mu_old
7753 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
7754 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: h_tendency, z_tendency
7756 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
7763 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
7767 REAL , INTENT(IN ) :: rdx, &
7773 INTEGER :: i, j, k, itf, jtf, ktf
7774 INTEGER :: i_start, i_end, j_start, j_end
7775 INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
7776 INTEGER :: jmin, jmax, jp, jm, imin, imax
7778 REAL :: mrdx, mrdy, ub, vb, uw, vw, mu
7779 REAL , DIMENSION(its:ite, kts:kte) :: vflux
7782 ! storage for high and low order fluxes
7784 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: fqx, fqy, fqz
7785 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: fqxl, fqyl, fqzl
7786 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: qmin, qmax
7787 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) :: scale_in, scale_out
7790 INTEGER :: horz_order, vert_order
7792 LOGICAL :: degrade_xs, degrade_ys
7793 LOGICAL :: degrade_xe, degrade_ye
7795 INTEGER :: jp1, jp0, jtmp
7797 REAL :: flux_out, ph_low, flux_in, ph_hi, scale
7798 REAL, PARAMETER :: eps=1.e-20
7801 ! definition of flux operators, 3rd, 4rth, 5th or 6th order
7803 REAL :: flux3, flux4, flux5, flux6, flux_upwind
7804 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr
7806 flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
7807 (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2)
7809 flux3(q_im2, q_im1, q_i, q_ip1, ua) = &
7810 flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
7811 sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1))
7813 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
7814 (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2) &
7815 +(1./60.)*(q_ip2+q_im3)
7817 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = &
7818 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) &
7819 -sign(1.,ua)*(1./60.)*( &
7820 (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )
7822 ! flux_upwind(q_im1, q_i, cr ) = 0.
7823 flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 &
7824 +0.5*(1.-sign(1.,cr))*q_i
7826 LOGICAL, PARAMETER :: mono_limit = .true.
7828 ! set order for the advection schemes
7831 horz_order = config_flags%h_sca_adv_order
7832 vert_order = config_flags%v_sca_adv_order
7837 qmin(i,k,j) = field_old(i,k,j)
7838 qmax(i,k,j) = field_old(i,k,j)
7839 scale_in(i,k,j) = 1.
7840 scale_out(i,k,j) = 1.
7851 ! begin with horizontal flux divergence
7852 ! here is the choice of flux operators
7855 horizontal_order_test : IF( horz_order == 5 ) THEN
7857 ! determine boundary mods for flux operators
7858 ! We degrade the flux operators from 3rd/4rth order
7859 ! to second order one gridpoint in from the boundaries for
7860 ! all boundary conditions except periodic and symmetry - these
7861 ! conditions have boundary zone data fill for correct application
7862 ! of the higher order flux stencils
7869 IF( config_flags%periodic_x .or. &
7870 config_flags%symmetric_xs .or. &
7871 (its > ids+3) ) degrade_xs = .false.
7872 IF( config_flags%periodic_x .or. &
7873 config_flags%symmetric_xe .or. &
7874 (ite < ide-4) ) degrade_xe = .false.
7875 IF( config_flags%periodic_y .or. &
7876 config_flags%symmetric_ys .or. &
7877 (jts > jds+3) ) degrade_ys = .false.
7878 IF( config_flags%periodic_y .or. &
7879 config_flags%symmetric_ye .or. &
7880 (jte < jde-4) ) degrade_ye = .false.
7882 !--------------- y - advection first
7884 !-- y flux compute; these bounds are for periodic and sym b.c.
7888 i_end = MIN(ite,ide-1)+1
7890 j_end = MIN(jte,jde-1)+1
7894 !-- modify loop bounds if open or specified
7897 ! IF(degrade_xs) i_start = its
7898 ! IF(degrade_xe) i_end = MIN(ite,ide-1)
7899 IF(degrade_xs) i_start = MAX(its-1,ids)
7900 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
7903 ! IF(degrade_ys) then
7904 ! j_start = MAX(jts,jds+1)
7908 ! IF(degrade_ye) then
7909 ! j_end = MIN(jte,jde-2)
7914 j_start = MAX(jts-1,jds+1)
7919 j_end = MIN(jte+1,jde-2)
7923 ! compute fluxes, 5th order
7925 j_loop_y_flux_5 : DO j = j_start, j_end+1
7927 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
7930 DO i = i_start, i_end
7934 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), vel)
7936 fqy( i, k, j ) = vel*flux5( &
7937 field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), &
7938 field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel )
7940 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7943 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1))
7944 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1))
7946 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j))
7947 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j))
7953 ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary
7956 DO i = i_start, i_end
7960 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
7962 fqy(i,k, j) = 0.5*rv(i,k,j)* &
7963 (field(i,k,j)+field(i,k,j-1))
7965 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7968 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1))
7969 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1))
7971 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j))
7972 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j))
7978 ELSE IF ( j == jds+2 ) THEN ! third of 4rth order flux 2 in from south boundary
7981 DO i = i_start, i_end
7985 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
7987 fqy( i, k, j ) = vel*flux3( &
7988 field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
7989 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
7992 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1))
7993 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1))
7995 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j))
7996 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j))
8002 ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary
8005 DO i = i_start, i_end
8009 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
8011 fqy(i, k, j ) = 0.5*rv(i,k,j)* &
8012 (field(i,k,j)+field(i,k,j-1))
8013 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
8016 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1))
8017 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1))
8019 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j))
8020 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j))
8026 ELSE IF ( j == jde-2 ) THEN ! 3rd or 4rth order flux 2 in from north boundary
8029 DO i = i_start, i_end
8033 fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr)
8035 fqy( i, k, j) = vel*flux3( &
8036 field(i,k,j-2),field(i,k,j-1), &
8037 field(i,k,j),field(i,k,j+1),vel )
8038 fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
8041 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k,j-1))
8042 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k,j-1))
8044 qmax(i,k,j-1) = amax1(qmax(i,k,j-1),field_old(i,k,j))
8045 qmin(i,k,j-1) = amin1(qmin(i,k,j-1),field_old(i,k,j))
8053 ENDDO j_loop_y_flux_5
8057 !-- these bounds are for periodic and sym conditions
8060 i_end = MIN(ite,ide-1)+1
8065 j_end = MIN(jte,jde-1)+1
8067 !-- modify loop bounds for open and specified b.c
8070 ! IF(degrade_ys) j_start = jts
8071 ! IF(degrade_ye) j_end = MIN(jte,jde-1)
8072 IF(degrade_ys) j_start = MAX(jts-1,jds)
8073 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
8076 ! IF(degrade_xs) then
8077 ! i_start = MAX(ids+1,its)
8078 ! i_start_f = i_start+2
8081 ! IF(degrade_xe) then
8082 ! i_end = MIN(ide-2,ite)
8087 i_start = MAX(ids+1,its-1)
8092 i_end = MIN(ide-2,ite+1)
8098 DO j = j_start, j_end
8100 ! 5th or 6th order flux
8103 DO i = i_start_f, i_end_f
8107 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
8109 fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), &
8110 field(i-1,k,j), field(i ,k,j), &
8111 field(i+1,k,j), field(i+2,k,j), &
8113 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
8116 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j))
8117 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j))
8119 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j))
8120 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j))
8126 ! lower order fluxes close to boundaries (if not periodic or symmetric)
8128 ! WCS 20090218 degrade_xs and xe recoded
8130 IF( degrade_xs ) THEN
8132 DO i=i_start,i_start_f-1
8134 IF(i == ids+1) THEN ! second order
8138 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
8140 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
8141 *(field(i,k,j)+field(i-1,k,j))
8143 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
8146 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j))
8147 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j))
8149 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j))
8150 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j))
8155 IF(i == ids+2) THEN ! third order
8159 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
8160 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
8161 field(i ,k,j), field(i+1,k,j), &
8163 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
8166 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j))
8167 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j))
8169 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j))
8170 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j))
8179 IF( degrade_xe ) THEN
8181 DO i = i_end_f+1, i_end+1
8183 IF( i == ide-1 ) THEN ! second order flux next to the boundary
8187 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
8188 fqx(i,k,j) = 0.5*(ru(i,k,j)) &
8189 *(field(i,k,j)+field(i-1,k,j))
8190 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
8193 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j))
8194 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j))
8196 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j))
8197 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j))
8202 IF( i == ide-2 ) THEN ! third order flux one in from the boundary
8206 fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr)
8207 fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), &
8208 field(i ,k,j), field(i+1,k,j), &
8210 fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
8213 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i-1,k,j))
8214 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i-1,k,j))
8216 qmax(i-1,k,j) = amax1(qmax(i-1,k,j),field_old(i,k,j))
8217 qmin(i-1,k,j) = amin1(qmin(i-1,k,j),field_old(i,k,j))
8224 ENDDO ! enddo for outer J loop
8228 WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_mono, h_order not known ',horz_order
8229 CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
8231 ENDIF horizontal_order_test
8233 ! pick up the rest of the horizontal radiation boundary conditions.
8234 ! (these are the computations that don't require 'cb'.
8235 ! first, set to index ranges
8238 i_end = MIN(ite,ide-1)
8240 j_end = MIN(jte,jde-1)
8242 ! compute x (u) conditions for v, w, or scalar
8244 IF( (config_flags%open_xs) .and. (its == ids) ) THEN
8246 DO j = j_start, j_end
8248 ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
8249 tendency(its,k,j) = tendency(its,k,j) &
8251 ub*( field_old(its+1,k,j) &
8252 - field_old(its ,k,j) ) + &
8253 field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) &
8260 IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
8262 DO j = j_start, j_end
8264 ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
8265 tendency(i_end,k,j) = tendency(i_end,k,j) &
8267 ub*( field_old(i_end ,k,j) &
8268 - field_old(i_end-1,k,j) ) + &
8269 field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) &
8276 IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
8278 DO i = i_start, i_end
8280 vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
8281 tendency(i,k,jts) = tendency(i,k,jts) &
8283 vb*( field_old(i,k,jts+1) &
8284 - field_old(i,k,jts ) ) + &
8285 field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) &
8292 IF( (config_flags%open_ye) .and. (jte == jde)) THEN
8294 DO i = i_start, i_end
8296 vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
8297 tendency(i,k,j_end) = tendency(i,k,j_end) &
8299 vb*( field_old(i,k,j_end ) &
8300 - field_old(i,k,j_end-1) ) + &
8301 field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) &
8308 !-------------------- vertical advection
8310 !-- loop bounds for periodic or sym conditions
8313 i_end = MIN(ite,ide-1)+1
8315 j_end = MIN(jte,jde-1)+1
8317 !-- loop bounds for open or specified conditions
8320 ! IF(degrade_xs) i_start = its
8321 ! IF(degrade_xe) i_end = MIN(ite,ide-1)
8322 ! IF(degrade_ys) j_start = jts
8323 ! IF(degrade_ye) j_end = MIN(jte,jde-1)
8325 IF(degrade_xs) i_start = MAX(its-1,ids)
8326 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
8327 IF(degrade_ys) j_start = MAX(jts-1,jds)
8328 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
8331 vert_order_test : IF (vert_order == 3) THEN
8333 DO j = j_start, j_end
8335 DO i = i_start, i_end
8343 DO i = i_start, i_end
8347 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
8349 fqz(i,k,j) = vel*flux3( &
8350 field(i,k-2,j), field(i,k-1,j), &
8351 field(i,k ,j), field(i,k+1,j), -vel )
8352 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
8355 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j))
8356 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j))
8358 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j))
8359 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j))
8365 DO i = i_start, i_end
8370 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
8371 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
8372 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
8375 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j))
8376 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j))
8378 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j))
8379 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j))
8385 fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr)
8386 fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
8387 fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
8390 qmax(i,k,j) = amax1(qmax(i,k,j),field_old(i,k-1,j))
8391 qmin(i,k,j) = amin1(qmin(i,k,j),field_old(i,k-1,j))
8393 qmax(i,k-1,j) = amax1(qmax(i,k-1,j),field_old(i,k,j))
8394 qmin(i,k-1,j) = amin1(qmin(i,k-1,j),field_old(i,k,j))
8402 WRITE (wrf_err_message,*) ' advect_scalar_mono, v_order not known ',vert_order
8403 CALL wrf_error_fatal ( wrf_err_message )
8405 ENDIF vert_order_test
8407 IF (mono_limit) THEN
8412 i_end = MIN(ite,ide-1)+1
8414 j_end = MIN(jte,jde-1)+1
8418 !-- loop bounds for open or specified conditions
8420 ! IF(degrade_xs) i_start = its
8421 ! IF(degrade_xe) i_end = MIN(ite,ide-1)
8422 ! IF(degrade_ys) j_start = jts
8423 ! IF(degrade_ye) j_end = MIN(jte,jde-1)
8425 ! IF(config_flags%specified .or. config_flags%nested) THEN
8426 ! IF (degrade_xs) i_start = MAX(its,ids+1)
8427 ! IF (degrade_xe) i_end = MIN(ite,ide-2)
8428 ! IF (degrade_ys) j_start = MAX(jts,jds+1)
8429 ! IF (degrade_ye) j_end = MIN(jte,jde-2)
8432 ! IF(config_flags%open_xs) THEN
8433 ! IF (degrade_xs) i_start = MAX(its,ids+1)
8435 ! IF(config_flags%open_xe) THEN
8436 ! IF (degrade_xe) i_end = MIN(ite,ide-2)
8438 ! IF(config_flags%open_ys) THEN
8439 ! IF (degrade_ys) j_start = MAX(jts,jds+1)
8441 ! IF(config_flags%open_ye) THEN
8442 ! IF (degrade_ye) j_end = MIN(jte,jde-2)
8445 IF(degrade_xs) i_start = MAX(its-1,ids)
8446 IF(degrade_xe) i_end = MIN(ite+1,ide-1)
8447 IF(degrade_ys) j_start = MAX(jts-1,jds)
8448 IF(degrade_ye) j_end = MIN(jte+1,jde-1)
8450 IF(config_flags%specified .or. config_flags%nested) THEN
8451 IF (degrade_xs) i_start = MAX(its-1,ids+1)
8452 IF (degrade_xe) i_end = MIN(ite+1,ide-2)
8453 IF (degrade_ys) j_start = MAX(jts-1,jds+1)
8454 IF (degrade_ye) j_end = MIN(jte+1,jde-2)
8457 IF(config_flags%open_xs) THEN
8458 IF (degrade_xs) i_start = MAX(its-1,ids+1)
8460 IF(config_flags%open_xe) THEN
8461 IF (degrade_xe) i_end = MIN(ite+1,ide-2)
8463 IF(config_flags%open_ys) THEN
8464 IF (degrade_ys) j_start = MAX(jts-1,jds+1)
8466 IF(config_flags%open_ye) THEN
8467 IF (degrade_ye) j_end = MIN(jte+1,jde-2)
8470 !-- here is the limiter...
8476 ph_upwind = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) &
8477 - dt*( msftx(i,j)*msfty(i,j)*( &
8478 rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) + &
8479 rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) &
8480 +msfty(i,j)*rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) )
8482 flux_in = -dt*( (msftx(i,j)*msfty(i,j))*( &
8483 rdx*( min(0.,fqx (i+1,k,j)) &
8484 -max(0.,fqx (i ,k,j)) ) &
8485 +rdy*( min(0.,fqy (i,k,j+1)) &
8486 -max(0.,fqy (i,k,j )) ) ) &
8487 +msfty(i,j)*rdzw(k)*( max(0.,fqz (i,k+1,j)) &
8488 -min(0.,fqz (i,k ,j)) ) )
8490 ph_hi = mut(i,j)*qmax(i,k,j) - ph_upwind
8491 IF( flux_in .gt. ph_hi ) scale_in(i,k,j) = max(0.,ph_hi/(flux_in+eps))
8494 flux_out = dt*( (msftx(i,j)*msfty(i,j))*( &
8495 rdx*( max(0.,fqx (i+1,k,j)) &
8496 -min(0.,fqx (i ,k,j)) ) &
8497 +rdy*( max(0.,fqy (i,k,j+1)) &
8498 -min(0.,fqy (i,k,j )) ) ) &
8499 +msfty(i,j)*rdzw(k)*( min(0.,fqz (i,k+1,j)) &
8500 -max(0.,fqz (i,k ,j)) ) )
8502 ph_low = ph_upwind - mut(i,j)*qmin(i,k,j)
8503 IF( flux_out .gt. ph_low ) scale_out(i,k,j) = max(0.,ph_low/(flux_out+eps))
8511 DO i=i_start, i_end+1
8512 IF( fqx (i,k,j) .gt. 0.) then
8513 fqx(i,k,j) = min(scale_in(i,k,j),scale_out(i-1,k,j))*fqx(i,k,j)
8515 fqx(i,k,j) = min(scale_out(i,k,j),scale_in(i-1,k,j))*fqx(i,k,j)
8521 DO j=j_start, j_end+1
8524 IF( fqy (i,k,j) .gt. 0.) then
8525 fqy(i,k,j) = min(scale_in(i,k,j),scale_out(i,k,j-1))*fqy(i,k,j)
8527 fqy(i,k,j) = min(scale_out(i,k,j),scale_in(i,k,j-1))*fqy(i,k,j)
8536 IF( fqz (i,k,j) .lt. 0.) then
8537 fqz(i,k,j) = min(scale_in(i,k,j),scale_out(i,k-1,j))*fqz(i,k,j)
8539 fqz(i,k,j) = min(scale_out(i,k,j),scale_in(i,k-1,j))*fqz(i,k,j)
8547 ! add in the mono-limited flux divergence
8548 ! we need to fix this for open b.c set ***********
8551 i_end = MIN(ite,ide-1)
8553 j_end = MIN(jte,jde-1)
8555 DO j = j_start, j_end
8557 DO i = i_start, i_end
8559 tendency (i,k,j) = tendency(i,k,j) &
8560 -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) &
8561 +fqzl(i,k+1,j)-fqzl(i,k,j))
8568 DO j = j_start, j_end
8570 DO i = i_start, i_end
8572 z_tendency (i,k,j) = 0. -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) &
8573 +fqzl(i,k+1,j)-fqzl(i,k,j))
8584 ! IF(degrade_xs) i_start = i_start + 1
8585 ! IF(degrade_xe) i_end = i_end - 1
8587 IF(degrade_xs) i_start = MAX(its,ids+1)
8588 IF(degrade_xe) i_end = MIN(ite,ide-2)
8590 DO j = j_start, j_end
8592 DO i = i_start, i_end
8594 ! Un-"canceled" map scale factor, ADT Eq. 48
8595 tendency (i,k,j) = tendency(i,k,j) &
8596 - msftx(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) &
8597 +fqxl(i+1,k,j)-fqxl(i,k,j)) )
8604 DO j = j_start, j_end
8606 DO i = i_start, i_end
8608 h_tendency (i,k,j) = 0. &
8609 - msftx(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) &
8610 +fqxl(i+1,k,j)-fqxl(i,k,j)) )
8620 i_end = MIN(ite,ide-1)
8623 ! IF(degrade_ys) j_start = j_start + 1
8624 ! IF(degrade_ye) j_end = j_end - 1
8626 IF(degrade_ys) j_start = MAX(jts,jds+1)
8627 IF(degrade_ye) j_end = MIN(jte,jde-2)
8629 DO j = j_start, j_end
8631 DO i = i_start, i_end
8633 ! Un-"canceled" map scale factor, ADT Eq. 48
8634 tendency (i,k,j) = tendency(i,k,j) &
8635 - msftx(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) &
8636 +fqyl(i,k,j+1)-fqyl(i,k,j)) )
8643 DO j = j_start, j_end
8645 DO i = i_start, i_end
8647 h_tendency (i,k,j) = h_tendency (i,k,j) &
8648 - msftx(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) &
8649 +fqyl(i,k,j+1)-fqyl(i,k,j)) )
8656 END SUBROUTINE advect_scalar_mono
8658 !-----------------------------------------------------------
8660 END MODULE module_advect_em