Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_advect_em.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_advect_em.F	(revision 195)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_advect_em.F	(revision 196)
@@ -124,6 +124,7 @@
 SUBROUTINE advect_u   ( u, u_old, tendency,            &
                         ru, rv, rom,                   &
-                        mut, config_flags,             &
-                        msfu, msfv, msft,              &
+                        mut, time_step, config_flags,  &
+                        msfu ,        msfv ,           &
+                        msft ,                         &
                         fzm, fzp,                      &
                         rdx, rdy, rdzw,                &
@@ -161,4 +162,5 @@
    REAL ,                                        INTENT(IN   ) :: rdx,  &
                                                                   rdy
+   INTEGER ,                                     INTENT(IN   ) :: time_step
 
    ! Local data
@@ -192,5 +194,5 @@
    flux3(q_im2, q_im1, q_i, q_ip1, ua) =                         &
             flux4(q_im2, q_im1, q_i, q_ip1, ua) +                &
-            sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+            sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
 
    flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =           &
@@ -200,5 +202,5 @@
    flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =           &
            flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)     &
-            -sign(1.,ua)*(                              &
+            -sign(1,time_step)*sign(1.,ua)*(                     &
               (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
 
@@ -227,21 +229,21 @@
 !   of the higher order flux stencils
 
-      degrade_xs = .true.
-      degrade_xe = .true.
-      degrade_ys = .true.
-      degrade_ye = .true.
-
-      IF( config_flags%periodic_x   .or. &
-          config_flags%symmetric_xs .or. &
-          (its > ids+2)                ) degrade_xs = .false.
-      IF( config_flags%periodic_x   .or. &
-          config_flags%symmetric_xe .or. &
-          (ite < ide-2)                ) degrade_xe = .false.
-      IF( config_flags%periodic_y   .or. &
-          config_flags%symmetric_ys .or. &
-          (jts > jds+2)                ) degrade_ys = .false.
-      IF( config_flags%periodic_y   .or. &
-          config_flags%symmetric_ye .or. &
-          (jte < jde-3)                ) degrade_ye = .false.
+   degrade_xs = .true.
+   degrade_xe = .true.
+   degrade_ys = .true.
+   degrade_ye = .true.
+
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xs .or. &
+       (its > ids+3)                ) degrade_xs = .false.
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xe .or. &
+       (ite < ide-2)                ) degrade_xe = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ys .or. &
+       (jts > jds+3)                ) degrade_ys = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ye .or. &
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -273,4 +275,5 @@
       ENDIF
 
+
 !  compute fluxes, 5th or 6th order
 
@@ -280,59 +283,57 @@
      j_loop_y_flux_6 : DO j = j_start, j_end+1
 
-        IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN  ! use full stencil
-
-           DO k=kts,ktf
-           DO i = i_start, i_end
-              vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
-              fqy( i, k, jp1 ) = vel*flux6(                                &
-                                 u(i,k,j-3), u(i,k,j-2), u(i,k,j-1),       &
-                                 u(i,k,j  ), u(i,k,j+1), u(i,k,j+2),  vel )
-           ENDDO
-           ENDDO
+      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN  ! use full stencil
+
+        DO k=kts,ktf
+        DO i = i_start, i_end
+          vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
+          fqy( i, k, jp1 ) = vel*flux6(               &
+                  u(i,k,j-3), u(i,k,j-2), u(i,k,j-1),       &
+                  u(i,k,j  ), u(i,k,j+1), u(i,k,j+2),  vel )
+        ENDDO
+        ENDDO
 
 !  we must be close to some boundary where we need to reduce the order of the stencil
 
-        ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
-
-           DO k=kts,ktf
-           DO i = i_start, i_end
+      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
               fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j))  &
-                                    *(u(i,k,j)+u(i,k,j-1))
-           ENDDO
-           ENDDO
-
-        ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
-
-           DO k=kts,ktf
-           DO i = i_start, i_end
+                                     *(u(i,k,j)+u(i,k,j-1))
+            ENDDO
+            ENDDO
+
+     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
               vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
               fqy( i, k, jp1 ) = vel*flux4(      &
                    u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel )
-           ENDDO
-           ENDDO
-
-        ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
-
-           DO k=kts,ktf
-           DO i = i_start, i_end
+            ENDDO
+            ENDDO
+
+     ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
               fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j))    &
                      *(u(i,k,j)+u(i,k,j-1))
-           ENDDO
-           ENDDO
-
-        ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
-
-           DO k=kts,ktf
-           DO i = i_start, i_end
+            ENDDO
+            ENDDO
+
+     ELSE IF ( j == jde-2 ) THEN  ! 3rd order flux 2 in from north boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
               vel = 0.5*(rv(i,k,j)+rv(i-1,k,j))
               fqy( i, k, jp1 ) = vel*flux4(     &
                    u(i,k,j-2),u(i,k,j-1),    &
                    u(i,k,j),u(i,k,j+1),vel )
-           ENDDO
-           ENDDO
-
-        END IF
-
-!stopped
+            ENDDO
+            ENDDO
+
+      END IF
 
 !  y flux-divergence into tendency
@@ -348,4 +349,5 @@
 
         ENDIF
+
 
 
@@ -475,5 +477,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -481,8 +483,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -514,4 +516,5 @@
       ENDIF
 
+
 !  compute fluxes, 5th or 6th order
 
@@ -562,5 +565,5 @@
             ENDDO
 
-     ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
+     ELSE IF ( j == jde-2 ) THEN  ! 3rd order flux 2 in from north boundary
 
             DO k=kts,ktf
@@ -587,4 +590,5 @@
 
         ENDIF
+
 
 
@@ -709,5 +713,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -715,8 +719,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !--------------- x - advection first
@@ -817,4 +821,5 @@
         j_end_f = jde-2
       ENDIF
+
 
 !  j flux loop for v flux of u momentum
@@ -835,6 +840,9 @@
        DO k = kts, ktf
        DO i = i_start, i_end
-         fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1))    &
-                *(u(i,k,j_end+1)+u(i,k,j_end))
+         ! Assumes j>j_end_f is ONLY j_end+1 ...
+!         fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1))    &
+!                *(u(i,k,j_end+1)+u(i,k,j_end))
+         fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j))    &
+                *(u(i,k,j)+u(i,k,j-1))
        ENDDO
        ENDDO
@@ -852,7 +860,7 @@
      END IF
 
+!  y flux-divergence into tendency
+
      IF (j > j_start) THEN
-
-!  y flux-divergence into tendency
 
        DO k=kts,ktf
@@ -865,4 +873,5 @@
      END IF
 
+
      jtmp = jp1
      jp1 = jp0
@@ -890,5 +899,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -896,8 +905,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !--------------- x - advection first
@@ -997,4 +1006,5 @@
         j_end_f = jde-2
       ENDIF
+
 
 !  j flux loop for v flux of u momentum
@@ -1015,6 +1025,9 @@
        DO k = kts, ktf
        DO i = i_start, i_end
-         fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1))    &
-                *(u(i,k,j_end+1)+u(i,k,j_end))
+         ! Assumes j>j_end_f is ONLY j_end+1 ...
+!         fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1))    &
+!                *(u(i,k,j_end+1)+u(i,k,j_end))
+         fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j))    &
+                *(u(i,k,j)+u(i,k,j-1))
        ENDDO
        ENDDO
@@ -1032,7 +1045,7 @@
      END IF
 
+!  y flux-divergence into tendency
+
      IF (j > j_start) THEN
-
-!  y flux-divergence into tendency
 
        DO k=kts,ktf
@@ -1044,4 +1057,5 @@
 
      END IF
+
 
      jtmp = jp1
@@ -1116,4 +1130,8 @@
       ENDDO
       ENDDO
+
+   ELSE IF ( horz_order == 0 ) THEN
+
+      ! Just in case we want to turn horizontal advection off, we can do it
 
    ELSE
@@ -1410,5 +1428,5 @@
 SUBROUTINE advect_v   ( v, v_old, tendency,            &
                         ru, rv, rom,                   &
-                        mut, config_flags,             &
+                        mut, time_step, config_flags,  &
                         msfu, msfv, msft,              &
                         fzm, fzp,                      &
@@ -1447,4 +1465,6 @@
    REAL ,                                        INTENT(IN   ) :: rdx,  &
                                                                   rdy
+   INTEGER ,                                     INTENT(IN   ) :: time_step
+
 
    ! Local data
@@ -1481,5 +1501,5 @@
    flux3(q_im2, q_im1, q_i, q_ip1, ua) =                     &
            flux4(q_im2, q_im1, q_i, q_ip1, ua) +                &
-           sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+           sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
 
    flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
@@ -1489,5 +1509,5 @@
    flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
            flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)    &
-            -sign(1.,ua)*(                             &
+            -sign(1,time_step)*sign(1.,ua)*(                    &
               (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
 
@@ -1518,25 +1538,23 @@
 !   of the higher order flux stencils
 
-      degrade_xs = .true.
-      degrade_xe = .true.
-      degrade_ys = .true.
-      degrade_ye = .true.
-
-      IF( config_flags%periodic_x   .or. &
-          config_flags%symmetric_xs .or. &
-          (its > ids+2)                ) degrade_xs = .false.
-      IF( config_flags%periodic_x   .or. &
-          config_flags%symmetric_xe .or. &
-          (ite < ide-3)                ) degrade_xe = .false.
-      IF( config_flags%periodic_y   .or. &
-          config_flags%symmetric_ys .or. &
-          (jts > jds+2)                ) degrade_ys = .false.
-      IF( config_flags%periodic_y   .or. &
-          config_flags%symmetric_ye .or. &
-          (jte < jde-2)                ) degrade_ye = .false.
+   degrade_xs = .true.
+   degrade_xe = .true.
+   degrade_ys = .true.
+   degrade_ye = .true.
+
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xs .or. &
+       (its > ids+3)                ) degrade_xs = .false.
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xe .or. &
+       (ite < ide-3)                ) degrade_xe = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ys .or. &
+       (jts > jds+3)                ) degrade_ys = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ye .or. &
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !--------------- y - advection first
-
-      ktf=MIN(kte,kde-1)
 
       i_start = its
@@ -1552,35 +1570,35 @@
 
       IF(degrade_ys) then
-         j_start = MAX(jts,jds+1)
-         j_start_f = jds+3
+        j_start = MAX(jts,jds+1)
+        j_start_f = jds+3
       ENDIF
 
       IF(degrade_ye) then
-         j_end = MIN(jte,jde-1)
-         j_end_f = jde-2
+        j_end = MIN(jte,jde-1)
+        j_end_f = jde-2
       ENDIF
 
 !  compute fluxes, 5th or 6th order
 
-      jp1 = 2
-      jp0 = 1
-
-      j_loop_y_flux_6 : DO j = j_start, j_end+1
-
-         IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
-
-            DO k=kts,ktf
-            DO i = i_start, i_end
-               vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
-               fqy( i, k, jp1 ) = vel*flux6(                                &
-                                  v(i,k,j-3), v(i,k,j-2), v(i,k,j-1),       &
-                                  v(i,k,j  ), v(i,k,j+1), v(i,k,j+2),  vel )
-            ENDDO
-            ENDDO
+     jp1 = 2
+     jp0 = 1
+
+     j_loop_y_flux_6 : DO j = j_start, j_end+1
+
+      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
+
+        DO k=kts,ktf
+        DO i = i_start, i_end
+          vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
+          fqy( i, k, jp1 ) = vel*flux6(               &
+                  v(i,k,j-3), v(i,k,j-2), v(i,k,j-1),       &
+                  v(i,k,j  ), v(i,k,j+1), v(i,k,j+2),  vel )
+        ENDDO
+        ENDDO
 
 !  we must be close to some boundary where we need to reduce the order of the stencil
 !  specified uses upstream normal wind at boundaries
 
-         ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
+      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
 
             DO k=kts,ktf
@@ -1593,16 +1611,16 @@
             ENDDO
 
-         ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
+     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
 
             DO k=kts,ktf
             DO i = i_start, i_end
