#include "cppdefs.h"
      MODULE set_avg_mod
#if defined AVERAGES && (!defined ADJOINT && defined NONLINEAR)
!
!svn $Id: set_avg.F 42 2007-05-04 20:58:11Z arango $
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2007 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This subroutine accumulates and computes output time-averaged       !
!  fields.  Due to synchronization, the time-averaged fields are       !
!  computed in delayed mode. All averages are accumulated at the       !
!  beggining of the next time-step.                                    !
!                                                                      !
!=======================================================================
!
      implicit none

      PRIVATE
      PUBLIC :: set_avg

      CONTAINS
!
!***********************************************************************
      SUBROUTINE set_avg (ng, tile)
!***********************************************************************
!
      USE mod_param
      USE mod_average
      USE mod_forces
# ifdef SOLVE3D
      USE mod_grid
      USE mod_mixing
# endif
      USE mod_ocean
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, iNLM, 5)
# endif
      CALL set_avg_tile (ng, Istr, Iend, Jstr, Jend,                    &
     &                   LBi, UBi, LBj, UBj,                            &
     &                   KOUT,                                          &
# ifdef SOLVE3D
     &                   NOUT,                                          &
     &                   GRID(ng) % pm,                                 &
     &                   GRID(ng) % pn,                                 &
#  ifdef AVERAGES_QUADRATIC
     &                   GRID(ng) % Huon,                               &
     &                   GRID(ng) % Hvom,                               &
#  endif
     &                   OCEAN(ng) % u,                                 &
     &                   OCEAN(ng) % v,                                 &
     &                   OCEAN(ng) % W,                                 &
     &                   OCEAN(ng) % wvel,                              &
     &                   OCEAN(ng) % t,                                 &
     &                   OCEAN(ng) % rho,                               &
#  ifdef AVERAGES_NEARSHORE
     &                   OCEAN(ng) % u_stokes,                          &
     &                   OCEAN(ng) % v_stokes,                          &
     &                   MIXING(ng) % rustr3d,                          &
     &                   MIXING(ng) % rvstr3d,                          &
     &                   MIXING(ng) % Sxx,                              &
     &                   MIXING(ng) % Sxy,                              &
     &                   MIXING(ng) % Syy,                              &
     &                   MIXING(ng) % Szx,                              &
     &                   MIXING(ng) % Szy,                              &
#  endif
#  ifdef LMD_SKPP
     &                   MIXING(ng) % hsbl,                             &
#  endif
#  ifdef LMD_BKPP
     &                   MIXING(ng) % hbbl,                             &
#  endif
#  ifdef AVERAGES_AKV
     &                   MIXING(ng) % Akv,                              &
#  endif
#  if defined AVERAGES_AKT || defined AVERAGES_AKS
     &                   MIXING(ng) % Akt,                              &
#  endif
#  ifdef AVERAGES_FLUXES
     &                   FORCES(ng) % stflx,                            &
#   ifdef BULK_FLUXES
     &                   FORCES(ng) % lhflx,                            &
     &                   FORCES(ng) % shflx,                            &
     &                   FORCES(ng) % lrflx,                            &
#    ifdef EMINUSP
     &                   FORCES(ng) % evap,                             &
     &                   FORCES(ng) % rain,                             &
#    endif
#   endif
#   ifdef SHORTWAVE
     &                   FORCES(ng) % srflx,                            &
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
     &                   FORCES(ng) % sustr,                            &
     &                   FORCES(ng) % svstr,                            &
# endif
     &                   OCEAN(ng) % ubar,                              &
     &                   OCEAN(ng) % vbar,                              &
     &                   OCEAN(ng) % zeta,                              &
# ifdef AVERAGES_NEARSHORE
     &                   OCEAN(ng) % ubar_stokes,                       &
     &                   OCEAN(ng) % vbar_stokes,                       &
     &                   MIXING(ng) % rustr2d,                          &
     &                   MIXING(ng) % rvstr2d,                          &
     &                   MIXING(ng) % Sxx_bar,                          &
     &                   MIXING(ng) % Sxy_bar,                          &
     &                   MIXING(ng) % Syy_bar,                          &
# endif
# ifdef SOLVE3D
     &                   AVERAGE(ng) % avgu3d,                          &
     &                   AVERAGE(ng) % avgv3d,                          &
     &                   AVERAGE(ng) % avgw3d,                          &
     &                   AVERAGE(ng) % avgwvel,                         &
     &                   AVERAGE(ng) % avgt,                            &
     &                   AVERAGE(ng) % avgrho,                          &
#  ifdef AVERAGES_NEARSHORE
     &                   AVERAGE(ng) % avgu3Sd,                         &
     &                   AVERAGE(ng) % avgv3Sd,                         &
     &                   AVERAGE(ng) % avgu3RS,                         &
     &                   AVERAGE(ng) % avgv3RS,                         &
     &                   AVERAGE(ng) % avgSxx3d,                        &
     &                   AVERAGE(ng) % avgSxy3d,                        &
     &                   AVERAGE(ng) % avgSyy3d,                        &
     &                   AVERAGE(ng) % avgSzx3d,                        &
     &                   AVERAGE(ng) % avgSzy3d,                        &
#  endif
#  ifdef LMD_SKPP
     &                   AVERAGE(ng) % avghsbl,                         &
#  endif
#  ifdef LMD_BKPP
     &                   AVERAGE(ng) % avghbbl,                         &
#  endif
#  ifdef AVERAGES_AKV
     &                   AVERAGE(ng) % avgAKv,                          &
#  endif
#  ifdef AVERAGES_AKT
     &                   AVERAGE(ng) % avgAKt,                          &
#  endif
#  ifdef AVERAGES_AKS
     &                   AVERAGE(ng) % avgAKs,                          &
#  endif
#  ifdef AVERAGES_FLUXES
     &                   AVERAGE(ng) % avgstf,                          &
     &                   AVERAGE(ng) % avgswf,                          &
#   ifdef BULK_FLUXES
     &                   AVERAGE(ng) % avglhf,                          &
     &                   AVERAGE(ng) % avgshf,                          &
     &                   AVERAGE(ng) % avglrf,                          &
#    ifdef EMINUSP
     &                   AVERAGE(ng) % avgevap,                         &
     &                   AVERAGE(ng) % avgrain,                         &
#    endif
#   endif
#   ifdef SHORTWAVE
     &                   AVERAGE(ng) % avgsrf,                          &
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
     &                   AVERAGE(ng) % avgsus,                          &
     &                   AVERAGE(ng) % avgsvs,                          &
# endif
# ifdef AVERAGES_QUADRATIC
#  ifdef SOLVE3D
     &                   AVERAGE(ng) % avgHuon,                         &
     &                   AVERAGE(ng) % avgHvom,                         &
     &                   AVERAGE(ng) % avgUU,                           &
     &                   AVERAGE(ng) % avgUV,                           &
     &                   AVERAGE(ng) % avgVV,                           &
     &                   AVERAGE(ng) % avgHuonT,                        &
     &                   AVERAGE(ng) % avgHvomT,                        &
     &                   AVERAGE(ng) % avgUT,                           &
     &                   AVERAGE(ng) % avgVT,                           &
     &                   AVERAGE(ng) % avgTT,                           &
#  endif     
     &                   AVERAGE(ng) % avgU2,                           &
     &                   AVERAGE(ng) % avgV2,                           &
     &                   AVERAGE(ng) % avgZZ,                           &
# endif
# ifdef AVERAGES_NEARSHORE
     &                   AVERAGE(ng) % avgu2Sd,                         &
     &                   AVERAGE(ng) % avgv2Sd,                         &
     &                   AVERAGE(ng) % avgu2RS,                         &
     &                   AVERAGE(ng) % avgv2RS,                         &
     &                   AVERAGE(ng) % avgSxx2d,                        &
     &                   AVERAGE(ng) % avgSxy2d,                        &
     &                   AVERAGE(ng) % avgSyy2d,                        &
# endif
     &                   AVERAGE(ng) % avgu2d,                          &
     &                   AVERAGE(ng) % avgv2d,                          &
     &                   AVERAGE(ng) % avgzeta)

