!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#ifdef AVERAGES
                                     ! Write time-averaged
      subroutine wrt_avg             ! fields into averages
                                     ! netCDF file.
!     implicit none
      integer ierr, record, lstr, lvar, lenstr
     &  , start(2), count(2), ibuff(4), nf_fwrite, type
# ifdef SOLVE3D
     &            , itrc
#  ifdef SEDIMENT
     &            , indxWrk
#  endif
# endif
# if defined MPI & !defined PARALLEL_FILES
      include 'mpif.h'
      integer status(MPI_STATUS_SIZE), blank
# endif
# include "param.h"
# include "scalars.h"
# include "averages.h"
# include "ncscrum.h"
# include "grid.h"
# include "ocean2d.h"
# include "ocean3d.h"
# include "mpi_cpl.h"
# include "work.h"
# include "netcdf.inc"
# ifdef SEDIMENT
#  include "sediment.h"
# endif
# ifdef NBQ
#  include "nbq.h"
# endif
# ifdef MASKING
#  define SWITCH *
# else
#  define SWITCH !
# endif
# if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
      if (mynode.gt.0) then
        call MPI_Recv (blank, 1, MPI_INTEGER, mynode-1,
     &                 1, MPI_COMM_WORLD, status, ierr)
      endif
# endif
!
! Create/open averages file; write grid arrays, if so needed.
!
      call def_avg (ncidavg, nrecavg, ierr)
      if (ierr .ne. nf_noerr) goto 99
      lstr=lenstr(avgname)
!                                            !!! WARNING: Here it is
! Set record within the file.                !!! assumed that global
!                                            !!! restart record index
      nrecavg=max(nrecavg,1)                 !!! nrecavg is already
      if (nrpfavg.eq.0) then                 !!! advanced by main.
        record=nrecavg
      else
        record=1+mod(nrecavg-1, nrpfavg)
      endif

!#define CR
CR      write(*,*) 'wrt_avg: Entry  ' MYID
!
! Write out time-averaged variables:
! ----- --- ------------- ----------
!
!     Time step and record indices.
!
      type=filetype_avg

      ibuff(1)=iic
      ibuff(2)=nrecrst
      ibuff(3)=nrechis
      ibuff(4)=nrecavg

      start(1)=1
      start(2)=record
      count(1)=4
      count(2)=1
      ierr=nf_put_vara_int (ncidavg, avgTstep, start, count, ibuff)
      if (ierr .ne. nf_noerr) then
         write(stdout,1) 'time_step', record,ierr MYID
         goto 99                !--> ERROR
      endif
CR      write(*,*) 'wrt_avg:  time  ' MYID
!
!     Averaged time
!
      ierr=nf_put_var1_FTYPE (ncidavg, avgTime, record, time_avg
# ifdef USE_CALENDAR
     &  - origin_date_in_sec
# endif
     &  )
      if (ierr .ne. nf_noerr) then
         lvar=lenstr(vname(1,indxTime))
         write(stdout,1) vname(1,indxTime)(1:lvar), record, ierr
     &        MYID
         goto 99                !--> ERROR
      endif
!
!     Averaged time2
!
      ierr=nf_put_var1_FTYPE (ncidavg, avgTime2, record, time_avg
# ifdef USE_CALENDAR
     &  - origin_date_in_sec
# endif
     &  )
      if (ierr .ne. nf_noerr) then
         lvar=lenstr(vname(1,indxTime2))
         write(stdout,1) vname(1,indxTime2)(1:lvar), record, ierr
     &        MYID
         goto 99                !--> ERROR
      endif
CR      write(*,*) 'wrt_avg:  time  ' MYID

!
!     Barotropic mode variables: free-surface and 2D momentum
!     components in XI-,ETA-directions.
!
      if (wrtavg(indxZ)) then
!$acc update if(compute_on_device) host(zeta_avg)
         work2d=zeta_avg
# ifdef WET_DRY
         do j=0,Mm
            do i=0,Lm
               if (h(i,j) .le. Dcrit(i,j)) then
                 work2d(i,j)=work2d(i,j)+h(i,j)
               endif
            enddo
         enddo
# endif
         call fillvalue2d(work2d,ncidavg,avgZ,vname(1,indxZ),
     &                   record,r2dvar,type)
      endif

      if (wrtavg(indxUb)) then
