!
!****************************************************************************************
!                                                                                       *
!       Read binary format v5.0 (hybrid S-Z models) of *hvel.67 (sidecenter vel.) and
!       compute element volumes and fluxes between elements.
!       WARNING: does not work for lat/lon!
!
!       Usage:
!       %proc MODEL_DATA_PATH OUTPUT_PATH FIRST_STACK LAST_STACK
!
!       Inputs: *_elev.61, *hvel.67, *temp.70, *salt.70, flux_regions.gr3 (depth=0,1,...M;
!               element flags based on max. of 3 node depths); screen inputs
!       Outputs: fluxes.out
!                volume.out
!       History:
!****************************************************************************************
! COMPILATION
!     ifort -Bstatic -O2 -assume byterecl -o compute_dual_fluxes3d_uvel compute_dual_fluxes3d_uvel.f90 selfe_geometry.f90
!  gfortran -Bstatic -Wall -O2 -o compute_dual_fluxes3d_uvel compute_dual_fluxes3d_uvel.f90 selfe_geometry.f90

      program read_out
      implicit none
      integer, parameter :: nbyte=4
      integer, allocatable :: idry(:),idry_s(:),idry_e(:)
      integer, allocatable :: kbp00(:), kbs(:),kbe(:),kbp2(:),kbs2(:),kbe2(:)
      integer, allocatable :: nm(:,:),openbndnodes(:),isopenbnd(:)
      real, allocatable :: zs(:,:), ze(:,:), ztot(:), ztmp(:,:),zeOld(:,:)
      integer, allocatable :: ic3(:,:),js(:,:),is(:,:),isidenode(:,:)
      real, allocatable :: sigma(:),cs(:),dp(:),dps(:),dpe(:)
      real, allocatable :: x(:),y(:),xcj(:),ycj(:)
      real, allocatable :: eta2(:),suv2(:,:,:),tsel(:,:,:),we_fv(:,:)
      real, allocatable :: area_el(:), area_e(:), sideuv(:,:), area_cell(:)
      real, allocatable :: sne(:,:), dihv(:,:),centeruv(:,:)
      double precision, allocatable :: hfluxes(:,:),vfluxes(:,:),nodal_flux(:,:)
      real, allocatable :: flx_adv(:,:,:) ! original horizontal flux (1:  the local x-driection) and vertical flux (2: positive upward)
      ! Transformation tensor for side frame: sframe(i,j,isd)
      real, allocatable :: sframe(:,:,:)
      real, allocatable :: distj(:)        ! Side lengths
      real, allocatable :: volume(:,:),volumeOld(:,:),volumeError(:,:),cum_influx(:,:)
      integer nopenbnds, nopenbndnodes, nopn, imaxmag
      integer first_stack, last_stack, nstacks
      real :: xel(3), yel(3)
      character*80   arg
      character*256  PATH
      character*256  OUTPATH
      real snx,sny,mag,flux,midx,midy,edgelen
      real dot1,vol, el_depth, uint, vint
      integer isd1,isd2,isd3,khh2,neW
      real dtout
      real h0, h_c, h_s, hmod2
      real theta_b, theta_f, thetan, time
      real vnor1, vnor2
      real xcon, ycon, zcon
      integer i, i23d, i34, ie,ie1,ie2, ics, istat, itmp, ivs, isd
      integer irec, irec_tmp, irec_TS, irec_u, irec_w, irec_dihv, itsflux
      integer irec0, irec0_TS, irec0_u, irec0_w, irec0_dihv
      integer iblock, istack, it
      integer j,k, kbpl, kin, kz
      integer n1, n2, n3, nargs
      integer ne, ne2, ne_e, np, np2, ns, ns2, nrec, ns_e, nspool, nvrt
      integer neTS, ne_eTS, ne_eW
      integer l, m, mm
      character*60 out_prefix
      character*12 it_char
      character*3 stack_char,stack_char2
      character*48 start_time,version,variable_nm,variable_dim
      character*48 data_format

!*********delwaq**********
      integer save_time
      save_time = -999
      open ( 20, file='delwaq.out'   , recl=300 )
      open ( 21, file='volumes.dwq'  , form='binary' )
      open ( 22, file='flows.dwq'    , form='binary' )
      open ( 23, file='exchanges.poi', form='binary' )
