program mrgl_archv
      use mod_plot  ! HYCOM plot array interface
      use mod_za    ! HYCOM array I/O interface
      implicit none
c
c --- merge layers in a HYCOM 2.0 archive file.
c
      character label*81,text*18,flnm_i*240,flnm_o*240
      logical initl,trcout,icegln

      logical strict
c
      integer          artype,iexpt,iversn,yrflag,kpalet,mxlflg
      integer          flxflg
      integer          i,ibad,j,k,kb,kl,kt,kkin,kkout,l,laybot(0:99)
      real             sigma(99),thbase,dpoij,dpuij,dpvij,
     &                 dij,dpthin,onem,pij
      real             puij,uij,pvij,vij
      double precision time3(3),time,year
c
      real, allocatable :: dpo(:,:,:),dpu(:,:),dpv(:,:)
c
      data trcout/.false./  ! must be .false. (no tracer remapping)
      data initl /.true. /
      data strict/.true./
c

      call xcspmd
      call zaiost
      lp=6
c
      onem=9806.0
c
c
c --- 'flnm_i' = name of original archive file
c --- 'flnm_o' = name of target   archive file
c --- 'iexpt ' = experiment number x10  (000=from archive file)
c --- 'yrflag' = days in year flag (0=360J16,1=366J16,2=366J01,3=actual)
c --- 'idm   ' = longitudinal array size
c --- 'jdm   ' = latitudinal  array size
c --- 'kdmold' = original number of layers
c --- 'kdmnew' = target   number of layers
c
      read (*,'(a)') flnm_i
      write (lp,'(2a)') ' input file: ',flnm_i(1:len_trim(flnm_i))
      call flush(lp)
      read (*,'(a)') flnm_o
      write (lp,'(2a)') 'output file: ',flnm_o(1:len_trim(flnm_o))
      call flush(lp)
      call blkini(iexpt, 'iexpt ')
      call blkini(yrflag,'yrflag')
      call blkini(ii,    'idm   ')
      call blkini(jj,    'jdm   ')
      call blkini(kkin,  'kdmold')
      call blkini(kkout, 'kdmnew')
      call blkini(flxflg,  'flxflg')
      if     (ii.ne.idm .or. jj.ne.jdm) then
        write(lp,*)
        write(lp,*) 'error - wrong idm or jdm (should be:',
     &                                         idm,jdm,')'
        write(lp,*)
        call flush(lp)
        stop
      endif
      iorign = 1
      jorign = 1
c
c --- 'thbase' = reference density (sigma units)
c
      call blkinr(thbase,
     &           'thbase','("blkinr: ",a6," =",f11.4," sig")')
c
c --- new layer conbinations and densities (sigma units)
c
      write(lp,*)
      laybot(0) = 0
      do k=1,kkout
c
c ---   'laybot' = last layer in next combination of layers
c
        call blkini(laybot(k),  'laybot')
        if     (laybot(k).le.laybot(k-1)) then
          write(lp,'(a)') 'error - laybot must be ascending'
          stop
        elseif (k.eq.kkout .and. laybot(k).ne.kkin) then
          write(lp,'(a)') 'error - last laybot must be kkin'
          stop
        elseif (k.ne.kkout .and. laybot(k).ge.kkin) then
          write(lp,'(a)') 'error - laybot reached kkin too soon'
          stop
        endif
        if     (laybot(k).gt.laybot(k-1)+1) then
          call blkinr(sigma(k),
     &                'sigma ','("blkinr: ",a6," =",f11.4," sig")')
        else
          sigma(k) = -1.0  ! get from input archive
        endif
      enddo
c
c --- array allocation
c
      kk    = 0
      kkmax = max(kkin,kkout)
      call plot_alloc
c
      dpthfil = 'regional.depth'
c
      do j=1,jdm
        do i=1,idm
          p(i,j,1)=0.
        enddo
      enddo
c
c --- read the archive file, from "*.[ab]".
c
      kk = kkin
      call getdatb(flnm_i,time3,artype,initl,icegln,trcout,
     &             iexpt,iversn,yrflag,kkin,flxflg)        ! hycom input
      time = time3(3)
c
      if     (artype.eq.3) then
        write(lp,'(/ a /)')
     &    'error - cannot merge std.dev. archive'
        call flush(lp)
        stop
      endif
c
      do k= 1,kkout
        if     (sigma(k).lt.0.0) then
          sigma(k) = theta(laybot(k))
        endif