-                vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
-                fqy( i, k, jp1 ) = vel*flux4(      &
-                                   v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-
-         ELSE IF ( j == jde ) THEN  ! 2nd order flux next to north boundary
+              vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
+              fqy( i, k, jp1 ) = vel*flux4(      &
+                   v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+
+     ELSE IF ( j == jde ) THEN  ! 2nd order flux next to north boundary
 
             DO k=kts,ktf
@@ -1615,34 +1633,37 @@
             ENDDO
 
-         ELSE IF ( j == jde-1 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
+     ELSE IF ( j == jde-1 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
 
             DO k=kts,ktf
             DO i = i_start, i_end
-                vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
-                fqy( i, k, jp1 ) = vel*flux4(     &
-                                   v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-         END IF
+              vel = 0.5*(rv(i,k,j)+rv(i,k,j-1))
+              fqy( i, k, jp1 ) = vel*flux4(     &
+                   v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+      END IF
 
 !  y flux-divergence into tendency
 
-         IF(j > j_start) THEN
-
-            DO k=kts,ktf
-            DO i = i_start, i_end
-               mrdy=msfv(i,j-1)*rdy
-               tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
-            ENDDO
-            ENDDO
-
-         ENDIF
-
-         jtmp = jp1
-         jp1 = jp0
-         jp0 = jtmp
-
-      ENDDO j_loop_y_flux_6
+
+        IF(j > j_start) THEN
+
+          DO k=kts,ktf
+          DO i = i_start, i_end
+            mrdy=msfv(i,j-1)*rdy
+            tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
+          ENDDO
+          ENDDO
+
+        ENDIF
+
+
+
+        jtmp = jp1
+        jp1 = jp0
+        jp0 = jtmp
+
+   ENDDO j_loop_y_flux_6
 
 !  next, x - flux divergence
@@ -1653,4 +1674,5 @@
       j_start = jts
       j_end   = jte
+
       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
@@ -1663,11 +1685,12 @@
 
       IF(degrade_xs) then
-         i_start = MAX(ids+1,its)
-         i_start_f = i_start+2
+        i_start = MAX(ids+1,its)
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
       ENDIF
 
       IF(degrade_xe) then
-         i_end = MIN(ide-2,ite)
-         i_end_f = ide-3
+        i_end = MIN(ide-2,ite)
+        i_end_f = ide-3
       ENDIF
 
@@ -1678,66 +1701,74 @@
 !  5th or 6th order flux
 
-         DO k=kts,ktf
-         DO i = i_start_f, i_end_f
-            vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
-            fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j),  &
-                                     v(i-1,k,j), v(i  ,k,j),  &
-                                     v(i+1,k,j), v(i+2,k,j),  &
-                                     vel                     )
-         ENDDO
-         ENDDO
+        DO k=kts,ktf
+        DO i = i_start_f, i_end_f
+          vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
+          fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j),  &
+                                         v(i-1,k,j), v(i  ,k,j),  &
+                                         v(i+1,k,j), v(i+2,k,j),  &
+                                         vel                     )
+        ENDDO
+        ENDDO
 
 !  lower order fluxes close to boundaries (if not periodic or symmetric)
 
-         IF( degrade_xs ) THEN
-
-            IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-               i = ids+1
-               DO k=kts,ktf
-                  fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
-                                 *(v(i,k,j)+v(i-1,k,j))
-               ENDDO
+        IF( degrade_xs ) THEN
+
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
+                                *(v(i,k,j)+v(i-1,k,j))
+              ENDDO
             ENDIF
 
-            i = ids+2
-            DO k=kts,ktf
-               vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
-               fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j),  &
-                                       v(i  ,k,j), v(i+1,k,j),  &
-                                       vel                     )
-            ENDDO
-
-         ENDIF
-
-         IF( degrade_xe ) THEN
-
-            IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-               i = ide-1
-               DO k=kts,ktf
-                  fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1))      &
-                                 *(v(i_end+1,k,j)+v(i_end,k,j))
-               ENDDO
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
+                fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j),  &
+                                        v(i  ,k,j), v(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
             ENDIF
 
-            i = ide-2
-            DO k=kts,ktf
-               vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
-               fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j),  &
-                                       v(i  ,k,j), v(i+1,k,j),  &
-                                       vel                     )
-            ENDDO
-
-         ENDIF
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1))      &
+                                *(v(i_end+1,k,j)+v(i_end,k,j))
+              ENDDO
+            ENDIF
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts,ktf
+                vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
+                fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j),  &
+                                        v(i  ,k,j), v(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+            ENDIF
+
+          ENDDO
+
+        ENDIF
 
 !  x flux-divergence into tendency
 
-         DO k=kts,ktf
-            DO i = i_start, i_end
+        DO k=kts,ktf
+          DO i = i_start, i_end
             mrdx=msfv(i,j)*rdx
             tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
-         ENDDO
-      ENDDO
-
-   ENDDO
+          ENDDO
+        ENDDO
+
+      ENDDO
 
    ELSE IF( horz_order == 5 ) THEN
@@ -1762,5 +1793,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -1768,8 +1799,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -1875,4 +1906,5 @@
         ENDIF
 
+
         jtmp = jp1
         jp1 = jp0
@@ -1899,5 +1931,6 @@
       IF(degrade_xs) then
         i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
       ENDIF
 
@@ -1927,38 +1960,46 @@
         IF( degrade_xs ) THEN
 
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts,ktf
-            fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
-                   *(v(i,k,j)+v(i-1,k,j))
-            ENDDO
-         ENDIF
-
-          i = ids+2
-          DO k=kts,ktf
-            vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
-            fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j),  &
-                                          v(i  ,k,j), v(i+1,k,j),  &
-                                          vel                     )
-          ENDDO
-
-        ENDIF
-
-        IF( degrade_xe ) THEN
-
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts,ktf
-              fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1))      &
-                              *(v(i_end+1,k,j)+v(i_end,k,j))
-            ENDDO
-          ENDIF
-
-          i = ide-2
-          DO k=kts,ktf
-          vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
-          fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j),  &
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) &
+                                *(v(i,k,j)+v(i-1,k,j))
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
+                fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j),  &
                                         v(i  ,k,j), v(i+1,k,j),  &
                                         vel                     )
+              ENDDO
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1))      &
+                                *(v(i_end+1,k,j)+v(i_end,k,j))
+              ENDDO
+            ENDIF
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts,ktf
+                vel = 0.5*(ru(i,k,j)+ru(i,k,j-1))
+                fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j),  &
+                                        v(i  ,k,j), v(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+            ENDIF
+
           ENDDO
 
@@ -1992,5 +2033,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -1998,8 +2039,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-1)                ) degrade_ye = .false.
+       (jte < jde-2)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -2074,4 +2115,6 @@
         ENDDO
         ENDDO
+
+
       END IF
 
@@ -2167,5 +2210,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -2173,8 +2216,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-1)                ) degrade_ye = .false.
+       (jte < jde-2)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -2249,4 +2292,6 @@
         ENDDO
         ENDDO
+
+
       END IF
 
@@ -2405,4 +2450,8 @@
       ENDDO
       ENDDO
+
+   ELSE IF ( horz_order == 0 ) THEN
+
+      ! Just in case we want to turn horizontal advection off, we can do it
 
   ELSE
@@ -2530,4 +2579,6 @@
       ENDDO
 
+      !
+      !
       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
@@ -2570,4 +2621,6 @@
          DO k=kts,ktf
          DO i = i_start, i_end
+            !
+            !
             tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
          ENDDO
@@ -2613,4 +2666,6 @@
          DO k=kts,ktf
          DO i = i_start, i_end
+            !
+            !
             tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
          ENDDO
@@ -2646,4 +2701,6 @@
          DO k=kts,ktf
          DO i = i_start, i_end
+            !
+            !
             tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
          ENDDO
@@ -2679,4 +2736,6 @@
          DO k=kts,ktf
          DO i = i_start, i_end
+            !
+            !
             tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
          ENDDO
@@ -2699,6 +2758,7 @@
       DO k=kts,ktf
       DO i = i_start, i_end
+            !
+            !
             tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k))
-
       ENDDO
       ENDDO
@@ -2718,5 +2778,5 @@
 SUBROUTINE advect_scalar   ( field, field_old, tendency,    &
                              ru, rv, rom,                   &
-                             mut, config_flags,             &
+                             mut, time_step, config_flags,  &
                              msfu, msfv, msft,              &
                              fzm, fzp,                      &
@@ -2755,4 +2815,6 @@
    REAL ,                                        INTENT(IN   ) :: rdx,  &
                                                                   rdy
+   INTEGER ,                                     INTENT(IN   ) :: time_step
+
 
    ! Local data
@@ -2788,5 +2850,5 @@
       flux3(q_im2, q_im1, q_i, q_ip1, ua) =                     &
            flux4(q_im2, q_im1, q_i, q_ip1, ua) +                &
-           sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+           sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
 
       flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
@@ -2796,5 +2858,5 @@
       flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
            flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)    &
-            -sign(1.,ua)*(                             &
+            -sign(1,time_step)*sign(1.,ua)*(                    &
               (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
 
@@ -2831,5 +2893,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -2837,8 +2899,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -2865,4 +2927,5 @@
         j_end_f = jde-3
       ENDIF
+
 
 !  compute fluxes, 5th or 6th order
@@ -2884,4 +2947,5 @@
         ENDDO
 
+
       ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
 
@@ -2894,5 +2958,5 @@
             ENDDO
 
-     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
+     ELSE IF  ( j == jds+2 ) THEN  ! 4th order flux 2 in from south boundary
 
             DO k=kts,ktf
@@ -2939,4 +3003,5 @@
         ENDIF
 
+
         jtmp = jp1
         jp1 = jp0
@@ -2961,5 +3026,6 @@
       IF(degrade_xs) then
         i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
       ENDIF
 
@@ -2989,19 +3055,22 @@
         IF( degrade_xs ) THEN
 
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts,ktf
-              fqx(i,k) = 0.5*(ru(i,k,j)) &
-                     *(field(i,k,j)+field(i-1,k,j))
-
-            ENDDO
-          ENDIF
-
-          i = ids+2
-          DO k=kts,ktf
-            vel = ru(i,k,j)
-            fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                     )
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                fqx(i,k) = 0.5*(ru(i,k,j)) &
+                       *(field(i,k,j)+field(i-1,k,j))
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
+                                              field(i  ,k,j), field(i+1,k,j),  &
+                                              vel                     )
+              ENDDO
+            END IF
+
           ENDDO
 
@@ -3010,21 +3079,25 @@
         IF( degrade_xe ) THEN
 
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts,ktf
-              fqx(i,k) = 0.5*(ru(i,k,j))      &
-                     *(field(i,k,j)+field(i-1,k,j))
-            ENDDO
-         ENDIF
-
-          i = ide-2
-          DO k=kts,ktf
-            vel = ru(i,k,j)
-            fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                             )
-          ENDDO
-
-        ENDIF
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                fqx(i,k) = 0.5*(ru(i,k,j))      &
+                       *(field(i,k,j)+field(i-1,k,j))
+              ENDDO
+           ENDIF
+
+           IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+             DO k=kts,ktf
+               vel = ru(i,k,j)
+               fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
+                                       field(i  ,k,j), field(i+1,k,j),  &
+                                       vel                             )
+             ENDDO
+           ENDIF
+
+         ENDDO
+
+       ENDIF
 
 !  x flux-divergence into tendency
@@ -3055,5 +3128,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -3061,8 +3134,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -3089,4 +3162,5 @@
         j_end_f = jde-3
       ENDIF
+
 
 !  compute fluxes, 5th or 6th order
@@ -3108,4 +3182,5 @@
         ENDDO
 
+
       ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
 
@@ -3186,5 +3261,6 @@
       IF(degrade_xs) then
         i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
       ENDIF
 
@@ -3214,19 +3290,22 @@
         IF( degrade_xs ) THEN
 
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts,ktf
-              fqx(i,k) = 0.5*(ru(i,k,j)) &
-                     *(field(i,k,j)+field(i-1,k,j))
-
-            ENDDO
-          ENDIF
-
-          i = ids+2
-          DO k=kts,ktf
-            vel = ru(i,k,j)
-            fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                     )
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                fqx(i,k) = 0.5*(ru(i,k,j)) &
+                       *(field(i,k,j)+field(i-1,k,j))
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+                                              field(i  ,k,j), field(i+1,k,j),  &
+                                              vel                     )
+              ENDDO
+            END IF
+
           ENDDO
 
@@ -3235,21 +3314,25 @@
         IF( degrade_xe ) THEN
 
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts,ktf
-              fqx(i,k) = 0.5*(ru(i,k,j))      &
-                     *(field(i,k,j)+field(i-1,k,j))
-            ENDDO
-         ENDIF
-
-          i = ide-2
-          DO k=kts,ktf
-            vel = ru(i,k,j)
-            fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                             )
-          ENDDO
-
-        ENDIF
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                fqx(i,k) = 0.5*(ru(i,k,j))      &
+                       *(field(i,k,j)+field(i-1,k,j))
+              ENDDO
+           ENDIF
+
+           IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+             DO k=kts,ktf
+               vel = ru(i,k,j)
+               fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+                                       field(i  ,k,j), field(i+1,k,j),  &
+                                       vel                             )
+             ENDDO
+           ENDIF
+
+         ENDDO
+
+       ENDIF
 
 !  x flux-divergence into tendency
@@ -3274,5 +3357,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -3280,8 +3363,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !  begin flux computations
@@ -3376,4 +3459,5 @@
         j_end_f = jde-2
       ENDIF
+
 
     jp1 = 2
@@ -3392,6 +3476,9 @@
       DO k = kts, ktf
       DO i = i_start, i_end
-         fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1)          &
-                *(field(i,k,j_end+1)+field(i,k,j_end))
+         ! Assumes j>j_end_f is ONLY j_end+1 ...
+!         fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1)          &
+!                *(field(i,k,j_end+1)+field(i,k,j_end))
+         fqy(i,k,jp1) = 0.5*rv(i,k,j)          &
+                *(field(i,k,j)+field(i,k,j-1))
       ENDDO
       ENDDO
@@ -3407,6 +3494,6 @@
     END IF
 
+!  y flux-divergence into tendency
     IF ( j > j_start ) THEN
-!  y flux-divergence into tendency
 
       DO k=kts,ktf
@@ -3416,4 +3503,6 @@
       ENDDO
       ENDDO
+
+
     END IF
 
@@ -3434,5 +3523,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -3440,8 +3529,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !  begin flux computations
@@ -3536,4 +3625,5 @@
         j_end_f = jde-2
       ENDIF
+
 
     jp1 = 2
@@ -3552,6 +3642,9 @@
       DO k = kts, ktf
       DO i = i_start, i_end
-         fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1)          &
-                *(field(i,k,j_end+1)+field(i,k,j_end))
+         ! Assumes j>j_end_f is ONLY j_end+1 ...
+!         fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1)          &
+!                *(field(i,k,j_end+1)+field(i,k,j_end))
+         fqy(i,k,jp1) = 0.5*rv(i,k,j)          &
+                *(field(i,k,j)+field(i,k,j-1))
       ENDDO
       ENDDO
@@ -3567,6 +3660,6 @@
     END IF
 
+!  y flux-divergence into tendency
     IF ( j > j_start ) THEN
-!  y flux-divergence into tendency
 
       DO k=kts,ktf
@@ -3576,4 +3669,6 @@
       ENDDO
       ENDDO
+
+
     END IF
 
@@ -3610,4 +3705,5 @@
       i_end   = MIN(ite,ide-1)
 
+      !
       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
@@ -3624,4 +3720,9 @@
       ENDDO
    
+
+   ELSE IF ( horz_order == 0 ) THEN
+
+      ! Just in case we want to turn horizontal advection off, we can do it
+
    ELSE
 
@@ -3859,5 +3960,4 @@
       ENDDO
 
-
    ELSE IF (vert_order == 2) THEN    
 
@@ -3890,5 +3990,5 @@
 SUBROUTINE advect_w    ( w, w_old, tendency,            &
                          ru, rv, rom,                   &
-                         mut, config_flags,             &
+                         mut, time_step, config_flags,  &
                          msfu, msfv, msft,              &
                          fzm, fzp,                      &
@@ -3927,4 +4027,6 @@
    REAL ,                                        INTENT(IN   ) :: rdx,  &
                                                                   rdy
+   INTEGER ,                                     INTENT(IN   ) :: time_step
+
 
    ! Local data
@@ -3958,5 +4060,5 @@
       flux3(q_im2, q_im1, q_i, q_ip1, ua) =                     &
            flux4(q_im2, q_im1, q_i, q_ip1, ua) +                &
-           sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+           sign(1,time_step)*sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
 
       flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
@@ -3966,5 +4068,5 @@
       flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
            flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)    &