!$acc update if(compute_on_device) host(ubar_avg)
         work2d=ubar_avg
         call fillvalue2d(work2d,ncidavg,avgUb,vname(1,indxUb),
     &                      record,u2dvar,type)
      endif

      if (wrtavg(indxVb)) then
!$acc update if(compute_on_device) host(vbar_avg)
         work2d=vbar_avg
         call fillvalue2d(work2d,ncidavg,avgVb,vname(1,indxVb),
     &                   record,v2dvar,type)
      endif

# ifdef MORPHODYN
      if (wrtavg(indxHm)) then
!$acc update if(compute_on_device) host(h_avg)
         work2d=h_avg
         call fillvalue2d(work2d,ncidavg,avgHm,vname(1,indxHm),
     &                    record,r2dvar,type)
      endif
# endif
!
!     Write out kinematic bottom stress (N/m2).
!
      if (wrtavg(indxBostr)) then
!$acc update if(compute_on_device) host(bostr_avg)
         work2d=bostr_avg
         call fillvalue2d(work2d,ncidavg,avgBostr,vname(1,indxBostr),
     &        record,r2dvar,type)
      endif

!
!     Write out kinematic U  bottom stress (N/m2).
!
      if (wrtavg(indxBustr)) then
!$acc update if(compute_on_device) host(bustr_avg)
         work2d=bustr_avg
         call fillvalue2d(work2d,ncidavg,avgBustr,vname(1,indxBustr),
     &        record,u2dvar,type)
      endif


!
!     Write out kinematic V  bottom stress (N/m2).
!
      if (wrtavg(indxBvstr)) then
!$acc update if(compute_on_device) host(bvstr_avg)
         work2d=bvstr_avg
         call fillvalue2d(work2d,ncidavg,avgBvstr,vname(1,indxBvstr),
     &        record,v2dvar,type)
      endif



!
!--   Atmospheric forcing : no mask, no fill value
!
!
!     Write out kinematic surface stress (N/m2).
!
      if (wrtavg(indxWstr)) then
!$acc update if(compute_on_device) host(wstr_avg)
         ierr=nf_fwrite(wstr_avg, ncidavg, avgWstr, record, r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxWstr))
            write(stdout,1) vname(1,indxWstr)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg:  Wstr' MYID
      endif
!
!     Write out kinematic U surface stress component (N/m2).
!
      if (wrtavg(indxUWstr)) then
!$acc update if(compute_on_device) host(sustr_avg)
# ifdef STRESS_AT_RHO_POINTS
         ierr=nf_fwrite(sustr_avg, ncidavg, avgUWstr, record, r2dvar)
# else
         ierr=nf_fwrite(sustr_avg, ncidavg, avgUWstr, record, u2dvar)
# endif
        if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxUWstr))
            write(stdout,1) vname(1,indxUWstr)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg: UWstr' MYID
      endif
!
!     Write out kinematic V surface stress component (N/m2).
!
      if (wrtavg(indxVWstr)) then
!$acc update if(compute_on_device) host(svstr_avg)
# ifdef STRESS_AT_RHO_POINTS
         ierr=nf_fwrite(svstr_avg, ncidavg, avgVWstr, record, r2dvar)
# else
         ierr=nf_fwrite(svstr_avg, ncidavg, avgVWstr, record, v2dvar)
# endif
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxVWstr))
            write(stdout,1) vname(1,indxVWstr)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg: VWstr' MYID
      endif
!
!--
!
# ifdef SOLVE3D
!
!     3D momentum components in XI- and ETA-directions.
!
      if (wrtavg(indxU)) then