c
        if     (k.gt.1) then
          if      (sigma(k).le.sigma(k-1)) then
            write(lp,'(/ a,i3,2f11.4 /)')
     &        'error - sigma is not stabally stratified',
     &        k,sigma(k-1),sigma(k)
            call flush(lp)
            stop
          endif
        endif
      enddo
c
c --- land masks.
c
      call bigrid(depths)
c
c --- check that bathymetry is consistent with this archive.
c
      ibad = 0
      do j= 1,jj
        do i= 1,ii
          if     (ip(i,j).eq.1) then
            if     (srfht(i,j).gt.2.0**99) then
              ibad = ibad + 1   ! topo sea, srfht land
            endif
          else
            if     (srfht(i,j).lt.2.0**99) then
              ibad = ibad + 1   ! topo land, srfht sea
            endif
          endif
        enddo
      enddo
      if     (ibad.ne.0) then
        write(lp,*)
        write(lp,*) 'error - wrong bathymetry for this archive file'
        write(lp,*) 'number of mismatches = ',ibad
        write(lp,*)
        call flush(lp)
        stop
      endif
c
c     remap layers.
c
c     fill massless layers on sea floor with fluid from above
c
      dpthin = 0.0001*onem
      do j= 1,jj
        do i= 1,ii
          if     (ip(i,j).eq.1) then
            dij = depths(i,j)*onem - dpthin
            pij = 0.0
            do k= 1,kkin
              pij = pij + dp(i,j,k)
              if     (pij.gt.dij .and. dp(i,j,k).le.dpthin) then
                temp(i,j,2*k) = temp(i,j,2*k-2)
                saln(i,j,2*k) = saln(i,j,2*k-2)
                th3d(i,j,2*k) = th3d(i,j,2*k-2)
              endif
            enddo !k
          endif !ip.eq.1
        enddo !i
      enddo !j
c
c
      allocate( dpo(idm,jdm,kkin), dpu(idm,jdm), dpv(idm,jdm) )
c
      dpo(:,:,1:kkin) = dp(:,:,1:kkin)
c
      do k= 1,kkout
        l  = 2*k-1
        kt = laybot(k-1)+1
        kb = laybot(k)
        write(lp,'(a,i3,i4,i3)') 'k,kt,kb = ',k,kt,kb
        call flush(lp)
c
        if     (kb.eq.kt) then
c
c ---     existing layer.
c
          do j= 1,jdm
            do i= 1,idm
                 u(i,j,l) =    u(i,j,2*kb)
                 v(i,j,l) =    v(i,j,2*kb)
                dp(i,j,k) =  dpo(i,j,  kb)
              temp(i,j,l) = temp(i,j,2*kb)
              saln(i,j,l) = saln(i,j,2*kb)
              th3d(i,j,l) = th3d(i,j,2*kb)
              if     (artype.eq.2) then
                ke(i,j,l) =   ke(i,j,2*kb)
              endif
            enddo
          enddo
        else
