
module bgrid_horiz_diff_mod

!=======================================================================
!
!                    linear horizontal mixing
!               with option for any order accuracy
!
!=======================================================================

use mpp_mod, only: input_nml_file 
use bgrid_horiz_mod      , only: horiz_grid_type
use bgrid_masks_mod      , only: grid_mask_type
use bgrid_prog_var_mod   , only: prog_var_type, var_init
use bgrid_halo_mod       , only: update_halo, vel_flux_boundary, &
                                 EAST, NORTH, NOPOLE, POLEONLY,  &
                                 TEMP, UWND, VWND
use bgrid_change_grid_mod, only: change_grid, TEMP_GRID, WIND_GRID

use         fms_mod, only:  error_mesg, FATAL, write_version_number,         &
                            open_namelist_file, check_nml_error, close_file, &
                            mpp_clock_id, mpp_clock_begin, mpp_clock_end,    &
                            MPP_CLOCK_SYNC, CLOCK_MODULE, file_exist,        &
                            mpp_pe, mpp_root_pe, uppercase, stdlog
use   constants_mod, only:  RADIUS

use  field_manager_mod, only: MODEL_ATMOS, parse
use tracer_manager_mod, only: query_method, get_tracer_names, get_number_tracers

implicit none
private

!--------- public interfaces ----------

public   horiz_diff, horiz_diff_init

!-----------------------------------------------------------------------
!  namelist bgrid_horiz_diff_nml

!  damp_scheme_wind   Determines how horizontal damping coefficients
!  damp_scheme_temp    vary with latitude.  Possible values are: 0,1,2,3.
!                  [0 = no damping; 1 = coeff are uniform with latitude;
!                       2 = coeffs vary as 1/(dx**2+dy**2);
!                       3 = coeffs vary as 1/(dx**2);
!                       4 = same as scheme 1 but with scheme 2 applied
!                       poleward of reflat]
!                      Notes: Schemes 2-4 provide increased damping in
!                      higher latitudes. Temperature and tracers are 
!                      controlled using "damp_scheme_temp".

   integer :: damp_scheme_temp = 1
   integer :: damp_scheme_wind = 1

!  reflat       Latitude cutoff (in degrees) at which increased
!               high latitude damping is applied.  Equatorward of
!               this latitude uniform damping (scheme=1) is applied;
!               poleward of this latitude enhanced damping (scheme=2)
!               is applied.  This variable is only used when
!               damp_scheme_wind=4 or damp_scheme_temp=4.

   real :: reflat = 85.
 
!  damp_order_wind      The horizontal damping order for momentum,
!  damp_order_temp       temperature, and default order for all 
!  damp_order_tracer     prognostic tracers. Only even numbers are allowed.

   integer :: damp_order_wind   = 4   ! use: 0, 2, 4,....
   integer :: damp_order_temp   = 4
   integer :: damp_order_tracer = 4

!  damp_coeff_wind      The horizontal damping coefficients for
!  damp_coeff_temp        momentum, temperature, and default value for
!  damp_coeff_tracer      all prognostic tracers. Only positive values
!                         are allowed.

   real    :: damp_coeff_wind   = 0.35  
   real    :: damp_coeff_temp   = 0.35  
   real    :: damp_coeff_tracer = 0.35

!  slope_corr_wind      The topography slope correction for horizontal
!  slope_corr_temp        damping of momentum, temperature, and default
!  slope_corr_tracer      for all prognostic tracers. 

   real, dimension(4) :: slope_corr_wind   = (/0.,0.,0.,0./)
   real, dimension(4) :: slope_corr_temp   = (/0.,0.,0.,0./)
   real, dimension(4) :: slope_corr_tracer = (/0.,0.,0.,0./)

namelist /bgrid_horiz_diff_nml/ damp_scheme_temp, damp_scheme_wind,      &
                    damp_order_wind, damp_order_temp, damp_order_tracer, &
                    damp_coeff_wind, damp_coeff_temp, damp_coeff_tracer, &
                    slope_corr_wind, slope_corr_temp, slope_corr_tracer, &
                    reflat

!-----------------------------------------------------------------------
!--------- private data ----------

! control parameters for horizontal damping
! set via the namelist or field table
 type hdiff_control_type
     integer, pointer :: order(:) =>NULL()
     real   , pointer :: coeff(:) =>NULL()
     real   , pointer :: slope(:,:) =>NULL()
     logical, pointer :: do_slope_adj(:) =>NULL()
     logical       :: do_damping, do_slope_adj_temp
     integer       :: damping_scheme_wind, damping_scheme_temp
     real, dimension(:),   pointer :: areahx =>NULL(), &
                                      areahy =>NULL(), &
                                      areavx =>NULL(), &
                                      areavy => NULL()
     real, dimension(:,:), pointer :: wth => NULL(), &
                                      wtv =>NULL()
 end type hdiff_control_type

 type(hdiff_control_type), save :: Control

 integer  :: nlev ! number of model levels

 character(len=128) :: version='$Id: bgrid_horiz_diff.F90,v 19.0 2012/01/06 19:54:01 fms Exp $'
 character(len=128) :: tagname='$Name: tikal $'
 logical :: do_log = .true.

