module m_archive_select
      !from m_fields_to_plot.F90 in MSCPROGS/src/Hyc2proj
#if defined (ARCHIVE_SELECT)

      implicit none

      private

      type fields
         character(len=8) fextract
         integer          laya
         integer          layb
         logical          option
         logical          vecflag
         character(len=8) vecpost
      end type fields

      type(fields)   :: archv_fields(1000)
      integer        :: n_archv_fields

      public :: archiv_init,archiv

      contains

      subroutine archiv_init()

      use mod_xc
      implicit none

      !set archv_fields,n_archv_fields for use by archiv routine
      call extract_fields(archv_fields,n_archv_fields)
      if (mnproc==1) print*,'archiv_init: setting fields to extract'
      end subroutine archiv_init

      subroutine extract_fields(fld,nfld)!,ftype)
      !determine which fields to dump to archv files
      !- save disk space by not dumping all variables to archive files
      !- could do the same for daily/weekly average files perhaps
      use mod_xc

      implicit none

      integer, parameter            :: iversion=4
      type(fields), intent(out)     :: fld(1000)
      integer, intent(out)          :: nfld
      !character(len=*),intent(in)   :: ftype
      character(len=*),parameter    :: ftype = 'archv'
      integer ndim,i,ivers,lay1,lay2,ifld
      character(len=8) char8, uname, vname, pre
      character(len=80) extrfile
      logical option,ex, vecflag
      
      if (mnproc==1) print*,'extract_fields: setting fields to extract'
      if (trim(ftype)=='archv') then
         extrfile='archv.extract'
      !else if (trim(ftype)=='nersc_daily') then
      !   extrfile='daily.extract'
      !else if (trim(ftype)=='nersc_weekly') then
      !   extrfile='weekly.extract'
      !else if (trim(ftype)=='archv_wav') then
      !   extrfile='archv_wav.extract'
      else
          if(mnproc==1)print *,'extract_fields> Unknown filetype '
     &                           //trim(ftype)
         call xcstop('(extract_fields)')
         stop '(extract_fields)'
      end if

      if (mnproc==1) print*,'Using:',trim(extrfile)
      inquire(exist=ex,file=trim(extrfile)) 

      if (.not.ex) then
         if(mnproc==1) print *,'Can not find extract file '
     &      //trim(extrfile)
         call xcstop('(extract_fields)')
         stop '(extract_fields)'
      end if

      open(21,file=trim(extrfile),status='old')
      read(21,*) ivers  
      if (ivers /= iversion) then
         close(21)
         if (mnproc==1) then
            print *,'Version of extract file should be:',iversion
            print *,'From 1 to 2 just add a F or T in column 3,5,7'//
     &         'of line 3'
            print *,'This makes it possible to'//
     &         'plot solution on a sphere,'
            print *,'to calculate and plot rotated velocities,'//
     &         'and to switch of'
            print *,'the index velocities. Use the line:'
            print *,'2 T T T     # 2-D, sphere, rotate,'//
     &         'index velocities'
            print *
            print *,'From version 3 to 4, field names are'//
     &         'b characters long, '
            print *,'starting on first column'
         end if
         call xcstop('(extract_fields)')
         stop 'wrong version of extract file'
      endif
      read(21,*) ! ignored

      read(21,*) ! ignored
      read(21,*)i         ! not used here
      read(21,*) lay1,lay2 ! not used here

      if (mnproc==1)print *,'The following variables'//
     &   '(and corresponding layers) will be printed:'

      nfld=0
      do i=1,1000
         !read(21,'(a8,t9,i2,t12,i2,t15,l1)',end=888) char8,lay1,lay2,option  
         read(21,*,end=888) char8,lay1,lay2,option  
         if (option) then
            nfld=nfld+1
            fld(nfld)%fextract=adjustl(char8)
            fld(nfld)%laya=lay1
            fld(nfld)%layb=lay2
            fld(nfld)%option=option
            !write(6,'(a)',advance='no') fld(nfld)%fextract//" "
            if (mnproc==1)print*,fld(nfld)%fextract
         endif
      enddo
      !write(6,'(a)',advance='yes') ' '