!*********delwaq**********

      nargs = iargc()
      if ( nargs /= 4 ) then
        write(*,*) 'Usage: '
        call getarg(0, arg)
        write(*,*) trim(arg)//' datadir outputdir first_stack last_stack'
        stop
      endif
      call getarg(1, PATH)
      call getarg(2, OUTPATH)
      call getarg(3, arg)
      read(arg,*) first_stack
      call getarg(4, arg)
      read(arg,*) last_stack

      write(*,*) ''
      write(*,*) 'Path:        ',trim(PATH)
      write(*,*) 'Output Dir:  ', trim(OUTPATH)
      write(*,*) 'First stack: ', first_stack
      write(*,*) 'Last stack:  ', last_stack

      nstacks = last_stack-first_stack+1
      itsflux = 0 ! 1 -> compute T,S
      ics = 1 ! Cartesian coordinates

      write(stack_char,'(i3)') first_stack
      stack_char = adjustl(stack_char)
      write(stack_char2,'(i3)') last_stack
      stack_char2 = adjustl(stack_char2)
      ! Prefix for output files
      write(out_prefix,'(a60)') trim(OUTPATH)//'/stacks_'//trim(stack_char)//'-'//trim(stack_char2)//'_'
      out_prefix= adjustl(out_prefix)