!  timing data
 integer :: id_total
 logical :: do_clock_init = .true.

!-----------------------------------------------------------------------

contains

!#######################################################################

 subroutine horiz_diff ( Hgrid, Masks, nplev, dt, dpde, pres, Var, Var_dt )

!-----------------------------------------------------------------------
!
!   Hgrid  = horizontal grid constants
!   Masks  = grid masking constants for eta coordinate
!   nplev  = number of "pure" pressure levels at the top of the model
!   dt     = adjustment time step
!   dpde   = pressure weight for model layers
!   pres   = pressure at full model levels
!   Var    = prognostic variables at the last updated time level
!   Var_dt = tendency of prognostic variables since the last
!            updated time level
!
!-----------------------------------------------------------------------

type(horiz_grid_type), intent(inout)  :: Hgrid
type (grid_mask_type),    intent(in)  :: Masks

integer,                intent(in)    :: nplev
  real,                 intent(in)    :: dt
  real,                 intent(in)    :: dpde(Hgrid%ilb:,Hgrid%jlb:,:),&
                                         pres(Hgrid%ilb:,Hgrid%jlb:,:)
  type (prog_var_type), intent(in)    :: Var
  type (prog_var_type), intent(inout) :: Var_dt

!-----------------------------------------------------------------------

  real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub,       &
                  size(dpde,3)) :: hkew3, hkns3, hkew, hkns, &
                                   hcew, hcns, dat, vdat, hdac

  real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub, &
                     size(dpde,3),3) :: pterms

  real    :: dt_inv, hsign, xcoeff, ycoeff
  integer :: i, j, k, n, is, ie, js, je, ntp

!=======================================================================
!  --- should damping be done ??? ----

   if ( .not. Control%do_damping    ) return
   if ( Control%damping_scheme_wind+Control%damping_scheme_temp == 0 ) return
   call mpp_clock_begin (id_total)

!-----------------------------------------------------------------------
!       --- check the horizontal dimensions of the input array ---

      nlev = size(dpde,3)

    if (size(dpde,1) /= Hgrid%isize .or.               &
        size(dpde,2) /= Hgrid%jsize ) call error_mesg  &
       ('bgrid_horiz_diff', 'input array has the wrong dimensions.', FATAL)

!   ---- time step related values ----

    dt_inv  = 1./dt

!   ---- initialize flux weights ----

    hkew3 = 0.0;  hkns3 = 0.0

!-----------------------------------------------------------------------
!-------------setup temperature and tracer damping ---------------------
!-----------------------------------------------------------------------
   ntp = count( Control%order(1:Var_dt%ntrace) > 0 )
   if (Control%order(0) > 0 .or. ntp > 0) then

      ! compute grid indiecs
      is = Hgrid%Tmp%is;  ie = Hgrid%Tmp%ie
      js = Hgrid%Tmp%js;  je = Hgrid%Tmp%je

      do j = js, je
         hdac(:,j,:) = Hgrid%Tmp%rarea(j)
      enddo
      do k = nplev+1, nlev
      do j = js, je
         hdac(:,j,k) = hdac(:,j,k) / (2.0*dpde(:,j,k))
      enddo
      enddo

!     ----- mass-weighted fluxes -----
      do j = js-1, je
         hkew3(:,j,1:nplev) = Control%areahx(j)
         hkns3(:,j,1:nplev) = Control%areahy(j)
      enddo

      do k = nplev+1, nlev
      do j = js-1, je
      do i = is-1, ie
         hkew3(i,j,k) = Control%areahx(j) * (dpde(i,j,k)+dpde(i+1,j,k))
         hkns3(i,j,k) = Control%areahy(j) * (dpde(i,j,k)+dpde(i,j+1,k))
      enddo
      enddo
      enddo

!     ----- mask out fluxes below step-mountain ----

    if (.not.Masks%sigma) then
       do k = 1, nlev
       do j = js-1, je
       do i = is-1, ie
          hkew3(i,j,k) = hkew3(i,j,k)*Masks%Tmp%mask(i,j,k)*Masks%Tmp%mask(i+1,j,k)
          hkns3(i,j,k) = hkns3(i,j,k)*Masks%Tmp%mask(i,j,k)*Masks%Tmp%mask(i,j+1,k)
       enddo
       enddo
       enddo
          hdac(:,js:je,:) = hdac(:,js:je,:) * Masks%Tmp%mask(:,js:je,:)
    endif

!-----------------------------------------------------------------------
!----- slope adjustment ------------------------------------------------

      if ( Control % do_slope_adj_temp ) then
         call slope_correction_init ( Hgrid, Masks, nplev, pres, pterms )
      else
         hcew = 0.0
         hcns = 0.0
      endif

!-----------------------------------------------------------------------
!------------------temperature damping----------------------------------

      if (Control%order(0) > 0) then

         ! setting up flux weights - limit to 1/8
         do k = 1, nlev
         do j = Hgrid%jlb, Hgrid%jub
            xcoeff = min(0.125,Control%coeff(0)*Control%wth(j,1))
            ycoeff = min(0.125,Control%coeff(0)*Control%wth(j,2))
            hkew(:,j,k) = hkew3(:,j,k) * xcoeff
            hkns(:,j,k) = hkns3(:,j,k) * ycoeff
         enddo
         enddo

         dat  = Var%t + dt*Var_dt%t