!$acc update if(compute_on_device) host(u_avg)
         workr=u_avg
         call fillvalue3d(workr,ncidavg,avgU,
     &        vname(1,indxU),record,u3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtU' MYID
      endif
!
      if (wrtavg(indxV)) then
!$acc update if(compute_on_device) host(v_avg)
         workr=v_avg
         call fillvalue3d(workr,ncidavg,avgV,
     &        vname(1,indxV),record,v3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtV' MYID
      endif
!
!     Tracer variables.
!
#  ifdef TRACERS
      do itrc=1,NT
         if (wrtavg(indxV+itrc)) then
!$acc update if(compute_on_device) host(t_avg(:,:,:,itrc))
            workr=t_avg(:,:,:,itrc)
            call fillvalue3d(workr,ncidavg,avgT(itrc),
     &        vname(1,indxV+itrc), record,r3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtT  ' MYID
         endif
      enddo
#  endif
!
!     Density anomaly.
!
      if (wrtavg(indxR)) then
!$acc update if(compute_on_device) host(rho_avg)
         workr=rho_avg+rho0-1000.
         call fillvalue3d(workr,ncidavg,avgR,
     &        vname(1,indxR), record,r3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtRHO' MYID
      endif
#  if defined ANA_VMIX || defined BVF_MIXING \
  || defined LMD_MIXING || defined LMD_SKPP || defined LMD_BKPP \
  || defined GLS_MIXING
!
!     Brunt vaisala frequency.
!
      if (wrtavg(indxbvf)) then
!$acc update if(compute_on_device) host(bvf_avg)
         work=bvf_avg
         call fillvalue3d_w(work,ncidavg,avgbvf,
     &        vname(1,indxbvf), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtbvf' MYID
      endif
#  endif
!
!     Write out S-coordinate omega vertical velocity (m/s).
!
      if (wrtavg(indxO)) then
!$acc update if(compute_on_device) host(omega_avg)
         work=omega_avg
         call fillvalue3d_w(work,ncidavg,avgO,
     &        vname(1,indxO), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtO  ' MYID
      endif
!
!     Write out true vertical velocity (m/s).
!
      if (wrtavg(indxW)) then
!$acc update if(compute_on_device) host(w_avg)
#  ifdef NBQ
         work=w_avg
         call fillvalue3d(work ,ncidavg,avgW,
     &        vname(1,indxW), record,w3dvar,type)
#  else
         workr=w_avg
         call fillvalue3d(workr,ncidavg,avgW,
     &        vname(1,indxW), record,r3dvar,type)
#  endif
CR      write(*,*) 'wrt_avg:  wrtW' MYID
      endif
!
#  if defined SST_SKIN && defined TEMPERATURE
!
!     Write out skin temperature (degC)
!
      if (wrtavg(indxT)) then
!$acc update if(compute_on_device) host(sst_skin_avg)
         work2d=sst_skin_avg
         call fillvalue2d(work2d,ncidavg,avgSST_skin,
     &        vname(1,indxSST_skin), record,r2dvar,type)
      endif
#  endif
!
!
#  ifdef VIS_COEF_3D
!
!     Write out Horizontal viscosity coefficient.
!
      if (wrtavg(indxVisc)) then
!$acc update if(compute_on_device) host(visc3d_avg)
         workr=visc3d_avg
         call fillvalue3d(workr,ncidavg,avgVisc,
     &        vname(1,indxVisc), record,r3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtVisc' MYID
      endif
#  endif
#  ifdef DIF_COEF_3D
!
!     Write out Horizontal diffusivity coefficient.
!
      if (wrtavg(indxDiff)) then
!$acc update if(compute_on_device) host(diff3d_avg)
         workr=diff3d_avg
         call fillvalue3d(workr,ncidavg,avgDiff,
     &        vname(1,indxDiff), record,r3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtDiff' MYID
      endif
#  endif
#  ifdef AVERAGES_K
!
!     Write out vertical viscosity coefficient.
!
      if (wrtavg(indxAkv)) then
!$acc update if(compute_on_device) host(Akv_avg)
         work=Akv_avg
         call fillvalue3d_w(work,ncidavg,avgAkv,
     &        vname(1,indxAkv), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtAkv' MYID
      endif
!
!     Write out vertical diffusion coefficient for potential temperature.
!
#   ifdef TEMPERATURE
      if (wrtavg(indxAkt)) then
!$acc update if(compute_on_device) host(Akt_avg(:,:,:,itemp))
         work=Akt_avg(:,:,:,itemp)
         call fillvalue3d_w(work,ncidavg,avgAkt,
     &        vname(1,indxAkt), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtAkt' MYID
      endif
#   endif /* TEMPERATURE */
#   ifdef SALINITY
!
!     Write out vertical diffusion coefficient for salinity.
!
      if (wrtavg(indxAks)) then
!$acc update if(compute_on_device) host(Akt_avg(:,:,:,isalt))
         work=Akt_avg(:,:,:,isalt)
         call fillvalue3d_w(work,ncidavg,avgAks,
     &        vname(1,indxAks), record,w3dvar,type)
CR       write(*,*) 'wrt_avg:  wrtAks' MYID
      endif
#   endif /* SALINITY */
#  endif  /* AVERAGES_K */

#  if defined LMD_SKPP || defined GLS_MIXING
!
!     Write out depth of planetary boundary layer (m).
!
      if (wrtavg(indxHbl)) then
!$acc update if(compute_on_device) host(hbl_avg)
         work2d=hbl_avg
         call fillvalue2d(work2d,ncidavg,avgHbl,
     &        vname(1,indxHbl), record,r2dvar,type)
CR      write(*,*) 'wrt_avg:  wrtHBL' MYID
      endif
#  endif
#  ifdef LMD_BKPP
!
!     Write out depth of planetary boundary layer (m).
!
      if (wrtavg(indxHbbl)) then
!$acc update if(compute_on_device) host(hbbl_avg)
         work2d=hbbl_avg
         call fillvalue2d(work2d,ncidavg,avgHbbl,
     &        vname(1,indxHbbl), record,r2dvar,type)
      endif
CR      write(*,*) 'wrt_avg:  wrtHbbl' MYID
#  endif

#  ifdef GLS_MIXING
!
!     Write out turbulent kinetic energy.
!
      if (wrtavg(indxTke)) then
!$acc update if(compute_on_device) host(tke_avg)
         work=tke_avg
         call fillvalue3d_w(work,ncidavg,avgTke,
     &        vname(1,indxTke), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtTke' MYID
      endif
!
!     Write out generic length scale
!
      if (wrtavg(indxGls)) then
!$acc update if(compute_on_device) host(gls_avg)
         work=gls_avg
         call fillvalue3d_w(work,ncidavg,avgGls,
     &        vname(1,indxGls), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtGls' MYID
      endif
!
!     Write out vertical mixing length scale
!
      if (wrtavg(indxLsc)) then
!$acc update if(compute_on_device) host(Lscale_avg)
         work=Lscale_avg
         call fillvalue3d_w(work,ncidavg,avgLsc,
     &        vname(1,indxLsc), record,w3dvar,type)
CR      write(*,*) 'wrt_avg:  wrtLsc' MYID
      endif
#  endif

#  ifdef TEMPERATURE
!
!     Write out total heat flux [W/m2]
!
      if (wrtavg(indxShflx)) then
!$acc update if(compute_on_device) host(stflx_avg(:,:,itemp))
         work2d=stflx_avg(:,:,itemp)
     &          SWITCH rmask
!
!        Mult by mask to avoid erroneous data
!        over land if Qcorrection cppkeys applied
!
         ierr=nf_fwrite(work2d, ncidavg, avgShflx, record, r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxShflx))
            write(stdout,1) vname(1,indxShflx)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg:  stflx(:,:,itemp)' MYID
       endif
#  endif
#  ifdef BHFLUX
!
!     Write out bottom heat flux
!
      if (wrtavg(indxBhflx)) then
!$acc update if(compute_on_device) host(btflx_avg(:,:,itemp))
        work2d=btflx_avg(:,:,itemp)
     &         SWITCH rmask
!
!     Mult by mask to avoid erroneous data
!     over land if Qcorrection cppkeys applied
!
        ierr=nf_fwrite(work2d, ncidavg, avgBhflx, record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxBhflx))
          write(stdout,1) vname(1,indxBhflx)(1:lvar), record, ierr
     &         MYID
          goto 99               !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  btflx(:,:,itemp)' MYID
      endif
#  endif
#  ifdef SALINITY
!
!     Write out surface E-P flux [cm/d]
!
      if (wrtavg(indxSwflx)) then
!$acc update if(compute_on_device) host(stflx_avg(:,:,isalt))
         work2d=stflx_avg(:,:,isalt)
     &         SWITCH rmask
!
!        Mult by mask to avoid erroneous data over land
!        if SFlux cppkeys correction applied and also
!        to remove enormous wrong value diciding by eps.
!
         ierr=nf_fwrite(work2d, ncidavg, avgSwflx, record, r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxSwflx))
            write(stdout,1) vname(1,indxSwflx)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg:  stflx(:,:,isalt)' MYID
       endif
#   ifdef BWFLUX
!
!     Write out bottom fresh water flux [cm/d]
!
      if (wrtavg(indxBwflx)) then
!$acc update if(compute_on_device) host(btflx_avg(:,:,isalt))
         work2d=btflx_avg(:,:,isalt)
     &         SWITCH rmask
!
!        Mult by mask to avoid erroneous data over land
!        if SFlux cppkeys correction applied and also
!        to remove enormous wrong value diciding by eps.
!
         ierr=nf_fwrite(work2d, ncidavg, avgBwflx, record, r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxBwflx))
            write(stdout,1) vname(1,indxBwflx)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg:  btflx(:,:,isalt)' MYID
       endif
#   endif
#  endif
!     Write out solar radiation flux [W/m2]
!
      if (wrtavg(indxShflx_rsw)) then
#  ifdef BULK_FLUX
!$acc update if(compute_on_device) host(shflx_rsw_avg)
         ierr=nf_fwrite(shflx_rsw_avg, ncidavg, avgShflx_rsw,
     &                  record, r2dvar)
#  else
!$acc update if(compute_on_device) host(srflx_avg)
         ierr=nf_fwrite(srflx_avg, ncidavg, avgShflx_rsw,
     &                  record, r2dvar)
#  endif
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxShflx_rsw))
            write(stdout,1) vname(1,indxShflx_rsw)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
CR      write(*,*) 'wrt_avg:  srflx' MYID
      endif
#  ifdef BULK_FLUX
      if (wrtavg(indxShflx_rlw)) then
!$acc update if(compute_on_device) host(shflx_rlw_avg)
         ierr=nf_fwrite(shflx_rlw_avg, ncidavg, avgShflx_rlw, record,
     &        r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxShflx_rlw))
            write(stdout,1) vname(1,indxShflx_rlw)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
      endif

      if (wrtavg(indxShflx_lat)) then
!$acc update if(compute_on_device) host(shflx_lat_avg)
         ierr=nf_fwrite(shflx_lat_avg, ncidavg, avgShflx_lat, record,
     &        r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxShflx_lat))
            write(stdout,1) vname(1,indxShflx_lat)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
      endif

      if (wrtavg(indxShflx_sen)) then
!$acc update if(compute_on_device) host(shflx_sen_avg)
         ierr=nf_fwrite(shflx_sen_avg, ncidavg, avgShflx_sen, record,
     &        r2dvar)
         if (ierr .ne. nf_noerr) then
            lvar=lenstr(vname(1,indxShflx_sen))
            write(stdout,1) vname(1,indxShflx_sen)(1:lvar), record, ierr
     &           MYID
            goto 99             !--> ERROR
         endif
      endif
#  endif

#  if defined BIOLOGY && !defined PISCES
!
!     Write out depth of the euphotic layer (m).
!
      if (wrtavg(indxHel)) then
!$acc update if(compute_on_device) host(hel_avg)
         work2d=hel_avg
         call fillvalue2d(work2d,ncidavg,avgHel,
     &        vname(1,indxHel), record,r2dvar,type)
      endif
!
!     Write out Chlorophyll to Carbon ratio (m).
!
#   ifdef BIO_NChlPZD
      if (wrtavg(indxChC)) then
!$acc update if(compute_on_device) host(theta_avg)
         workr=theta_avg
         call fillvalue3d(workr,ncidavg,avgChC,
     &        vname(1,indxChC), record,r3dvar,type)
      endif
#    ifdef OXYGEN
      if (wrtavg(indxU10)) then
!$acc update if(compute_on_device) host(u10_avg)
         work2d=u10_avg
         call fillvalue2d(work2d,ncidavg,avgU10,
     &        vname(1,indxU10), record,r2dvar,type)
      endif
      if (wrtavg(indxKvO2)) then
!$acc update if(compute_on_device) host(Kv_O2_avg)
         work2d=Kv_O2_avg
         call fillvalue2d(work2d,ncidavg,avgKvO2,
     &        vname(1,indxKvO2), record,r2dvar,type)
      endif
      if (wrtavg(indxO2sat)) then
!$acc update if(compute_on_device) host(O2satu_avg)
         work2d=O2satu_avg
         call fillvalue2d(work2d,ncidavg,avgO2sat,
     &        vname(1,indxO2sat), record,r2dvar,type)
      endif
#    endif /* OXYGEN */
!
#   elif defined BIO_BioEBUS
      if (wrtavg(indxAOU)) then
!$acc update if(compute_on_device) host(AOU_avg)
         workr=AOU_avg
         call fillvalue3d(workr,ncidavg,avgAOU,
     &        vname(1,indxAOU), record,r3dvar,type)
      endif
       if (wrtavg(indxWIND10)) then
!$acc update if(compute_on_device) host(wind10_avg)
         work2d=wind10_avg
         call fillvalue2d(work2d,ncidavg,avgwind10,
     &        vname(1,indxWIND10), record,r2dvar,type)
      endif
#   endif
#  endif /* BIOLOGY */
!
#  ifdef SEDIMENT
!
!     Write out sediment bed layer thickness, porosity, volume
!     fraction of size class in sediment bed (2+2*NST b3dgrd variables)
!
      if (wrtavg(indxBTHK)) then
!$acc update if(compute_on_device) host(bed_thick_avg)
         worksed_bed=bed_thick_avg
         call fillvalue3d(worksed_bed,ncidavg,avgSed(2),
     &        vname(1,indxBTHK), record,b3dvar,type)
      endif
!
      if (wrtavg(indxBPOR)) then
!$acc update if(compute_on_device) host(bed_poros_avg)
         worksed_bed=bed_poros_avg
         call fillvalue3d(worksed_bed,ncidavg,avgSed(3),
     &        vname(1,indxBPOR), record,b3dvar,type)
      endif

      do itrc=1,NST
         indxWrk=indxBFRA(1)+itrc-1
         if (wrtavg(indxWrk)) then
!$acc update if(compute_on_device) host(bed_frac_avg(:,:,:,itrc))
            worksed_frac=bed_frac_avg(:,:,:,itrc)
            call fillvalue3d(worksed_frac,ncidavg,avgSed(itrc+3),
     &        vname(1,indxWrk), record,b3dvar,type)
         endif
      enddo

#   ifdef SUSPLOAD
      do itrc=1,NST
         indxWrk=indxDFLX(1)+itrc-1
         if (wrtavg(indxWrk)) then
!$acc update if(compute_on_device) host(settling_flux_avg(:,:,itrc))
            work2d=settling_flux_avg(:,:,itrc)/dt
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3+NST),
     &        vname(1,indxWrk), record,r2dvar,type)
         endif
      enddo
      do itrc=1,NST
         indxWrk=indxEFLX(1)+itrc-1
         if (wrtavg(indxWrk)) then
!$acc update if(compute_on_device) host(ero_flux_avg(:,:,itrc))
            work2d=ero_flux_avg(:,:,itrc)/dt
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3+2*NST),
     &        vname(1,indxWrk),record,r2dvar,type)
         endif
      enddo