c
c ---     combination of layers.
c
          do j= 1,jdm
            do i= 1,idm
              if     (iu(i,j).eq.1) then
                 dpuij          = dpo(i,j,kt) + dpo(i+1,j,kt)
                 if(strict) then

                    dpuij=
     $                   max(0.5*(dpo(i-1,j,kt)+dpo(i,j,kt)),0.0)
                 endif
                 dpu(i,j)       = dpuij
                   u(i,j,l)     = dpuij*   u(i,j,2*kt)
              endif
              if     (iv(i,j).eq.1) then
                 dpvij          = dpo(i,j,kt) + dpo(i,j+1,kt)
                 if(strict) then

                    dpvij=
     $                   max(0.5*(dpo(i,j-1,kt)+dpo(i,j,kt)),0.0)
                 endif
                 dpv(i,j)       = dpvij
                   v(i,j,l)     = dpvij*   v(i,j,2*kt)
              endif
              if     (ip(i,j).eq.1) then
                 dpoij          = dpo(i,j,kt)
                  dp(i,j,k)     = dpoij
                temp(i,j,l)     = dpoij*temp(i,j,2*kt)
                saln(i,j,l)     = dpoij*saln(i,j,2*kt)
                th3d(i,j,l)     = dpoij*th3d(i,j,2*kt)
                if     (artype.eq.2) then
                  ke(i,j,l)     = dpoij*ke(  i,j,2*kt)
                endif
              endif
              do kl= kt+1,kb
                if     (iu(i,j).eq.1) then
                   dpuij          = dpo(i,j,kl) + dpo(i+1,j,kl)
                   if(strict) then
                      dpuij=
     $                     max(0.5*(dpo(i-1,j,kl)+dpo(i,j,kl)),0.0)
                   endif
                   dpu(i,j)       = dpuij                + dpu(i,j)
                     u(i,j,l)     = dpuij*   u(i,j,2*kl) +   u(i,j,l)
                endif
                if     (iv(i,j).eq.1) then
                   dpvij          = dpo(i,j,kl) + dpo(i,j+1,kl)
                   if(strict) then
                      dpvij=
     $                     max(0.5*(dpo(i,j-1,kl)+dpo(i,j,kl)),0.0)
                   endif
                   dpv(i,j)       = dpvij                + dpv(i,j)
                     v(i,j,l)     = dpvij*   v(i,j,2*kl) +   v(i,j,l)
                endif
                if     (ip(i,j).eq.1) then
                   dpoij          = dpo(i,j,kl)
                    dp(i,j,k)     = dpoij                +   dp(i,j,k)
                  temp(i,j,l)     = dpoij*temp(i,j,2*kl) + temp(i,j,l)
                  saln(i,j,l)     = dpoij*saln(i,j,2*kl) + saln(i,j,l)
                  th3d(i,j,l)     = dpoij*th3d(i,j,2*kl) + th3d(i,j,l)
                  if     (artype.eq.2) then
                    ke(i,j,l)     = dpoij*ke(  i,j,2*kl) +   ke(i,j,l)
                  endif
                endif
              enddo
              if     (iu(i,j).eq.1 .and. dpu(i,j).gt.0.0) then
                   u(i,j,l)     = u(i,j,l)/dpu(i,j)
              else  ! allow for land mask and zero thickness layer
                   u(i,j,l)     = u(i,j,2*kb)
              endif
              if     (iv(i,j).eq.1 .and. dpv(i,j).gt.0.0) then
                   v(i,j,l)     = u(i,j,l)/dpv(i,j)
              else  ! allow for land mask and zero thickness layer
                   v(i,j,l)     = v(i,j,2*kb)
              endif
              if     (ip(i,j).eq.1 .and. dp(i,j,k).gt.0.0) then
                temp(i,j,l)     = temp(i,j,l)/dp(i,j,k)
                saln(i,j,l)     = saln(i,j,l)/dp(i,j,k)
                th3d(i,j,l)     = th3d(i,j,l)/dp(i,j,k)
              else  ! allow for land mask and zero thickness layer
                temp(i,j,l)     = temp(i,j,2*kb)
                saln(i,j,l)     = saln(i,j,2*kb)
                th3d(i,j,l)     = th3d(i,j,2*kb)
              endif
            enddo
          enddo
        endif
      enddo  ! k=1,kkout
c
      theta(1:kkout) = sigma(1:kkout)

      if(strict) then
c
c -- internal velocities average to zero
c
         do j= 1,jdm
            do i= 1,idm
               puij=0.0
               uij=0.0
               pvij=0.0
               vij=0.0
               do k=1,kkout
                  if     (iu(i,j).eq.1) then
                     dpuij=
     $                    max(0.5*(dp(i-1,j,k)+dp(i,j,k)),0.0)
                     uij     = dpuij* u(i,j,k)+uij
                     puij=puij+dpuij
                  endif
                  if     (iv(i,j).eq.1) then
                     dpvij=
     $                    max(0.5*(dp(i,j-1,k)+dp(i,j,k)),0.0)
                     vij     = dpvij*   v(i,j,k)+vij
                     pvij=pvij+dpvij
                  endif
               enddo
               if     (iu(i,j).eq.1) then
                  uij=uij/puij
                  do k=1,kkout
                     u(i,j,k)=u(i,j,k)-uij
                  enddo
               endif
               if     (iv(i,j).eq.1) then
                  vij=vij/pvij
                  do k=1,kkout
                     v(i,j,k)=v(i,j,k)-vij
                  enddo
               endif
            enddo
         enddo

      endif
c
c --- write the archive file.
c
      l = len_trim(flnm_o)
      if     (flnm_o(l-1:l).eq.'.a' .or. flnm_o(l-1:l).eq.'.b') then
        flnm_o(l-1:l) = '  '  ! to prevent putdat from using '*.[AB]'
      endif
      kk = kkout
      call putdat(flnm_o,artype,time3,icegln,trcout,
     &            iexpt,iversn,yrflag,kkout,-1,flxflg)
      end