-            -sign(1.,ua)*(                             &
+            -sign(1,time_step)*sign(1.,ua)*(                    &
               (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
 
@@ -4001,8 +4103,585 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
+       (its > ids+3)                ) degrade_xs = .false.
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xe .or. &
+       (ite < ide-3)                ) degrade_xe = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ys .or. &
+       (jts > jds+3)                ) degrade_ys = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ye .or. &
+       (jte < jde-4)                ) degrade_ye = .false.
+
+!--------------- y - advection first
+
+      i_start = its
+      i_end   = MIN(ite,ide-1)
+      j_start = jts
+      j_end   = MIN(jte,jde-1)
+
+!  higher order flux has a 5 or 7 point stencil, so compute
+!  bounds so we can switch to second order flux close to the boundary
+
+      j_start_f = j_start
+      j_end_f   = j_end+1
+
+      IF(degrade_ys) then
+        j_start = MAX(jts,jds+1)
+        j_start_f = jds+3
+      ENDIF
+
+      IF(degrade_ye) then
+        j_end = MIN(jte,jde-2)
+        j_end_f = jde-3
+      ENDIF
+
+
+!  compute fluxes, 5th or 6th order
+
+     jp1 = 2
+     jp0 = 1
+
+     j_loop_y_flux_6 : DO j = j_start, j_end+1
+
+      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
+
+        DO k=kts+1,ktf
+        DO i = i_start, i_end
+          vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+          fqy( i, k, jp1 ) = vel*flux6(                     &
+                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
+                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
+        ENDDO
+        ENDDO
+
+        k = ktf+1
+        DO i = i_start, i_end
+          vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+          fqy( i, k, jp1 ) = vel*flux6(                     &
+                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
+                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
+        ENDDO
+
+      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*          &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*          &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+
+     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+              fqy( i, k, jp1 ) = vel*flux4(              &
+                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+              fqy( i, k, jp1 ) = vel*flux4(              &
+                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+
+     ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*      &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*      &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+
+     ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+              fqy( i, k, jp1 ) = vel*flux4(             &
+                   w(i,k,j-2),w(i,k,j-1),    &
+                   w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+              fqy( i, k, jp1 ) = vel*flux4(             &
+                   w(i,k,j-2),w(i,k,j-1),    &
+                   w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+
+     ENDIF
+
+!  y flux-divergence into tendency
+
+        IF(j > j_start) THEN
+
+          DO k=kts+1,ktf+1
+          DO i = i_start, i_end
+            mrdy=msft(i,j-1)*rdy
+            tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
+          ENDDO
+          ENDDO
+
+       ENDIF
+
+
+        jtmp = jp1
+        jp1 = jp0
+        jp0 = jtmp
+
+      ENDDO j_loop_y_flux_6
+
+!  next, x - flux divergence
+
+      i_start = its
+      i_end   = MIN(ite,ide-1)
+
+      j_start = jts
+      j_end   = MIN(jte,jde-1)
+
+!  higher order flux has a 5 or 7 point stencil, so compute
+!  bounds so we can switch to second order flux close to the boundary
+
+      i_start_f = i_start
+      i_end_f   = i_end+1
+
+      IF(degrade_xs) then
+        i_start = MAX(ids+1,its)
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
+      ENDIF
+
+      IF(degrade_xe) then
+        i_end = MIN(ide-2,ite)
+        i_end_f = ide-3
+      ENDIF
+
+!  compute fluxes
+
+      DO j = j_start, j_end
+
+!  5th or 6th order flux
+
+        DO k=kts+1,ktf
+        DO i = i_start_f, i_end_f
+          vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+          fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j),  &
+                                  w(i-1,k,j), w(i  ,k,j),  &
+                                  w(i+1,k,j), w(i+2,k,j),  &
+                                  vel                     )
+        ENDDO
+        ENDDO
+
+        k = ktf+1
+        DO i = i_start_f, i_end_f
+          vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+          fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j),  &
+                                  w(i-1,k,j), w(i  ,k,j),  &
+                                  w(i+1,k,j), w(i+2,k,j),  &
+                                  vel                     )
+        ENDDO
+
+!  lower order fluxes close to boundaries (if not periodic or symmetric)
+
+        IF( degrade_xs ) THEN
+
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts+1,ktf
+                fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
+                                *(w(i,k,j)+w(i-1,k,j))
+              ENDDO
+              k = ktf+1
+              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
+                     *(w(i,k,j)+w(i-1,k,j))
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts+1,ktf
+                vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+                fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
+                                        w(i  ,k,j), w(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+              k = ktf+1
+              vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+              fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
+                                      w(i  ,k,j), w(i+1,k,j),  &
+                                      vel                     )
+            END IF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts+1,ktf
+                fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j))      &
+                                  *(w(i,k,j)+w(i-1,k,j))
+              ENDDO
+              k = ktf+1
+              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j))      &
+                     *(w(i,k,j)+w(i-1,k,j))
+            ENDIF
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts+1,ktf
+                vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+                fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
+                                        w(i  ,k,j), w(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+              k = ktf+1
+              vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+              fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
+                                      w(i  ,k,j), w(i+1,k,j),  &
+                                      vel                     )
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+!  x flux-divergence into tendency
+
+        DO k=kts+1,ktf+1
+          DO i = i_start, i_end
+            mrdx=msft(i,j)*rdx
+            tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
+          ENDDO
+        ENDDO
+
+      ENDDO
+
+ELSE IF (horz_order == 5 ) THEN
+
+!  determine boundary mods for flux operators
+!  We degrade the flux operators from 3rd/4th order
+!   to second order one gridpoint in from the boundaries for
+!   all boundary conditions except periodic and symmetry - these
+!   conditions have boundary zone data fill for correct application
+!   of the higher order flux stencils
+
+   degrade_xs = .true.
+   degrade_xe = .true.
+   degrade_ys = .true.
+   degrade_ye = .true.
+
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xs .or. &
+       (its > ids+3)                ) degrade_xs = .false.
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xe .or. &
+       (ite < ide-3)                ) degrade_xe = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ys .or. &
+       (jts > jds+3)                ) degrade_ys = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ye .or. &
+       (jte < jde-4)                ) degrade_ye = .false.
+
+!--------------- y - advection first
+
+      i_start = its
+      i_end   = MIN(ite,ide-1)
+      j_start = jts
+      j_end   = MIN(jte,jde-1)
+
+!  higher order flux has a 5 or 7 point stencil, so compute
+!  bounds so we can switch to second order flux close to the boundary
+
+      j_start_f = j_start
+      j_end_f   = j_end+1
+
+      IF(degrade_ys) then
+        j_start = MAX(jts,jds+1)
+        j_start_f = jds+3
+      ENDIF
+
+      IF(degrade_ye) then
+        j_end = MIN(jte,jde-2)
+        j_end_f = jde-3
+      ENDIF
+
+
+!  compute fluxes, 5th or 6th order
+
+     jp1 = 2
+     jp0 = 1
+
+     j_loop_y_flux_5 : DO j = j_start, j_end+1
+
+      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
+
+        DO k=kts+1,ktf
+        DO i = i_start, i_end
+          vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+          fqy( i, k, jp1 ) = vel*flux5(                     &
+                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
+                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
+        ENDDO
+        ENDDO
+
+        k = ktf+1
+        DO i = i_start, i_end
+          vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+          fqy( i, k, jp1 ) = vel*flux5(                     &
+                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
+                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
+        ENDDO
+
+      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*          &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*          &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+
+     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+              fqy( i, k, jp1 ) = vel*flux3(              &
+                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+              fqy( i, k, jp1 ) = vel*flux3(              &
+                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+
+     ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*      &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*      &
+                     (w(i,k,j)+w(i,k,j-1))
+            ENDDO
+
+     ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
+
+            DO k=kts+1,ktf
+            DO i = i_start, i_end
+              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
+              fqy( i, k, jp1 ) = vel*flux3(             &
+                   w(i,k,j-2),w(i,k,j-1),    &
+                   w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+            ENDDO
+
+            k = ktf+1
+            DO i = i_start, i_end
+              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
+              fqy( i, k, jp1 ) = vel*flux3(             &
+                   w(i,k,j-2),w(i,k,j-1),    &
+                   w(i,k,j),w(i,k,j+1),vel )
+            ENDDO
+
+     ENDIF
+
+!  y flux-divergence into tendency
+
+        IF(j > j_start) THEN
+
+          DO k=kts+1,ktf+1
+          DO i = i_start, i_end
+            mrdy=msft(i,j-1)*rdy
+            tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
+          ENDDO
+          ENDDO
+
+       ENDIF
+
+
+        jtmp = jp1
+        jp1 = jp0
+        jp0 = jtmp
+
+      ENDDO j_loop_y_flux_5
+
+!  next, x - flux divergence
+
+      i_start = its
+      i_end   = MIN(ite,ide-1)
+
+      j_start = jts
+      j_end   = MIN(jte,jde-1)
+
+!  higher order flux has a 5 or 7 point stencil, so compute
+!  bounds so we can switch to second order flux close to the boundary
+
+      i_start_f = i_start
+      i_end_f   = i_end+1
+
+      IF(degrade_xs) then
+        i_start = MAX(ids+1,its)
+!        i_start_f = i_start+2
+        i_start_f = MIN(i_start+2,ids+3)
+      ENDIF
+
+      IF(degrade_xe) then
+        i_end = MIN(ide-2,ite)
+        i_end_f = ide-3
+      ENDIF
+
+!  compute fluxes
+
+      DO j = j_start, j_end
+
+!  5th or 6th order flux
+
+        DO k=kts+1,ktf
+        DO i = i_start_f, i_end_f
+          vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+          fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j),  &
+                                  w(i-1,k,j), w(i  ,k,j),  &
+                                  w(i+1,k,j), w(i+2,k,j),  &
+                                  vel                     )
+        ENDDO
+        ENDDO
+
+        k = ktf+1
+        DO i = i_start_f, i_end_f
+          vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+          fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j),  &
+                                  w(i-1,k,j), w(i  ,k,j),  &
+                                  w(i+1,k,j), w(i+2,k,j),  &
+                                  vel                     )
+        ENDDO
+
+!  lower order fluxes close to boundaries (if not periodic or symmetric)
+
+        IF( degrade_xs ) THEN
+
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts+1,ktf
+                fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
+                                *(w(i,k,j)+w(i-1,k,j))
+              ENDDO
+              k = ktf+1
+              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
+                     *(w(i,k,j)+w(i-1,k,j))
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts+1,ktf
+                vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+                fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
+                                        w(i  ,k,j), w(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+              k = ktf+1
+              vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+              fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
+                                      w(i  ,k,j), w(i+1,k,j),  &
+                                      vel                     )
+            END IF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts+1,ktf
+                fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j))      &
+                                  *(w(i,k,j)+w(i-1,k,j))
+              ENDDO
+              k = ktf+1
+              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j))      &
+                     *(w(i,k,j)+w(i-1,k,j))
+            ENDIF
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts+1,ktf
+                vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
+                fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
+                                        w(i  ,k,j), w(i+1,k,j),  &
+                                        vel                     )
+              ENDDO
+              k = ktf+1
+              vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
+              fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
+                                      w(i  ,k,j), w(i+1,k,j),  &
+                                      vel                     )
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+!  x flux-divergence into tendency
+
+        DO k=kts+1,ktf+1
+          DO i = i_start, i_end
+            mrdx=msft(i,j)*rdx
+            tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
+          ENDDO
+        ENDDO
+
+      ENDDO
+
+ELSE IF ( horz_order == 4 ) THEN
+
+   degrade_xs = .true.
+   degrade_xe = .true.
+   degrade_ys = .true.
+   degrade_ye = .true.
+
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xs .or. &
        (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
-       (ite < ide-3)                ) degrade_xe = .false.
+       (ite < ide-2)                ) degrade_xe = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
@@ -4012,559 +4691,4 @@
        (jte < jde-3)                ) degrade_ye = .false.
 