#   endif
#   ifdef BEDLOAD
      do itrc=1,NST
         indxWrk=indxBDLU(1)+itrc-1
         if (wrtavg(indxWrk)) then
!$acc update if(compute_on_device) host(bedldu_avg(:,:,itrc))
            work2d=bedldu_avg(:,:,itrc)
#    ifdef SUSPLOAD
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3+3*NST),
     &        vname(1,indxWrk),record,r2dvar,type)
#    else
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3*NST),
     &        vname(1,indxWrk),record,r2dvar,type)
#    endif
         endif
      enddo
      do itrc=1,NST
         indxWrk=indxBDLV(1)+itrc-1
         if (wrtavg(indxWrk)) then
!$acc update if(compute_on_device) host(bedldv_avg(:,:,itrc))
            work2d=bedldv_avg(:,:,itrc)
#    ifdef SUSPLOAD
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3+4*NST),
     &        vname(1,indxWrk),record,r2dvar,type)
#    else
            call fillvalue2d(work2d,ncidavg,avgSed(itrc+3+2*NST),
     &        vname(1,indxWrk),record,r2dvar,type)
#    endif
         endif
      enddo
#   endif
#  endif /* SEDIMENT */
!
# endif /* SOLVE3D */

# ifdef WAVE_IO
      if (wrtavg(indxHRM)) then