!        --- slope adjustent ---
         if ( Control%do_slope_adj(0) ) &
         call slope_correction ( Hgrid, Masks, nplev, Control%slope(:,0), &
                                 pterms, dat, hcew, hcns )
         call diff_mass (Hgrid, dat, hcew, hcns, hkew, hkns, hdac,  &
                         Control%order(0))
         hsign = 1.; if (mod(Control%order(0),4) == 0) hsign = -1.
         Var_dt%t = Var_dt%t + hsign * dt_inv * dat

      endif

!-----------------------------------------------------------------------
!------------ tracer damping (prognostic tracers only) ---------------

      if (ntp > 0) then

         do n = 1, Var_dt%ntrace

            ! setting up flux weights - limit to 1/8
            do k = 1, nlev
            do j = Hgrid%jlb, Hgrid%jub
               xcoeff = min(0.125,Control%coeff(n)*Control%wth(j,1))
               ycoeff = min(0.125,Control%coeff(n)*Control%wth(j,2))
               hkew(:,j,k) = hkew3(:,j,k) * xcoeff
               hkns(:,j,k) = hkns3(:,j,k) * ycoeff
            enddo
            enddo

            if (Control%order(n) == 0) cycle
            dat  = Var%r(:,:,:,n) + dt*Var_dt%r(:,:,:,n)
           !--- slope adjustent ---
            if ( Control%do_slope_adj(n) ) &
            call slope_correction ( Hgrid, Masks, nplev, Control%slope(:,n), &
                                    pterms, dat, hcew, hcns )
            call diff_mass (Hgrid, dat, hcew, hcns, hkew, hkns, &
                            hdac, Control%order(n))
            hsign = 1.; if (mod(Control%order(n),4) == 0) hsign = -1.
            Var_dt%r(:,:,:,n) = Var_dt%r(:,:,:,n) + hsign * dt_inv * dat
         enddo

      endif

   endif
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------setup momentum damping -------------------------------

   if (Control%order(-1) > 0) then

      is = Hgrid%Vel%is;  ie = Hgrid%Vel%ie
      js = Hgrid%Vel%js;  je = Hgrid%Vel%je

      call change_grid (Hgrid, TEMP_GRID, WIND_GRID, &
                        dpde(:,:,nplev+1:nlev), vdat(:,:,nplev+1:nlev))
      call update_halo (Hgrid, UWND, vdat(:,:,nplev+1:nlev), &
                        halos=EAST+NORTH, flags=NOPOLE)

      do j = js, je
         hdac(:,j,:) = Hgrid%Vel%rarea(j)
      enddo
      do k = nplev+1, nlev
         hdac(:,js:je,k) = hdac(:,js:je,k) / (2.*vdat(:,js:je,k))
      enddo

!     ----- mass-weighted fluxes -----
      do j = js, je+1
         hkew3(:,j,1:nplev) = Control%areavx(j)
         hkns3(:,j,1:nplev) = Control%areavy(j)
      enddo

      do k = nplev+1, nlev
      do j = js, je+1
      do i = is, ie+1
         hkew3(i,j,k) = Control%areavx(j)*(vdat(i,j,k)+vdat(i-1,j,k))
         hkns3(i,j,k) = Control%areavy(j)*(vdat(i,j,k)+vdat(i,j-1,k))
      enddo
      enddo
      enddo

!     ----- mask out fluxes below step-mountain ----

    if (.not.Masks%sigma) then
       do k = 1, nlev
       do j = js, je+1
       do i = is, ie+1
          hkew3(i,j,k) = hkew3(i,j,k)*Masks%Vel%mask(i,j,k)*Masks%Vel%mask(i-1,j,k)
          hkns3(i,j,k) = hkns3(i,j,k)*Masks%Vel%mask(i,j,k)*Masks%Vel%mask(i,j-1,k)
       enddo
       enddo
       enddo
          hdac(:,js:je,:) = hdac(:,js:je,:) * Masks%Tmp%mask(:,js:je,:)
    endif

!-----------------------------------------------------------------------
!----- slope adjustment setup ------

    if ( Control % do_slope_adj(-1) ) then
!      ---- pressure at velocity points ----
       vdat(:,:,1:nplev) = pres(:,:,1:nplev)
       call change_grid ( Hgrid, TEMP_GRID, WIND_GRID, &
                          pres(:,:,nplev+1:nlev), vdat(:,:,nplev+1:nlev) )
       call update_halo (Hgrid, UWND, vdat(:,:,nplev+1:nlev), &
                                        halos=EAST+NORTH, flags=NOPOLE)

       call vel_slope_correction_init ( Hgrid, Masks, nplev, &
                                        Control%slope(:,-1), &
                                        vdat, pterms         )
    endif