-!--------------- y - advection first
-
-      i_start = its
-      i_end   = MIN(ite,ide-1)
-      j_start = jts
-      j_end   = MIN(jte,jde-1)
-
-!  higher order flux has a 5 or 7 point stencil, so compute
-!  bounds so we can switch to second order flux close to the boundary
-
-      j_start_f = j_start
-      j_end_f   = j_end+1
-
-      IF(degrade_ys) then
-        j_start = MAX(jts,jds+1)
-        j_start_f = jds+3
-      ENDIF
-
-      IF(degrade_ye) then
-        j_end = MIN(jte,jde-2)
-        j_end_f = jde-3
-      ENDIF
-
-!  compute fluxes, 5th or 6th order
-
-     jp1 = 2
-     jp0 = 1
-
-     j_loop_y_flux_6 : DO j = j_start, j_end+1
-
-      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
-
-        DO k=kts+1,ktf
-        DO i = i_start, i_end
-          vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-          fqy( i, k, jp1 ) = vel*flux6(                     &
-                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
-                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
-        ENDDO
-        ENDDO
-
-        k = ktf+1
-        DO i = i_start, i_end
-          vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-          fqy( i, k, jp1 ) = vel*flux6(                     &
-                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
-                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
-        ENDDO
-
-      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*          &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-
-     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-              fqy( i, k, jp1 ) = vel*flux4(              &
-                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-              fqy( i, k, jp1 ) = vel*flux4(              &
-                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-
-     ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*      &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*      &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-
-     ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-              fqy( i, k, jp1 ) = vel*flux4(             &
-                   w(i,k,j-2),w(i,k,j-1),    &
-                   w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-              fqy( i, k, jp1 ) = vel*flux4(             &
-                   w(i,k,j-2),w(i,k,j-1),    &
-                   w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-
-     ENDIF
-
-!  y flux-divergence into tendency
-
-        IF(j > j_start) THEN
-
-          DO k=kts+1,ktf+1
-          DO i = i_start, i_end
-            mrdy=msft(i,j-1)*rdy
-            tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
-          ENDDO
-          ENDDO
-
-       ENDIF
-
-        jtmp = jp1
-        jp1 = jp0
-        jp0 = jtmp
-
-      ENDDO j_loop_y_flux_6
-
-!  next, x - flux divergence
-
-      i_start = its
-      i_end   = MIN(ite,ide-1)
-
-      j_start = jts
-      j_end   = MIN(jte,jde-1)
-
-!  higher order flux has a 5 or 7 point stencil, so compute
-!  bounds so we can switch to second order flux close to the boundary
-
-      i_start_f = i_start
-      i_end_f   = i_end+1
-
-      IF(degrade_xs) then
-        i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
-      ENDIF
-
-      IF(degrade_xe) then
-        i_end = MIN(ide-2,ite)
-        i_end_f = ide-3
-      ENDIF
-
-!  compute fluxes
-
-      DO j = j_start, j_end
-
-!  5th or 6th order flux
-
-        DO k=kts+1,ktf
-        DO i = i_start_f, i_end_f
-          vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-          fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j),  &
-                                         w(i-1,k,j), w(i  ,k,j),  &
-                                         w(i+1,k,j), w(i+2,k,j),  &
-                                         vel                             )
-        ENDDO
-        ENDDO
-
-        k = ktf+1
-        DO i = i_start_f, i_end_f
-          vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-          fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j),  &
-                                         w(i-1,k,j), w(i  ,k,j),  &
-                                         w(i+1,k,j), w(i+2,k,j),  &
-                                         vel                             )
-        ENDDO
-
-!  lower order fluxes close to boundaries (if not periodic or symmetric)
-
-        IF( degrade_xs ) THEN
-
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts+1,ktf
-              fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
-                     *(w(i,k,j)+w(i-1,k,j))
-            ENDDO
-              k = ktf+1
-              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
-                     *(w(i,k,j)+w(i-1,k,j))
-          ENDIF
-
-          DO k=kts+1,ktf
-            i = i_start+1
-            vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-            fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                     )
-          ENDDO
-
-            k = ktf+1
-            vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-            fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                     )
-        ENDIF
-
-        IF( degrade_xe ) THEN
-
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts+1,ktf
-              fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j))      &
-                     *(w(i,k,j)+w(i-1,k,j))
-            ENDDO
-              k = ktf+1
-              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j))      &
-                     *(w(i,k,j)+w(i-1,k,j))
-          ENDIF
-
-          i = ide-2
-          DO k=kts+1,ktf
-            vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-            fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                             )
-          ENDDO
-
-            k = ktf+1
-            vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-            fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                             )
-        ENDIF
-
-!  x flux-divergence into tendency
-
-        DO k=kts+1,ktf+1
-          DO i = i_start, i_end
-            mrdx=msft(i,j)*rdx
-            tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
-          ENDDO
-        ENDDO
-
-      ENDDO
-
-
-ELSE IF (horz_order == 5 ) THEN
-
-!  determine boundary mods for flux operators
-!  We degrade the flux operators from 3rd/4th order
-!   to second order one gridpoint in from the boundaries for
-!   all boundary conditions except periodic and symmetry - these
-!   conditions have boundary zone data fill for correct application
-!   of the higher order flux stencils
-
-   degrade_xs = .true.
-   degrade_xe = .true.
-   degrade_ys = .true.
-   degrade_ye = .true.
-
-   IF( config_flags%periodic_x   .or. &
-       config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
-   IF( config_flags%periodic_x   .or. &
-       config_flags%symmetric_xe .or. &
-       (ite < ide-3)                ) degrade_xe = .false.
-   IF( config_flags%periodic_y   .or. &
-       config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
-   IF( config_flags%periodic_y   .or. &
-       config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
-
-!--------------- y - advection first
-
-      i_start = its
-      i_end   = MIN(ite,ide-1)
-      j_start = jts
-      j_end   = MIN(jte,jde-1)
-
-!  higher order flux has a 5 or 7 point stencil, so compute
-!  bounds so we can switch to second order flux close to the boundary
-
-      j_start_f = j_start
-      j_end_f   = j_end+1
-
-      IF(degrade_ys) then
-        j_start = MAX(jts,jds+1)
-        j_start_f = jds+3
-      ENDIF
-
-      IF(degrade_ye) then
-        j_end = MIN(jte,jde-2)
-        j_end_f = jde-3
-      ENDIF
-
-!  compute fluxes, 5th or 6th order
-
-     jp1 = 2
-     jp0 = 1
-
-     j_loop_y_flux_5 : DO j = j_start, j_end+1
-
-      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN
-
-        DO k=kts+1,ktf
-        DO i = i_start, i_end
-          vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-          fqy( i, k, jp1 ) = vel*flux5(                     &
-                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
-                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
-        ENDDO
-        ENDDO
-
-        k = ktf+1
-        DO i = i_start, i_end
-          vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-          fqy( i, k, jp1 ) = vel*flux5(                     &
-                  w(i,k,j-3), w(i,k,j-2), w(i,k,j-1),       &
-                  w(i,k,j  ), w(i,k,j+1), w(i,k,j+2),  vel )
-        ENDDO
-
-      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*          &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*          &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-
-     ELSE IF  ( j == jds+2 ) THEN  ! third of 4th order flux 2 in from south boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-              fqy( i, k, jp1 ) = vel*flux3(              &
-                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-              fqy( i, k, jp1 ) = vel*flux3(              &
-                   w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-
-     ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))*      &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))*      &
-                     (w(i,k,j)+w(i,k,j-1))
-            ENDDO
-
-     ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
-
-            DO k=kts+1,ktf
-            DO i = i_start, i_end
-              vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)
-              fqy( i, k, jp1 ) = vel*flux3(             &
-                   w(i,k,j-2),w(i,k,j-1),    &
-                   w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-            ENDDO
-
-            k = ktf+1
-            DO i = i_start, i_end
-              vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)
-              fqy( i, k, jp1 ) = vel*flux3(             &
-                   w(i,k,j-2),w(i,k,j-1),    &
-                   w(i,k,j),w(i,k,j+1),vel )
-            ENDDO
-
-     ENDIF
-
-!  y flux-divergence into tendency
-
-        IF(j > j_start) THEN
-
-          DO k=kts+1,ktf+1
-          DO i = i_start, i_end
-            mrdy=msft(i,j-1)*rdy
-            tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0))
-          ENDDO
-          ENDDO
-
-       ENDIF
-
-        jtmp = jp1
-        jp1 = jp0
-        jp0 = jtmp
-
-      ENDDO j_loop_y_flux_5
-
-!  next, x - flux divergence
-
-      i_start = its
-      i_end   = MIN(ite,ide-1)
-
-      j_start = jts
-      j_end   = MIN(jte,jde-1)
-
-!  higher order flux has a 5 or 7 point stencil, so compute
-!  bounds so we can switch to second order flux close to the boundary
-
-      i_start_f = i_start
-      i_end_f   = i_end+1
-
-      IF(degrade_xs) then
-        i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
-      ENDIF
-
-      IF(degrade_xe) then
-        i_end = MIN(ide-2,ite)
-        i_end_f = ide-3
-      ENDIF
-
-!  compute fluxes
-
-      DO j = j_start, j_end
-
-!  5th or 6th order flux
-
-        DO k=kts+1,ktf
-        DO i = i_start_f, i_end_f
-          vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-          fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j),  &
-                                  w(i-1,k,j), w(i  ,k,j),  &
-                                  w(i+1,k,j), w(i+2,k,j),  &
-                          vel                             )
-        ENDDO
-        ENDDO
-
-        k = ktf+1
-        DO i = i_start_f, i_end_f
-          vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-          fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j),  &
-                                  w(i-1,k,j), w(i  ,k,j),  &
-                                  w(i+1,k,j), w(i+2,k,j),  &
-                          vel                             )
-        ENDDO
-
-!  lower order fluxes close to boundaries (if not periodic or symmetric)
-
-        IF( degrade_xs ) THEN
-
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts+1,ktf
-              fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) &
-                     *(w(i,k,j)+w(i-1,k,j))
-            ENDDO
-              k = ktf+1
-              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) &
-                     *(w(i,k,j)+w(i-1,k,j))
-          ENDIF
-
-          i = i_start+1
-          DO k=kts+1,ktf
-            vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-            fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
-                                    w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                     )
-          ENDDO
-            k = ktf+1
-            vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-            fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
-                                    w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                     )
-
-        ENDIF
-
-        IF( degrade_xe ) THEN
-
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts+1,ktf
-              fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j))      &
-                     *(w(i,k,j)+w(i-1,k,j))
-            ENDDO
-              k = ktf+1
-              fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j))      &
-                     *(w(i,k,j)+w(i-1,k,j))
-          ENDIF
-
-          i = ide-2
-          DO k=kts+1,ktf
-            vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)
-            fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                             )
-          ENDDO
-            k = ktf+1
-            vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)
-            fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j),  &
-                                          w(i  ,k,j), w(i+1,k,j),  &
-                                          vel                             )
-        ENDIF
-
-!  x flux-divergence into tendency
-
-        DO k=kts+1,ktf+1
-          DO i = i_start, i_end
-            mrdx=msft(i,j)*rdx
-            tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k))
-          ENDDO
-        ENDDO
-
-      ENDDO
-
-ELSE IF ( horz_order == 4 ) THEN
-
-   degrade_xs = .true.
-   degrade_xe = .true.
-   degrade_ys = .true.
-   degrade_ye = .true.
-
-   IF( config_flags%periodic_x   .or. &
-       config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
-   IF( config_flags%periodic_x   .or. &
-       config_flags%symmetric_xe .or. &
-       (ite < ide-2)                ) degrade_xe = .false.
-   IF( config_flags%periodic_y   .or. &
-       config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
-   IF( config_flags%periodic_y   .or. &
-       config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
-
 !  begin flux computations
 !  start with x flux divergence
@@ -4676,4 +4800,5 @@
       ENDIF
 
+
         jp1 = 2
         jp0 = 1
@@ -4698,14 +4823,22 @@
           DO k = kts+1, ktf
           DO i = i_start, i_end
+            ! Assumes j>j_end_f is ONLY j_end+1 ...
+!            fqy(i, k, jp1) =                             &
+!               0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1))     &
+!                   *(w(i,k,j_end+1)+w(i,k,j_end))
             fqy(i, k, jp1) =                             &
-               0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1))     &
-                   *(w(i,k,j_end+1)+w(i,k,j_end))
+               0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))     &
+                   *(w(i,k,j)+w(i,k,j-1))
           ENDDO
           ENDDO
           k = ktf+1
           DO i = i_start, i_end
+            ! Assumes j>j_end_f is ONLY j_end+1 ...
+!            fqy(i, k, jp1) =                                         &
+!               0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1))     &
+!                   *(w(i,k,j_end+1)+w(i,k,j_end))
             fqy(i, k, jp1) =                                         &
-               0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1))     &
-                   *(w(i,k,j_end+1)+w(i,k,j_end))
+               0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))     &
+                   *(w(i,k,j)+w(i,k,j-1))
           ENDDO
        ELSE
@@ -4728,6 +4861,7 @@
        END IF
 
+!  y flux-divergence into tendency
        IF( j > j_start ) THEN
-!  y flux-divergence into tendency
+
           DO k = kts+1, ktf+1
           DO i = i_start, i_end
@@ -4736,4 +4870,6 @@
           ENDDO
           ENDDO
+
+
        END IF
 
@@ -4753,5 +4889,5 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
@@ -4759,8 +4895,8 @@
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-3)                ) degrade_ye = .false.
 
 !  begin flux computations
@@ -4873,4 +5009,5 @@
       ENDIF
 
+
         jp1 = 2
         jp0 = 1
@@ -4895,14 +5032,22 @@
           DO k = kts+1, ktf
           DO i = i_start, i_end
+            ! Assumes j>j_end_f is ONLY j_end+1 ...
+!            fqy(i, k, jp1) =                             &
+!               0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1))     &
+!                   *(w(i,k,j_end+1)+w(i,k,j_end))
             fqy(i, k, jp1) =                             &
-               0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1))     &
-                   *(w(i,k,j_end+1)+w(i,k,j_end))
+               0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))     &
+                   *(w(i,k,j)+w(i,k,j-1))
           ENDDO
           ENDDO
           k = ktf+1
           DO i = i_start, i_end
+            ! Assumes j>j_end_f is ONLY j_end+1 ...
+!            fqy(i, k, jp1) =                             &
+!               0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1))     &
+!                   *(w(i,k,j_end+1)+w(i,k,j_end))
             fqy(i, k, jp1) =                             &
-               0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1))     &
-                   *(w(i,k,j_end+1)+w(i,k,j_end))
+               0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))     &
+                   *(w(i,k,j)+w(i,k,j-1))
           ENDDO
        ELSE
@@ -4925,6 +5070,7 @@
        END IF
 
+!  y flux-divergence into tendency
        IF( j > j_start ) THEN
-!  y flux-divergence into tendency
+
           DO k = kts+1, ktf+1
           DO i = i_start, i_end
@@ -4933,4 +5079,6 @@
           ENDDO
           ENDDO
+
+
        END IF
 
@@ -4985,4 +5133,5 @@
       i_start = its
       i_end   = MIN(ite,ide-1)
+      !
       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
@@ -5017,4 +5166,8 @@
 
       ENDDO
+
+   ELSE IF ( horz_order == 0 ) THEN
+
+      ! Just in case we want to turn horizontal advection off, we can do it
 
    ELSE
@@ -5501,4 +5654,7 @@
       flux_upwind(q_im1, q_i, cr ) = 0.5*min( 1.0,(cr+abs(cr)))*q_im1 &
                                     +0.5*max(-1.0,(cr-abs(cr)))*q_i
+
+!      flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 &
+!                                    +0.5*(1.-sign(1.,cr))*q_i
 !      flux_upwind(q_im1, q_i, cr ) = 0.
 
@@ -5535,14 +5691,14 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
-       (ite < ide-3)                ) degrade_xe = .false.
+       (ite < ide-4)                ) degrade_xe = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -5560,14 +5716,16 @@
 !--  modify loop bounds if open or specified
 
-      IF(degrade_xs) i_start = its
-      IF(degrade_xe) i_end   = MIN(ite,ide-1)
+!      IF(degrade_xs) i_start = MAX(its-1,ids-1)
+!      IF(degrade_xe) i_end   = MIN(ite+1,ide-2)
+      IF(degrade_xs) i_start = MAX(its-1,ids)
+      IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
 
       IF(degrade_ys) then