888   close(21)



      ! Special case (and rules) for vectors -
      ! 1) No scalar variables should begin with u,v,taux or tauy
      ! 2) v(tauy)-component must come immediately after u(taux)-component
      fld(1:nfld)%vecflag=.false.
      do ifld=1,nfld-1
         if (fld(ifld)%fextract(1:1)=='u' .or.
     &       fld(ifld)%fextract(1:4)=='taux') then
         if (ifld==nfld) then
            if (mnproc==1) then
               print *,'For vectors, the v-component must come'//
     &            'immediately after the u-component;'
               print *,'for stresses, the tauy-component must come'//
     &            'immediately after the taux-component'
               print *,fld(ifld)
               print *,'- this is the last field'
            end if
            call xcstop('(extract_fields)')
         else if (fld(ifld)%fextract(1:1)=='u' .and.
     &      (fld(ifld+1)%fextract(1:1)/='v' .or.
     &       fld(ifld)%fextract(2:8)/=fld(ifld+1)%fextract(2:8) )) then
            if (mnproc==1) then
               print *,'For vectors, the v-component must come'//
     &            'immediately after the u-component - 1'
               print *,fld(ifld)
               print *,fld(ifld+1)
            end if
            call xcstop('(extract_fields)')
         else if (fld(ifld)%fextract(1:4)=='taux' .and.
     &      ( fld(ifld)%fextract(5:8)/=fld(ifld+1)%fextract(5:8) .or.
     &        fld(ifld+1)%fextract(1:4)/='tauy' )) then
            if (mnproc==1) then
               print *,'For stresses, the tauy-component must come'//
     &            'immediately after the taux-component'
               print *,fld(ifld)%fextract(1:4)
               print *,fld(ifld  )%fextract
               print *,fld(ifld+1)%fextract
            end if
            call xcstop('(extract_fields)')
         end if
         end if

         
         uname=fld(ifld  )%fextract
         vname=fld(ifld+1)%fextract ; 
         ! Flag indicating vector component
         vecflag=
     &      (uname(1:1) == 'u'   .and.vname(1:1) == 'v') .or.
     &      (uname(1:4) == 'taux'.and.vname(1:4) == 'tauy')  
         ! Now test that the rest of their names are equal
         if (vecflag) then
            if (uname(1:4) == 'taux') then
               vecflag=uname(5:8) == vname(5:8)
               pre='tau'//uname(5:8)
            elseif (uname(1:1) == 'u' .and. len_trim(uname)>1) then
               vecflag=uname(2:8) == vname(2:8)
               pre=uname(2:8)
            end if
         end if

         if (vecflag) then
            fld(ifld  )%vecflag = .true.
            fld(ifld+1)%vecflag = .true.
            fld(ifld  )%vecpost = pre
            fld(ifld+1)%vecpost = pre
            fld(ifld+1)%option = .false. ! This switches off scalar processing of 2nd comp
            if (mnproc==1) then
               write(6,'(a)') 'Found vector pair: '//uname//vname
            end if
         end if

      end do

      end subroutine extract_fields

      subroutine archiv(n, kkout, iyear,iday,ihour, intvl,dtime0)
!     Routine for dumping instantaneous fields into binary.
!     Set frequency in blkdat.input (dsurfq)      
!     TW: Use archv.extract file to set archv_fields
      use mod_xc  ! HYCOM communication interface
      use mod_za  ! HYCOM I/O interface
#if defined (NERSC_VERSION)
      use mod_year_info
      use mod_hycom_nersc
      use mod_common_ice
      use m_depthave
#if defined (WAVES)
      use mod_common_wavesice, only: dfloe,swh_out,mwp_out,mwd_out
#endif
#endif
#if defined (THERM_DIAG)
      ! diagnostics from m_thermf_nersc
      use m_thermf_nersc, only: flx_cool,flx_ow,flx_itop                &
     &   ,ficem_therm_init,hicem_therm_init                             &
     &   ,hsnwm_therm_init,hsnwm_therm_int
#endif
      use mod_forcing_nersc, only:
     &   clouds
     &  ,uwind
     &  ,vwind
     &  ,slp
     &  ,dewtmp
     &  ,relhum
      !NB airtmp is in common_blocks.h
      implicit none
c
      include 'common_blocks.h'
c
      !inputs
      integer     :: n, kkout, iyear,iday,ihour
      character   :: intvl*3
      real*8      :: dtime0
      !!
      real  :: dtime1
      real  :: sssc,sstc
c
      include 'stmt_fns.h'
c
c --- write an archive file.
c
      character(len=80) :: flnm_arc
      type(year_info) tt0,tt
      character*80 cformat
      character*6  ctime,ctime0
      character*8  cdate,cdate0
      character*2  cdm0,cmin0,csec0
      character*2  cdm ,cmin ,csec
      integer      imin,isec
      integer      year0,mon0,dd0,year1,mon1,dd1
      integer      i,j,k,ktr,l,ldot,nop,nopa
      real         coord,xmin,xmax

      integer  :: ifld

      integer,parameter :: uvsurf_opt  = 2   !1: avg over upper 'surf_depth' meters;
                                             !2: 1st layer
                                             !TODO link to hyc2proj? (opt 2) read dp00 from blkdat.input?

      real,parameter    :: surf_depth  = 3.0 !depth to average over to get surface current
                                             !TODO make this a parameter in an infile?
                                             !TODO how can hyc2proj get this value?

#if defined (NERSC_VERSION)
      real fld(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)
#if defined (WAVES)
      real out_mask(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)
#endif
      character*3 cmem
C --- JB if nersc version iday=iday-1
      iday = iday-1