!-----------------------------------------------------------------------

      ! setting up flux weights - limit to 1/8
      do k = 1, nlev
      do j = js, je+1
         xcoeff = min(0.125,Control%coeff(-1)*Control%wtv(j,1))
         ycoeff = min(0.125,Control%coeff(-1)*Control%wtv(j,2))
         hkew(:,j,k) = hkew3(:,j,k) * xcoeff
         hkns(:,j,k) = hkns3(:,j,k) * ycoeff
      enddo
      enddo
      ! zero-out cross polar fluxes
      call vel_flux_boundary (Hgrid, hkns)

!-----------------------------------------------------------------------
!-------------------------momentum damping------------------------------

       dat = Var%u + dt*Var_dt%u
      vdat = Var%v + dt*Var_dt%v
      if ( Control % do_slope_adj(-1) ) then
         if (Masks%sigma) then
             call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac,  &
                            Control%order(-1), pterms)
         else
             call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac,  &
                            Control%order(-1), pterms, Masks%Vel%kbot)
         endif
      else
         call diff_vel (Hgrid, dat,vdat, hkew,hkns, hdac,  &
                        Control%order(-1))
      endif
      hsign = 1.; if (mod(Control%order(-1),4) == 0) hsign = -1.
      Var_dt%u = Var_dt%u + hsign * dt_inv *  dat
      Var_dt%v = Var_dt%v + hsign * dt_inv * vdat

   endif

   call mpp_clock_end (id_total)

!-----------------------------------------------------------------------

 end subroutine horiz_diff

!#######################################################################
!#######################################################################

 subroutine diff_mass (Hgrid, rdat, hcew, hcns, hkew, hkns, hdac, &
                       order)

!----------------------------------------------------------------------
!
!        diff_mass is a private interface that performs multiple
!        2nd order lapacians.
!
!----------------------------------------------------------------------

   type(horiz_grid_type), intent(inout) :: Hgrid
   integer, intent(in)    :: order
   real   , intent(inout), dimension(Hgrid%ilb:Hgrid%iub, &
                                     Hgrid%jlb:Hgrid%jub, &
                                     nlev) :: rdat
   real   , intent(in),    dimension(Hgrid%ilb:Hgrid%iub, &
                                     Hgrid%jlb:Hgrid%jub, &
                                     nlev) :: hcew, hcns, &
                                              hkew, hkns, hdac

!----------------------------------------------------------------------
!  hcew, hcns are weighted corrections to the diffusive fluxes for
!  sloping sigma surfaces
!----------------------------------------------------------------------

   real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) :: &
              rew, rns
   integer :: i, j, k, n, is, ie, js, je, nordr

!-----------------------------------------------------------------------

   is = Hgrid%Tmp%is;  ie = Hgrid%Tmp%ie
   js = Hgrid%Tmp%js;  je = Hgrid%Tmp%je

!-----------------------------------------------------------------------
!------------loop for order of damping scheme-------------------------

   do n = 1, order/2

        do k = 1, nlev

!------------------contributions (fluxes) ------------------------------

     if ( n == 1 ) then
        do j = js-1, je
        do i = is-1, ie
           rew(i,j) = (rdat(i+1,j,k)-rdat(i,j,k)+hcew(i,j,k))*hkew(i,j,k)
           rns(i,j) = (rdat(i,j+1,k)-rdat(i,j,k)+hcns(i,j,k))*hkns(i,j,k)
        enddo
        enddo
     else
        do j = js-1, je
        do i = is-1, ie
           rew(i,j) = (rdat(i+1,j,k)-rdat(i,j,k))*hkew(i,j,k)
           rns(i,j) = (rdat(i,j+1,k)-rdat(i,j,k))*hkns(i,j,k)
        enddo
        enddo
     endif

!-----------------------------------------------------------------------

     do j = js, je
     do i = is, ie
        rdat(i,j,k)=(rew(i,j)-rew(i-1,j)+rns(i,j)-rns(i,j-1)) &
                    *hdac(i,j,k)
     enddo
     enddo

     enddo

!-----------------------------------------------------------------------
!---- update all halo rows ? ----
!   do not update on last pass, halos will updated in the main program

     if ( n < order/2 ) then
         call update_halo (Hgrid, TEMP, rdat)
     endif

!-----------------------------------------------------------------------

   enddo

!-----------------------------------------------------------------------

 end subroutine diff_mass

!#######################################################################

 subroutine diff_vel (Hgrid, udat, vdat, vkew, vkns, hdac, order, terms, kbot)

!----------------------------------------------------------------------
!
!        diff_vel is a private interface that performs multiple
!        2nd order lapacians for the momentum components.
!
!----------------------------------------------------------------------

   type(horiz_grid_type), intent(inout) :: Hgrid
   integer, intent(in)    :: order
   real   , intent(inout), dimension(Hgrid%ilb:Hgrid%iub, &
                                     Hgrid%jlb:Hgrid%jub, &
                                     nlev) :: udat, vdat
   real   , intent(in),    dimension(Hgrid%ilb:Hgrid%iub, &
                                     Hgrid%jlb:Hgrid%jub, &
                                     nlev) :: vkew, vkns, hdac
   real   , intent(in),    dimension(Hgrid%ilb:Hgrid%iub, &
                                     Hgrid%jlb:Hgrid%jub, &
                                     nlev, 3), optional :: terms
   integer, intent(in),    dimension(Hgrid%ilb:Hgrid%iub,  &
                                     Hgrid%jlb:Hgrid%jub), &
                                               optional :: kbot