-        j_start = MAX(jts,jds+1)
+        j_start = MAX(jts-1,jds+1)
         j_start_f = jds+3
       ENDIF
 
       IF(degrade_ye) then
-        j_end = MIN(jte,jde-2)
+        j_end = MIN(jte+1,jde-2)
         j_end_f = jde-3
       ENDIF
@@ -5652,5 +5810,5 @@
             ENDDO
 
-      ELSE IF ( j == jde-2 ) THEN  ! 4th order flux 2 in from north boundary
+      ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4th order flux 2 in from north boundary
 
             DO k=kts,ktf
@@ -5673,5 +5831,5 @@
       ENDIF
 
-    ENDDO j_loop_y_flux_6
+   ENDDO j_loop_y_flux_6
 
 !  next, x flux
@@ -5689,14 +5847,16 @@
 !--  modify loop bounds for open and specified b.c
 
-      IF(degrade_ys) j_start = jts
-      IF(degrade_ye) j_end   = MIN(jte,jde-1)
+!      IF(degrade_ys) j_start = MAX(jts-1,jds+1)
+!      IF(degrade_ye) j_end   = MIN(jte+1,jde-2)
+      IF(degrade_ys) j_start = MAX(jts-1,jds)
+      IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
 
       IF(degrade_xs) then
-        i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
+        i_start = MAX(ids+1,its-1)
+        i_start_f = ids+3
       ENDIF
 
       IF(degrade_xe) then
-        i_end = MIN(ide-2,ite)
+        i_end = MIN(ide-2,ite+1)
         i_end_f = ide-3
       ENDIF
@@ -5706,5 +5866,5 @@
       DO j = j_start, j_end
 
-!  6th order flux
+!  5th order flux
 
         DO k=kts,ktf
@@ -5730,67 +5890,68 @@
         IF( degrade_xs ) THEN
 
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts,ktf
-
-              dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-              mu = 0.5*(mut(i,j)+mut(i-1,j))
-              vel = ru(i,k,j)/mu
-              cr = vel*dt/dx
-              fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-
-              fqx(i,k,j) = 0.5*(ru(i,k,j)) &
-                     *(field(i,k,j)+field(i-1,k,j))
-
-              fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-            ENDDO
-          ENDIF
-
-          i = ids+2
-          DO k=kts,ktf
-            dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-            mu = 0.5*(mut(i,j)+mut(i-1,j))
-            vel = ru(i,k,j)
-            cr = vel*dt/dx/mu
-            fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-            fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                     )
-            fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-          ENDDO
-
-        ENDIF
-
-        IF( degrade_xe ) THEN
-
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts,ktf
-              dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-              mu = 0.5*(mut(i,j)+mut(i-1,j))
-              vel = ru(i,k,j)
-              cr = vel*dt/dx/mu
-              fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-              fqx(i,k,j) = 0.5*(ru(i,k,j))      &
-                     *(field(i,k,j)+field(i-1,k,j))
-              fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-            ENDDO
-          ENDIF
-
-          i = ide-2
-          DO k=kts,ktf
-
-            dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-            mu = 0.5*(mut(i,j)+mut(i-1,j))
-            vel = ru(i,k,j)
-            cr = vel*dt/dx/mu
-            fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-            fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)/mu
+                cr = vel*dt/dx
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx(i,k,j) = 0.5*(ru(i,k,j)) &
+                       *(field(i,k,j)+field(i-1,k,j))
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! fourth order
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
                                           field(i  ,k,j), field(i+1,k,j),  &
                                           vel                             )
-            fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx(i,k,j) = 0.5*(ru(i,k,j))      &
+                       *(field(i,k,j)+field(i-1,k,j))
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+
+            IF( i == ide-2 ) THEN ! fourth order flux one in from the boundary
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j),  &
+                                          field(i  ,k,j), field(i+1,k,j),  &
+                                          vel                             )
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
 
           ENDDO
@@ -5806,14 +5967,14 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+2)                ) degrade_xs = .false.
+       (its > ids+3)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
-       (ite < ide-3)                ) degrade_xe = .false.
+       (ite < ide-4)                ) degrade_xe = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+2)                ) degrade_ys = .false.
+       (jts > jds+3)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-3)                ) degrade_ye = .false.
+       (jte < jde-4)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -5831,14 +5992,16 @@
 !--  modify loop bounds if open or specified
 
-      IF(degrade_xs) i_start = its
-      IF(degrade_xe) i_end   = MIN(ite,ide-1)
+!      IF(degrade_xs) i_start = MAX(its-1,ids-1)
+!      IF(degrade_xe) i_end   = MIN(ite+1,ide-2)
+      IF(degrade_xs) i_start = MAX(its-1,ids)
+      IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
 
       IF(degrade_ys) then
-        j_start = MAX(jts,jds+1)
+        j_start = MAX(jts-1,jds+1)
         j_start_f = jds+3
       ENDIF
 
       IF(degrade_ye) then
-        j_end = MIN(jte,jde-2)
+        j_end = MIN(jte+1,jde-2)
         j_end_f = jde-3
       ENDIF
@@ -5960,14 +6123,16 @@
 !--  modify loop bounds for open and specified b.c
 
-      IF(degrade_ys) j_start = jts
-      IF(degrade_ye) j_end   = MIN(jte,jde-1)
+!      IF(degrade_ys) j_start = MAX(jts-1,jds+1)
+!      IF(degrade_ye) j_end   = MIN(jte+1,jde-2)
+      IF(degrade_ys) j_start = MAX(jts-1,jds)
+      IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
 
       IF(degrade_xs) then
-        i_start = MAX(ids+1,its)
-        i_start_f = i_start+2
+        i_start = MAX(ids+1,its-1)
+        i_start_f = ids+3
       ENDIF
 
       IF(degrade_xe) then
-        i_end = MIN(ide-2,ite)
+        i_end = MIN(ide-2,ite+1)
         i_end_f = ide-3
       ENDIF
@@ -6001,67 +6166,68 @@
         IF( degrade_xs ) THEN
 
-          IF( i_start == ids+1 ) THEN ! second order flux next to the boundary
-            i = ids+1
-            DO k=kts,ktf
-
-              dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-              mu = 0.5*(mut(i,j)+mut(i-1,j))
-              vel = ru(i,k,j)/mu
-              cr = vel*dt/dx
-              fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-
-              fqx(i,k,j) = 0.5*(ru(i,k,j)) &
-                     *(field(i,k,j)+field(i-1,k,j))
-
-              fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-            ENDDO
-          ENDIF
-
-          i = ids+2
-          DO k=kts,ktf
-            dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-            mu = 0.5*(mut(i,j)+mut(i-1,j))
-            vel = ru(i,k,j)
-            cr = vel*dt/dx/mu
-            fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-            fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                     )
-            fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-          ENDDO
-
-        ENDIF
-
-        IF( degrade_xe ) THEN
-
-          IF( i_end == ide-2 ) THEN ! second order flux next to the boundary
-            i = ide-1
-            DO k=kts,ktf
-              dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-              mu = 0.5*(mut(i,j)+mut(i-1,j))
-              vel = ru(i,k,j)
-              cr = vel*dt/dx/mu
-              fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-              fqx(i,k,j) = 0.5*(ru(i,k,j))      &
-                     *(field(i,k,j)+field(i-1,k,j))
-              fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
-
-            ENDDO
-          ENDIF
-
-          i = ide-2
-          DO k=kts,ktf
-
-            dx = 2./(msft(i,j)+msft(i-1,j))/rdx
-            mu = 0.5*(mut(i,j)+mut(i-1,j))
-            vel = ru(i,k,j)
-            cr = vel*dt/dx/mu
-            fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-            fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)/mu
+                cr = vel*dt/dx
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx(i,k,j) = 0.5*(ru(i,k,j)) &
+                       *(field(i,k,j)+field(i-1,k,j))
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
                                           field(i  ,k,j), field(i+1,k,j),  &
                                           vel                             )
-            fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx(i,k,j) = 0.5*(ru(i,k,j))      &
+                       *(field(i,k,j)+field(i-1,k,j))
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
+
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts,ktf
+                dx = 2./(msft(i,j)+msft(i-1,j))/rdx
+                mu = 0.5*(mut(i,j)+mut(i-1,j))
+                vel = ru(i,k,j)
+                cr = vel*dt/dx/mu
+                fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+                                          field(i  ,k,j), field(i+1,k,j),  &
+                                          vel                             )
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+              ENDDO
+            ENDIF
 
           ENDDO
@@ -6275,14 +6441,14 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids+1)                ) degrade_xs = .false.
+       (its > ids+2)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
-       (ite < ide-2)                ) degrade_xe = .false.
+       (ite < ide-1)                ) degrade_xe = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds+1)                ) degrade_ys = .false.
+       (jts > jds+2)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-2)                ) degrade_ye = .false.
+       (jte < jde-1)                ) degrade_ye = .false.
 
 !--------------- y - advection first
@@ -6475,14 +6641,14 @@
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xs .or. &
-       (its > ids)                ) degrade_xs = .false.
+       (its > ids+1)                ) degrade_xs = .false.
    IF( config_flags%periodic_x   .or. &
        config_flags%symmetric_xe .or. &
-       (ite < ide-1)                ) degrade_xe = .false.
+       (ite < ide-2)                ) degrade_xe = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ys .or. &
-       (jts > jds)                ) degrade_ys = .false.
+       (jts > jds+1)                ) degrade_ys = .false.
    IF( config_flags%periodic_y   .or. &
        config_flags%symmetric_ye .or. &
-       (jte < jde-1)                ) degrade_ye = .false.
+       (jte < jde-2)                ) degrade_ye = .false.
 
 !--  y flux compute; these bounds are for periodic and sym b.c.
@@ -6530,7 +6696,7 @@
             cr = vel*dt/dx/mu
             fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
-            fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
-                                          field(i  ,k,j), field(i+1,k,j),  &
-                                          vel                             )
+            fqx( i,k,j ) = 0.5*ru(i,k,j)*          &
+                  (field(i,k,j)+field(i-1,k,j))
+
             fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
         ENDDO
@@ -6538,5 +6704,5 @@
       ENDDO
 
-!--- end of 3nd order horizontal flux calculation
+!--- end of 2nd order horizontal flux calculation
 
    ELSE
@@ -6633,8 +6799,8 @@
 !-- loop bounds for open or specified conditions
 
-    IF(degrade_xs) i_start = its
-    IF(degrade_xe) i_end   = MIN(ite,ide-1)
-    IF(degrade_ys) j_start = jts
-    IF(degrade_ye) j_end   = MIN(jte,jde-1)
+    IF(degrade_xs) i_start = MAX(its-1,ids)
+    IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
+    IF(degrade_ys) j_start = MAX(jts-1,jds)
+    IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
 
     vert_order_test : IF (vert_order == 6) THEN    
@@ -6931,27 +7097,27 @@
 !-- loop bounds for open or specified conditions
 
-   IF(degrade_xs) i_start = its
-   IF(degrade_xe) i_end   = MIN(ite,ide-1)
-   IF(degrade_ys) j_start = jts
-   IF(degrade_ye) j_end   = MIN(jte,jde-1)
+   IF(degrade_xs) i_start = MAX(its-1,ids)
+   IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
+   IF(degrade_ys) j_start = MAX(jts-1,jds)
+   IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
 
    IF(config_flags%specified .or. config_flags%nested) THEN
-     IF (degrade_xs) i_start = MAX(its,ids+1)
-     IF (degrade_xe) i_end   = MIN(ite,ide-2)
-     IF (degrade_ys) j_start = MAX(jts,jds+1)
-     IF (degrade_ye) j_end   = MIN(jte,jde-2)
+     IF (degrade_xs) i_start = MAX(its-1,ids+1)
+     IF (degrade_xe) i_end   = MIN(ite+1,ide-2)
+     IF (degrade_ys) j_start = MAX(jts-1,jds+1)
+     IF (degrade_ye) j_end   = MIN(jte+1,jde-2)
    END IF
 
    IF(config_flags%open_xs) THEN
-     IF (degrade_xs) i_start = MAX(its,ids+1)
+     IF (degrade_xs) i_start = MAX(its-1,ids+1)
    END IF
    IF(config_flags%open_xe) THEN
-     IF (degrade_xe) i_end   = MIN(ite,ide-2)
+     IF (degrade_xe) i_end   = MIN(ite+1,ide-2)
    END IF
    IF(config_flags%open_ys) THEN
-     IF (degrade_ys) j_start = MAX(jts,jds+1)
+     IF (degrade_ys) j_start = MAX(jts-1,jds+1)
    END IF
    IF(config_flags%open_ye) THEN
-     IF (degrade_ye) j_end   = MIN(jte,jde-2)
+     IF (degrade_ye) j_end   = MIN(jte+1,jde-2)
    END IF
 
@@ -6962,14 +7128,16 @@
    DO i=i_start, i_end
 
-     ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j)                   &
-                - dt*( msft(i,j)*( rdx*(fqxl(i+1,k,j)-fqxl(i,k,j))      &
-                                  +rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) )    &
-                              +rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j))   )  
-
-     flux_out = dt*(msft(i,j)*( rdx*(  max(0.,fqx (i+1,k,j))      &
+     ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j)        &
+                - dt*( msft(i,j)*(                           & 
+                       rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) +     &
+                       rdy*(fqyl(i,k,j+1)-fqyl(i,k,j))  )    &
+                      +         rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) )
+
+     flux_out = dt*(  msft(i,j)             *(                    &
+                                rdx*(  max(0.,fqx (i+1,k,j))      &
                                       -min(0.,fqx (i  ,k,j)) )    &
                                +rdy*(  max(0.,fqy (i,k,j+1))      &
                                       -min(0.,fqy (i,k,j  )) ) )  &
-                           +rdzw(k)*(  min(0.,fqz (i,k+1,j))      &
+                +           rdzw(k)*(  min(0.,fqz (i,k+1,j))      &
                                       -max(0.,fqz (i,k  ,j)) )   )
 
@@ -7015,6 +7183,6 @@
 ! x flux divergence
 !
-  IF(degrade_xs) i_start = i_start + 1
-  IF(degrade_xe) i_end   = i_end   - 1
+  IF(degrade_xs) i_start = MAX(its,ids+1)
+  IF(degrade_xe) i_end   = MIN(ite,ide-2)
 
   DO j = j_start, j_end
@@ -7023,5 +7191,5 @@
 
      tendency (i,k,j) = tendency(i,k,j)                           &
-                - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j)     &
+               - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j)     &
                                    +fqxl(i+1,k,j)-fqxl(i,k,j))   )
 