# ifdef PROFILE
      CALL wclock_off (ng, iNLM, 5)
# endif
      RETURN
      END SUBROUTINE set_avg
!
!***********************************************************************
      SUBROUTINE set_avg_tile (ng, Istr, Iend, Jstr, Jend,              &
     &                         LBi, UBi, LBj, UBj,                      &
     &                         Kout,                                    &
# ifdef SOLVE3D
     &                         Nout,                                    &
     &                         pm, pn,                                  &
#  ifdef AVERAGES_QUADRATIC
     &                         Huon, Hvom,                              &
#  endif
     &                         u, v,                                    &
     &                         W, wvel, t, rho,                         &
#  ifdef AVERAGES_NEARSHORE
     &                         u_stokes, v_stokes,                      &
     &                         rustr3d, rvstr3d,                        &
     &                         Sxx, Sxy, Syy, Szx, Szy,                 &
#  endif
#  ifdef LMD_SKPP
     &                         hsbl,                                    &
#  endif
#  ifdef LMD_BKPP
     &                         hbbl,                                    &
#  endif
#  ifdef AVERAGES_AKV
     &                         Akv,                                     &
#  endif
#  if defined AVERAGES_AKT || defined AVERAGES_AKS
     &                         Akt,                                     &
#  endif
#  ifdef AVERAGES_FLUXES
     &                         stflx,                                   &
#   ifdef BULK_FLUXES
     &                         lhflx, shflx, lrflx,                     &
#    ifdef EMINUSP
     &                         evap, rain,                              &
#    endif
#   endif
#   ifdef SHORTWAVE
     &                         srflx,                                   &
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
     &                         sustr, svstr,                            &
# endif
     &                         ubar, vbar,                              &
     &                         zeta,                                    &
#  ifdef AVERAGES_NEARSHORE
     &                         ubar_stokes, vbar_stokes,                &
     &                         rustr2d, rvstr2d,                        &
     &                         Sxx_bar, Sxy_bar, Syy_bar,               &
#  endif
# ifdef SOLVE3D
     &                         avgu3d, avgv3d,                          &
     &                         avgw3d, avgwvel,                         &
     &                         avgt, avgrho,                            &
#  ifdef AVERAGES_NEARSHORE
     &                         avgu3Sd, avgv3Sd,                        &
     &                         avgu3RS, avgv3RS,                        &
     &                         avgSxx3d, avgSxy3d, avgSyy3d,            &
     &                         avgSzx3d, avgSzy3d,                      &
#  endif
#  ifdef LMD_SKPP
     &                         avghsbl,                                 &
#  endif
#  ifdef LMD_BKPP
     &                         avghbbl,                                 &
#  endif
#  ifdef AVERAGES_AKV
     &                         avgAKv,                                  &
#  endif
#  ifdef AVERAGES_AKT
     &                         avgAKt,                                  &
#  endif
#  ifdef AVERAGES_AKS
     &                         avgAKs,                                  &
#  endif
#  ifdef AVERAGES_FLUXES
     &                         avgstf, avgswf,                          &
#   ifdef BULK_FLUXES
     &                         avglhf, avgshf, avglrf,                  &
#    ifdef EMINUSP
     &                         avgevap, avgrain,                        &
#    endif
#   endif
#   ifdef SHORTWAVE
     &                         avgsrf,                                  &
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
     &                         avgsus, avgsvs,                          &
# endif
# ifdef AVERAGES_QUADRATIC
#  ifdef SOLVE3D
     &                         avgHuon, avgHvom,                        &
     &                         avgUU, avgUV, avgVV,                     &
     &                         avgHuonT, avgHvomT,                      &
     &                         avgUT, avgVT, avgTT,                     &
#  endif     
     &                         avgU2, avgV2, avgZZ,                     &
# endif
# ifdef AVERAGES_NEARSHORE
     &                         avgu2Sd, avgv2Sd,                        &
     &                         avgu2RS, avgv2RS,                        &
     &                         avgSxx2d, avgSxy2d, avgSyy2d,            &
# endif
     &                         avgu2d, avgv2d,                          &
     &                         avgzeta)
