!   Copyright 2014 College of William and Mary
!
!   Licensed under the Apache License, Version 2.0 (the "License");
!   you may not use this file except in compliance with the License.
!   You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
!   Unless required by applicable law or agreed to in writing, software
!   distributed under the License is distributed on an "AS IS" BASIS,
!   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!   See the License for the specific language governing permissions and
!   limitations under the License.

!
!****************************************************************************************
!											*
!	Read binary format v5.0 (hybrid S-Z models) of *hvel.67 (sidecenter vel.) and
!       compute fluxes at the borders between multiple regions (excluding open/land boundaries).
!       WARNING: does not work for lat/lon!
!       Works with mixed grids and LSC2
!    
!       Inputs: vgrid.in, *_elev.61, *hvel.67, fluxflag.prop (depth=-1,0,1,...);
!               screen inputs. If itsflux/=0, *temp.70, *salt.70 (not completed)
!       Outputs: fluxes.out (if the difference between region #s=1, and neith=-1)
!       History: 
!****************************************************************************************
!     ifort -cpp -O2 -mcmodel=medium -CB -Bstatic -o compute_fluxes_ns.WW ../UtilLib/schism_geometry.f90 ../UtilLib/compute_zcor.f90 compute_fluxes_ns.f90
      program read_out
      use compute_zcor
      use schism_geometry_mod
      parameter(nbyte=4)
      character*30 file63
      character*12 it_char
      character*48 start_time,version,variable_nm,variable_dim
      character*48 data_format
      integer,allocatable :: elnode(:,:),elside(:,:),i34(:)
      allocatable :: sigma(:),cs(:),ztot(:),x(:),y(:),dp(:),kbp(:),kfp(:)
      allocatable :: icum(:,:,:),eta2(:),ztmp(:,:)
      allocatable :: xcj(:),ycj(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:)
      allocatable :: zs(:,:),ze(:,:),idry(:),idry_s(:),idry_e(:)
      allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),nwild(:),iflux_e(:)
      allocatable :: fluxes(:),suv2(:,:,:),tsel(:,:,:),kbp2(:),kbs2(:),kbe2(:),sigma_lcl(:,:)
      
!      print*, 'Outputs in Cartesian (1) or lat/lon (2)?'
!      read(*,*) ics
      ics=1
      print*, 'Input start and end stack #:'
      read(*,*)iday1,iday2

      print*, 'Do you want to compute T,S fluxes? 0: no; 1: yes'
      read(*,*)itsflux