@@ -7034,6 +7202,6 @@
   i_start = its
   i_end   = MIN(ite,ide-1)
-  IF(degrade_ys) j_start = j_start + 1
-  IF(degrade_ye) j_end   = j_end   - 1
+  IF(degrade_ys) j_start = MAX(jts,jds+1)
+  IF(degrade_ye) j_end   = MIN(jte,jde-2)
 
   DO j = j_start, j_end
@@ -7042,5 +7210,5 @@
 
      tendency (i,k,j) = tendency(i,k,j)                           &
-                - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j)     &
+               - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j)     &
                                    +fqyl(i,k,j+1)-fqyl(i,k,j))   )
 
@@ -7053,4 +7221,904 @@
 !----------------------------------------------------------------
 
+SUBROUTINE advect_scalar_mono   ( field, field_old, tendency,    &
+                                  ru, rv, rom,                   &
+                                  mut, mub, mu_old,              &
+                                  config_flags,                  &
+                                  msfu, msfv,                    &
+                                  msft,                          &
+                                  fzm, fzp,                      &
+                                  rdx, rdy, rdzw, dt,            &
+                                  ids, ide, jds, jde, kds, kde,  &
+                                  ims, ime, jms, jme, kms, kme,  &
+                                  its, ite, jts, jte, kts, kte  )
+
+!  monotonic advection option
+!  for scalars in WRF RK3 advection.  This version is memory intensive ->
+!  we save 3d arrays of x, y and z both high and low order fluxes
+!  (six in all).  Alternatively, we could sweep in a direction 
+!  and lower the cost considerably.
+
+!  uses the Smolarkiewicz MWR 1989 approach, with addition of first-order
+!  fluxes initially
+
+   IMPLICIT NONE
+   
+   ! Input data
+   
+   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
+
+   INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
+                                              ims, ime, jms, jme, kms, kme, &
+                                              its, ite, jts, jte, kts, kte
+
+   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field,     &
+                                                                      field_old, &
+                                                                      ru,    &
+                                                                      rv,    &
+                                                                      rom
+
+   REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mut, mub, mu_old
+   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
+
+   REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfu,  &
+                                                                    msfv,  &
+                                                                    msft
+
+   REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
+                                                                  fzp,  &
+                                                                  rdzw
+
+   REAL ,                                        INTENT(IN   ) :: rdx,  &
+                                                                  rdy,  &
+                                                                  dt
+
+   ! Local data
+   
+   INTEGER :: i, j, k, itf, jtf, ktf
+   INTEGER :: i_start, i_end, j_start, j_end
+   INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
+   INTEGER :: jmin, jmax, jp, jm, imin, imax
+
+   REAL    :: mrdx, mrdy, ub, vb, uw, vw, mu
+   REAL , DIMENSION(its:ite, kts:kte) :: vflux
+
+
+!  storage for high and low order fluxes
+
+   REAL,  DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2  ) :: fqx, fqy, fqz
+   REAL,  DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2  ) :: fqxl, fqyl, fqzl
+   REAL,  DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2  ) :: qmin, qmax
+   REAL,  DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2  ) :: scale_in, scale_out
+   REAL :: ph_upwind
+
+   INTEGER :: horz_order, vert_order
+   
+   LOGICAL :: degrade_xs, degrade_ys
+   LOGICAL :: degrade_xe, degrade_ye
+
+   INTEGER :: jp1, jp0, jtmp
+
+   REAL :: flux_out, ph_low, flux_in, ph_hi, scale
+   REAL, PARAMETER :: eps=1.e-20
+
+
+! definition of flux operators, 3rd, 4rth, 5th or 6th order
+
+   REAL    :: flux3, flux4, flux5, flux6, flux_upwind
+   REAL    :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr
+
+      flux4(q_im2, q_im1, q_i, q_ip1, ua) =                     &
+            (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2)
+
+      flux3(q_im2, q_im1, q_i, q_ip1, ua) =                     &
+           flux4(q_im2, q_im1, q_i, q_ip1, ua) +                &
+           sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1))
+
+      flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
+            (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2)      &
+            +(1./60.)*(q_ip2+q_im3)
+
+      flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) =       &
+           flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)    &
+            -sign(1.,ua)*(1./60.)*(                             &
+              (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )
+
+!      flux_upwind(q_im1, q_i, cr ) = 0.
+      flux_upwind(q_im1, q_i, cr ) = 0.5*(1.+sign(1.,cr))*q_im1 &
+                                    +0.5*(1.-sign(1.,cr))*q_i
+
+    LOGICAL, PARAMETER :: mono_limit = .true.
+
+! set order for the advection schemes
+
+  ktf=MIN(kte,kde-1)
+  horz_order = config_flags%h_sca_adv_order
+  vert_order = config_flags%v_sca_adv_order
+
+  do j=jts-2,jte+2
+  do k=kts,kte
+  do i=its-2,ite+2
+    qmin(i,k,j) = field_old(i,k,j)
+    qmax(i,k,j) = field_old(i,k,j)
+    scale_in(i,k,j) = 1.
+    scale_out(i,k,j) = 1.
+    fqx(i,k,j) = 0.
+    fqy(i,k,j) = 0.
+    fqz(i,k,j) = 0.
+    fqxl(i,k,j) = 0.
+    fqyl(i,k,j) = 0.
+    fqzl(i,k,j) = 0.
+  enddo
+  enddo
+  enddo
+
+!  begin with horizontal flux divergence
+!  here is the choice of flux operators
+
+
+  horizontal_order_test : IF( horz_order == 5 ) THEN
+
+!  determine boundary mods for flux operators
+!  We degrade the flux operators from 3rd/4rth order
+!   to second order one gridpoint in from the boundaries for
+!   all boundary conditions except periodic and symmetry - these
+!   conditions have boundary zone data fill for correct application
+!   of the higher order flux stencils
+
+   degrade_xs = .true.
+   degrade_xe = .true.
+   degrade_ys = .true.
+   degrade_ye = .true.
+
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xs .or. &
+       (its > ids+3)                ) degrade_xs = .false.
+   IF( config_flags%periodic_x   .or. &
+       config_flags%symmetric_xe .or. &
+       (ite < ide-4)                ) degrade_xe = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ys .or. &
+       (jts > jds+3)                ) degrade_ys = .false.
+   IF( config_flags%periodic_y   .or. &
+       config_flags%symmetric_ye .or. &
+       (jte < jde-4)                ) degrade_ye = .false.
+
+!--------------- y - advection first
+
+!--  y flux compute; these bounds are for periodic and sym b.c.
+
+      ktf=MIN(kte,kde-1)
+      i_start = its-1
+      i_end   = MIN(ite,ide-1)+1
+      j_start = jts-1
+      j_end   = MIN(jte,jde-1)+1
+      j_start_f = j_start
+      j_end_f   = j_end+1
+
+!--  modify loop bounds if open or specified
+
+!  WCS 20090218
+!      IF(degrade_xs) i_start = its
+!      IF(degrade_xe) i_end   = MIN(ite,ide-1)
+      IF(degrade_xs) i_start = MAX(its-1,ids)
+      IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
+
+!  WCS 20090218
+!      IF(degrade_ys) then
+!        j_start = MAX(jts,jds+1)
+!        j_start_f = jds+3
+!      ENDIF
+!
+!      IF(degrade_ye) then
+!        j_end = MIN(jte,jde-2)
+!        j_end_f = jde-3
+!      ENDIF
+
+      IF(degrade_ys) then
+        j_start = MAX(jts-1,jds+1)
+        j_start_f = jds+3
+      ENDIF
+
+      IF(degrade_ye) then
+        j_end = MIN(jte+1,jde-2)
+        j_end_f = jde-3
+      ENDIF
+
+!  compute fluxes, 5th order
+
+      j_loop_y_flux_5 : DO j = j_start, j_end+1
+
+      IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil
+
+        DO k=kts,ktf
+        DO i = i_start, i_end
+
+          vel = rv(i,k,j)
+          cr = vel
+          fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j  ), vel)
+
+          fqy( i, k, j  ) = vel*flux5(                                  &
+                  field(i,k,j-3), field(i,k,j-2), field(i,k,j-1),       &
+                  field(i,k,j  ), field(i,k,j+1), field(i,k,j+2),  vel )
+
+          fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k,j-1))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k,j-1))
+          else
+             qmax(i,k,j-1)  = amax1(qmax(i,k,j-1),field_old(i,k,j))
+             qmin(i,k,j-1)  = amin1(qmin(i,k,j-1),field_old(i,k,j))
+          end if
+
+        ENDDO
+        ENDDO
+
+      ELSE IF ( j == jds+1 ) THEN   ! 2nd order flux next to south boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
+
+              vel = rv(i,k,j)
+              cr = vel
+              fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j  ), cr)
+
+              fqy(i,k, j) = 0.5*rv(i,k,j)*          &
+                     (field(i,k,j)+field(i,k,j-1))
+
+              fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k,j-1))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k,j-1))
+          else
+             qmax(i,k,j-1)  = amax1(qmax(i,k,j-1),field_old(i,k,j))
+             qmin(i,k,j-1)  = amin1(qmin(i,k,j-1),field_old(i,k,j))
+          end if
+
+            ENDDO
+            ENDDO
+
+      ELSE IF  ( j == jds+2 ) THEN  ! third of 4rth order flux 2 in from south boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
+
+              vel = rv(i,k,j)
+              cr = vel
+              fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j  ), cr)
+
+              fqy( i, k, j ) = vel*flux3(              &
+                   field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel )
+              fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k,j-1))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k,j-1))
+          else
+             qmax(i,k,j-1)  = amax1(qmax(i,k,j-1),field_old(i,k,j))
+             qmin(i,k,j-1)  = amin1(qmin(i,k,j-1),field_old(i,k,j))
+          end if
+
+            ENDDO
+            ENDDO
+
+      ELSE IF ( j == jde-1 ) THEN  ! 2nd order flux next to north boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
+
+              vel = rv(i,k,j)
+              cr = vel
+              fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j  ), cr)
+
+              fqy(i, k, j ) = 0.5*rv(i,k,j)*      &
+                     (field(i,k,j)+field(i,k,j-1))
+              fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k,j-1))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k,j-1))
+          else
+             qmax(i,k,j-1)  = amax1(qmax(i,k,j-1),field_old(i,k,j))
+             qmin(i,k,j-1)  = amin1(qmin(i,k,j-1),field_old(i,k,j))
+          end if
+
+            ENDDO
+            ENDDO
+
+      ELSE IF ( j == jde-2 ) THEN  ! 3rd or 4rth order flux 2 in from north boundary
+
+            DO k=kts,ktf
+            DO i = i_start, i_end
+
+              vel = rv(i,k,j)
+              cr = vel
+              fqyl(i,k,j) = vel*flux_upwind(field_old(i,k,j-1), field_old(i,k,j  ), cr)
+
+              fqy( i, k, j) = vel*flux3(             &
+                   field(i,k,j-2),field(i,k,j-1),    &
+                   field(i,k,j),field(i,k,j+1),vel )
+              fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k,j-1))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k,j-1))
+          else
+             qmax(i,k,j-1)  = amax1(qmax(i,k,j-1),field_old(i,k,j))
+             qmin(i,k,j-1)  = amin1(qmin(i,k,j-1),field_old(i,k,j))
+          end if
+
+            ENDDO
+            ENDDO
+
+      ENDIF
+
+   ENDDO j_loop_y_flux_5
+
+!  next, x flux
+
+!--  these bounds are for periodic and sym conditions
+
+      i_start = its-1
+      i_end   = MIN(ite,ide-1)+1
+      i_start_f = i_start
+      i_end_f   = i_end+1
+
+      j_start = jts-1
+      j_end   = MIN(jte,jde-1)+1
+
+!--  modify loop bounds for open and specified b.c
+
+!  WCS 20090218
+!      IF(degrade_ys) j_start = jts
+!      IF(degrade_ye) j_end   = MIN(jte,jde-1)
+      IF(degrade_ys) j_start = MAX(jts-1,jds)
+      IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
+
+!  WCS 20090218
+!      IF(degrade_xs) then
+!        i_start = MAX(ids+1,its)
+!        i_start_f = i_start+2
+!      ENDIF
+
+!      IF(degrade_xe) then
+!        i_end = MIN(ide-2,ite)
+!        i_end_f = ide-3
+!      ENDIF
+
+      IF(degrade_xs) then
+        i_start = MAX(ids+1,its-1)
+        i_start_f = ids+3
+      ENDIF
+
+      IF(degrade_xe) then
+        i_end = MIN(ide-2,ite+1)
+        i_end_f = ide-3
+      ENDIF
+
+!  compute fluxes
+
+      DO j = j_start, j_end
+
+!  5th or 6th order flux
+
+        DO k=kts,ktf
+        DO i = i_start_f, i_end_f
+
+          vel = ru(i,k,j)
+          cr = vel
+          fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+
+          fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j),  &
+                                         field(i-1,k,j), field(i  ,k,j),  &
+                                         field(i+1,k,j), field(i+2,k,j),  &
+                                         vel                             )
+          fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i-1,k,j))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i-1,k,j))
+          else
+             qmax(i-1,k,j)  = amax1(qmax(i-1,k,j),field_old(i,k,j))
+             qmin(i-1,k,j)  = amin1(qmin(i-1,k,j),field_old(i,k,j))
+          end if
+
+        ENDDO
+        ENDDO
+
+!  lower order fluxes close to boundaries (if not periodic or symmetric)
+
+!  WCS 20090218 degrade_xs and xe recoded
+
+        IF( degrade_xs ) THEN
+
+          DO i=i_start,i_start_f-1
+
+            IF(i == ids+1) THEN ! second order
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                cr = vel
+                fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+
+                fqx(i,k,j) = 0.5*(ru(i,k,j)) &
+                       *(field(i,k,j)+field(i-1,k,j))
+
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+
+                if(cr.gt. 0) then
+                  qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i-1,k,j))
+                  qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i-1,k,j))
+                else
+                  qmax(i-1,k,j)  = amax1(qmax(i-1,k,j),field_old(i,k,j))
+                  qmin(i-1,k,j)  = amin1(qmin(i-1,k,j),field_old(i,k,j))
+                end if
+              ENDDO
+            ENDIF
+
+            IF(i == ids+2) THEN  ! third order
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                cr = vel
+                fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+                                          field(i  ,k,j), field(i+1,k,j),  &
+                                          vel                             )
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+
+                if(cr.gt. 0) then
+                  qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i-1,k,j))
+                  qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i-1,k,j))
+                else
+                  qmax(i-1,k,j)  = amax1(qmax(i-1,k,j),field_old(i,k,j))
+                  qmin(i-1,k,j)  = amin1(qmin(i-1,k,j),field_old(i,k,j))
+                end if
+              ENDDO
+            ENDIF
+
+          ENDDO
+
+        ENDIF
+
+        IF( degrade_xe ) THEN
+
+          DO i = i_end_f+1, i_end+1
+
+            IF( i == ide-1 ) THEN ! second order flux next to the boundary
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                cr = vel
+                fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx(i,k,j) = 0.5*(ru(i,k,j))      &
+                       *(field(i,k,j)+field(i-1,k,j))
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+
+                if(cr.gt. 0) then
+                  qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i-1,k,j))
+                  qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i-1,k,j))
+                else
+                  qmax(i-1,k,j)  = amax1(qmax(i-1,k,j),field_old(i,k,j))
+                  qmin(i-1,k,j)  = amin1(qmin(i-1,k,j),field_old(i,k,j))
+                end if
+              ENDDO
+            ENDIF
+
+            IF( i == ide-2 ) THEN ! third order flux one in from the boundary
+              DO k=kts,ktf
+                vel = ru(i,k,j)
+                cr = vel
+                fqxl(i,k,j) = vel*flux_upwind(field_old(i-1,k,j), field_old(i,k,j  ), cr)
+                fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j),  &
+                                          field(i  ,k,j), field(i+1,k,j),  &
+                                          vel                             )
+                fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j)
+
+                if(cr.gt. 0) then
+                  qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i-1,k,j))
+                  qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i-1,k,j))
+                else
+                  qmax(i-1,k,j)  = amax1(qmax(i-1,k,j),field_old(i,k,j))
+                  qmin(i-1,k,j)  = amin1(qmin(i-1,k,j),field_old(i,k,j))
+                end if
+              ENDDO
+            ENDIF
+          ENDDO
+        ENDIF
+
+      ENDDO  ! enddo for outer J loop
+
+   ELSE
+
+      WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_mono, h_order not known ',horz_order
+      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
+
+   ENDIF horizontal_order_test
+
+!  pick up the rest of the horizontal radiation boundary conditions.
+!  (these are the computations that don't require 'cb'.
+!  first, set to index ranges
+
+      i_start = its
+      i_end   = MIN(ite,ide-1)
+      j_start = jts
+      j_end   = MIN(jte,jde-1)
+
+!  compute x (u) conditions for v, w, or scalar
+
+   IF( (config_flags%open_xs) .and. (its == ids) ) THEN
+
+       DO j = j_start, j_end
+       DO k = kts, ktf
+         ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. )
+         tendency(its,k,j) = tendency(its,k,j)                     &
+               - rdx*(                                             &
+                       ub*(   field_old(its+1,k,j)                 &
+                            - field_old(its  ,k,j)   ) +           &
+                       field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j))  &
+                                                                )
+       ENDDO
+       ENDDO
+
+   ENDIF
+
+   IF( (config_flags%open_xe) .and. (ite == ide) ) THEN
+
+       DO j = j_start, j_end
+       DO k = kts, ktf
+         ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. )
+         tendency(i_end,k,j) = tendency(i_end,k,j)                   &
+               - rdx*(                                               &
+                       ub*(  field_old(i_end  ,k,j)                  &
+                           - field_old(i_end-1,k,j) ) +              &
+                       field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j))  &
+                                                                    )
+       ENDDO
+       ENDDO
+
+   ENDIF
+
+   IF( (config_flags%open_ys) .and. (jts == jds) ) THEN
+
+       DO i = i_start, i_end
+       DO k = kts, ktf
+         vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. )
+         tendency(i,k,jts) = tendency(i,k,jts)                     &
+               - rdy*(                                             &
+                       vb*(  field_old(i,k,jts+1)                  &
+                           - field_old(i,k,jts  ) ) +              &
+                       field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts))  &
+                                                                )
+       ENDDO
+       ENDDO
+
+   ENDIF
+
+   IF( (config_flags%open_ye) .and. (jte == jde)) THEN
+
+       DO i = i_start, i_end
+       DO k = kts, ktf
+         vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. )
+         tendency(i,k,j_end) = tendency(i,k,j_end)                   &
+               - rdy*(                                               &
+                       vb*(   field_old(i,k,j_end  )                 &
+                            - field_old(i,k,j_end-1) ) +             &
+                       field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1))  &
+                                                                    )
+       ENDDO
+       ENDDO
+
+   ENDIF
+
+!-------------------- vertical advection
+
+!-- loop bounds for periodic or sym conditions
+
+      i_start = its-1
+      i_end   = MIN(ite,ide-1)+1
+      j_start = jts-1
+      j_end   = MIN(jte,jde-1)+1
+
+!-- loop bounds for open or specified conditions
+
+!  WCS 20090218
+!    IF(degrade_xs) i_start = its
+!    IF(degrade_xe) i_end   = MIN(ite,ide-1)
+!    IF(degrade_ys) j_start = jts
+!    IF(degrade_ye) j_end   = MIN(jte,jde-1)
+
+    IF(degrade_xs) i_start = MAX(its-1,ids)
+    IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
+    IF(degrade_ys) j_start = MAX(jts-1,jds)
+    IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
+
+
+    vert_order_test : IF (vert_order == 3) THEN    
+
+      DO j = j_start, j_end
+
+         DO i = i_start, i_end
+           fqz(i,1,j)  = 0.
+           fqzl(i,1,j) = 0.
+           fqz(i,kde,j)  = 0.
+           fqzl(i,kde,j) = 0.
+         ENDDO
+
+         DO k=kts+2,ktf-1
+         DO i = i_start, i_end
+
+           vel = rom(i,k,j)
+           cr = -vel
+           fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j  ), cr)
+
+           fqz(i,k,j) = vel*flux3(                      &
+                   field(i,k-2,j), field(i,k-1,j),      &
+                   field(i,k  ,j), field(i,k+1,j),  -vel )
+           fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k-1,j))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k-1,j))
+          else
+             qmax(i,k-1,j)  = amax1(qmax(i,k-1,j),field_old(i,k,j))
+             qmin(i,k-1,j)  = amin1(qmin(i,k-1,j),field_old(i,k,j))
+          end if
+
+         ENDDO
+         ENDDO
+
+         DO i = i_start, i_end
+
+           k=kts+1
+           vel = rom(i,k,j)
+           cr = -vel
+           fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j  ), cr)
+           fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
+           fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k-1,j))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k-1,j))
+          else
+             qmax(i,k-1,j)  = amax1(qmax(i,k-1,j),field_old(i,k,j))
+             qmin(i,k-1,j)  = amin1(qmin(i,k-1,j),field_old(i,k,j))
+          end if
+
+           k=ktf
+           vel = rom(i,k,j)
+           cr = -vel
+           fqzl(i,k,j) = vel*flux_upwind(field_old(i,k-1,j), field_old(i,k,j  ), cr)
+           fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j))
+           fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j)
+
+          if(cr.gt. 0) then
+             qmax(i,k,j)  = amax1(qmax(i,k,j),field_old(i,k-1,j))
+             qmin(i,k,j)  = amin1(qmin(i,k,j),field_old(i,k-1,j))
+          else
+             qmax(i,k-1,j)  = amax1(qmax(i,k-1,j),field_old(i,k,j))
+             qmin(i,k-1,j)  = amin1(qmin(i,k-1,j),field_old(i,k,j))
+          end if
+         ENDDO
+
+      ENDDO
+
+   ELSE
+
+      WRITE (wrf_err_message,*) ' advect_scalar_mono, v_order not known ',vert_order
+      CALL wrf_error_fatal ( wrf_err_message )
+
+   ENDIF vert_order_test
+
+   IF (mono_limit) THEN
+
+! montonic filter
+
+   i_start = its-1
+   i_end   = MIN(ite,ide-1)+1
+   j_start = jts-1
+   j_end   = MIN(jte,jde-1)+1
+
+! WCS 20090218
+
+!-- loop bounds for open or specified conditions
+!
+!   IF(degrade_xs) i_start = its
+!   IF(degrade_xe) i_end   = MIN(ite,ide-1)
+!   IF(degrade_ys) j_start = jts
+!   IF(degrade_ye) j_end   = MIN(jte,jde-1)
+!
+!   IF(config_flags%specified .or. config_flags%nested) THEN
+!     IF (degrade_xs) i_start = MAX(its,ids+1)
+!     IF (degrade_xe) i_end   = MIN(ite,ide-2)
+!     IF (degrade_ys) j_start = MAX(jts,jds+1)
+!     IF (degrade_ye) j_end   = MIN(jte,jde-2)
+!   END IF
+!
+!   IF(config_flags%open_xs) THEN
+!     IF (degrade_xs) i_start = MAX(its,ids+1)
+!   END IF
+!   IF(config_flags%open_xe) THEN
+!     IF (degrade_xe) i_end   = MIN(ite,ide-2)
+!   END IF
+!   IF(config_flags%open_ys) THEN
+!     IF (degrade_ys) j_start = MAX(jts,jds+1)
+!   END IF
+!   IF(config_flags%open_ye) THEN
+!     IF (degrade_ye) j_end   = MIN(jte,jde-2)
+!   END IF
+
+   IF(degrade_xs) i_start = MAX(its-1,ids)
+   IF(degrade_xe) i_end   = MIN(ite+1,ide-1)
+   IF(degrade_ys) j_start = MAX(jts-1,jds)
+   IF(degrade_ye) j_end   = MIN(jte+1,jde-1)
+
+   IF(config_flags%specified .or. config_flags%nested) THEN
+     IF (degrade_xs) i_start = MAX(its-1,ids+1)
+     IF (degrade_xe) i_end   = MIN(ite+1,ide-2)
+     IF (degrade_ys) j_start = MAX(jts-1,jds+1)
+     IF (degrade_ye) j_end   = MIN(jte+1,jde-2)
+   END IF
+
+   IF(config_flags%open_xs) THEN
+     IF (degrade_xs) i_start = MAX(its-1,ids+1)
+   END IF
+   IF(config_flags%open_xe) THEN
+     IF (degrade_xe) i_end   = MIN(ite+1,ide-2)
+   END IF
+   IF(config_flags%open_ys) THEN
+     IF (degrade_ys) j_start = MAX(jts-1,jds+1)
+   END IF
+   IF(config_flags%open_ye) THEN
+     IF (degrade_ye) j_end   = MIN(jte+1,jde-2)
+   END IF
+
+!-- here is the limiter...
+
+   DO j=j_start, j_end
+   DO k=kts, ktf
+   DO i=i_start, i_end
+
+     ph_upwind = (mub(i,j)+mu_old(i,j))*field_old(i,k,j)        &
+                   - dt*( msft(i,j)            *(               &
+                          rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) +     &
+                          rdy*(fqyl(i,k,j+1)-fqyl(i,k,j))  )    &
+                         +           rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) )
+
+     flux_in = -dt*(  msft(i,j)             *(                   &
+                               rdx*(  min(0.,fqx (i+1,k,j))      &
+                                     -max(0.,fqx (i  ,k,j)) )    &
+                              +rdy*(  min(0.,fqy (i,k,j+1))      &
+                                     -max(0.,fqy (i,k,j  )) ) )  &
+               +           rdzw(k)*(  max(0.,fqz (i,k+1,j))      &
+                                     -min(0.,fqz (i,k  ,j)) )   )
+
+     ph_hi = mut(i,j)*qmax(i,k,j) - ph_upwind
+     IF( flux_in .gt. ph_hi ) scale_in(i,k,j) = max(0.,ph_hi/(flux_in+eps))
+
+
+     flux_out = dt*( msft(i,j)*(                    &
+                                rdx*(  max(0.,fqx (i+1,k,j))      &
+                                      -min(0.,fqx (i  ,k,j)) )    &
+                               +rdy*(  max(0.,fqy (i,k,j+1))      &
+                                      -min(0.,fqy (i,k,j  )) ) )  &
+                +           rdzw(k)*(  min(0.,fqz (i,k+1,j))      &
+                                      -max(0.,fqz (i,k  ,j)) )   )
+
+     ph_low = ph_upwind - mut(i,j)*qmin(i,k,j)
+     IF( flux_out .gt. ph_low ) scale_out(i,k,j) = max(0.,ph_low/(flux_out+eps))
+
+   ENDDO
+   ENDDO
+   ENDDO
+
+   DO j=j_start, j_end
+   DO k=kts, ktf
+   DO i=i_start, i_end+1
+       IF( fqx (i,k,j) .gt. 0.) then
+         fqx(i,k,j) = min(scale_in(i,k,j),scale_out(i-1,k,j))*fqx(i,k,j)
+       ELSE
+         fqx(i,k,j) = min(scale_out(i,k,j),scale_in(i-1,k,j))*fqx(i,k,j)
+       ENDIF
+   ENDDO
+   ENDDO
+   ENDDO
+
+   DO j=j_start, j_end+1
+   DO k=kts, ktf
+   DO i=i_start, i_end
+       IF( fqy (i,k,j) .gt. 0.) then
+         fqy(i,k,j) = min(scale_in(i,k,j),scale_out(i,k,j-1))*fqy(i,k,j)
+       ELSE
+         fqy(i,k,j) = min(scale_out(i,k,j),scale_in(i,k,j-1))*fqy(i,k,j)
+       ENDIF
+   ENDDO
+   ENDDO
+   ENDDO
+
+   DO j=j_start, j_end
+   DO k=kts+1, ktf
+   DO i=i_start, i_end
+       IF( fqz (i,k,j) .lt. 0.) then
+         fqz(i,k,j) = min(scale_in(i,k,j),scale_out(i,k-1,j))*fqz(i,k,j)
+       ELSE
+         fqz(i,k,j) = min(scale_out(i,k,j),scale_in(i,k-1,j))*fqz(i,k,j)
+       ENDIF
+   ENDDO
+   ENDDO
+   ENDDO
+
+   END IF
+
+! add in the mono-limited flux divergence
+! we need to fix this for open b.c set ***********
+
+  i_start = its
+  i_end   = MIN(ite,ide-1)
+  j_start = jts
+  j_end   = MIN(jte,jde-1)
+
+  DO j = j_start, j_end
+  DO k = kts, ktf
+  DO i = i_start, i_end
+
+     tendency (i,k,j) = tendency(i,k,j)                           &
+                            -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j)  &
+                                      +fqzl(i,k+1,j)-fqzl(i,k,j))
+
+  ENDDO
+  ENDDO
+  ENDDO
+
+! x flux divergence
+!
+
+! WCS 20090218
+!  IF(degrade_xs) i_start = i_start + 1
+!  IF(degrade_xe) i_end   = i_end   - 1
+
+  IF(degrade_xs) i_start = MAX(its,ids+1)
+  IF(degrade_xe) i_end   = MIN(ite,ide-2)
+
+  DO j = j_start, j_end
+  DO k = kts, ktf
+  DO i = i_start, i_end
+
+     ! Un-"canceled" map scale factor, ADT Eq. 48
+     tendency (i,k,j) = tendency(i,k,j)                           &
+               - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j)     &
+                                   +fqxl(i+1,k,j)-fqxl(i,k,j))   )
+
+  ENDDO
+  ENDDO
+  ENDDO
+
+! y flux divergence
+!
+  i_start = its
+  i_end   = MIN(ite,ide-1)
+
+! WCS 20090218
+!  IF(degrade_ys) j_start = j_start + 1
+!  IF(degrade_ye) j_end   = j_end   - 1
+
+  IF(degrade_ys) j_start = MAX(jts,jds+1)
+  IF(degrade_ye) j_end   = MIN(jte,jde-2)
+
+  DO j = j_start, j_end
+  DO k = kts, ktf
+  DO i = i_start, i_end
+
+     ! Un-"canceled" map scale factor, ADT Eq. 48
+     tendency (i,k,j) = tendency(i,k,j)                           &
+               - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j)     &
+                                   +fqyl(i,k,j+1)-fqyl(i,k,j))   )
+
+  ENDDO
+  ENDDO
+  ENDDO
+
+END SUBROUTINE advect_scalar_mono
+
+!-----------------------------------------------------------
+
 END MODULE module_advect_em
 
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_em.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_em.F	(revision 195)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_em.F	(revision 196)
@@ -16,5 +16,6 @@
                            u, v, w, t, ph, mu,              &
                            moist,                           &