!----------------------------------------------------------------------

   real, dimension(Hgrid%ilb:Hgrid%iub, Hgrid%jlb:Hgrid%jub) ::  &
              uew, uns, vew, vns, dudp, dvdp, ucew, ucns, vcew, vcns
   integer :: i, j, k, n, is, ie, js, je, k1, k2

!-----------------------------------------------------------------------

   is = Hgrid%Vel%is;  ie = Hgrid%Vel%ie
   js = Hgrid%Vel%js;  je = Hgrid%Vel%je

!------------loop for order of damping scheme---------------------------

   do n = 1, order/2

      do k = 1, nlev

!------------------contributions (fluxes) ------------------------------

      if ( n == 1 .and. present(terms) ) then

        k1 = max(k-1,1)
        if (present(kbot)) then
           do j = Hgrid%jlb, Hgrid%jub
           do i = Hgrid%ilb, Hgrid%iub
              k2 = min(k+1,kbot(i,j))
              dudp(i,j) = (udat(i,j,k2)-udat(i,j,k1))*terms(i,j,k,1)
              dvdp(i,j) = (vdat(i,j,k2)-vdat(i,j,k1))*terms(i,j,k,1)
           enddo
           enddo
        else
           k2 = min(k+1,nlev)
           dudp(:,:) = (udat(:,:,k2)-udat(:,:,k1))*terms(:,:,k,1)
           dvdp(:,:) = (vdat(:,:,k2)-vdat(:,:,k1))*terms(:,:,k,1)
        endif

        do j = js, je+1
        do i = is, ie+1
!          ---- slope correction terms ----
           ucew(i,j) = (dudp(i,j)+dudp(i-1,j))*terms(i,j,k,2)
           ucns(i,j) = (dudp(i,j)+dudp(i,j-1))*terms(i,j,k,3)
           vcew(i,j) = (dvdp(i,j)+dvdp(i-1,j))*terms(i,j,k,2)
           vcns(i,j) = (dvdp(i,j)+dvdp(i,j-1))*terms(i,j,k,3)

           uew(i,j) = (udat(i,j,k)-udat(i-1,j  ,k)+ucew(i,j))*vkew(i,j,k)
           uns(i,j) = (udat(i,j,k)-udat(i  ,j-1,k)+ucns(i,j))*vkns(i,j,k)
           vew(i,j) = (vdat(i,j,k)-vdat(i-1,j  ,k)+vcew(i,j))*vkew(i,j,k)
           vns(i,j) = (vdat(i,j,k)-vdat(i  ,j-1,k)+vcns(i,j))*vkns(i,j,k)
        enddo
        enddo

      else

        do j = js, je+1
        do i = is, ie+1
           uew(i,j) = (udat(i,j,k)-udat(i-1,j  ,k))*vkew(i,j,k)
           uns(i,j) = (udat(i,j,k)-udat(i  ,j-1,k))*vkns(i,j,k)
           vew(i,j) = (vdat(i,j,k)-vdat(i-1,j  ,k))*vkew(i,j,k)
           vns(i,j) = (vdat(i,j,k)-vdat(i  ,j-1,k))*vkns(i,j,k)
        enddo
        enddo

      endif

!-----------------------------------------------------------------------
      do j = js, je
      do i = is, ie
        udat(i,j,k) = (uew(i+1,j  )-uew(i,j)+ &
                       uns(i  ,j+1)-uns(i,j))*hdac(i,j,k)
        vdat(i,j,k) = (vew(i+1,j  )-vew(i,j)+ &
                       vns(i  ,j+1)-vns(i,j))*hdac(i,j,k)
      enddo
      enddo

      enddo
!-----------------------------------------------------------------------
!---- update all halo rows ? ----
!   do not update on last pass, halos will updated in the main program

      if ( n < order/2 ) then
         call update_halo (Hgrid, UWND, udat)
         call update_halo (Hgrid, VWND, vdat)
      endif

!-----------------------------------------------------------------------

   enddo

!-----------------------------------------------------------------------

 end subroutine diff_vel

!#######################################################################

 subroutine horiz_diff_init ( Hgrid )

!-----------------------------------------------------------------------
!       initialization of horizontal damping coefficients
!-----------------------------------------------------------------------

 type(horiz_grid_type), intent(in) :: Hgrid ! horizontal grid constants

!-----------------------------------------------------------------------

real    :: eps = 1.e-6
integer :: j

 integer            :: n, nv, order, ntrace, unit, ierr, io, logunit
 real               :: coeff, slope(4)
 character(len=128) :: scheme, params, tname

!-----------------------------------------------------------------------
! read namelist
   if (file_exist('input.nml')) then
#ifdef INTERNAL_FILE_NML
       read (input_nml_file, nml=bgrid_horiz_diff_nml, iostat=io)
       ierr = check_nml_error(io,'bgrid_horiz_diff_nml')