!$acc update if(compute_on_device) host(whrm_avg)
        ierr=nf_fwrite (whrm_avg(START_2D_ARRAY), ncidavg, avgWAVE(1),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxHRM))
          write(stdout,1) vname(1,indxHRM)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtHrm' MYID
      endif
      if (wrtavg(indxFRQ)) then
!$acc update if(compute_on_device) host(wfrq_avg)
        ierr=nf_fwrite (wfrq_avg(START_2D_ARRAY), ncidavg, avgWAVE(2),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxFRQ))
          write(stdout,1) vname(1,indxFRQ)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtFrq' MYID
      endif
#  ifdef WKB_WWAVE
      if (wrtavg(indxWAC)) then
!$acc update if(compute_on_device) host(wac_avg)
        ierr=nf_fwrite (wac_avg(START_2D_ARRAY), ncidavg, avgWAVE(3),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxWAC))
          write(stdout,1) vname(1,indxWAC)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtWac' MYID
      endif
#  endif
      if (wrtavg(indxWKX)) then
!$acc update if(compute_on_device) host(wkx_avg)
        ierr=nf_fwrite (wkx_avg(START_2D_ARRAY), ncidavg, avgWAVE(4),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxWKX))
          write(stdout,1) vname(1,indxWKX)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtWkx' MYID
      endif
      if (wrtavg(indxWKE)) then