!      Open output file
      write(*,*) 'unit 90: ',trim(out_prefix)//'horz_flux.out'
      open(90,file=trim(out_prefix)//'horz_flux.out')
      write(*,*) 'unit 91: ',trim(out_prefix)//'vert_flux.out'
      open(91,file=trim(out_prefix)//'vert_flux.out')
      write(*,*) 'unit 98: ',trim(out_prefix)//'volume.out'
      open(98,file=trim(out_prefix)//'volume.out')
      write(*,*) 'unit 92: ',trim(OUTPATH)//'sides.out'
      open(92,file=trim(OUTPATH)//'sides.out')
      write(*,*) 'unit 93: ',trim(out_prefix)//'drymask.out'
      open(93,file=trim(out_prefix)//'drymask.out')
      write(*,*) 'unit 94: ',trim(OUTPATH)//'bottom_index.out'
      open(94,file=trim(OUTPATH)//'bottom_index.out')
      write(*,*) 'unit 95: ',trim(OUTPATH)//'open_bnd_nodes.out'
      open(95,file=trim(OUTPATH)//'open_bnd_nodes.out')

!      Open input files

      write(it_char,'(i12)') first_stack
      it_char = adjustl(it_char)
      write(*,*) 'unit 63: ',trim(PATH)//trim(it_char)//'_elev.61'
      open(63,file=trim(PATH)//trim(it_char)//'_elev.61',status='old',access='direct',recl=nbyte)
      write(*,*) 'unit 69: ',trim(PATH)//trim(it_char)//'_uvel.67'
      open(69,file=trim(PATH)//trim(it_char)//'_uvel.67',status='old',access='direct',recl=nbyte)
      write(*,*) 'unit 70: ',trim(PATH)//trim(it_char)//'_vvel.67'
      open(70,file=trim(PATH)//trim(it_char)//'_vvel.67',status='old',access='direct',recl=nbyte)
      write(*,*) 'unit 67: ',trim(PATH)//trim(it_char)//'_dihv.66'
      open(67,file=trim(PATH)//trim(it_char)//'_dihv.66',status='old',access='direct',recl=nbyte)
!       open(68,file=trim(PATH)//'1_vert.69',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file=trim(PATH)//trim(it_char)//'_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file=trim(PATH)//trim(it_char)//'_salt.70',status='old',access='direct',recl=nbyte)
      endif !itsflux

!..... Elevation file
! -    Header
      irec=0
      do m=1,48/nbyte
        read(63,rec=irec+m) data_format(nbyte*(m-1)+1:nbyte*m)
      enddo
      if(data_format.ne.'DataFormat v5.0') then
        write(*,*) 'This code reads only v5.0:  ',data_format
        stop
      endif
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) version(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) start_time(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) variable_nm(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) variable_dim(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte

      write(*,'(a48)')data_format
      write(*,'(a48)')version
      write(*,'(a48)')start_time
      write(*,'(a48)')variable_nm
      write(*,'(a48)')variable_dim

      read(63,rec=irec+1) nrec   ! # of output time steps (int)
      read(63,rec=irec+2) dtout  ! output time step (real)
      read(63,rec=irec+3) nspool ! skip (int)
      read(63,rec=irec+4) ivs    ! ivs (=1 or 2 for scalar or vector)
      read(63,rec=irec+5) i23d   ! i23d (=2 or 3 for 2D or 3D)
      irec=irec+5

!     Vertical grid
      read(63,rec=irec+1) nvrt    ! total # of vertical levels
      read(63,rec=irec+2) kz      ! # of Z-levels
      read(63,rec=irec+3) h0      ! dry threshold
      read(63,rec=irec+4) h_s     ! Z-sigma transition depth
      read(63,rec=irec+5) h_c     ! ~thickness of finely resolved surf layer
      read(63,rec=irec+6) theta_b ! sigma zoom param: 1: resolve surf+bottom, 0: only surf
      read(63,rec=irec+7) theta_f ! sigma zoom param: ->0: uniform sigma, >>1: coarser in middle
      irec=irec+7
      allocate(ztot(nvrt),sigma(nvrt),cs(nvrt),stat=istat)
      if(istat/=0) stop 'Failed to allocate (1)'

      do k=1,kz-1
        read(63,rec=irec+k) ztot(k) ! z-coordinates of each z-level
        ! deepest point ztot(1)
        ! actual bottom located at index kbp00 (or kbp or kbp2)
        ! kz : last index in z grid
      enddo
      do k=kz,nvrt
        kin=k-kz+1
        read(63,rec=irec+k) sigma(kin) ! sigma-coordinates of each S-level
        cs(kin)=(1-theta_b)*sinh(theta_f*sigma(kin))/sinh(theta_f)+ &
     &theta_b*(tanh(theta_f*(sigma(kin)+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
      enddo
      irec=irec+nvrt

      irec_tmp=irec !save for non-standard outputs

      !     Horizontal grid
      read(63,rec=irec+1) np ! # of nodes
      read(63,rec=irec+2) ne ! # of elements
      irec=irec+2
      write(*,*) ' np, ne: ',np,ne

      allocate(x(np),y(np),dp(np),kbp00(np),nm(ne,4),idry(np),idry_e(ne), &
     &eta2(np),dpe(ne),kbe(ne),ztmp(nvrt,np),ze(nvrt,ne),zeOld(nvrt,ne), &
     &kbp2(np),kbe2(ne),stat=istat)
      allocate(area_el(ne),stat=istat)
      if(istat/=0) stop 'Failed to allocate 2D arrays'
      allocate(area_cell(np),stat=istat)
      if(istat/=0) stop 'Failed to allocate 2D arrays'

      do m=1,np
        read(63,rec=irec+1) x(m)     ! x
        read(63,rec=irec+2) y(m)     ! y
        read(63,rec=irec+3) dp(m)    ! depth
        read(63,rec=irec+4) kbp00(m) ! initial bottom indices
        irec=irec+4
      enddo !m=1,np

      do m=1,ne
        read(63,rec=irec+1) i34      ! element type (currently must be 3)
        irec=irec+1
        do mm=1,i34
          read(63,rec=irec+1) nm(m,mm) ! node # (connectivity table)
          irec=irec+1
        enddo !mm
      enddo !m
      irec0=irec

      ! read hgrid.gr3 to get open boundary nodes
      open(14,file=trim(PATH)//'../hgrid.gr3',status='old')
      read(14,*); read(14,*) ne2,np2
      if(ne/=ne2) then
        write(*,*) np,ne
        write(*,*) np2,ne2
        stop 'mismatch in hgrid nodes'
      endif
      do i=1,np2; read(14,*); enddo ! skip node coordinates
      do i=1,ne2; read(14,*); enddo ! skip connectivity table
      read(14,*) nopenbnds
      read(14,*) nopenbndnodes
      write(*,*) nopenbnds,nopenbndnodes
      allocate(openbndnodes(nopenbndnodes),isopenbnd(np),stat=istat)
      if(istat/=0) stop 'Failed to allocate open bnd arrays'
      isopenbnd = 0
      l = 1
      do i=1,nopenbnds
        read(14,*) nopn
        do j=1,nopn
          read(14,*) openbndnodes(l)
          isopenbnd( openbndnodes(l) ) = 1
          l=l+1
        enddo
      enddo
      close(14)

!     Write out open boundary nodes
      write(95,*) nopenbndnodes, '# total number of nodes'
      write(95,'(i9)') openbndnodes
      close(95)

!     dihv.66
      read(67,rec=irec_tmp+1) ne2
      read(67,rec=irec_tmp+2) ne_e
!       if(np/=np2) stop 'mismatch in nodes'
      if(ne/=ne2) stop 'mismatch in dihv.66 nodes'
      allocate(dihv(2,ne),stat=istat)
      if(istat/=0) stop 'Failed to allocate hvel arrays'
      ! irec0 should work for dihv.66 too
!       irec0_dihv = irec0
      irec0_dihv = irec_tmp + 2 + 4*ne + 4*ne_e

!     vert.69
!       read(68,rec=irec_tmp+1) neW
!       read(68,rec=irec_tmp+2) ne_eW
!       if(neW/=ne) stop 'mismatch in elem sides'
!       allocate(we_fv(nvrt,ne),stat=istat)
!       if(istat/=0) stop 'Failed to allocate w arrays'
!       irec0_w = irec_tmp + 2 + 4*ne + 4*ne_eW

!     hvel.67
      read(69,rec=irec_tmp+1) ns

!*********delwaq**********
      write ( 20, * ) 'number of volumes          = ', np * (nvrt-1)
      write ( 20, * ) 'number of horizontal flows = ', ns * (nvrt-1)
      write ( 20, * ) 'total number of flows      = ', ns * (nvrt-1) + np * (nvrt-2)
!*********delwaq**********

      read(69,rec=irec_tmp+2) ns_e
      write(*,*) ' ns, nse:',ns,ns_e
      allocate(dps(ns),kbs(ns),kbs2(ns),zs(nvrt,ns),idry_s(ns),suv2(2,nvrt,ns),stat=istat)
      if(istat/=0) stop 'Failed to allocate hvel arrays'
      allocate(sideuv(nvrt,2),centeruv(np,2),stat=istat)
      irec0_u=irec_tmp+2
      do m=1,ns
        read(69,rec=irec0_u+3) dps(m)
        read(69,rec=irec0_u+4) kbs(m)
        irec0_u=irec0_u+4
      enddo !m
      irec0_u=irec0_u+4*ns_e

      ! allocate fluxes arrays
      allocate(volume(np,nvrt),volumeOld(np,nvrt),volumeError(np,nvrt),cum_influx(np,nvrt),stat=istat)
!       allocate(volume(nvrt,ne),stat=istat)
      if(istat/=0) stop 'Failed to allocate volumes'
      allocate(flx_adv(2,nvrt,ns),stat=istat)
      if(istat/=0) stop 'Transport: fail to allocate'
      allocate(hfluxes(ns,nvrt),vfluxes(np,nvrt),nodal_flux(3,nvrt),stat=istat)
      if(istat/=0) stop 'Transport: fail to allocate (2)'


!     *.70
      if(itsflux==1) then
        read(65,rec=irec_tmp+1) neTS
        read(65,rec=irec_tmp+2) ne_eTS
        allocate(tsel(2,nvrt,ne),stat=istat)
        if(istat/=0) stop 'Failed to allocate T,S arrays'
        if(neTS/=ne) then
          stop 'mismatch in elem'
        endif
        irec0_TS=irec_tmp+2
        do m=1,ne
          read(65,rec=irec0_TS+3) dpe(m)
          read(65,rec=irec0_TS+4) kbe(m)
          irec0_TS=irec0_TS+4
        enddo !m
        irec0_TS=irec0_TS+4*ne_eTS
      endif !itsflux

!     Compute geometry
      call compute_nside(np,ne,nm,ns2)
      if(ns/=ns2) stop 'mismatch in sides'
      write(*,*) 'done compute_nside'
      allocate(ic3(ne,3),js(ne,3),is(ns2,2),isidenode(ns2,2),xcj(ns2),ycj(ns2),stat=istat)
      if(istat/=0) stop 'Failed to allocate side arrays'
      !  ns: # of sides
      !  ic3(ne,3): 3 neighboring elements of an element
      !  js(ne,3): 3 sides of an element
      !  is(ns,2): 2 adjacent elements of a side
      !  isidenode(ns,2): 2 end nodes of a side
      !  xcj(ns),ycj(ns): x,y of each side center
      call selfe_geometry(np,ne,ns2,x,y,nm,ic3,js,is,isidenode,xcj,ycj)
      write(*,*) 'done selfe_geometry'
      write(*,*) 'last element',(nm(ne,j),j=1,3)


      ! Allocate and compute side frame tensor for ics=1 or 2
      ! sframe(i,j,isd): where j is the axis id, i is the component id, isd is the local side id
      ! For ics=1, only sframe(1:2,1:2,isd) are used
      allocate(sframe(3,3,ns),stat=istat);
      if(istat/=0) stop 'AQUIRE_HGRID: sframe allocation failure'
      allocate(distj(ns),sne(nvrt,3),area_e(ne),stat=istat);
      if(istat/=0) stop 'AQUIRE_HGRID: distj allocation failure'

!   ! Build side data for augmented subdomain
!   do ie=1,ne
!     do j=1,3
!       jsj=js(ie,j)
!       n1=nm(ie,nx(j,1)) ! isidenode(jsj,1)
!       n2=nm(ie,nx(j,2))
! !       dps(jsj)=(dp(n1)+dp(n2))/2d0
!       distj(jsj)=sqrt((x(n2)-x(n1))**2+(y(n2)-y(n1))**2+(z(n2)-z(n1))**2)
!       if(distj(jsj)==0) then
!         write(*,*) 'AQUIRE_HGRID: Zero side',jsj
!         stop
!       endif
!     enddo !j=1,3
!   enddo !ie=1,nea

      ! Build side data for augmented subdomain
      do i=1,ns    ! for all sides
        !Internal wet sides
        ie1=is(i,1); ie2=is(i,2) ! the elements sharing the side
        n1=isidenode(i,1); n2=isidenode(i,2) ! side end nodes
        distj(i)=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2) ! side length
        ! znd zero for cartesian
!         distj(jsj)=sqrt((x(n2)-x(n1))**2+(y(n2)-y(n1))**2+(znd(n2)-znd(n1))**2)
      enddo

      thetan=-1.e10 !max. deviation between ze and zs axes
      sframe=0 !for ics=1
      do j=1,ns
        n1=isidenode(j,1)
        n2=isidenode(j,2)
        if(ics==1) then ! cartesian
          thetan=atan2(x(n1)-x(n2),y(n2)-y(n1))
          sframe(1,1,j)=cos(thetan) !old snx
          sframe(2,1,j)=sin(thetan) !old sny
          sframe(1,2,j)=-sframe(2,1,j)
          sframe(2,2,j)=sframe(1,1,j)
        else !lat/lon
          ! omitted
        endif !ics
      enddo !j=1,ns

!...  Element area
      area_el=0
      area_cell=0
      do i=1,ne
        xel(1:3)=x(nm(i,1:3))
        yel(1:3)=y(nm(i,1:3))
        area_el(i)=((xel(1)-xel(3))*(yel(2)-yel(3))-(xel(2)-xel(3))*(yel(1)-yel(3)))/2d0
        do j=1,3 ! node
          k = nm(i,j)
          area_cell(k) = area_cell(k) + area_el(i)/3
        enddo
      enddo !i=1,nea

!     Write side definition
      write(92,'(a)') 'iSide    cell1    cell2'
      write(92,'(i9)') ns
      do j=1,ns
        write(92,'(i9,i9,i9)') j,isidenode(j,1),isidenode(j,2)
      enddo !i

!*********delwaq**********
      write ( 23 ) ( ( isidenode(isd,1)+k*np, isidenode(isd,2)+ k   *np, 0, 0, isd = 1,ns ) , k=0,nvrt-2 )
      write ( 23 ) ( ( i               +k*np, i               +(k+1)*np, 0, 0, i   = 1,np ) , k=0,nvrt-3 )
!*********delwaq**********

!...  Time iteration
!...
      do istack=first_stack,last_stack
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')istack
      it_char = adjustl(it_char)
      open(63,file=trim(PATH)//trim(it_char)//'_elev.61',status='old',access='direct',recl=nbyte)
      open(69,file=trim(PATH)//trim(it_char)//'_uvel.67',status='old',access='direct',recl=nbyte)
      open(70,file=trim(PATH)//trim(it_char)//'_vvel.67',status='old',access='direct',recl=nbyte)
      open(67,file=trim(PATH)//trim(it_char)//'_dihv.66',status='old',access='direct',recl=nbyte)
!       open(68,file=trim(PATH)//trim(it_char)//'_vert.69',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file=trim(PATH)//trim(it_char)//'_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file=trim(PATH)//trim(it_char)//'_salt.70',status='old',access='direct',recl=nbyte)
      endif !itsflux

      irec=irec0
      irec_u=irec0_u
      irec_dihv=irec0_dihv
      irec_TS=irec0_TS
      irec_w=irec0_w
      do iblock=1,nrec
        read(63,rec=irec+1) time
        read(63,rec=irec+2) it
        irec=irec+2

        if (iblock == 1 .or. iblock == nrec) then
          write(*,*) 'i=',trim(it_char),' iblock=',iblock,' time (days) =',time/86400
        endif

        do i=1,np
          read(63,rec=irec+i) eta2(i)
        enddo !i
        irec=irec+np+np !done with (63 for this step
        irec_dihv=irec_dihv+ne+2 ! skip eta field
        irec_u=irec_u+ns+2
        irec_TS=irec_TS+ne+2
        irec_w=irec_w+ne+2

!       Compute z coordinates
        ztmp=-1.e22
        do i=1,np
          if(eta2(i)+dp(i)<h0) then !node dry
            idry(i)=1
            kbp2(i)=-1
          else !wet
            idry(i)=0
!           Compute z-coordinates
            do k=kz,nvrt
              kin=k-kz+1
              hmod2=min(dp(i),h_s)
              if(hmod2<=h_c) then
                ztmp(k,i)=sigma(kin)*(hmod2+eta2(i))+eta2(i)
              else if(eta2(i)<=-h_c-(hmod2-h_c)*theta_f/sinh(theta_f)) then
                write(*,*)'Pls choose a larger h_c (2):',eta2(i),h_c
                stop
              else
                ztmp(k,i)=eta2(i)*(1+sigma(kin))+h_c*sigma(kin)+(hmod2-h_c)*cs(kin)
              endif

!             Following to prevent underflow
              if(k==kz) ztmp(k,i)=-hmod2
              if(k==nvrt) ztmp(k,i)=eta2(i)
            enddo !k

            if(dp(i)<=h_s) then !If node depth smaller (+) than transition depth between S and Z
              kbpl=kz           !bottom point id equals number of Z levels
            else !z levels
!             Find bottom index
              kbpl=0
              do k=1,kz-1
                if(-dp(i)>=ztot(k).and.-dp(i)<ztot(k+1)) then
                  kbpl=k
                  exit
                endif
              enddo !k
              if(kbpl==0) then
                write(*,*)'Cannot find a bottom level:',dp(i)
                stop
              endif
              ztmp(kbpl,i)=-dp(i)
              do k=kbpl+1,kz-1
                ztmp(k,i)=ztot(k)
              enddo !k
            endif
            !if(kbpl/=kbp00(i)) then
            !  write(*,*)'Bottom index wrong:',i,kbpl,kbp00(i),eta2(i),dp(i)
            !  stop
            !endif

            do k=kbpl+1,nvrt
              if(ztmp(k,i)-ztmp(k-1,i)<=0) then
                write(*,*)'Inverted z-level:',eta2(i),dp(i),ztmp(k,i)-ztmp(k-1,i)
                stop
              endif
            enddo !k
            kbp2(i)=kbpl
          endif !dry/wet
        enddo !i=1,np

!       idry_e, idry_s
        do i=1,ne
          idry_e(i)=maxval(idry(nm(i,1:3)))
          kbe2(i)=maxval(kbp2(nm(i,1:3)))
        enddo !i
        idry_s=1 !init.
        do i=1,ns
          kbs2(i)=maxval(kbp2(isidenode(i,1:2)))
          do j=1,2
            ie=is(i,j)
            if(ie/=0.and.idry_e(max(1,ie))==0) idry_s(i)=0
          enddo !j
          !Check
          if(idry_s(i)==0) then
            itmp=maxval(idry(isidenode(i,1:2)))
            if(itmp/=0) then
              write(*,*)'Wet side has dry nodes:',i,isidenode(i,1:2)
              stop
            endif
          endif
        enddo !i
        kbe = kbe2

        ! Write bottom indices (only once)
        if(istack==first_stack.and.iblock==1) then
          write(94,'(a)') 'iCell    kBottom'
          write(94,'(i9)') ne
          do ie=1,ne
            write(94,'(i9,i9)') ie,kbe(ie)
          enddo !i
        endif


!       dihv.66
        dihv=0 !for below bottom
        do i=1,ne
            do m=1,2
              read(67,rec=irec_dihv+1) dihv(m,i)
              irec_dihv=irec_dihv+1
            enddo !m
          !Debug
          !if(istack==nstacks.and.iblock==nrec) write(97,*)xcj(i),ycj(i),suv2(1:2,nvrt,i)
        enddo !i

!       hvel.67
        suv2=0 !for below bottom
        do i=1,ns
          do k=max0(1,kbs(i)),nvrt
!             do m=1,2
!               read(64,rec=irec_u+1) suv2(m,k,i)
!               irec_u=irec_u+1
!             enddo !m
            read(69,rec=irec_u+1) suv2(1,k,i)
            read(70,rec=irec_u+1) suv2(2,k,i)
            irec_u=irec_u+1
          enddo !k

          !Debug
          !if(istack==nstacks.and.iblock==nrec) write(97,*)xcj(i),ycj(i),suv2(1:2,nvrt,i)
        enddo !i

        do i=1,ns
          if(idry_s(i)==1) cycle

          !Compute zcor for wet sides (deal with some inconsistency)
          n1=isidenode(i,1); n2=isidenode(i,2)
          do k=kbs2(i),nvrt
            zs(k,i)=(ztmp(k,n1)+ztmp(k,n2))/2
            if(abs(zs(k,i))>1.e18) then
              write(*,*)'Dry side (1):',i,k,n1,n2,ztmp(k,n1),ztmp(k,n2)
              stop
            endif
          enddo !k
        enddo !i=1,ns

!       T,S
        if(itsflux==1) then
          tsel=-99
          do i=1,ne
            do k=max0(1,kbe(i)),nvrt
              read(65,rec=irec_TS+1) tsel(1,k,i)
              read(66,rec=irec_TS+1) tsel(2,k,i)
              irec_TS=irec_TS+1
            enddo !k

            !Debug
            !if(istack==ndays.and.iblock==nrec) write(99,*)i,tsel(2,nvrt,i)
          enddo !i
        endif !itsflux

!       vert.69
!         we_fv=0 !for below bottom
!         do i=1,ne
!           do k=max0(1,kbe2(i)),nvrt
!             read(68,rec=irec_w+1) we_fv(k,i)
!             irec_w=irec_w+1
!           enddo !k
!         enddo !i

        zeOld=ze
        do i=1,ne
          if(idry_e(i)==1) cycle

          !Compute zcor for wet elem.
          n1=nm(i,1); n2=nm(i,2); n3=nm(i,3)
          do k=kbe2(i),nvrt
            ze(k,i)=(ztmp(k,n1)+ztmp(k,n2)+ztmp(k,n3))/3
            if(abs(ze(k,i))>1.e18) then
              write(*,*)'Dry elem:',i,ztmp(k,n1),ztmp(k,n2),ztmp(k,n3)
              stop
            endif
          enddo !k

        enddo !i=1,ne

        hfluxes=0
        vfluxes=0
        ! loop over elems
        do ie=1,ne
          if(idry_e(ie)==1) cycle
          el_depth = ze(nvrt,ie)-ze(kbe2(ie),ie)
          ! compute hvel at triangle center
          centeruv = 0
          uint = 0 ! vertical integral
          vint = 0 ! vertical integral
          do k=kbe2(ie),nvrt-1
            do j=1,3 ! sides
              isd = js(ie,j)
              centeruv(k,1) = centeruv(k,1)+ (suv2(1,k+1,isd)+suv2(1,k,isd))/2/3
              centeruv(k,2) = centeruv(k,2)+ (suv2(2,k+1,isd)+suv2(2,k,isd))/2/3
            enddo
            uint=uint+centeruv(k,1)*(ze(k+1,ie)-ze(k,ie))
            vint=vint+centeruv(k,2)*(ze(k+1,ie)-ze(k,ie))
          enddo
          ! correct 3D velocity field
          if( abs(dihv(1,ie)/uint) < 1.6 ) then
            ! correct by multiplying the profile
            do k=kbe2(ie),nvrt-1
              centeruv(k,1) = centeruv(k,1)/uint*dihv(1,ie)
            enddo
          else
            ! correct by substracting the error
            uint  = uint - dihv(1,ie)
            ! correction
            do k=kbe2(ie),nvrt-1
              centeruv(k,1) = centeruv(k,1) - uint/el_depth
            enddo
          endif
          if( abs(dihv(2,ie)/vint) < 1.6 ) then
            ! correct by multiplying the profile
            do k=kbe2(ie),nvrt-1
              centeruv(k,2) = centeruv(k,2)/vint*dihv(2,ie)
            enddo
          else
            ! correct by substracting the error
            vint  = vint - dihv(2,ie)
            ! correction
            do k=kbe2(ie),nvrt-1
              centeruv(k,2) = centeruv(k,2) - vint/el_depth
            enddo
          endif
          ! compute the 3 nodal fluxes ( flux towards node nm(ie,i) )
          do i=1,3 ! nodes
            n1 = nm(ie,i)
            ! other nodes
            n2 = nm(ie,modulo(i+0,3)+1)
            n3 = nm(ie,modulo(i+1,3)+1)
            ! opposite edge mid point
            midx = 0.5*(x(n2)+x(n3))
            midy = 0.5*(y(n2)+y(n3))
            ! length of opposite edge
            edgelen = sqrt( (x(n2)-x(n3))**2 + (y(n2)-y(n3))**2 )
            ! vector normal to opposite edge (towards node n1)
            snx = -( y(n2)-y(n3) ) ! normal
            sny = +( x(n2)-x(n3) )
            if ( (x(n1)-midx)*snx+(y(n1)-midy)*sny < 0 ) then
              snx=-snx
              sny=-sny
            endif
            ! flux towards node n1 from element ie
            do k=kbe2(ie),nvrt-1
              nodal_flux(i,k) = 1/2.0*( centeruv(k,1)*snx + centeruv(k,2)*sny )*(ze(k+1,ie)-ze(k,ie))
              if ( isnan(nodal_flux(i,k)) .or. abs(nodal_flux(i,k)) > huge(nodal_flux(i,k)) ) then
                write(*,*)'flux is nan/inf',nodal_flux(i,k),i,k,centeruv(k,1)*snx, centeruv(k,2)*sny, (ze(k+1,ie)-ze(k,ie))
                stop
              endif
            enddo

          enddo
          ! assemble fluxes on sides
          do j=1,3 ! sides
            isd = js(ie,j)
            ! end nodes
            n1 = isidenode(isd,1)
            n2 = isidenode(isd,2)
            do i=1,3 ! nodes
              if(nm(ie,i)==n2) then
              ! end node is target node
                do k=kbe2(ie),nvrt-1
                  hfluxes(isd,k) = hfluxes(isd,k) + nodal_flux(i,k)/3
                enddo
              elseif (nm(ie,i)==n1) then
              ! start node is target node
                do k=kbe2(ie),nvrt-1
                  hfluxes(isd,k) = hfluxes(isd,k) - nodal_flux(i,k)/3
                enddo
              endif
            enddo
          enddo
        enddo

        ! element volume
        volumeOld=volume
        volume=0

        ! compute volume change based on horizontal fluxes
        cum_influx = 0
        do isd=1,ns
!           if(idry_s(isd)==1) cycle
!         choose the end node
          n1 = isidenode(isd,1)
          n2 = isidenode(isd,2)
          do k=kbs2(isd),nvrt-1
            cum_influx(n1,k) = cum_influx(n1,k) - hfluxes(isd,k)
            cum_influx(n2,k) = cum_influx(n2,k) + hfluxes(isd,k)
          enddo
        enddo
        ! compute new volume from z coords
        ! compute vertical flux from geometry and h. fluxes
        do i=1,np
!           if(idry(i)==1) cycle
          ! bottom to top
          do k=kbp2(i),nvrt-1
            volume(i,k) = area_cell(i)*(ztmp(k+1,i) - ztmp(k,i))
            ! recursion on vertical fluxes, vfluxes(i,kbp2(i)) = 0
            vfluxes(i,k+1) = -(volume(i,k)-volumeOld(i,k))/dtout +cum_influx(i,k) +vfluxes(i,k)
          enddo
!           vfluxes(i,nvrt) = 0 ! force impermeable surface
!           ! top to bottom
!           do k=nvrt-1,kbp2(i),-1
!             volume(i,k) = area_cell(i)*(ztmp(k+1,i) - ztmp(k,i))
!             ! recursion on vertical fluxes, vfluxes(i,kbp2(i)) = 0
!             vfluxes(i,k) = +(volume(i,k)-volumeOld(i,k))/dtout -cum_influx(i,k) +vfluxes(i,k+1)
!           enddo
!           ! force impermeable surface ( vfluxes(i,nvrt)=0 )
!           do k=kbp2(i),nvrt-1
!             vfluxes(i,k+1) = vfluxes(i,k+1) - vfluxes(i,nvrt)*(ztmp(k+1,i) - ztmp(k,i))/((ztmp(nvrt,i) - ztmp(k,i)))
!           enddo
        enddo

!*********delwaq**********
        write ( 20, * ) 'time: ',istack, int(time), nvrt
        write ( 21 ) int(time), volume ( : , 1:nvrt-1 )
        if ( save_time .ne. -999 ) write ( 22 ) save_time, sngl(hfluxes(:,1:nvrt-1)), sngl(vfluxes(:,2:nvrt-1))
        save_time = int(time)
!*********delwaq**********

        ! check volume conservation
        mag = 0 ! max abs error
        imaxmag = 0
        do i=1,np
          if(idry(i)==1) cycle
          do k=kbp2(i),nvrt-1
            vol = volume(i,k)-volumeOld(i,k)
            volumeError(i,k) = ( vol-cum_influx(i,k)*dtout - (vfluxes(i,k)-vfluxes(i,k+1))*dtout )/(volume(i,k))
            if(isopenbnd(i)==0) then
!               mag = max( mag, abs( volumeError(i,k) ) )
              if (abs( volumeError(i,k) ) > mag ) then
                mag = abs( volumeError(i,k) )
                imaxmag = i
              endif
            endif
!             write(*,*) i, k, volumeError(i,k), vfluxes(i,k+1)
          enddo
!           do k=kbp2(i),nvrt
!             write(*,*) i, k, vfluxes(i,k)
!           enddo
        enddo
        write(*,*) iblock, 'max error',mag,imaxmag
!       stop

!#####   Output

!       Horizontal fluxes (from elem1 to elem2 across the side)
        write(90,'(a,f10.7)') 'Time(days) = ',time/86400
        write(90,'(i9,i9)') ns, nvrt-1
        do j=1,ns
          write(90,'(i9,6000(1x,f9.3))') j,hfluxes(j,1:nvrt-1)
        enddo !i

!       Vertical fluxes (positive upwards)
        write(91,'(a,f10.7)') 'Time(days) = ',time/86400
        write(91,'(i9,i9)') np, nvrt
        do i=1,np
          write(91,'(i9,6000(1x,f9.3))') i,vfluxes(i,:)
!           write(90,'(i9,6000(1x,e25.16))') ie,flx_adv(1,2:nvrt,i)
        enddo !i

!       Volume
        write(98,'(a,f10.7)') 'Time(days) = ',time/86400
        write(98,'(i9,i9)') np, nvrt-1
        do i =1,np
          write(98,'(i9,6000(1x,f9.3))') i, volume(i,1:nvrt-1)
!           write(98,'(i9,6000(1x,e25.16))') ie, volume(2:nvrt,ie)
        enddo
!         write(98,'(i9,i9)') np
!         do i =1,np
!           write(98,'(i9,6000(1x,e20.11))') i, volume(i)
! !           write(98,'(i9,6000(1x,e25.16))') ie, volume(2:nvrt,ie)
!         enddo

!       Dry mask (0 - active element, 1 - dry element)
        write(93,'(a,f10.7)') 'Time(days) = ',time/86400
        write(93,'(i9)') np
        do i=1,np
          write(93,'(i9,i9)') i,idry(i)
        enddo !i

      enddo !iblock=1,nrec

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      enddo !istack=1,nstacks

!*********delwaq**********
      hfluxes = 0.0
      vfluxes = 0.0
      write ( 22 ) int(time), sngl(hfluxes(:,1:nvrt-1)),sngl(vfluxes(:,2:nvrt-1))
!*********delwaq**********

      close(90)
      close(91)
      close(92)
      close(93)
      close(94)
      close(98)

      close(63)
      close(69)
      close(70)


      stop
      end