#else
       unit = open_namelist_file ( )
       ierr=1; do while (ierr /= 0)
          read (unit, nml=bgrid_horiz_diff_nml, iostat=io, end=5)
          ierr = check_nml_error (io, 'bgrid_horiz_diff_nml')
       enddo
 5     call close_file (unit)
#endif
   endif

   logunit = stdlog()
   if (do_log) then
      call write_version_number (version,tagname)
      if (mpp_pe() == mpp_root_pe()) write (logunit, nml=bgrid_horiz_diff_nml)
      do_log = .false.
   endif

   call get_number_tracers ( MODEL_ATMOS, num_prog=ntrace )

   allocate ( Control%order        (-1:ntrace), &
              Control%coeff        (-1:ntrace), &
              Control%do_slope_adj (-1:ntrace), &
              Control%slope      (4,-1:ntrace)  )

 ! defaults
   Control%order = 4
   Control%coeff = .35
   Control%slope =  0.

 ! namelist arguments
   Control%order  (-1) = damp_order_wind
   Control%coeff  (-1) = damp_coeff_wind
   Control%slope(:,-1) = slope_corr_wind

   Control%order  (0) = damp_order_temp
   Control%coeff  (0) = damp_coeff_temp
   Control%slope(:,0) = slope_corr_temp

   Control%order  (1:ntrace) = damp_order_tracer
   Control%coeff  (1:ntrace) = damp_coeff_tracer
   Control%slope(:,1:ntrace) = spread(slope_corr_tracer,2,ntrace)

 ! set the damping scheme (check for errors below)
   Control % damping_scheme_wind = damp_scheme_wind
   Control % damping_scheme_temp = damp_scheme_temp

 ! process tracer table information for horizontal damping methods
   do n = 1, ntrace
      if (query_method('diff_horiz', MODEL_ATMOS, n, scheme, params)) then
          if (uppercase(trim(scheme)) /= 'LINEAR' .and. uppercase(trim(scheme)) /= 'NONE') &
              call error_mesg  ('bgrid_horiz_diff_mod',  &
                            'invalid damping method, '//uppercase(trim(scheme)), FATAL)
          if (parse(params,'order', order) == 1) Control%order(n) = order
          if (uppercase(trim(scheme)) == 'NONE') Control%order(n) = 0
          if (parse(params,'coeff', coeff) == 1) Control%coeff(n) = coeff
          nv = parse(params,'slope', slope)
          Control%slope (1:nv,n) = slope(1:nv)
      endif
      if (mpp_pe() == mpp_root_pe()) then
         call get_tracer_names (MODEL_ATMOS, n, tname)
         write (logunit,10) n, trim(tname), Control%order(n), Control%coeff(n), Control%slope(:,n)
      10 format (i3, a24, ', Order=',i2, ', Coeff=',f10.5, ', Slope=',4f10.5)
      endif
   enddo

 ! error checking
   do n = -1, ntrace
     !-- order
      if (Control%order(n) < 0 .or. mod(Control%order(n),2) /= 0) &
          call error_mesg ('bgrid_horiz_diff_mod', 'invalid damping order', FATAL)
     !-- non-dimension, normalized coefficient
      if (Control%coeff(n) < 0. .or. Control%coeff(n) > 1.) call error_mesg &
                   ('bgrid_horiz_diff_mod', 'invalid damping coeff', FATAL)
     !-- slope correction weights
      if (minval(Control%slope(:,n)) < 0. .or. maxval(Control%slope(:,n)) > 1.) &
            call error_mesg ('bgrid_horiz_diff_mod', 'invalid slope correction coeff', FATAL)
   enddo

 ! set flags
   Control%do_damping        = maxval(Control%order) > 0
   Control%do_slope_adj_temp = .false.
   do n = -1, ntrace
     Control%do_slope_adj(n) = maxval(Control%slope(:, n)) > 1.e-6
     ! set flag for temp OR tracer slope adjustment
     if (n >= 0 .and. Control%do_slope_adj(n)) Control%do_slope_adj_temp = .true.
   enddo

!-----------------------------------------------------------------------
!----- pre-compute metric terms ------

     allocate ( Control%areahx (Hgrid%jlb:Hgrid%jub), &
                Control%areahy (Hgrid%jlb:Hgrid%jub), &
                Control%areavx (Hgrid%jlb:Hgrid%jub), &
                Control%areavy (Hgrid%jlb:Hgrid%jub), &
                Control%wth    (Hgrid%jlb:Hgrid%jub,2), &
                Control%wtv    (Hgrid%jlb:Hgrid%jub,2)  )

!    ---- areas averaged along axes ----

     do j = Hgrid%Tmp%jsd, Hgrid%Tmp%jed-1
        Control%areahx(j) = 0.5*(Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j))
        Control%areahy(j) = 0.5*(Hgrid%Tmp%area(j)+Hgrid%Tmp%area(j+1))
     enddo

     do j = Hgrid%Vel%jsd+1, Hgrid%Vel%jed
        Control%areavx(j) = 0.5*(Hgrid%Vel%area(j)+Hgrid%Vel%area(j))
        Control%areavy(j) = 0.5*(Hgrid%Vel%area(j)+Hgrid%Vel%area(j-1))
     enddo