!$acc update if(compute_on_device) host(wke_avg)
        ierr=nf_fwrite (wke_avg(START_2D_ARRAY), ncidavg, avgWAVE(5),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxWKE))
          write(stdout,1) vname(1,indxWKE)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtWke' MYID
      endif
      if (wrtavg(indxEPB)) then
!$acc update if(compute_on_device) host(wepb_avg)
        ierr=nf_fwrite (wepb_avg(START_2D_ARRAY), ncidavg, avgWAVE(6),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxEPB))
          write(stdout,1) vname(1,indxEPB)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtEpb' MYID
      endif
      if (wrtavg(indxEPD)) then
!$acc update if(compute_on_device) host(wepd_avg)
        ierr=nf_fwrite (wepd_avg(START_2D_ARRAY), ncidavg, avgWAVE(7),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxEPD))
          write(stdout,1) vname(1,indxEPD)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtEpd' MYID
      endif
#  ifdef WAVE_ROLLER
#   ifdef WKB_WWAVE
      if (wrtavg(indxWAR)) then
!$acc update if(compute_on_device) host(war_avg)
        ierr=nf_fwrite (war_avg(START_2D_ARRAY), ncidavg, avgWAVE(8),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxWAR))
          write(stdout,1) vname(1,indxWAR)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtWar' MYID
      endif