!...  Header
      open(63,file='1_elev.61',status='old',access='direct',recl=nbyte)
      open(64,file='1_hvel.67',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file='1_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file='1_salt.70',status='old',access='direct',recl=nbyte)
      endif
      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
        print*, '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
      read(63,rec=irec+2) dtout
      read(63,rec=irec+3) nspool
      read(63,rec=irec+4) ivs
      read(63,rec=irec+5) i23d
      irec=irec+5

      print*, 'ivs=',ivs,'; i23d=',i23d,' ;nrec= ',nrec

!     Vertical grid (obsolete; will be overwritten later)
      read(63,rec=irec+1) nvrt
      read(63,rec=irec+2) kz
      read(63,rec=irec+3) h0
      read(63,rec=irec+4) h_s
      read(63,rec=irec+5) h_c
      read(63,rec=irec+6) theta_b
      read(63,rec=irec+7) theta_f
      irec=irec+7

!      do k=1,kz-1
!        read(63,rec=irec+k) ztot(k)
!      enddo
!      do k=kz,nvrt
!        kin=k-kz+1
!        read(63,rec=irec+k) sigma(kin)
!        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 
      read(63,rec=irec+2) ne
      irec=irec+2

      allocate(x(np),y(np),dp(np),kbp(np),kfp(np),elnode(4,ne),idry(np),idry_e(ne), &
     &eta2(np),xctr(ne),yctr(ne),dpe(ne),kbe(ne),ztmp(nvrt,np),ze(nvrt,ne), &
     &tsel(2,nvrt,ne),kbp2(np),kbe2(ne),i34(ne),stat=istat)
      if(istat/=0) stop 'Falied to allocate (2)'

      do m=1,np
        read(63,rec=irec+1)x(m)
        read(63,rec=irec+2)y(m)
        read(63,rec=irec+3)dp(m)
        read(63,rec=irec+4)kbp(m)
        irec=irec+4
      enddo !m=1,np

      do m=1,ne
        read(63,rec=irec+1)i34(m)
        irec=irec+1
        do mm=1,i34(m)
          read(63,rec=irec+1)elnode(mm,m)
          irec=irec+1
        enddo !mm
      enddo !m
      irec0=irec

!     hvel.67
      read(64,rec=irec_tmp+1) ns 
      read(64,rec=irec_tmp+2) ns_e
      print*, 'side grid:',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 'Allocation error: side(8)'
      irec0_u=irec_tmp+2
      do m=1,ns
        read(64,rec=irec0_u+3)dps(m)
        read(64,rec=irec0_u+4)kbs(m)
        irec0_u=irec0_u+4
      enddo !m

!      do m=1,ns_e
!        read(64,rec=irec0_u+1)i34
!        irec0_u=irec0_u+1
!        do mm=1,i34
!          read(64,rec=irec0_u+1)itmp
!          irec0_u=irec0_u+1
!          if(m==ns_e) print*, 'last elem in side grid:',itmp
!        enddo !mm
!      enddo !m

      irec0_u=irec0_u+4*ns_e !pure tri for sidecenters.gr3

!     *.70
      if(itsflux==1) then
        read(65,rec=irec_tmp+1) ne1
        read(65,rec=irec_tmp+2) ne_e
        if(ne1/=ne) stop 'mismatch in elem.'
        print*, 'Elem. grid:',ne1,ne_e
        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_e
      endif !itsflux

!     Compute geometry
      call compute_nside(np,ne,i34,elnode,ns2)
      if(ns/=ns2) stop 'mismatch in sides'
      print*, 'done compute_nside'
      allocate(ic3(4,ne),elside(4,ne),isdel(2,ns),isidenode(2,ns),xcj(ns),ycj(ns),stat=istat)
      if(istat/=0) stop 'Allocation error: side(7)'
      call schism_geometry_single(np,ne,ns2,x,y,i34,elnode,ic3,elside,isdel,isidenode,xcj,ycj)
      print*, 'done schism_geometry'

      print*, 'last element',(elnode(j,ne),j=1,i34(ne))

!     Read vgrid.in
      allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp(np),stat=istat)
      if(istat/=0) stop 'Falied to allocate (1)'
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp)

!     Read in fluxflag.prop
      allocate(iflux_e(ne),stat=istat)
      open(32,file='fluxflag.prop',status='old')
      do i=1,ne
        read(32,*)j,tmp1
        itmp=tmp1
        if(itmp<-1) stop 'fluxflag.prop wrong'
        iflux_e(i)=itmp
      enddo
      close(32)