!***********************************************************************
!
      USE mod_param
      USE mod_scalars
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, Iend, Istr, Jend, Jstr
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: Kout
# ifdef SOLVE3D
      integer, intent(in) :: Nout
# endif
!
# ifdef ASSUMED_SHAPE
#  ifdef SOLVE3D
      real(r8), intent(in) :: pm(LBi:,LBj:)
      real(r8), intent(in) :: pn(LBi:,LBj:)
#   ifdef AVERAGES_QUADRATIC
      real(r8), intent(in) :: Huon(LBi:,LBj:,:)
      real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
#   endif
      real(r8), intent(in) :: u(LBi:,LBj:,:,:)
      real(r8), intent(in) :: v(LBi:,LBj:,:,:)
      real(r8), intent(in) :: W(LBi:,LBj:,0:)
      real(r8), intent(in) :: wvel(LBi:,LBj:,0:)
      real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
      real(r8), intent(in) :: rho(LBi:,LBj:,:)
#   ifdef AVERAGES_NEARSHORE
      real(r8), intent(in) :: u_stokes(LBi:,LBj:,:)
      real(r8), intent(in) :: v_stokes(LBi:,LBj:,:)
      real(r8), intent(in) :: rustr3d(LBi:,LBj:,:)
      real(r8), intent(in) :: rvstr3d(LBi:,LBj:,:)
      real(r8), intent(in) :: Sxx(LBi:,LBj:,:)
      real(r8), intent(in) :: Sxy(LBi:,LBj:,:)
      real(r8), intent(in) :: Syy(LBi:,LBj:,:)
      real(r8), intent(in) :: Szx(LBi:,LBj:,:)
      real(r8), intent(in) :: Szy(LBi:,LBj:,:)
#   endif
#   ifdef LMD_SKPP
      real(r8), intent(in) :: hsbl(LBi:,LBj:)
#   endif
#   ifdef LMD_BKPP
      real(r8), intent(in) :: hbbl(LBi:,LBj:)
#   endif
#   ifdef AVERAGES_AKV
      real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
#   endif
#   if defined AVERAGES_AKT || defined AVERAGES_AKS
      real(r8), intent(in) :: Akt(LBi:,LBj:,0:,:)
#   endif
#   ifdef AVERAGES_FLUXES
      real(r8), intent(in) :: stflx(LBi:,LBj:,:)
#    ifdef BULK_FLUXES
      real(r8), intent(in) :: lhflx(LBi:,LBj:)
      real(r8), intent(in) :: shflx(LBi:,LBj:)
      real(r8), intent(in) :: lrflx(LBi:,LBj:)
#     ifdef EMINUSP
      real(r8), intent(in) :: evap(LBi:,LBj:)
      real(r8), intent(in) :: rain(LBi:,LBj:)
#     endif
#    endif
#    ifdef SHORTWAVE
      real(r8), intent(in) :: srflx(LBi:,LBj:)
#    endif
#   endif
#  endif
#  ifdef AVERAGES_FLUXES
      real(r8), intent(in) :: sustr(LBi:,LBj:)
      real(r8), intent(in) :: svstr(LBi:,LBj:)
#  endif
      real(r8), intent(in) :: ubar(LBi:,LBj:,:)
      real(r8), intent(in) :: vbar(LBi:,LBj:,:)
      real(r8), intent(in) :: zeta(LBi:,LBj:,:)
#  ifdef AVERAGES_NEARSHORE
      real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
      real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
      real(r8), intent(in) :: rustr2d(LBi:,LBj:)
      real(r8), intent(in) :: rvstr2d(LBi:,LBj:)
      real(r8), intent(in) :: Sxx_bar(LBi:,LBj:)
      real(r8), intent(in) :: Sxy_bar(LBi:,LBj:)
      real(r8), intent(in) :: Syy_bar(LBi:,LBj:)