!    ---- damping weight of x and y axis varies ----
!    ---- depending which damping scheme is used ----

     call damp_scheme_init ( Hgrid, WIND_GRID, &
                             Control%damping_scheme_wind, Control%wtv )
     call damp_scheme_init ( Hgrid, TEMP_GRID, &
                             Control%damping_scheme_temp, Control%wth )

! initialize code sections for performance timing 

  if (do_clock_init) then
     id_total = mpp_clock_id ('BGRID: horiz_diff (TOTAL)', &
                      flags=MPP_CLOCK_SYNC, grain=CLOCK_MODULE)
     do_clock_init = .false. 
  endif

!-----------------------------------------------------------------------

 end subroutine horiz_diff_init

!#######################################################################

 subroutine damp_scheme_init ( Hgrid, grid, scheme, weights )
 type(horiz_grid_type), intent(in) :: Hgrid
 integer, intent(in)  :: grid, scheme
 real,    intent(out) :: weights(Hgrid%jlb:Hgrid%jub,2)
 
 real, dimension(Hgrid%jlb:Hgrid%jub) :: dx2
 real    :: dxdy2eq, dx2eq, dy2, factor2
 integer :: j
 

!    ---- damping weight of x and y axis varies ----
!    ---- depending which damping scheme is used ----
!
!      scheme 1:   equal/constant
!      scheme 2:   function of diagonal grid distance
!      scheme 3:   function of x-axis grid distance
!

  select case (scheme)

!   ---- uniform damping ----
    case (1)
       weights = 0.125

    case (2:5)
       dx2=0.0
       if (grid == TEMP_GRID) then
       do j = Hgrid%jlb, Hgrid%jub-1
             dx2(j) = Hgrid%Tmp%dx(j)**2
       enddo
             dy2 = Hgrid%Tmp%dy**2
       else if (grid == WIND_GRID) then
       do j = Hgrid%jlb+1, Hgrid%jub
             dx2(j) = Hgrid%Vel%dx(j)**2
       enddo
             dy2 = Hgrid%Vel%dy**2
       else
          ! error condition needed
       endif

       ! function of diagonal distance (poleward of reflat)
       if (scheme == 2 .or. scheme == 4) then
           factor2 = 1./cos(reflat*acos(0.)/90.)**2
           dxdy2eq = RADIUS**2*(Hgrid%dlm**2+Hgrid%dph**2)
           weights(:,1) = 0.125*max(1.,dxdy2eq/(factor2*dx2+dy2))
           if (scheme == 2) then
               weights(:,2) = weights(:,1)
           else
               weights(:,2) = 0.125
           endif
       endif

       ! function of x-distance 
       if (scheme == 3 .or. scheme == 5) then
           dx2eq = (RADIUS*Hgrid%dlm*cos(reflat*acos(0.)/90.))**2
           where (dx2 /= 0.0) weights(:,1) = 0.125*max(1.,dx2eq/dx2)
           if (scheme == 3) then
               weights(:,2) = weights(:,1)
           else
               weights(:,2) = 0.125
           endif
       endif


    case default

       call error_mesg ('bgrid_horiz_diff_mod', 'invalid damping scheme', &
                         FATAL)

  end select

 end subroutine damp_scheme_init

!#######################################################################

 subroutine slope_correction_init ( Hgrid, Masks, nplev, pres, terms )

  type(horiz_grid_type), intent(in)  :: Hgrid
  type (grid_mask_type), intent(in)  :: Masks
  integer,               intent(in)  :: nplev
  real,                  intent(in)  :: pres (Hgrid%ilb:,Hgrid%jlb:,:)
  real,                  intent(out) :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:)
  integer :: i, j, k, k1, k2, nlev

!  initialization of pressure terms for the sigma slope correction
!  these pressure terms do not change between mass variables
!  may want make weight a function of variable and/or level
!      USE ONE-HALF OF SPECIFIED WEIGHT AT LOWEST LEVEL

   nlev = size(pres,3)

   do k = nplev+1, nlev
     !--- reciprocal of vert gradient ---
      k1 = max(k-1,1)
      if (Masks%sigma) then
         k2 = min(k+1,nlev)
         terms(:,:,k,1) = 1.0/(pres(:,:,k2)-pres(:,:,k1))
      else
         do j = Hgrid%jlb, Hgrid%jub
         do i = Hgrid%ilb, Hgrid%iub
            k2 = min(k+1,Masks%Tmp%kbot(i,j))
            if (k2 > k1) then
              terms(i,j,k,1) = 1.0/(pres(i,j,k2)-pres(i,j,k1))
            else
              terms(i,j,k,1) = 0.0
            endif
         enddo
         enddo
      endif

!     --- horiz gradients (flip sign) ---
      do j = Hgrid%jlb, Hgrid%jub-1
      do i = Hgrid%ilb, Hgrid%iub-1
           terms(i,j,k,2) = (pres(i,j,k)-pres(i+1,j,k))
           terms(i,j,k,3) = (pres(i,j,k)-pres(i,j+1,k))
      enddo
      enddo
   enddo

 end subroutine slope_correction_init