#   endif
      if (wrtavg(indxEPR)) then
!$acc update if(compute_on_device) host(wepr_avg)
        ierr=nf_fwrite (wepr_avg(START_2D_ARRAY), ncidavg, avgWAVE(9),
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxEPR))
          write(stdout,1) vname(1,indxEPR)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtEpr' MYID
      endif
#  endif
# endif

# ifdef MRL_WCI
      if (wrtavg(indxSUP)) then
!$acc update if(compute_on_device) host(sup_avg)
        ierr=nf_fwrite (sup_avg(START_2D_ARRAY), ncidavg, avgSUP,
     &                                            record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxSUP))
          write(stdout,1) vname(1,indxSUP)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtSup' MYID
      endif
      if (wrtavg(indxUST2D)) then
!$acc update if(compute_on_device) host(ust2d_avg)
        ierr=nf_fwrite (ust2d_avg(START_2D_ARRAY), ncidavg, avgUST2D,
     &                                                record, u2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxUST2D))
          write(stdout,1) vname(1,indxUST2D)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtUst2D' MYID
      endif
      if (wrtavg(indxVST2D)) then
!$acc update if(compute_on_device) host(vst2d_avg)
        ierr=nf_fwrite (vst2d_avg(START_2D_ARRAY), ncidavg, avgVST2D,
     &                                                record, v2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxVST2D))
          write(stdout,1) vname(1,indxVST2D)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtVst2D' MYID
      endif
#  ifdef SOLVE3D
      if (wrtavg(indxUST)) then
!$acc update if(compute_on_device) host(ust_avg)
        ierr=nf_fwrite (ust_avg(START_2D_ARRAY,1), ncidavg, avgUST,
     &                                              record, u3dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxUST))
          write(stdout,1) vname(1,indxUST)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtUst' MYID
      endif
      if (wrtavg(indxVST)) then
!$acc update if(compute_on_device) host(vst_avg)
        ierr=nf_fwrite (vst_avg(START_2D_ARRAY,1), ncidavg, avgVST,
     &                                              record, v3dvar)
        if (ierr .ne. nf_noerr) then
        lvar=lenstr(vname(1,indxVST))
          write(stdout,1) vname(1,indxVST)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtVst' MYID
      endif
      if (wrtavg(indxWST)) then