#endif

      if (time<dtime0) then
         ! time not set yet - just use start time of run
         dtime1   = real(dtime0)
      else
         ! time of dump is set in common_blocks.h
         dtime1   = time
      end if

      call juliantodate(floor(dtime0),year0,mon0,dd0,1900,1,1)
      call year_day(dtime0,year0,tt0,yrflag,.false.)    !refyear is current year - be careful with idd,cdd
      write(cdm0 ,'(i2.2)') 1+tt0%idm                   !add 1 to day in month
      write(cmin0,'(i2.2)') tt0%iss/60                  !minutes
      write(csec0,'(i2.2)') mod(tt0%iss,60)             !seconds
      cdate0   = tt0%cyy//tt0%cmm//cdm0
      ctime0   = tt0%chh//cmin0//csec0

      call juliantodate(floor(dtime1),year1,mon1,dd1,1900,1,1)
      call year_day(dtime1  ,year1,tt ,yrflag,.false.)    !refyear is current year
      write(cdm ,'(i2.2)') 1+tt%idm                     !add 1 to day in month
      write(cmin,'(i2.2)') tt%iss/60                    !minutes
      write(csec,'(i2.2)') mod(tt%iss,60)               !seconds
      cdate = tt%cyy//tt%cmm//cdm
      ctime = tt%chh//cmin//csec

      isec=nint(((dtime1-floor(dtime1))*24.0d0 - tt%ihh)*3600.0d0 )
      imin=isec/60
      isec=mod(isec,60)
      write(ctime,'(i2.2,i2.2,i2.2)') tt%ihh,imin,isec
c
#if defined (NERSC_VERSION)
      !!TW [20150401]: changing name of archv to [eg] TP4archv.[YYYY]_[DDD]_[HHMMSS]Z.[ab]
      if(ensflag) then 
         write(cmem,'(i3.3)') imem
         flnm_arc = rungen//'archv.'//tt%cyy//'_'
     &                //tt%cdd//'_'//ctime//'Z_mem'//cmem
      else
         flnm_arc = rungen//'archv.'//tt%cyy//'_'
     &                //tt%cdd//'_'//ctime//'Z'
      endif
      if (mnproc==1) then
         print*,'archiv (ARCHIVE_SELECT): Making archive files:'
         print*,trim(flnm_arc)//'.[ab]'
      end if
#else
      ldot = index(flnmarc,'.',back=.true.)
      if     (ldot.eq.0) then
        if     (mnproc.eq.1) then
        write (lp,*) 'need decimal point in flnmarc'
        write (lp,*) 'flnmarc = ',trim(flnmarc)
        endif
        call xcstop('(flnmarc)')
               stop '(flnmarc)'
      endif
      ldot = min(ldot,len(flnmarc)-11)  !need 11 characters for archive date
c
      if     ((kkout.eq.1 .and. dsurfq.ge.1.0/24.0) .or.
     &        (kkout.gt.1 .and. diagfq.ge.1.0/24.0)     ) then
c ---   indicate the archive date
        write(flnmarc(ldot+1:ldot+11),'(i4.4,a1,i3.3,a1,i2.2)') 
     &   iyear,'_',iday,'_',ihour
        ldot=ldot+11
      else
c ---   indicate the archive time step
        write(flnmarc(ldot+1:ldot+11),'(i11.11)') nstep
        ldot=ldot+11
      endif
      if (mnproc==1) then
         print*,'archiv (ARCHIVE_SELECT): Making archive files:'
         print*,trim(flnmarc)//'.[ab]'
      end if
#endif
      nopa=13
      nop =13+uoff

      if (mnproc==1) then
         print*,'-----------------------------:'
         print*,'m_archive_select: field maximums:'
         print*,'uice',maxval(iceu)
         print*,'vice',maxval(icev)
         print*,'u',maxval(u(:,:,1,1))
         print*,'v',maxval(v(:,:,1,1))
         print*,'u 2',maxval(u(:,:,1,2))
         print*,'v 2',maxval(v(:,:,1,2))
         print*,'-----------------------------:'
      end if
c
c --- no .[ab] files for 1-D cases (<=6x6) or for dsur1p surface cases.
c
      if     (max(itdm,jtdm).gt.6 .and.
     &        .not.(dsur1p .and. kkout.eq.1)) then  !not 1-D output