-                           ru, rv, rw, ww, php, alt, muu, muv,  &
+                           ru, rv, rw, ww, php, alt,        &
+                           muu, muv,                        &
                            mub, mut, phb, pb, p, al, alb,   &
                            cqu, cqv, cqw,                   &
@@ -280,4 +281,5 @@
    REAL    :: kdift, khdq, kvdq, cfn, cfn1, cf1, cf2, cf3
    INTEGER :: i,j,k
+   INTEGER :: time_step
 
 !<DESCRIPTION>
@@ -362,7 +364,8 @@
 
      !  advection tendencies
+     CALL nl_get_time_step ( 1, time_step )
 
      CALL advect_u ( u, u , ru_tend, ru, rv, ww,   &
-                     mut, config_flags,            &
+                     mut, time_step, config_flags, &
                      msfu, msfv, msft,             &
                      fnm, fnp, rdx, rdy, rdnw,     &
@@ -372,5 +375,5 @@
 
      CALL advect_v ( v, v , rv_tend, ru, rv, ww,   &
-                     mut, config_flags,            &
+                     mut, time_step, config_flags, &
                      msfu, msfv, msft,             &
                      fnm, fnp, rdx, rdy, rdnw,     &
@@ -381,5 +384,5 @@
      IF (non_hydrostatic)                          &
      CALL advect_w ( w, w, rw_tend, ru, rv, ww,    &
-                     mut, config_flags,            &
+                     mut, time_step, config_flags, &
                      msfu, msfv, msft,             &
                      fnm, fnp, rdx, rdy, rdn,      &
@@ -391,5 +394,5 @@
 
      CALL advect_scalar ( t, t, t_tend, ru, rv, ww,     &
-                          mut, config_flags,            &
+                          mut, time_step, config_flags, &
                           msfu, msfv, msft, fnm, fnp,   &
                           rdx, rdy, rdnw,               &
@@ -418,5 +421,4 @@
                   ims, ime, jms, jme, kms, kme,      &
                   its, ite, jts, jte, kts, kte      )
-
 
   CALL horizontal_pressure_gradient( ru_tend,rv_tend,                &
@@ -454,5 +456,4 @@
                                  its, ite, jts, jte, kts, kte )
   ELSE
-
     CALL coriolis ( ru, rv, rw,                   &
                     ru_tend,  rv_tend,  rw_tend,  &
@@ -463,6 +464,5 @@
                     its, ite, jts, jte, kts, kte )
 
-   END IF
-
+  END IF
 
   CALL curvature ( ru, rv, rw, u, v, w,            &
@@ -515,4 +515,5 @@
                                          its, ite, jts, jte, kts, kte )
 
+
 !!!****MARS: vertical diffusion is done in the physics (TODO: consider the nonhydrostatic case ?)
         pbl_test : IF ( (config_flags%bl_pbl_physics .eq. 0) &
@@ -548,5 +549,4 @@
 
         ENDIF pbl_test
-
 
    !  Theta tendency computations.
@@ -654,5 +654,4 @@
    INTEGER :: i, j, k
 
-
 !<DESCRIPTION>
 !
@@ -679,5 +678,7 @@
    DO k = kts,kte-1
    DO i = its,ite
+     ! multiply by m to uncouple u
      IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) +  u_save(i,k,j)*msfu(i,j)
+     ! divide by m to couple u
      ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfu(i,j)
    ENDDO
@@ -688,5 +689,7 @@
    DO k = kts,kte-1
    DO i = its,MIN(ite,ide-1)
+     ! multiply by m to uncouple v
      IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) +  v_save(i,k,j)*msfv(i,j)
+     ! divide by m to couple v
      rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)/msfv(i,j)
    ENDDO
@@ -697,7 +700,10 @@
    DO k = kts,kte
    DO i = its,MIN(ite,ide-1)
+     ! multiply by m to uncouple w
      IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) +  w_save(i,k,j)*msft(i,j)
