!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MODULE SMOOTH_MODULE
!
! This module provides routines for smoothing.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module smooth_module

   use parallel_module


   contains

 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: one_two_one
   !
   ! Purpose: Apply the 1-2-1 smoother from the MM5 program TERRAIN 
   !   (found in smth121.F) to array.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine one_two_one(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                          start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(inout) :: array 
  
      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch
  
      allocate(scratch(start_x+1:end_x-1, start_y:end_y, start_z:end_z))
  
      do ipass=1,npass

         do iy=start_y,end_y
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  scratch(ix,iy,iz) = 0.50*array(ix,iy,iz)+0.25*(array(ix-1,iy,iz)+array(ix+1,iy,iz))
               end do
            end do
         end do
   
         do iy=start_y+1,end_y-1
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  array(ix,iy,iz) = 0.50*scratch(ix,iy,iz)+0.25*(scratch(ix,iy-1,iz)+scratch(ix,iy+1,iz))
               end do
             end do
          end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

      end do
  
      deallocate(scratch)
 
   end subroutine one_two_one 
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: smth_desmth
   !
   ! Purpose: Apply the smoother-desmoother from the MM5 program TERRAIN 
   !   (found in smther.F) to array.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine smth_desmth(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                          start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(inout) :: array 
  
      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch
  
      allocate(scratch(start_x+1:end_x-1, start_y:end_y, start_z:end_z))
  
      do ipass=1,npass

         !
         ! Smoothing pass
         !
         do iy=start_y,end_y
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  scratch(ix,iy,iz) = 0.5*array(ix,iy,iz) + 0.25*(array(ix-1,iy,iz)+array(ix+1,iy,iz))
               end do
            end do
         end do
   
         do iy=start_y+1,end_y-1
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  array(ix,iy,iz) = 0.5*scratch(ix,iy,iz) + 0.25*(scratch(ix,iy-1,iz)+scratch(ix,iy+1,iz))
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)
   
         !
         ! Desmoothing pass
         !
         do iy=start_y,end_y
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  scratch(ix,iy,iz) = 1.52*array(ix,iy,iz) - 0.26*(array(ix-1,iy,iz)+array(ix+1,iy,iz))
               end do
            end do
         end do
   
         do iy=start_y+1,end_y-1
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  array(ix,iy,iz) = 1.52*scratch(ix,iy,iz) - 0.26*(scratch(ix,iy-1,iz)+scratch(ix,iy+1,iz))
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

      end do
  
      deallocate(scratch)
 
   end subroutine smth_desmth


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: smth_desmth_special
   !
   ! Purpose: Apply the smoother-desmoother from the MM5 program TERRAIN 
   !   (found in smther.F) to array; however, any grid points that were not
   !   originally negative but which have been smoothed to a negative value
   !   will be restored to their original values.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine smth_desmth_special(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                                  start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval)

      implicit none

      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(inout) :: array

      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch, orig_array

      allocate(scratch(start_x+1:end_x-1, start_y:end_y, start_z:end_z))
      allocate(orig_array(start_x:end_x, start_y:end_y, start_z:end_z))

      orig_array = array

      do ipass=1,npass

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Terrain Modification
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Initially modified by Greg West to remove Sierras from 2 domains
!Revised by Trevor Alcott 27 Apr 2010 to modify terrain in
!12-km West, 4-km IM, 1.33 km GSL regions
!
!we need to differentiate between d01 and d02 and d03
!based on number of east-west gridpoints
!domain 1 e_we = 181
!domain 2 e_we = 163
!domain 3 e_we = 256

!domain 1
      if ( end_x .GT. 186.00 ) then
      do iy=start_y,end_y
         do ix=start_x,end_x
            do iz=start_z,end_z
              if ((ix .GT. 20) .AND. ( ix .LT. 235) .AND. &
(iy .GT. 20) .AND. (iy .LT. 310)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-0.99*(array(ix,iy,iz)-1300.0)
                  END if
              END if
              if ((ix .GT. 10) .AND. ( ix .LE. 20) .AND. &
(iy .GT. 20) .AND. (iy .LT. 310)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((ix-10.0)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GE. 235) .AND. ( ix .LT. 245) .AND. &
(iy .GT. 20) .AND. (iy .LT. 310)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((245.0-ix)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GT. 20) .AND. ( ix .LT. 235) .AND. &
(iy .GT. 10) .AND. (iy .LE. 20)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((iy-10.0)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GT. 20) .AND. ( ix .LT. 235) .AND. &
(iy .GE. 310) .AND. (iy .LT. 320)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((320.0-iy)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GE. 10) .AND. ( ix .LE. 20) .AND. &
(iy .GE. 310) .AND. (iy .LE. 320)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((((ix-10.0)+(320.0-iy))*0.5)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GE. 235) .AND. ( ix .LE. 245) .AND. &
(iy .GE. 310) .AND. (iy .LE. 320)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((((245.0-ix)+(320.0-iy))*0.5)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GE. 10) .AND. ( ix .LE. 20) .AND. &
(iy .GE. 10) .AND. (iy .LE. 20)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((((ix-10.0)+(iy-10.0))*0.5)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
              if ((ix .GE. 235) .AND. ( ix .LE. 245) .AND. &
(iy .GE. 10) .AND. (iy .LE. 20)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-((((245.0-ix)+(iy-10.0))*0.5)*0.099*(array(ix,iy,iz)-1300.0))
                  END if
              END if
            end do
         end do
      end do
     else if ( end_x .GT. 180.00 ) then
      do iy=start_y,end_y
         do ix=start_x,end_x
            do iz=start_z,end_z
               if ((ix .GT. 74) .AND. ( ix .LT. 96) .AND. &
(iy .GT. 68) .AND. (iy .LT. 96)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-0.99*(array(ix,iy,iz)-1300)
                  END if
              END if
            end do
         end do
      end do
      else
        do iy=start_y,end_y
         do ix=start_x,end_x
            do iz=start_z,end_z
              if ((ix .GT. 45) .AND. ( ix .LT. 110) .AND. &
(iy .GT. 32) .AND. (iy .LT. 130)) then
                  if ( array(ix,iy,iz) .GT. 1300.00 ) then
                     array(ix,iy,iz)=array(ix,iy,iz)-0.99*(array(ix,iy,iz)-1300)
                  END if
              END if
            end do
         end do
      end do
      END if

!!!!!!!!!!!!!!!!!!!!!!!!!

         !
         ! Smoothing pass
         !
         do iy=start_y,end_y
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  scratch(ix,iy,iz) = 0.5*array(ix,iy,iz) + 0.25*(array(ix-1,iy,iz)+array(ix+1,iy,iz))
               end do
            end do
         end do

         do iy=start_y+1,end_y-1
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  array(ix,iy,iz) = 0.5*scratch(ix,iy,iz) + 0.25*(scratch(ix,iy-1,iz)+scratch(ix,iy+1,iz))
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

         !
         ! Desmoothing pass
         !
         do iy=start_y,end_y
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  scratch(ix,iy,iz) = 1.52*array(ix,iy,iz) - 0.26*(array(ix-1,iy,iz)+array(ix+1,iy,iz))
               end do
            end do
         end do

         do iy=start_y+1,end_y-1
            do ix=start_x+1,end_x-1
               do iz=start_z,end_z
                  array(ix,iy,iz) = 1.52*scratch(ix,iy,iz) - 0.26*(scratch(ix,iy-1,iz)+scratch(ix,iy+1,iz))
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

      end do

      ! Remove artificially negative values
      do iy=start_y,end_y
         do ix=start_x,end_x
            do iz=start_z,end_z
               if (array(ix,iy,iz) < 0. .and. orig_array(ix,iy,iz) >= 0.) then
                  array(ix,iy,iz) = orig_array(ix,iy,iz)
               end if
            end do
         end do
      end do

      deallocate(scratch)
      deallocate(orig_array)

   end subroutine smth_desmth_special


   !
   ! Smoothing routines for E-grid, contributed by Matthew Pyle
   !

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: one_two_one_egrid
   !
   ! Purpose: Apply the 1-2-1 smoother from the MM5 program TERRAIN 
   !   (found in smth121.F) to array.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine one_two_one_egrid(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                                start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval, hflag)

      implicit none

      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval, hflag
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(inout) :: array

      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch
      integer, dimension(start_y:end_y) :: ihe, ihw, istart

      allocate(scratch(start_x:end_x, start_y:end_y, start_z:end_z))

      do iy=start_y,end_y
         if (hflag == 1.0) then
            ihe(iy) = abs(mod(iy+1,2))
            ihw(iy) = ihe(iy)-1
         else
            ! assign ive,ivw equivs to ihe,ihw
            ihe(iy) = abs(mod(iy,2))
            ihw(iy) = ihe(iy)-1
         end if
      end do

      do iy=start_y,end_y
         if (hflag == 1.0) then
            if (mod(iy,2) == 0) then
               istart(iy) = start_x
            else
               istart(iy) = start_x+1
            end if
         else ! v points
            if (abs(mod(iy,2)) == 1) then
               istart(iy) = start_x
            else
               istart(iy) = start_x+1
            end if
         end if
      end do

      do ipass=1,npass

         do iy=start_y,end_y
            do ix=start_x,end_x
               scratch(ix,iy,1) = array(ix,iy,1) ! for points used in 2nd computation but not defined in 1st computation
            end do
         end do

         ! SW-NE direction
         do iy=start_y+1,end_y-1
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval == 1.0 .and. array(ix,iy,iz) /= 0.) .or. msgval /= 1.0) then
                     scratch(ix,iy,iz) = 0.50*array(ix,iy,iz)+ &
                                      0.25*(array(ix+ihw(iy),iy-1,iz)+array(ix+ihe(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         ! NW-SE direction
         do iy=start_y+1,end_y-1
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval == 1.0 .and. array(ix,iy,iz) /= 0.) .or. msgval /= 1.0) then
                     array(ix,iy,iz) = 0.50*scratch(ix,iy,iz)+ &
                                    0.25*(scratch(ix+ihe(iy),iy-1,iz)+scratch(ix+ihw(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

      end do

      deallocate(scratch)

   end subroutine one_two_one_egrid


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: smth_desmth_egrid_old
   !
   ! Purpose: Apply the smoother-desmoother for E grid
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine smth_desmth_egrid_old(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                                    start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval, hflag)

      implicit none

      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval, hflag
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
               intent(inout) :: array

      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch
      integer, dimension(start_y:end_y) :: ihe, ihw, istart
      real, parameter:: cenwgt = 1.52
      real, parameter:: endwgt = 0.13

      allocate(scratch(start_x:end_x, start_y:end_y, start_z:end_z))

      do iy=start_y,end_y
         if (hflag == 1.0) then
            ihe(iy) = abs(mod(iy+1,2))
            ihw(iy) = ihe(iy)-1
         else
            ! assign ive,ivw equivs to ihe,ihw
            ihe(iy) = abs(mod(iy,2))
            ihw(iy) = ihe(iy)-1
         end if
      end do

      do iy=start_y,end_y
         if (hflag == 1.0) then
            if (mod(iy,2) == 0) then
               istart(iy) = start_x
            else
               istart(iy) = start_x+1
            endif
         else ! v points
            if (abs(mod(iy,2)) == 1) then
               istart(iy) = start_x
            else
               istart(iy) = start_x+1
            endif
         endif
      end do

      do ipass=1,npass

         !
         ! Smoothing pass
         !

         do iy=start_y,end_y
            do ix=start_x,end_x
               scratch(ix,iy,1) = array(ix,iy,1) 
            end do
         end do

         do iy=start_y+1,end_y-1
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval == 1.0 .and. array(ix,iy,iz) /= 0.) .or. msgval /= 1.0) then
                     scratch(ix,iy,iz) = 0.50*array(ix,iy,iz)+ &
                                      0.125*(array(ix+ihw(iy),iy-1,iz)+array(ix+ihe(iy),iy+1,iz)+ &
                                             array(ix+ihw(iy),iy+1,iz)+array(ix+ihe(iy),iy-1,iz))
                  end if
               end do
            end do
         end do


         !
         ! Desmoothing pass
         !

         do iy=start_y+2,end_y-2
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval == 1.0 .and. scratch(ix,iy,iz) /= 0.) .or. msgval /= 1.0) then
                     array(ix,iy,iz) = cenwgt*scratch(ix,iy,iz) - &
                                      endwgt*(scratch(ix+ihw(iy),iy-1,iz)+scratch(ix+ihe(iy),iy+1,iz) + &
                                              scratch(ix+ihw(iy),iy+1,iz)+scratch(ix+ihe(iy),iy-1,iz))
                  end if
               end do
            end do
         end do

      end do

      deallocate(scratch)

   end subroutine smth_desmth_egrid_old


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: smth_desmth_egrid
   !
   ! Purpose: Apply the smoother-desmoother for E grid 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine smth_desmth_egrid(array, start_dom_x, end_dom_x, start_dom_y, end_dom_y, &
                                start_x, end_x, start_y, end_y, start_z, end_z, npass, msgval, hflag)

      implicit none

      ! Arguments
      integer, intent(in) :: start_dom_x, start_dom_y, start_x, start_y, start_z
      integer, intent(in) :: end_dom_x, end_dom_y, end_x, end_y, end_z
      integer, intent(in) :: npass
      real, intent(in) :: msgval, hflag
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
               intent(inout) :: array

      ! Local variables
      integer :: ix, iy, iz, ipass
      real, pointer, dimension(:,:,:) :: scratch
      integer, dimension(start_y:end_y) :: ihe, ihw, istart
      real, parameter :: cenwgt = 1.52
      real, parameter :: endwgt = 0.26

      allocate(scratch(start_x:end_x, start_y:end_y, start_z:end_z))

      do iy=start_y,end_y

         if (hflag .eq. 1.0) then
            ihe(iy)=abs(mod(iy+1,2))
            ihw(iy)=ihe(iy)-1

         ! assign ive,ivw equivs to ihe,ihw
         else
            ihe(iy)=abs(mod(iy,2))
            ihw(iy)=ihe(iy)-1

         end if

      end do

      do iy=start_y,end_y

         if (hflag .eq. 1.0) then
            if (mod(iy,2) .eq. 0) then
               istart(iy)=start_x
            else
               istart(iy)=start_x+1
            endif

         else ! v points
            if (abs(mod(iy,2)) .eq. 1) then
               istart(iy)=start_x
            else
               istart(iy)=start_x+1
            end if

         end if

      end do


      do ipass=1,npass

         !
         ! Smoothing pass
         !

         do iy=start_y,end_y
         do ix=start_x,end_x
            scratch(ix,iy,1)=array(ix,iy,1) ! for points used in 2nd computation but 
                                            !    not defined in 1st
         end do
         end do

         ! SW-NE direction
         do iy=start_y+1,end_y-1
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval .eq. 1.0 .and. array(ix,iy,iz) .ne. 0.) .or. msgval .ne. 1.0) then
                     scratch(ix,iy,iz) = 0.50*array(ix,iy,iz)+ &
                     0.25*(array(ix+ihw(iy),iy-1,iz)+array(ix+ihe(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         ! NW-SE direction
         do iy=start_y+1,end_y-1
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval .eq. 1.0 .and. array(ix,iy,iz) .ne. 0.) .or. msgval .ne. 1.0) then
                     array(ix,iy,iz) = 0.50*scratch(ix,iy,iz)+ &
                     0.25*(scratch(ix+ihe(iy),iy-1,iz)+scratch(ix+ihw(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)



         !
         ! Desmoothing pass
         !

         ! SW-NE direction
         do iy=start_y+2,end_y-2
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval .eq. 1.0 .and. array(ix,iy,iz) .ne. 0.) .or. msgval .ne. 1.0) then
                     scratch(ix,iy,iz) = cenwgt*array(ix,iy,iz) - &
                       endwgt*(array(ix+ihw(iy),iy-1,iz)+array(ix+ihe(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         ! NW-SE direction
         do iy=start_y+2,end_y-2
            do ix=istart(iy),end_x-1
               do iz=start_z,end_z
                  if ( (msgval .eq. 1.0 .and. array(ix,iy,iz) .ne. 0.) .or. msgval .ne. 1.0) then
                     array(ix,iy,iz) = cenwgt*scratch(ix,iy,iz) - &
                       endwgt*(scratch(ix+ihe(iy),iy-1,iz)+scratch(ix+ihw(iy),iy+1,iz))
                  end if
               end do
            end do
         end do

         call exchange_halo_r(array, &
                              start_x, end_x, start_y, end_y, start_z, end_z, &
                              start_dom_x, end_dom_x, start_dom_y, end_dom_y, start_z, end_z)

      end do

      deallocate(scratch)

   end subroutine smth_desmth_egrid

end module smooth_module