!$acc update if(compute_on_device) host(wst_avg)
        ierr=nf_fwrite (wst_avg(START_2D_ARRAY,1), ncidavg, avgWST,
     &                                              record, r3dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxWST))
          write(stdout,1) vname(1,indxWST)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtWst' MYID
      endif
      if (wrtavg(indxAkb)) then
!$acc update if(compute_on_device) host(akb_avg)
        ierr=nf_fwrite (akb_avg(START_2D_ARRAY,0), ncidavg,
     &                                avgAkb, record, w3dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxAkb))
          write(stdout,1) vname(1,indxAkb)(1:lvar), record, ierr,
     &                    nf_strerror(ierr) MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtAkb' MYID
      endif
      if (wrtavg(indxAkw)) then
!$acc update if(compute_on_device) host(akw_avg)
        ierr=nf_fwrite (akw_avg(START_2D_ARRAY,0), ncidavg,
     &                                avgAkw, record, w3dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxAkw))
          write(stdout,1) vname(1,indxAkw)(1:lvar), record, ierr,
     &                    nf_strerror(ierr) MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtAkw' MYID
      endif
      if (wrtavg(indxKVF)) then
!$acc update if(compute_on_device) host(kvf_avg)
        ierr=nf_fwrite (kvf_avg(START_2D_ARRAY,1), ncidavg, avgKVF,
     &                                              record, r3dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxKVF))
          write(stdout,1) vname(1,indxKVF)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtKvf' MYID
      endif
      if (wrtavg(indxCALP)) then
!$acc update if(compute_on_device) host(calp_avg)
        ierr=nf_fwrite (calp_avg(START_2D_ARRAY), ncidavg, avgCALP,
     &                                              record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxCALP))
          write(stdout,1) vname(1,indxCALP)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtCalp' MYID
      endif
      if (wrtavg(indxKAPS)) then
!$acc update if(compute_on_device) host(kaps_avg)
        ierr=nf_fwrite (kaps_avg(START_2D_ARRAY), ncidavg, avgKAPS,
     &                                              record, r2dvar)
        if (ierr .ne. nf_noerr) then
          lvar=lenstr(vname(1,indxKAPS))
          write(stdout,1) vname(1,indxKAPS)(1:lvar), record, ierr
     &                    MYID
          goto 99                                         !--> ERROR
        endif
CR      write(*,*) 'wrt_avg:  wrtKaps' MYID
      endif
#  endif  /* SOLVE3D */
# endif   /* MRL_WCI */
!
 1    format(/' WRT_AVG - ERROR while writing variable(',1x,a,1x,
     &     ')into averages file.',/,11x,'Time record:',
     &     i6,3x,'netCDF error code',i4,3x,a,i4)
      goto 100
 99   may_day_flag=3
 100  continue

!
!     Synchronize netCDF file to disk to allow other processes
!     to access data immediately after it is written.
!
# if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
      ierr=nf_close(ncidavg)
      if (nrpfavg.gt.0 .and. record.ge.nrpfavg) ncidavg=-1
# else
      if (nrpfavg.gt.0 .and. record.ge.nrpfavg) then
         ierr=nf_close(ncidavg)
         ncidavg=-1
      else
         ierr=nf_sync(ncidavg)
      endif
# endif
      if (ierr .eq. nf_noerr) then
      MPI_master_only write(stdout,'(6x,A,2(A,I4,1x),A,I3)')
     &        'WRT_AVG -- wrote ',
     &        'averaged fields into time record =', record, '/',
     &        nrecavg  MYID
      else
       MPI_master_only write(stdout,'(/1x,2A/)')
     &        'WRT_AVG ERROR: Cannot ',
     &        'synchronize/close averages netCDF file.'
         may_day_flag=3
      endif

# if defined MPI & !defined PARALLEL_FILES & !defined NC4PAR
      if (mynode .lt. NNODES-1) then
         call MPI_Send (blank, 1, MPI_INTEGER, mynode+1,
     &        1, MPI_COMM_WORLD,  ierr)
      endif
# endif
      return
      end
#else
      subroutine wrt_avg_empty
      end
#endif /* AVERAGES */