!     Only compute if diff(region #)=1; fluxes(i) is from i to i-1
      max_reg=maxval(iflux_e)
      allocate(fluxes(max_reg),stat=istat)
      if(istat/=0) stop 'Falied to allocate (6)'
      open(60,file='fluxes.out')
      
!...  Compute relative record # for a node and level for 3D outputs
!...
!      icount=0
!      do i=1,np
!        do k=max0(1,kbp(i)),nvrt
!          do m=1,ivs
!            icount=icount+1
!            icum(i,k,m)=icount
!          enddo !m
!        enddo !k
!      enddo !i=1,np

!...  Time iteration
!...
      it_tot=0
      do iday=iday1,iday2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')iday
      open(63,file=it_char//'_elev.61',status='old',access='direct',recl=nbyte)
      open(64,file=it_char//'_hvel.67',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file=it_char//'_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file=it_char//'_salt.70',status='old',access='direct',recl=nbyte)
      endif

      irec=irec0
      irec_u=irec0_u
      irec_TS=irec0_TS
      do it1=1,nrec
        read(63,rec=irec+1) time
        read(63,rec=irec+2) it
        irec=irec+2
        it_tot=it_tot+1
!        time=it_tot*dtout

        print*, 'time in days =',time/86400

        do i=1,np
          read(63,rec=irec+i) eta2(i)
        enddo !i
        irec=irec+np+np !done with (63 for this step

        !Compute z coordinates
        do i=1,np
          if(ivcor==1) then !localized
            if(dp(i)+eta2(i)<=h0) then
              idry(i)=1
            else !wet
              idry(i)=0
              do k=kbp(i),nvrt
                ztmp(k,i)=(eta2(i)+dp(i))*sigma_lcl(k,i)+eta2(i)
              enddo !k
            endif !wet/dry
          else if(ivcor==2) then !SZ
            call zcor_SZ_single(dp(i),eta2(i),h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot, &
     &sigma,ztmp(:,i),idry(i),kbpl)
            if(kbpl/=kbp(i)) stop 'mismatch (2)'
          endif
        enddo !i=1,np
          
!       idry_e, idry_s
        do i=1,ne      
          idry_e(i)=maxval(idry(elnode(1:i34(i),i)))
          kbe2(i)=minval(kbp(elnode(1:i34(i),i)))
        enddo !i
        idry_s=1 !init.
        do i=1,ns
          kbs2(i)=minval(kbp(isidenode(1:2,i)))
          do j=1,2
            ie=isdel(j,i)
            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(1:2,i)))
            if(itmp/=0) then
              write(*,*)'Wet side has dry nodes:',i,isidenode(1:2,i)
              stop
            endif
          endif
        enddo !i=1,ns

!       hvel.67
        irec_u=irec_u+ns+2
        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
          enddo !k

          !Debug
          !if(iday==ndays.and.it1==nrec) write(97,*)xcj(i),ycj(i),suv2(1:2,nvrt,i)
        enddo !i

        !Compute zcor for wet sides (deal with some inconsistency)
        do i=1,ns
          if(idry_s(i)==1) cycle

          n1=isidenode(1,i); n2=isidenode(2,i)
          do k=kbs2(i),nvrt
            k1=max(k,kbp(n1))
            k2=max(k,kbp(n2))
            zs(k,i)=(ztmp(k1,n1)+ztmp(k2,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
          irec_TS=irec_TS+ne+2
          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(iday==ndays.and.it1==nrec) write(99,*)i,tsel(2,nvrt,i)
          enddo !i

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

            !Compute zcor for wet elem.
            do k=kbe2(i),nvrt
              ze(k,i)=0
              do j=1,i34(i)
                nd=elnode(j,i)
                kl=max(k,kbp(nd))
                ze(k,i)=ze(k,i)+ztmp(kl,nd)/i34(i)
              enddo !j
              if(abs(ze(k,i))>1.e18) then
                write(*,*)'Dry elem:',i
                stop
              endif
            enddo !k

            !Extend valid values to bottom for wet elem.
            !After this all wet elem. have valid values from 1:nvrt
            kk0=-1
            do k=1,nvrt
              if(tsel(1,k,i)>=0) then
                kk0=k; exit
              endif
            enddo !k
            if(kk0<=0) then
              write(*,*)'Weird T:',i,kk0,tsel(1,:,i)
              stop
            else if(tsel(2,kk0,i)<0) then
              write(*,*)'Weird S:',i,kk0,tsel(2,:,i)
              stop
            endif
            do k=1,kk0-1
              tsel(1:2,k,i)=tsel(1:2,kk0,i)
            enddo !k
          enddo !i=1,ne
        endif !itsflux

!       Fluxes
        fluxes=0
        do i=1,ns
          if(isdel(2,i)==0.or.idry_s(i)==1) cycle 

          !Internal wet sides
          ie1=isdel(1,i); ie2=isdel(2,i)
          if(iflux_e(ie1)==-1.or.iflux_e(ie2)==-1) cycle
          if(iabs(iflux_e(ie1)-iflux_e(ie2))/=1) cycle

          n1=isidenode(1,i); n2=isidenode(2,i)
          distj=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2)
          thetan=atan2(x(n1)-x(n2),y(n2)-y(n1))
          snx=cos(thetan); sny=sin(thetan)
          ireg_hi=max(iflux_e(ie1),iflux_e(ie2))
          ireg_lo=min(iflux_e(ie1),iflux_e(ie2))
          if(ireg_hi==iflux_e(ie1)) then
            isgn=1 !'hi' to 'lo' is same as elem. '1' to '2'
          else
            isgn=-1
          endif

          do k=kbs2(i),nvrt-1
            !Normal vel. from 1 to 2
            vnn=(suv2(1,k+1,i)+suv2(1,k,i))/2*snx+(suv2(2,k+1,i)+suv2(2,k,i))/2*sny 
            fluxes(ireg_hi)=fluxes(ireg_hi)+isgn*vnn*distj*(zs(k+1,i)-zs(k,i)) !m^3/s
          enddo !k
        enddo !i=1,ns

!       Output. fluxes(i) is from region i to i-1
        write(60,'(f16.6,60000(1x,e14.4))')time/86400,fluxes(:)
      enddo !it1=1,nrec

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      enddo !iday

      stop
      end