!#######################################################################

 subroutine vel_slope_correction_init ( Hgrid, Masks, nplev, weights, pres, terms )

  type(horiz_grid_type), intent(in)  :: Hgrid
  type (grid_mask_type), intent(in)  :: Masks
  integer,               intent(in)  :: nplev
  real,                  intent(in)  :: weights(4)
  real,                  intent(in)  :: pres (Hgrid%ilb:,Hgrid%jlb:,:)
  real,                  intent(out) :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:)
  real    :: wt2
  integer :: i, j, k, k1, k2, ks, nlev, isd, ied, jsd, jed
  real :: dp(size(pres,1),size(pres,2))

!  initialization of pressure terms for the sigma slope correction
!  these pressure terms do not change between mass variables
!  may want make weight a function of variable and/or level
!      USE ONE-HALF OF SPECIFIED WEIGHT AT LOWEST LEVEL

   nlev = size(pres,3)
   isd = Hgrid%Vel%isd; ied = Hgrid%Vel%ied
   jsd = Hgrid%Vel%jsd; jed = Hgrid%Vel%jed

   do k = 1, nplev
      terms(:,:,k,:) = 0.0
   enddo

   do k = nplev+1, nlev
      k1 = max(k-1,1)
      if (Masks%sigma) then
          k2 = min(k+1,nlev)
          dp = (pres(:,:,k2)-pres(:,:,k1))
      else
          do j = Hgrid%jlb, Hgrid%jub
          do i = Hgrid%ilb, Hgrid%iub
             k2 = min(k+1,Masks%Vel%kbot(i,j))
             dp(i,j) = (pres(i,j,k2)-pres(i,j,k1))
          enddo
          enddo
      endif
     !--- reciprocal of vert gradient ---
      where (dp > 0.)
        terms(:,:,k,1) = 1.0/dp
      elsewhere
        terms(:,:,k,1) = 1.e30 ! these values should not be used where it counts
      endwhere

!     --- horiz gradients (flip sign) ---
      ks = max(1,k-nlev+4)
      wt2  = 0.5*weights(ks)
      do j = Hgrid%jlb+1, Hgrid%jub
      do i = Hgrid%ilb+1, Hgrid%iub
           terms(i,j,k,2) = wt2*(pres(i-1,j,k)-pres(i,j,k))
           terms(i,j,k,3) = wt2*(pres(i,j-1,k)-pres(i,j,k))
      enddo
      enddo
   enddo

 end subroutine vel_slope_correction_init

!#######################################################################

 subroutine slope_correction ( Hgrid, Masks, nplev, weights, terms, temp, cew, cns )

  type(horiz_grid_type), intent(in)  :: Hgrid
  type (grid_mask_type), intent(in)  :: Masks
  integer,               intent(in)  :: nplev
  real,                  intent(in)  :: weights(4)
  real,                  intent(in)  :: terms(Hgrid%ilb:,Hgrid%jlb:,:,:)
  real,                  intent(in)  :: temp (Hgrid%ilb:,Hgrid%jlb:,:)
  real,                  intent(out) :: cew  (Hgrid%ilb:,Hgrid%jlb:,:),&
                                        cns  (Hgrid%ilb:,Hgrid%jlb:,:)
  integer :: i, j, k, k1, k2, ks, nlev
  real :: wt2
  real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: dtdp

!  computes weighted-corrections for the slope of sigma surfaces
!  to east-west and north-south fluxes of field temp

 ! check for no correction
   do k = 1, 4
      if (weights(k) > 1.e-6) go to 10
   enddo
   cew = 0.0; cns = 0.0
   return

10 nlev = size(temp,3)

   do k = 1, nplev
      cew(:,:,k) = 0.0
      cns(:,:,k) = 0.0
   enddo

   do k = nplev+1, nlev
      k1 = max(k-1,1)
      if (Masks%sigma) then
         k2 = min(k+1,nlev) ! one-sided derivative at surface
         dtdp(:,:) = (temp(:,:,k2)-temp(:,:,k1))*terms(:,:,k,1)
      else
         do j = Hgrid%jlb, Hgrid%jub
         do i = Hgrid%ilb, Hgrid%iub
            k2 = min(k+1,Masks%Tmp%kbot(i,j))
            dtdp(i,j) = (temp(i,j,k2)-temp(i,j,k1))*terms(i,j,k,1)
         enddo
         enddo
      endif

      ks = max(1,k-nlev+4)
      wt2  = 0.5*weights(ks)
      do j = Hgrid%jlb, Hgrid%jub-1
      do i = Hgrid%ilb, Hgrid%iub-1
           ! note: terms are grouped for reproducibility with previous version
           cew(i,j,k) = (dtdp(i,j)+dtdp(i+1,j))*(terms(i,j,k,2)*wt2)
           cns(i,j,k) = (dtdp(i,j)+dtdp(i,j+1))*(terms(i,j,k,3)*wt2)
      enddo
      enddo
   enddo

 end subroutine slope_correction

!#######################################################################

end module bgrid_horiz_diff_mod