c
#if defined (NERSC_VERSION)
      call zaiopf(trim(flnm_arc)//'.a', 'replace', nopa) !TW: new name of files
      if     (mnproc.eq.1) then
      open (unit=nop,file=trim(flnm_arc)//'.b',status='replace') !uoff+13
      write(nop,116) ctitle,iversn,iexpt,yrflag,itdm,jtdm
     &   ,floor(dtime0),cdate0,ctime0 
     &   ,floor(dtime1),cdate,ctime 
      call flush(nop)
      endif !1st tile
#else
      call zaiopf(flnmarc(1:ldot)//'.a', 'replace', nopa)
      if     (mnproc.eq.1) then
      open (unit=nop,file=flnmarc(1:ldot)//'.b',status='replace') !uoff+13
      write(nop,116) ctitle,iversn,iexpt,yrflag,itdm,jtdm
      call flush(nop)
      endif !1st tile
#endif

#if defined (NERSC_VERSION)
      ! new format with extra info in .b files
 116  format (a80/a80/a80/a80/
     & i5,7x, '''iversn'' = hycom version number x10'/
     & i5,7x, '''iexpt '' = experiment number x10'/
     & i5,7x, '''yrflag'' = days in year flag'/
     & i5,7x, '''idm   '' = longitudinal array size'/
     & i5,7x, '''jdm   '' = latitudinal  array size'/
     & i7,5x, '''dtime0'' = model start day'/
     & a ,4x, '''date0 '' = start date'/
     & a ,6x, '''time0 '' = start time (HHMMSS)'/
     & i7,5x, '''dtime '' = model dump day'/
     & a ,4x, '''date  '' = dump  date'/
     & a ,6x, '''time  '' = dump  time (HHMMSS)'/
     & 'field       time step  model day',
     & '  k  dens        min              max')
#else
      ! old format
 116  format (a80/a80/a80/a80/
     & i5,4x, '''iversn'' = hycom version number x10'/
     & i5,4x, '''iexpt '' = experiment number x10'/
     & i5,4x, '''yrflag'' = days in year flag'/
     & i5,4x, '''idm   '' = longitudinal array size'/
     & i5,4x, '''jdm   '' = latitudinal  array size'/
     & 'field       time step  model day',
     & '  k  dens        min              max')

#endif

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !TW: loop over fields to dump
      if (mnproc==1) print*,'archiv (ARCHIVE_SELECT): '//
     &                        'starting loop over fields to dump'

      do ifld=1,n_archv_fields
c
c --- surface fields
c
      coord=0.
c
      if (mnproc==1) print*,'TW:',archv_fields(ifld)%fextract

      if (trim(archv_fields(ifld)%fextract)=='montg1') then
         call zaiowr(montg1,ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
c ---    identify the equation of state on the first record
         write (nop,117) 'montg1  ',nstep,dtime1,sigver,thbase,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if

      if (trim(archv_fields(ifld)%fextract)=='srfhgt') then
         call zaiowr(srfhgt,ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'srfhgt  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if

      if (trim(archv_fields(ifld)%fextract)=='srftmp') then
         fld   = temp(:,:,1,n)
         call zaiowr(fld,ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'srftmp  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if

      if (trim(archv_fields(ifld)%fextract)=='srfsal') then
         fld   = saln(:,:,1,n)
         call zaiowr(fld,ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'srfsal  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if

      if     (sshflg.ne.0) then
      if (trim(archv_fields(ifld)%fextract)=='steric') then
c ---   write out steric SSH.
        call zaiowr(steric,ip,.true.,
     &              xmin,xmax, nopa, .false.)
        if     (mnproc.eq.1) then
        write (nop,117) 'steric  ',nstep,dtime1,0,coord,xmin,xmax
        call flush(nop)
        endif !1st tile
      endif
      endif !sshflg
c
      if (trim(archv_fields(ifld)%fextract)=='surflx') then
         call zaiowr(surflx,ip,.true., xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'surflx  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif

      if (trim(archv_fields(ifld)%fextract)=='salflx') then
         call zaiowr(salflx,ip,.true., xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'salflx  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif
c
      if (trim(archv_fields(ifld)%fextract)=='bl_dpth') then
         call zaiowr(dpbl,ip,.true., xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'bl_dpth ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif

      if (trim(archv_fields(ifld)%fextract)=='mix_dpth') then
         call zaiowr(dpmixl(1-nbdy,1-nbdy,n),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'mix_dpth',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif
      if     (iceflg.ne.0) then
         if (trim(archv_fields(ifld)%fextract)=='covice') then
            call zaiowr(covice,ip,.true., xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 'covice  ',nstep,dtime1,0,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         endif
         if (trim(archv_fields(ifld)%fextract)=='thkice') then
            call zaiowr(thkice,ip,.true., xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 'thkice  ',nstep,dtime1,0,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         endif
         if (trim(archv_fields(ifld)%fextract)=='temice') then
            call zaiowr(temice,ip,.true., xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 'temice  ',nstep,dtime1,0,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         endif
      endif  !write ice fields
c
c --- depth averaged fields
c
      if (trim(archv_fields(ifld)%fextract)=='u_btrop') then
         call zaiowr(ubavg(1-nbdy,1-nbdy,n),iu,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'u_btrop ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif

      if (trim(archv_fields(ifld)%fextract)=='v_btrop') then
         call zaiowr(vbavg(1-nbdy,1-nbdy,n),iv,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'v_btrop ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      endif
c
c --- layer loop.
c
      do 75 k=1,kkout
      coord=sigma(k)

      if (trim(archv_fields(ifld)%fextract)=='u-vel') then
         call zaiowr(u(1-nbdy,1-nbdy,k,n),iu,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'u-vel.  ',nstep,dtime1,k,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if
      if (trim(archv_fields(ifld)%fextract)=='v-vel') then
         call zaiowr(v(1-nbdy,1-nbdy,k,n),iv,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'v-vel.  ',nstep,dtime1,k,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if
      if (trim(archv_fields(ifld)%fextract)=='thknss') then
         call zaiowr(dp(1-nbdy,1-nbdy,k,n),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'thknss  ',nstep,dtime1,k,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if
      if (trim(archv_fields(ifld)%fextract)=='temp') then
         call zaiowr(temp(1-nbdy,1-nbdy,k,n),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'temp    ',nstep,dtime1,k,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if
      if (trim(archv_fields(ifld)%fextract)=='salin') then
         call zaiowr(saln(1-nbdy,1-nbdy,k,n),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'salin   ',nstep,dtime1,k,coord,xmin,xmax
         call flush(nop)
         endif !1st tile
      end if
      if (trim(archv_fields(ifld)%fextract)=='tracer') then
         do ktr= 1,ntracr
           call zaiowr(tracer(1-nbdy,1-nbdy,k,n,ktr),ip,.true.,
     &                 xmin,xmax, nopa, .false.)
           if     (mnproc.eq.1) then
           write (nop,117) 'tracer  ',nstep,dtime1,k,coord,xmin,xmax
           call flush(nop)
           endif !1st tile
         enddo !ktr
      end if
      if     (difout) then
         if (trim(archv_fields(ifld)%fextract)=='viscty') then
            call zaiowr(vcty(1-nbdy,1-nbdy,k+1),ip,.true.,
     &                  xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 'viscty  ',nstep,dtime1,k,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         end if
         if (trim(archv_fields(ifld)%fextract)=='t-diff') then
            call zaiowr(dift(1-nbdy,1-nbdy,k+1),ip,.true.,
     &                  xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 't-diff  ',nstep,dtime1,k,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         end if
         if (trim(archv_fields(ifld)%fextract)=='s-diff') then
            call zaiowr(difs(1-nbdy,1-nbdy,k+1),ip,.true.,
     &                  xmin,xmax, nopa, .false.)
            if     (mnproc.eq.1) then
            write (nop,117) 's-diff  ',nstep,dtime1,k,coord,xmin,xmax
            call flush(nop)
            endif !1st tile
         end if
      endif


 75   continue !end layer loop

#if defined (NERSC_VERSION)
C --- Additional NERSC 2D fields

C --- Only save sea ice variables as surface archive fields 
C
      if (trim(archv_fields(ifld)%fextract)=='uice') then
         call zaiowr(iceu(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'uice    ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='vice') then
C
         call zaiowr(icev(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'vice    ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='ficem') then
         call zaiowr(ficem(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'ficem   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='hicem') then
         call zaiowr(hicem(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'hicem   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='hsnwm') then
         call zaiowr(hsnwm(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'hsnwm   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if

!     !!wind velocities
!     if (trim(archv_fields(ifld)%fextract)=='uwind') then
!        call zaiowr(uwind(1-nbdy,1-nbdy),ip,.true.,
!    &               xmin,xmax, nopa, .false.)
!        if     (mnproc.eq.1) then
!        write (nop,117) 'uwind   ',nstep,dtime1,0,coord,xmin,xmax
!        call flush(nop)
!        endif
!     end if
!     if (trim(archv_fields(ifld)%fextract)=='vwind') then
!        call zaiowr(vwind(1-nbdy,1-nbdy),ip,.true.,
!    &               xmin,xmax, nopa, .false.)
!        if     (mnproc.eq.1) then
!        write (nop,117) 'vwind   ',nstep,dtime1,0,coord,xmin,xmax
!        call flush(nop)
!        endif
!     end if

      !!surface currents
      if (trim(archv_fields(ifld)%fextract)=='usurf') then
         if (uvsurf_opt.eq.1) then
            !!average baroclinic currents over surf_depth
            fld   = depthave(u(:,:,1:kdm,n),dpu(:,:,1:kdm,n),
     &                        surf_depth,'u')
         else
            !!just use 1st layer - harder to find out the depth though
            fld   = u(:,:,1,n)
         end if

         !!add barotropic currents (const over depth)
         fld   = fld+ubavg(:,:,n)
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'usurf   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='vsurf') then
         if (uvsurf_opt.eq.1) then
            !!average baroclinic currents over surf_depth
            fld   = depthave(v(:,:,1:kdm,n),dpv(:,:,1:kdm,n),
     &                        surf_depth,'v')
         else
            !!just use 1st layer
            fld   = depthave(v(:,:,1:kdm,n),dpv(:,:,1:kdm,n),
     &                        surf_depth,'v')
         end if

         !!add barotropic currents (const over depth)
         fld   = fld+vbavg(:,:,n)
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'vsurf   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if

#if defined (THERM_DIAG)
      if (trim(archv_fields(ifld)%fextract)=='flx_itop') then
         !flux into ice from snow/atm above
         call zaiowr(flx_itop(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'flx_itop',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='flx_cool') then
         !flux to cool mixed (top) layer to freezing
         call zaiowr(flx_cool(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'flx_cool',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='flx_ow') then
         !flux into open water from atmosphere
         call zaiowr(flx_ow(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'flx_ow  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='hi_old') then
         !ice thickness at start of therm
         call zaiowr(hicem_therm_init(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'hi_old  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='fi_old') then
         !ice fraction at start of therm
         call zaiowr(ficem_therm_init(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'fi_old  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='hs_old') then
         !snow thickness at start of therm
         call zaiowr(hsnwm_therm_init(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'hs_old  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='hs_int') then
         !snow thickness at start of ice therm
         !(after snow melt/precip)
         call zaiowr(hsnwm_therm_int(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'hs_int  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='dfice') then
         !ice fraction change in thermo
         call zaiowr(delta_ficem(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'dfice   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='dhice') then
         !ice thickness change in thermo
         call zaiowr(delta_hicem(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'dhice   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
#endif

      if (trim(archv_fields(ifld)%fextract)=='taux') then
         !wind stress taux - treat like airtmp
         fld = taux(:,:,l0)*w0+taux(:,:,l1)*w1+
     &         taux(:,:,l2)*w2+taux(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'taux    ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if
      if (trim(archv_fields(ifld)%fextract)=='tauy') then
         !wind stress tauy - treat like airtmp
         fld = tauy(:,:,l0)*w0+tauy(:,:,l1)*w1+
     &         tauy(:,:,l2)*w2+tauy(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'tauy    ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      end if

C
#if defined(WAVES)
      !! Putting waves variables here may not be optimal as
      !! archiv is called after the waves run
      !! - ie 3/6h after current time so could be problem with synchronisation
      !!
      !! archv_wav is called before the waves run so has waves variables
      !! at current time (updated in the wave mask by the new external
      !! wave forcings
      if (trim(archv_fields(ifld)%fextract)=='dfloe') then
         call zaiowr(dfloe(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'dfloe   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      !!mask wave outputs if swh<1cm
      fld      = swh_out
      out_mask = 0.
      where (fld>.01) out_mask = 1.

      if (trim(archv_fields(ifld)%fextract)=='swh') then
         fld   = out_mask*swh_out
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'swh     ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='mwp') then
         fld   = out_mask*mwp_out
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'mwp     ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='mwd') then
         fld   = out_mask*mwd_out
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'mwd     ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
C
#endif/*WAVES*/
C
#if defined(ICE_DYN_DIAG)
C
      if (trim(archv_fields(ifld)%fextract)=='stressp') then
         call zaiowr(stressp(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'stressp ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='stressm') then
         call zaiowr(stressm(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'stressm ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='stress12') then
         call zaiowr(stress12(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'stress12',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='pice') then
         call zaiowr(pice(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'pice    ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='strainI') then
         call zaiowr(strainI(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'strainI ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='strainII') then
         call zaiowr(strainII(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'strainII',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
#endif/*ICE_DYN_DIAG*/
C
C
      if (trim(archv_fields(ifld)%fextract)=='sswflx') then
         call zaiowr(sswflx(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'sswflx  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='sssflx') then
         call zaiowr(sssflx(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'sssflx  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      if (trim(archv_fields(ifld)%fextract)=='sstflx') then
         call zaiowr(sstflx(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'sstflx  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! ==================================================================
      ! atmospheric forcings
      ! precip (climatology)
      if (trim(archv_fields(ifld)%fextract)=='precip') then
         fld = precip(:,:,l0)*w0+precip(:,:,l1)*w1+
     &         precip(:,:,l2)*w2+precip(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'precip  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! air temp (T2M)
      if (trim(archv_fields(ifld)%fextract)=='airtmp') then
         fld = airtmp(:,:,l0)*w0+airtmp(:,:,l1)*w1+
     &         airtmp(:,:,l2)*w2+airtmp(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'airtmp  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! clouds (TCC)
      if (trim(archv_fields(ifld)%fextract)=='clouds') then
         fld = clouds(:,:,l0)*w0+clouds(:,:,l1)*w1+
     &         clouds(:,:,l2)*w2+clouds(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'clouds  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! uwind (U10M)
      if (trim(archv_fields(ifld)%fextract)=='uwind') then
         fld = uwind(:,:,l0)*w0+uwind(:,:,l1)*w1+
     &         uwind(:,:,l2)*w2+uwind(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'uwind   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! vwind (V10M)
      if (trim(archv_fields(ifld)%fextract)=='vwind') then
         fld = vwind(:,:,l0)*w0+vwind(:,:,l1)*w1+
     &         vwind(:,:,l2)*w2+vwind(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'vwind   ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! slp (sea level pressure: MSL)
      if (trim(archv_fields(ifld)%fextract)=='slpres') then
         fld = slp(:,:,l0)*w0+slp(:,:,l1)*w1+
     &         slp(:,:,l2)*w2+slp(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'slpres  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! dew point temperature: D2M
      if (trim(archv_fields(ifld)%fextract)=='dewtmp') then
         fld = dewtmp(:,:,l0)*w0+dewtmp(:,:,l1)*w1+
     &         dewtmp(:,:,l2)*w2+dewtmp(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'dewtmp  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif

      ! relative humidity (from dew point temperature: D2M)
      if (trim(archv_fields(ifld)%fextract)=='relhum') then
         fld = relhum(:,:,l0)*w0+relhum(:,:,l1)*w1+
     &         relhum(:,:,l2)*w2+relhum(:,:,l3)*w3
         call zaiowr(fld(1-nbdy,1-nbdy),ip,.true.,
     &               xmin,xmax, nopa, .false.)
         if     (mnproc.eq.1) then
         write (nop,117) 'relhum  ',nstep,dtime1,0,coord,xmin,xmax
         call flush(nop)
         endif
      endif
      ! ==================================================================
#endif/*NERSC_VERSION*/

c
 117  format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7)
c
c --- output time-averaged mass fluxes, if required
c
      if (.not. (mxlkpp .or. mxlmy .or. mxlgiss) .and. kkout.eq.kk) then
         if (trim(archv_fields(ifld)%fextract)=='diafx') then
           do k=1,kk
             coord=sigma(k)
             call zaiowr(diaflx(1-nbdy,1-nbdy,k),ip,.true.,
     &                   xmin,xmax, nopa, .false.)
             if     (mnproc.eq.1) then
             write (nop,118) 'diafx',intvl,nstep,dtime1
     &                           ,k,coord,xmin,xmax
             call flush(nop)
             endif !1st tile
 118         format (a5,a3,' =',i11,f11.3,i3,f7.3,1p2e16.7)
           enddo
         endif
      endif !diaflx
c
      end do !TW: loop over fields to dump
c
c     if(mnproc==1)print*,'archiv: closing the files'
      close (unit=nop)
      call zaiocl(nopa)
      call xcsync(no_flush)
c     if(mnproc==1)print*,'archiv: files closed'
c
      endif  !not 1-D

c ----------------------------------------------------------
      !1-D results go into a text file
      if     (itest.gt.0 .and. jtest.gt.0) then
        if     (relaxf .and. sstflg.le.1) then
          sstc = twall(itest,jtest,1,lc0)*wc0+
     &           twall(itest,jtest,1,lc1)*wc1+
     &           twall(itest,jtest,1,lc2)*wc2+
     &           twall(itest,jtest,1,lc3)*wc3
        else !synoptic observed sst
          sstc = seatmp(itest,jtest,l0)*w0+
     &           seatmp(itest,jtest,l1)*w1+
     &           seatmp(itest,jtest,l2)*w2+
     &           seatmp(itest,jtest,l3)*w3
        endif
        sssc = swall(itest,jtest,1,lc0)*wc0+
     &         swall(itest,jtest,1,lc1)*wc1+
     &         swall(itest,jtest,1,lc2)*wc2+
     &         swall(itest,jtest,1,lc3)*wc3

#if defined (NERSC_VERSION)
        !new filenames
        print*,'archiv (ARCHIVE_SELECT): saving 1D diagnostics to: '
     &            //trim(flnm_arc)//'.txt'
        open (unit=nop,file=trim(flnm_arc)//'.txt',status='replace') !uoff+13
#else
        print*,'archiv (ARCHIVE_SELECT): saving 1D diagnostics to: '
     &            //flnmarc(1:ldot)//'.txt'
        open (unit=nop,file=flnmarc(1:ldot)//'.txt',status='replace') !uoff+13
#endif

        write (nop,'(3a / a,6i7,2f8.1,i7,i7.4,i7.3,i7.2)')
     &      '##   expt    idm    jdm    kdm',
     &        '  itest  jtest  lontst  lattst',
     &        ' yrflag   year    day     hr',
     &      '##',iexpt,  itdm,  jtdm,   kdm,
     &          ittest,jttest,
     &          mod(plon(itest,jtest),360.0),plat(itest,jtest),
     &          yrflag, iyear,  iday, ihour
        write (nop,'(7a / a,f10.3, f8.2,4f8.1, 2f9.2,2f8.4,
     &                    f9.5,4f9.3, 2f8.3, 3f8.3, 4f8.2)')
     &    '## model-day',
     &    '  srfhgt  sswflx  mixflx  surflx  sstflx',
     &    '      E-P   sssE-P  bhtflx  buoflx',
     &    '    ustar   hekman    dpbbl     dpbl   dpmixl',
     &    '   tclim   sclim',
     &    '    tmix    smix   thmix    umix    vmix',
     &    '   ubavg   vbavg',
     &    '#',time,                                              !model-day
     &    srfhgt(itest,jtest)*100.0/g,                           !cm
     &    sswflx(itest,jtest),                                   !W/m**2
     &    mixflx(itest,jtest),                                   !W/m**2
     &    surflx(itest,jtest),                                   !W/m**2
     &    sstflx(itest,jtest),                                   !W/m**2
     &    salflx(itest,jtest)*thref*8.64E7/saln(itest,jtest,1,n),!mm/day
     &    sssflx(itest,jtest)*thref*8.64E7/saln(itest,jtest,1,n),!mm/day
     &    bhtflx(itest,jtest)*1.e6,                         !1.e6*m**2/sec**3
     &    buoflx(itest,jtest)*1.e6,                         !1.e6*m**2/sec**3
     &     ustar(itest,jtest),                                   !m/s?
     &    min(hekman(itest,jtest),         9999.999),            !m
     &    min( dpbbl(itest,jtest)  *qonem, 9999.999),            !m
     &    min(  dpbl(itest,jtest)  *qonem, 9999.999),            !m
     &    min(dpmixl(itest,jtest,n)*qonem, 9999.999),            !m
     &      sstc,                                                !degC
     &      sssc,                                                !psu
     &      tmix(itest,jtest),                                   !degC
     &      smix(itest,jtest),                                   !psu
     &     thmix(itest,jtest)+thbase,                            !SigmaT
     &    max(-999.99,min(999.99,
     &        (umix(itest,jtest)+ubavg(itest,jtest,n))*100.0)),  !cm/s
     &    max(-999.99,min(999.99,
     &        (vmix(itest,jtest)+vbavg(itest,jtest,n))*100.0)),  !cm/s
     &    max(-999.99,min(999.99,
     &        ubavg(itest,jtest,n)*100.0)),                      !cm/s
     &    max(-999.99,min(999.99,
     &        vbavg(itest,jtest,n)*100.0))                       !cm/s
        if     (iceflg.ne.0) then
          write (nop,'(2a / a,f10.3, 3f8.2,2f8.1,f9.2)')
     &    '## model-day',
     &    '  covice  thkice  temice  flxice  fswice   iceE-P',
     &    '#',time,                                                !model-day
     &      covice(itest,jtest)*100.0,                             !%
     &      thkice(itest,jtest),                                   !m
     &      temice(itest,jtest),                                   !degC
     &      flxice(itest,jtest),                                   !W/m**2
     &      fswice(itest,jtest),                                   !W/m**2
     &      sflice(itest,jtest)*thref*8.64E7/saln(itest,jtest,1,n) !mm/day
        endif !iceflg
        if     (ntracr.eq.0) then
          write(cformat,'(a)')
     &      '(3a / (i4,2f8.2,3f8.3,f9.3,f10.3,2f8.2))'
        else
          write(cformat,'(a,i2,a,i2,a)')
     &      '(3a,', ntracr,
     &      'a / (i4,2f8.2,3f8.3,f9.3,f10.3,2f8.2,', ntracr,
     &      'f8.3))'
        endif
        write (nop,cformat)
     &      '#  k',
     &      '    utot    vtot    temp    saln    dens',
     &      '    thkns      dpth  viscty  t-diff',
     &      ('  tracer',ktr=1,ntracr),
     &      (k,
     &       max(-999.99,min(999.99,
     &           (u(itest,jtest,k,n)+ubavg(itest,jtest,n))*100.0)),  !cm/s
     &       max(-999.99,min(999.99,
     &           (v(itest,jtest,k,n)+vbavg(itest,jtest,n))*100.0)),  !cm/s
     &       temp(itest,jtest,k,n),                                  !degC
     &       saln(itest,jtest,k,n),                                  !psu
     &       th3d(itest,jtest,k,n)+thbase,                           !SigmaT
     &         dp(itest,jtest,k,n)*qonem,                            !m
     &         (p(itest,jtest,k+1)+p(itest,jtest,k))*0.5*qonem,      !m
     &       vcty(itest,jtest,k+1)*1.e4,                             !m**2/s*2
     &       dift(itest,jtest,k+1)*1.e4,                             !m**2/s*2
     &       (tracer(itest,jtest,k,n,ktr),ktr=1,ntracr),             !0-999?
     &       k=1,kk)
        close (unit=nop)
      endif !test point tile
      call xcsync(no_flush)
c ----------------------------------------------------------

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      return
      end subroutine archiv
#endif/*ARCHIVE_SELECT*/
      end module m_archive_select