#  endif
#  ifdef SOLVE3D
      real(r8), intent(inout) :: avgu3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgv3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgw3d(LBi:,LBj:,0:)
      real(r8), intent(inout) :: avgwvel(LBi:,LBj:,0:)
      real(r8), intent(inout) :: avgt(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: avgrho(LBi:,LBj:,:)
#   ifdef AVERAGES_NEARSHORE
      real(r8), intent(inout) :: avgu3Sd(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgv3Sd(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgu3RS(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgv3RS(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgSxx3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgSxy3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgSyy3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgSzx3d(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgSzy3d(LBi:,LBj:,:)
#   endif
#   ifdef LMD_SKPP
      real(r8), intent(inout) :: avghsbl(LBi:,LBj:)
#   endif
#   ifdef LMD_BKPP
      real(r8), intent(inout) :: avghbbl(LBi:,LBj:)
#   endif
#   ifdef AVERAGES_AKV
      real(r8), intent(inout) :: avgAKv(LBi:,LBj:,0:)
#   endif
#   ifdef AVERAGES_AKT
      real(r8), intent(inout) :: avgAKt(LBi:,LBj:,0:)
#   endif
#   ifdef AVERAGES_AKS
      real(r8), intent(inout) :: avgAKs(LBi:,LBj:,0:)
#   endif
#   ifdef AVERAGES_FLUXES
      real(r8), intent(inout) :: avgstf(LBi:,LBj:)
      real(r8), intent(inout) :: avgswf(LBi:,LBj:)
#    ifdef BULK_FLUXES
      real(r8), intent(inout) :: avglhf(LBi:,LBj:)
      real(r8), intent(inout) :: avgshf(LBi:,LBj:)
      real(r8), intent(inout) :: avglrf(LBi:,LBj:)
#     ifdef EMINUSP
      real(r8), intent(inout) :: avgevap(LBi:,LBj:)
      real(r8), intent(inout) :: avgrain(LBi:,LBj:)
#     endif
#    endif
#    ifdef SHORTWAVE
      real(r8), intent(inout) :: avgsrf(LBi:,LBj:)
#    endif
#   endif
#  endif
#  ifdef AVERAGES_FLUXES
      real(r8), intent(inout) :: avgsus(LBi:,LBj:)
      real(r8), intent(inout) :: avgsvs(LBi:,LBj:)
#  endif
#  ifdef AVERAGES_QUADRATIC
#   ifdef SOLVE3D
      real(r8), intent(inout) :: avgHuon(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgHvom(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgUU(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgUV(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgVV(LBi:,LBj:,:)
      real(r8), intent(inout) :: avgHuonT(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: avgHvomT(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: avgUT(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: avgVT(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: avgTT(LBi:,LBj:,:,:)
#   endif     
      real(r8), intent(inout) :: avgU2(LBi:,LBj:)
      real(r8), intent(inout) :: avgV2(LBi:,LBj:)
      real(r8), intent(inout) :: avgZZ(LBi:,LBj:)
#  endif
#  ifdef AVERAGES_NEARSHORE
      real(r8), intent(inout) :: avgu2Sd(LBi:,LBj:)
      real(r8), intent(inout) :: avgv2Sd(LBi:,LBj:)
      real(r8), intent(inout) :: avgu2RS(LBi:,LBj:)
      real(r8), intent(inout) :: avgv2RS(LBi:,LBj:)
      real(r8), intent(inout) :: avgSxx2d(LBi:,LBj:)
      real(r8), intent(inout) :: avgSxy2d(LBi:,LBj:)
      real(r8), intent(inout) :: avgSyy2d(LBi:,LBj:)
#  endif
      real(r8), intent(inout) :: avgu2d(LBi:,LBj:)
      real(r8), intent(inout) :: avgv2d(LBi:,LBj:)
      real(r8), intent(inout) :: avgzeta(LBi:,LBj:)

# else

#  ifdef SOLVE3D
      real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
#   ifdef AVERAGES_QUADRATIC
      real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
#   endif
      real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: wvel(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
      real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
#   ifdef AVERAGES_NEARSHORE
      real(r8), intent(in) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: rustr3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: rvstr3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Sxx(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Sxy(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Syy(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Szx(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Szy(LBi:UBi,LBj:UBj,N(ng))
#   endif
#   ifdef LMD_SKPP
      real(r8), intent(in) :: hsbl(LBi:UBi,LBj:UBj)
#   endif
#   ifdef LMD_BKPP
      real(r8), intent(in) :: hbbl(LBi:UBi,LBj:UBj)
#   endif
#   ifdef AVERAGES_AKV
      real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
#   endif
#   if defined AVERAGES_AKT || defined AVERAGES_AKS
      real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
#   endif
#   ifdef AVERAGES_FLUXES
      real(r8), intent(in) :: stflx(LBi:UBi,LBj:UBj,NT(ng))
#    ifdef BULK_FLUXES
      real(r8), intent(in) :: lhflx(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: shflx(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: lrflx(LBi:UBi,LBj:UBj)
#     ifdef EMINUSP
      real(r8), intent(in) :: evap(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: rain(LBi:UBi,LBj:UBj)
#     endif
#    endif
#    ifdef SHORTWAVE
      real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
#    endif
#   endif
#  endif
#  ifdef AVERAGES_FLUXES
      real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(in) :: ubar(LBi:UBi,LBj:UBj,3)
      real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,3)
      real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
#   ifdef AVERAGES_NEARSHORE
      real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: rustr2d(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: rvstr2d(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Sxx_bar(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Sxy_bar(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Syy_bar(LBi:UBi,LBj:UBj)
#   endif
#  ifdef SOLVE3D
      real(r8), intent(inout) :: avgu3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgv3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgw3d(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(inout) :: avgwvel(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(inout) :: avgt(LBi:UBi,LBj:UBj,N(ng),NT(ng))
      real(r8), intent(inout) :: avgrho(LBi:UBi,LBj:UBj,N(ng))
#   ifdef AVERAGES_NEARSHORE
      real(r8), intent(inout) :: avgu3Sd(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgv3Sd(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgu3RS(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgv3RS(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgSxx3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgSxy3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgSyy3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgSzx3d(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgSzy3d(LBi:UBi,LBj:UBj,N(ng))
#   endif
#   ifdef LMD_SKPP
      real(r8), intent(inout) :: avghsbl(LBi:UBi,LBj:UBj)
#   endif
#   ifdef LMD_BKPP
      real(r8), intent(inout) :: avghbbl(LBi:UBi,LBj:UBj)
#   endif
#   ifdef AVERAGES_AKV
      real(r8), intent(inout) :: avgAKv(LBi:UBi,LBj:UBj,0:N(ng))
#   endif
#   ifdef AVERAGES_AKT
      real(r8), intent(inout) :: avgAKt(LBi:UBi,LBj:UBj,0:N(ng))
#   endif
#   ifdef AVERAGES_AKS
      real(r8), intent(inout) :: avgAKs(LBi:UBi,LBj:UBj,0:N(ng))
#   endif
#   ifdef AVERAGES_FLUXES
      real(r8), intent(inout) :: avgstf(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgswf(LBi:UBi,LBj:UBj)
#    ifdef BULK_FLUXES
      real(r8), intent(inout) :: avglhf(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgshf(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avglrf(LBi:UBi,LBj:UBj)
#     ifdef EMINUSP
      real(r8), intent(inout) :: avgevap(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgrain(LBi:UBi,LBj:UBj)
#     endif
#    endif
#    ifdef SHORTWAVE
      real(r8), intent(inout) :: avgsrf(LBi:UBi,LBj:UBj)
#    endif
#   endif
#  endif
#  ifdef AVERAGES_FLUXES
      real(r8), intent(inout) :: avgsus(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgsvs(LBi:UBi,LBj:UBj)
#  endif
#  ifdef AVERAGES_QUADRATIC
#   ifdef SOLVE3D
      real(r8), intent(inout) :: avgHuon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgHvom(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgUU(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgUV(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgVV(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: avgHuonT(LBi:UBi,LBj:UBj,N(ng),NAT)
      real(r8), intent(inout) :: avgHvomT(LBi:UBi,LBj:UBj,N(ng),NAT)
      real(r8), intent(inout) :: avgUT(LBi:UBi,LBj:UBj,N(ng),NAT)
      real(r8), intent(inout) :: avgVT(LBi:UBi,LBj:UBj,N(ng),NAT)
      real(r8), intent(inout) :: avgTT(LBi:UBi,LBj:UBj,N(ng),NAT)
#   endif     
      real(r8), intent(inout) :: avgU2(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgV2(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgZZ(LBi:UBi,LBj:UBj)
#  endif
#  ifdef AVERAGES_NEARSHORE
      real(r8), intent(inout) :: avgu2Sd(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgv2Sd(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgu2RS(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgv2RS(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgSxx2d(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgSxy2d(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgSyy2d(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(inout) :: avgu2d(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgv2d(LBi:UBi,LBj:UBj)
      real(r8), intent(inout) :: avgzeta(LBi:UBi,LBj:UBj)
# endif
!
!  Local variable declarations.
!
      integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
      integer :: i, itrc, j, k

      real(r8) :: fac

# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Initialize time-averaged arrays when appropriate.  Notice that
! fields are initilized twice during re-start.  However, the time-
! averaged fields are computed correctly.
!-----------------------------------------------------------------------
!
      IF (((iic(ng).gt.ntsAVG(ng)).and.                                 &
     &     (MOD(iic(ng)-1,nAVG(ng)).eq.1)).or.                          &
     &    ((nrrec(ng).gt.0).and.(iic(ng).eq.ntstart(ng)))) THEN
!
!  Initialize 2D fields.
!
        DO j=JstrR,JendR
          DO i=IstrR,IendR
            avgzeta(i,j)=zeta(i,j,Kout)
            avgu2d (i,j)=ubar(i,j,Kout)
            avgv2d (i,j)=vbar(i,j,Kout)
# ifdef AVERAGES_NEARSHORE
            avgu2Sd(i,j)=ubar_stokes(i,j)
            avgv2Sd(i,j)=vbar_stokes(i,j)
            avgu2RS(i,j)=rustr2d(i,j)
            avgv2RS(i,j)=rvstr2d(i,j)
            avgSxx2d(i,j)=Sxx_bar(i,j)
            avgSxy2d(i,j)=Sxy_bar(i,j)
            avgSyy2d(i,j)=Syy_bar(i,j)
# endif
# ifdef AVERAGES_QUADRATIC
            avgZZ(i,j)=zeta(i,j,Kout)*zeta(i,j,Kout)
            avgU2(i,j)=ubar(i,j,Kout)*ubar(i,j,Kout)
            avgV2(i,j)=vbar(i,j,Kout)*vbar(i,j,Kout)
# endif
# ifdef SOLVE3D
#  ifdef LMD_SKPP
            avghsbl(i,j)=hsbl(i,j)
#  endif
#  ifdef LMD_BKPP
            avghbbl(i,j)=hbbl(i,j)
#  endif
#  ifdef AVERAGES_FLUXES
            avgstf(i,j)=stflx(i,j,itemp)
            avgswf(i,j)=stflx(i,j,isalt)
#   ifdef BULK_FLUXES
            avglhf(i,j)=lhflx(i,j)
            avgshf(i,j)=shflx(i,j)
            avglrf(i,j)=lrflx(i,j)
#    ifdef EMINUSP
            avgevap(i,j)=evap(i,j)
            avgrain(i,j)=rain(i,j)
#    endif
#   endif
#   ifdef SHORTWAVE
            avgsrf(i,j)=srflx(i,j)
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
            avgsus(i,j)=sustr(i,j)
            avgsvs(i,j)=svstr(i,j)
# endif
          END DO
        END DO

# ifdef SOLVE3D
!
!  Initialize fields associated with 3D horizontal momentum.
!
        DO k=1,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgu3d(i,j,k)=u(i,j,k,Nout)
              avgv3d(i,j,k)=v(i,j,k,Nout)
              avgrho(i,j,k)=rho(i,j,k)
#  ifdef AVERAGES_NEARSHORE
              avgu3Sd(i,j,k)=u_stokes(i,j,k)
              avgv3Sd(i,j,k)=v_stokes(i,j,k)
              avgu3RS(i,j,k)=rustr3d(i,j,k)
              avgv3RS(i,j,k)=rvstr3d(i,j,k)
              avgSxx3d(i,j,k)=Sxx(i,j,k)
              avgSxy3d(i,j,k)=Sxy(i,j,k)
              avgSyy3d(i,j,k)=Syy(i,j,k)
              avgSzx3d(i,j,k)=Szx(i,j,k)
              avgSzy3d(i,j,k)=Szy(i,j,k)
#  endif
#  ifdef AVERAGES_QUADRATIC
              avgHuon(i,j,k)=Huon(i,j,k)
              avgHvom(i,j,k)=Hvom(i,j,k)
              avgUU(i,j,k)=u(i,j,k,Nout)*u(i,j,k,Nout)
              avgVV(i,j,k)=v(i,j,k,Nout)*v(i,j,k,Nout)
#  endif
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          DO j=Jstr,Jend
            DO i=Istr,Iend
              avgUV(i,j,k)=0.25_r8*(u(i  ,j  ,k,Nout)+                  &
     &                              u(i+1,j  ,k,Nout))*                 &
     &                             (v(i  ,j  ,k,Nout)+                  &
     &                              v(i  ,j+1,k,Nout))
            END DO
          END DO
#  endif
        END DO
!
!  Initialized fields associated with tracers.
!
        DO itrc=1,NT(ng)
          DO k=1,N(ng)
            DO j=JstrR,JendR
              DO i=IstrR,IendR
                avgt(i,j,k,itrc)=t(i,j,k,Nout,itrc)
              END DO
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          IF (itrc.le.NAT) THEN
            DO k=1,N(ng)
              DO j=JstrR,JendR
                DO i=IstrR,IendR
                  avgTT(i,j,k,itrc)=t(i,j,k,Nout,itrc)*                 &
     &                              t(i,j,k,Nout,itrc)
                END DO
                DO i=Istr,Iend
                  avgUT(i,j,k,itrc)=0.5_r8*u(i,j,k,Nout)*               &
     &                              (t(i-1,j,k,Nout,itrc)+              &
     &                               t(i  ,j,k,Nout,itrc))
                  avgHuonT(i,j,k,itrc)=0.5_r8*Huon(i,j,k)*              &
     &                                 (t(i-1,j,k,Nout,itrc)+           &
     &                                  t(i  ,j,k,Nout,itrc))
                END DO
              END DO
              DO j=Jstr,Jend
                DO i=IstrR,IendR
                  avgVT(i,j,k,itrc)=0.5_r8*v(i,j,k,Nout)*               &
     &                              (t(i,j-1,k,Nout,itrc)+              &
     &                               t(i,j  ,k,Nout,itrc))
                  avgHvomT(i,j,k,itrc)=0.5_r8*Hvom(i,j,k)*              &
     &                                 (t(i,j-1,k,Nout,itrc)+           &
     &                                  t(i,j  ,k,Nout,itrc))
                END DO
              END DO
            END DO
          END IF
#  endif
        END DO
!
!  Initialize fields at W-points.
!
        DO k=0,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgw3d(i,j,k)=W(i,j,k)*pm(i,j)*pn(i,j)
              avgwvel(i,j,k)=wvel(i,j,k)
#  ifdef AVERAGES_AKV
              avgAKv(i,j,k)=Akv(i,j,k)
#  endif
#  ifdef AVERAGES_AKT
              avgAKt(i,j,k)=Akt(i,j,k,itemp)
#  endif
#  ifdef AVERAGES_AKS
              avgAKs(i,j,k)=Akt(i,j,k,isalt)
#  endif
            END DO
          END DO
        END DO
# endif
!
!-----------------------------------------------------------------------
!  Accumulate time-averaged fields.
!-----------------------------------------------------------------------
!
      ELSE IF (iic(ng).gt.ntsAVG(ng)) THEN
!
!  Accumulate 2D fields.
!
        DO j=JstrR,JendR
          DO i=IstrR,IendR
            avgzeta(i,j)=avgzeta(i,j)+zeta(i,j,Kout)
            avgu2d (i,j)=avgu2d (i,j)+ubar(i,j,Kout)
            avgv2d (i,j)=avgv2d (i,j)+vbar(i,j,Kout)
# ifdef AVERAGES_NEARSHORE
            avgu2Sd(i,j)=avgu2Sd(i,j)+ubar_stokes(i,j)
            avgv2Sd(i,j)=avgv2Sd(i,j)+vbar_stokes(i,j)
            avgu2RS(i,j)=avgu2RS(i,j)+rustr2d(i,j)
            avgv2RS(i,j)=avgv2RS(i,j)+rvstr2d(i,j)
            avgSxx2d(i,j)=avgSxx2d(i,j)+Sxx_bar(i,j)
            avgSxy2d(i,j)=avgSxy2d(i,j)+Sxy_bar(i,j)
            avgSyy2d(i,j)=avgSyy2d(i,j)+Syy_bar(i,j)
# endif
# ifdef AVERAGES_QUADRATIC
            avgZZ(i,j)=avgZZ(i,j)+zeta(i,j,Kout)*zeta(i,j,Kout)
            avgU2(i,j)=avgU2(i,j)+ubar(i,j,Kout)*ubar(i,j,Kout)
            avgV2(i,j)=avgU2(i,j)+vbar(i,j,Kout)*vbar(i,j,Kout)
# endif
# ifdef SOLVE3D
#  ifdef LMD_SKPP
            avghsbl(i,j)=avghsbl(i,j)+hsbl(i,j)
#  endif
#  ifdef LMD_BKPP
            avghbbl(i,j)=avghbbl(i,j)+hbbl(i,j)
#  endif
#  ifdef AVERAGES_FLUXES
            avgstf(i,j)=avgstf(i,j)+stflx(i,j,itemp)
            avgswf(i,j)=avgswf(i,j)+stflx(i,j,isalt)
#   ifdef BULK_FLUXES
            avglhf(i,j)=avglhf(i,j)+lhflx(i,j)
            avgshf(i,j)=avgshf(i,j)+shflx(i,j)
            avglrf(i,j)=avglrf(i,j)+lrflx(i,j)
#    ifdef EMINUSP
            avgevap(i,j)=avgevap(i,j)+evap(i,j)
            avgrain(i,j)=avgrain(i,j)+rain(i,j)
#    endif
#   endif
#   ifdef SHORTWAVE
            avgsrf(i,j)=avgsrf(i,j)+srflx(i,j)
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
            avgsus(i,j)=avgsus(i,j)+sustr(i,j)
            avgsvs(i,j)=avgsvs(i,j)+svstr(i,j)
# endif
          END DO
        END DO

# ifdef SOLVE3D
!
!  Accumulate fields associated with 3D horizontal momentum.
!
        DO k=1,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgu3d(i,j,k)=avgu3d(i,j,k)+u(i,j,k,Nout)
              avgv3d(i,j,k)=avgv3d(i,j,k)+v(i,j,k,Nout)
              avgrho(i,j,k)=avgrho(i,j,k)+rho(i,j,k)
#  ifdef AVERAGES_NEARSHORE
              avgu3Sd(i,j,k)=avgu3Sd(i,j,k)+u_stokes(i,j,k)
              avgv3Sd(i,j,k)=avgv3Sd(i,j,k)+v_stokes(i,j,k)
              avgu3RS(i,j,k)=avgu3RS(i,j,k)+rustr3d(i,j,k)
              avgv3RS(i,j,k)=avgv3RS(i,j,k)+rvstr3d(i,j,k)
              avgSxx3d(i,j,k)=avgSxx3d(i,j,k)+Sxx(i,j,k)
              avgSxy3d(i,j,k)=avgSxy3d(i,j,k)+Sxy(i,j,k)
              avgSyy3d(i,j,k)=avgSyy3d(i,j,k)+Syy(i,j,k)
              avgSzx3d(i,j,k)=avgSzx3d(i,j,k)+Szx(i,j,k)
              avgSzy3d(i,j,k)=avgSzy3d(i,j,k)+Szy(i,j,k)
#  endif
#  ifdef AVERAGES_QUADRATIC
              avgHuon(i,j,k)=avgHuon(i,j,k)+Huon(i,j,k)
              avgHvom(i,j,k)=avgHvom(i,j,k)+Hvom(i,j,k)
              avgUU(i,j,k)=avgUU(i,j,k)+u(i,j,k,Nout)*u(i,j,k,Nout)
              avgVV(i,j,k)=avgVV(i,j,k)+v(i,j,k,Nout)*v(i,j,k,Nout)
#  endif
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          DO j=Jstr,Jend
            DO i=Istr,Iend
              avgUV(i,j,k)=avgUV(i,j,k)+                                &
     &                     0.25_r8*(u(i  ,j  ,k,Nout)+                  &
     &                              u(i+1,j  ,k,Nout))*                 &
     &                             (v(i  ,j  ,k,Nout)+                  &
     &                              v(i  ,j+1,k,Nout))
            END DO
          END DO
#  endif
        END DO
!
!  Accumulate fields associated with tracers.
!
        DO itrc=1,NT(ng)
          DO k=1,N(ng)
            DO j=JstrR,JendR
              DO i=IstrR,IendR
                avgt(i,j,k,itrc)=avgt(i,j,k,itrc)+t(i,j,k,Nout,itrc)
              END DO
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          IF (itrc.le.NAT) THEN
            DO k=1,N(ng)
              DO j=JstrR,JendR
                DO i=IstrR,IendR
                  avgTT(i,j,k,itrc)=avgTT(i,j,k,itrc)+                  &
     &                              t(i,j,k,Nout,itrc)*                 &
     &                              t(i,j,k,Nout,itrc)
                END DO
                DO i=Istr,Iend
                  avgUT(i,j,k,itrc)=avgUT(i,j,k,itrc)+                  &
     &                              0.5_r8*u(i,j,k,Nout)*               &
     &                              (t(i-1,j,k,Nout,itrc)+              &
     &                               t(i  ,j,k,Nout,itrc))
                  avgHuonT(i,j,k,itrc)=avgHuonT(i,j,k,itrc)+            &
     &                                 0.5_r8*Huon(i,j,k)*              &
     &                                 (t(i-1,j,k,Nout,itrc)+           &
     &                                  t(i  ,j,k,Nout,itrc))
                END DO
              END DO
              DO j=Jstr,Jend
                DO i=IstrR,IendR
                  avgVT(i,j,k,itrc)=avgVT(i,j,k,itrc)+                  &
     &                              0.5_r8*v(i,j,k,Nout)*               &
     &                              (t(i,j-1,k,Nout,itrc)+              &
     &                               t(i,j  ,k,Nout,itrc))
                  avgHvomT(i,j,k,itrc)=avgVT(i,j,k,itrc)+               &
     &                                 0.5_r8*Hvom(i,j,k)*              &
     &                                 (t(i,j-1,k,Nout,itrc)+           &
     &                                  t(i,j  ,k,Nout,itrc))
                END DO
              END DO
            END DO
          END IF
#  endif
        END DO
!
!  Accumulate fields at W-points.
!
        DO k=0,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgw3d(i,j,k)=avgw3d(i,j,k)+W(i,j,k)*pm(i,j)*pn(i,j)
              avgwvel(i,j,k)=avgwvel(i,j,k)+wvel(i,j,k)
#  ifdef AVERAGES_AKV
              avgAKv(i,j,k)=avgAKv(i,j,k)+Akv(i,j,k)
#  endif
#  ifdef AVERAGES_AKT
              avgAKt(i,j,k)=avgAKt(i,j,k)+Akt(i,j,k,itemp)
#  endif
#  ifdef AVERAGES_AKS
              avgAKs(i,j,k)=avgAKs(i,j,k)+Akt(i,j,k,isalt)
#  endif
            END DO
          END DO
        END DO
# endif
      END IF
!
!-----------------------------------------------------------------------
!  Convert accumulated sums into time-averages, if appropriate.
!-----------------------------------------------------------------------
!
      IF ((iic(ng).gt.ntsAVG(ng)).and.                                  &
     &    (MOD(iic(ng)-1,nAVG(ng)).eq.0).and.                           &
     &    ((iic(ng).ne.ntstart(ng)).or.(nrrec(ng).eq.0))) THEN
        fac=1.0_r8/REAL(nAVG(ng),r8)
        IF (SOUTH_WEST_TEST) THEN
          AVGtime(ng)=AVGtime(ng)+REAL(nAVG(ng),r8)*dt(ng)
        END IF
!
!  Process 2D fields.
!
        DO j=JstrR,JendR
          DO i=IstrR,IendR
            avgzeta(i,j)=fac*avgzeta(i,j)
            avgu2d (i,j)=fac*avgu2d (i,j)
            avgv2d (i,j)=fac*avgv2d (i,j)
# ifdef AVERAGES_NEARSHORE
            avgu2Sd(i,j)=fac*avgu2Sd(i,j)
            avgv2Sd(i,j)=fac*avgv2Sd(i,j)
            avgu2RS(i,j)=fac*avgu2RS(i,j)
            avgv2RS(i,j)=fac*avgv2RS(i,j)
            avgSxx2d(i,j)=fac*avgSxx2d(i,j)
            avgSxy2d(i,j)=fac*avgSxy2d(i,j)
            avgSyy2d(i,j)=fac*avgSyy2d(i,j)
# endif
# ifdef AVERAGES_QUADRATIC
            avgZZ(i,j)=fac*avgZZ(i,j)
            avgU2(i,j)=fac*avgU2(i,j)
            avgV2(i,j)=fac*avgU2(i,j)
# endif
# ifdef SOLVE3D
#  ifdef LMD_SKPP
            avghsbl(i,j)=fac*avghsbl(i,j)
#  endif
#  ifdef LMD_BKPP
            avghbbl(i,j)=fac*avghbbl(i,j)
#  endif
#  ifdef AVERAGES_FLUXES
            avgstf(i,j)=fac*avgstf(i,j)
            avgswf(i,j)=fac*avgswf(i,j)
#   ifdef BULK_FLUXES
            avglhf(i,j)=fac*avglhf(i,j)
            avgshf(i,j)=fac*avgshf(i,j)
            avglrf(i,j)=fac*avglrf(i,j)
#    ifdef EMINUSP
            avgevap(i,j)=fac*avgevap(i,j)
            avgrain(i,j)=fac*avgrain(i,j)
#    endif
#   endif
#   ifdef SHORTWAVE
            avgsrf(i,j)=fac*avgsrf(i,j)
#   endif
#  endif
# endif
# ifdef AVERAGES_FLUXES
            avgsus(i,j)=fac*avgsus(i,j)
            avgsvs(i,j)=fac*avgsvs(i,j)
# endif
          END DO
        END DO

# ifdef SOLVE3D
!
!  Process fields associated with 3D horizontal momentum.
!
        DO k=1,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgu3d(i,j,k)=fac*avgu3d(i,j,k)
              avgv3d(i,j,k)=fac*avgv3d(i,j,k)
              avgrho(i,j,k)=fac*avgrho(i,j,k)
#  ifdef AVERAGES_NEARSHORE
              avgu3Sd(i,j,k)=fac*avgu3Sd(i,j,k)
              avgv3Sd(i,j,k)=fac*avgv3Sd(i,j,k)
              avgu3RS(i,j,k)=fac*avgu3RS(i,j,k)
              avgv3RS(i,j,k)=fac*avgv3RS(i,j,k)
              avgSxx3d(i,j,k)=fac*avgSxx3d(i,j,k)
              avgSxy3d(i,j,k)=fac*avgSxy3d(i,j,k)
              avgSyy3d(i,j,k)=fac*avgSyy3d(i,j,k)
              avgSzx3d(i,j,k)=fac*avgSzx3d(i,j,k)
              avgSzy3d(i,j,k)=fac*avgSzy3d(i,j,k)
#  endif
#  ifdef AVERAGES_QUADRATIC
              avgHuon(i,j,k)=fac*avgHuon(i,j,k)
              avgHvom(i,j,k)=fac*avgHvom(i,j,k)
              avgUU(i,j,k)=fac*avgUU(i,j,k)
              avgVV(i,j,k)=fac*avgVV(i,j,k)
#  endif
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          DO j=Jstr,Jend
            DO i=Istr,Iend
              avgUV(i,j,k)=fac*avgUV(i,j,k)
            END DO
          END DO
#  endif
        END DO
!
!  Process fields associated with tracers.
!
        DO itrc=1,NT(ng)
          DO k=1,N(ng)
            DO j=JstrR,JendR
              DO i=IstrR,IendR
                avgt(i,j,k,itrc)=fac*avgt(i,j,k,itrc)
              END DO
            END DO
          END DO
#  ifdef AVERAGES_QUADRATIC
          IF (itrc.le.NAT) THEN
            DO k=1,N(ng)
              DO j=JstrR,JendR
                DO i=IstrR,IendR
                  avgTT(i,j,k,itrc)=fac*avgTT(i,j,k,itrc)
                END DO
                DO i=Istr,Iend
                  avgUT(i,j,k,itrc)=fac*avgUT(i,j,k,itrc)
                  avgHuonT(i,j,k,itrc)=fac*avgHuonT(i,j,k,itrc)
                END DO
              END DO
              DO j=Jstr,Jend
                DO i=IstrR,IendR
                  avgVT(i,j,k,itrc)=fac*avgVT(i,j,k,itrc)
                  avgHvomT(i,j,k,itrc)=fac*avgHvomT(i,j,k,itrc)
                END DO
              END DO
            END DO
          END IF
#  endif
        END DO
!
!  Process fields at W-points.
!
        DO k=0,N(ng)
          DO j=JstrR,JendR
            DO i=IstrR,IendR
              avgw3d(i,j,k)=fac*avgw3d(i,j,k)
              avgwvel(i,j,k)=fac*avgwvel(i,j,k)
#  ifdef AVERAGES_AKV
              avgAKv(i,j,k)=fac*avgAKv(i,j,k)
#  endif
#  ifdef AVERAGES_AKT
              avgAKt(i,j,k)=fac*avgAKt(i,j,k)
#  endif
#  ifdef AVERAGES_AKS
              avgAKs(i,j,k)=fac*avgAKs(i,j,k)
#  endif
            END DO
          END DO
        END DO
# endif
      END IF

      RETURN
      END SUBROUTINE set_avg_tile
#endif
      END MODULE set_avg_mod