+     ! divide by m to couple w
      rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msft(i,j)
      IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) +  ph_save(i,k,j)
+     ! divide by m to couple scalar
      ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msft(i,j)
    ENDDO
@@ -709,6 +715,8 @@
    DO i = its,MIN(ite,ide-1)
      IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) +  t_save(i,k,j)
-     t_tend(i,k,j) =  t_tend(i,k,j) +  t_tendf(i,k,j)/msft(i,j)  &
+     ! divide by m to couple theta
+      t_tend(i,k,j) =  t_tend(i,k,j) +  t_tendf(i,k,j)/msft(i,j)  &
                                      +  mut(i,j)*h_diabatic(i,k,j)/msft(i,j)
+     ! divide by m to couple heating
    ENDDO
    ENDDO
@@ -757,10 +765,10 @@
 
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                &
-                                         INTENT(INOUT)  :: scalar, scalar_old
+                                         INTENT(IN   )  :: scalar, scalar_old
 
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                      &
                                          INTENT(INOUT)  :: scalar_tends
                                                     
-   REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ) :: advect_tend
+   REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(INOUT) :: advect_tend
 
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(OUT  ) :: RQVFTEN
@@ -801,4 +809,5 @@
   
    INTEGER :: im, i,j,k
+   INTEGER :: time_step
 
    REAL    :: khdq, kvdq, tendency
@@ -822,28 +831,43 @@
                       its, ite, jts, jte, kts, kte )
 
+     CALL nl_get_time_step ( 1, time_step )
+
       IF( (rk_step == 3) .and. pd_advection ) THEN
 
-       CALL advect_scalar_pd       ( scalar(ims,kms,jms,im),             &
-                                     scalar_old(ims,kms,jms,im),         &
-                                     advect_tend(ims,kms,jms),           &
-                                     ru, rv, ww, mut, mub, mu_old,       &
-                                     config_flags,                       &
-                                     msfu, msfv, msft, fnm, fnp,         &
-                                     rdx, rdy, rdnw,dt,                  &
-                                     ids, ide, jds, jde, kds, kde,       &
-                                     ims, ime, jms, jme, kms, kme,       &
-                                     its, ite, jts, jte, kts, kte     )
-
+        CALL advect_scalar_pd       ( scalar(ims,kms,jms,im),             &
+                                      scalar_old(ims,kms,jms,im),         &
+                                      advect_tend(ims,kms,jms),           &
+                                      ru, rv, ww, mut, mub, mu_old,       &
+                                      config_flags,                       &
+                                      msfu, msfv, msft, fnm, fnp,         &
+                                      rdx, rdy, rdnw,dt,                  &
+                                      ids, ide, jds, jde, kds, kde,       &
+                                      ims, ime, jms, jme, kms, kme,       &
+                                      its, ite, jts, jte, kts, kte     )
+
+
+
+        CALL advect_scalar_mono       ( scalar(ims,kms,jms,im),             &
+                                        scalar_old(ims,kms,jms,im),         &
+                                        advect_tend(ims,kms,jms),           &
+                                        ru, rv, ww, mut, mub, mu_old,       &
+                                        config_flags,                       &
+                                        msfu, msfv, msft, fnm, fnp,         &
+                                        rdx, rdy, rdnw,dt,                  &
+                                        ids, ide, jds, jde, kds, kde,       &
+                                        ims, ime, jms, jme, kms, kme,       &
+                                        its, ite, jts, jte, kts, kte     )
       ELSE
 
-       CALL advect_scalar     ( scalar(ims,kms,jms,im),        &
-                                scalar(ims,kms,jms,im),        &
-                                advect_tend(ims,kms,jms),      &
-                                ru, rv, ww, mut, config_flags, &
-                                msfu, msfv, msft, fnm, fnp,    &
-                                rdx, rdy, rdnw,                &
-                                ids, ide, jds, jde, kds, kde,  &
-                                ims, ime, jms, jme, kms, kme,  &
-                                its, ite, jts, jte, kts, kte  )
+        CALL advect_scalar     ( scalar(ims,kms,jms,im),        &
+                                 scalar(ims,kms,jms,im),        &
+                                 advect_tend(ims,kms,jms),      &
+                                 ru, rv, ww, mut, time_step,    &
+                                 config_flags,                  &
+                                 msfu, msfv, msft, fnm, fnp,    &
+                                 rdx, rdy, rdnw,                &
+                                 ids, ide, jds, jde, kds, kde,  &
+                                 ims, ime, jms, jme, kms, kme,  &
+                                 its, ite, jts, jte, kts, kte  )
       END IF
 
@@ -863,5 +887,6 @@
                                         scalar_tends(ims,kms,jms,im), mut, &
                                         config_flags,                      &
-                                        msfu, msfv, msft, khdq , xkmhd, rdx, rdy, &
+                                        msfu, msfv, msft,                  & 
+                                        khdq , xkmhd, rdx, rdy,            &
                                         ids, ide, jds, jde, kds, kde,      &
                                         ims, ime, jms, jme, kms, kme,      &
@@ -940,6 +965,8 @@
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
          INTENT(INOUT)                                  :: scalar_1,    &
-                                                           scalar_2,  &
-                                                           sc_tend
+                                                           scalar_2
+
+   REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
+         INTENT(IN)                                     :: sc_tend
 
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme ),                &
@@ -1014,4 +1041,5 @@
        DO  k = k_start,k_end
        DO  i = i_start,i_end
+          ! scalar was coupled with m
            tendency(i,k,j) = advect_tend(i,k,j) * msft(i,j)
        ENDDO
@@ -1064,4 +1092,5 @@
        DO  k = k_start,k_end
        DO  i = i_start,i_end
+           ! scalar was coupled with m
            tendency(i,k,j) = advect_tend(i,k,j) * msft(i,j)
        ENDDO
@@ -1104,5 +1133,4 @@
 SUBROUTINE rk_update_scalar_pd( scs, sce,                      &
                                 scalar, sc_tend,               &
-                                msft,                          &
                                 mu_old, mu_new, mu_base,       &
                                 rk_step, dt, spec_zone,        &
@@ -1131,6 +1159,5 @@
    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
                                                           mu_new,  &
-                                                          mu_base, &
-                                                          msft
+                                                          mu_base
 
    INTEGER :: i,j,k,im
@@ -1523,4 +1550,5 @@
 !!****MARS
    IF ( (config_flags%bl_pbl_physics .gt. 0) .OR. (config_flags%modif_wrf) ) THEN
+
       DO J=jts,jtf
       DO K=kts,ktf
@@ -1564,6 +1592,4 @@
 
     ENDIF
-
-
 
 ! fdda
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F	(revision 195)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F	(revision 196)
@@ -2046,5 +2046,4 @@
                                  moist_old(ims,kms,jms,im),                &
                                  moist_tend(ims,kms,jms,im),               &
-                                 grid%msft,                                &
                                  grid%em_mu_1, grid%em_mu_1, grid%em_mub,  &
                                  rk_step, dt_rk, grid%spec_zone,           &
@@ -2111,5 +2110,4 @@
                                  scalar_old(ims,kms,jms,im),              &
                                  scalar_tend(ims,kms,jms,im),             &
-                                 grid%msft,                               &
                                  grid%em_mu_1, grid%em_mu_1, grid%em_mub, &
                                  rk_step, dt_rk, grid%spec_zone,          &
@@ -2178,5 +2176,4 @@
                                  chem_old(ims,kms,jms,im),                &
                                  chem_tend(ims,kms,jms,im),               &
-                                 grid%msft,                               &
                                  grid%em_mu_1, grid%em_mu_1, grid%em_mub, &
                                  rk_step, dt_rk, grid%spec_zone,          &
@@ -2246,5 +2243,4 @@
                                  grid%em_tke_1,                           &
                                  tke_tend(ims,kms,jms),                   &
-                                 grid%msft,                               &
                                  grid%em_mu_1, grid%em_mu_1, grid%em_mub, &
                                  rk_step, dt_rk, grid%spec_zone,          &
