!   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 hvel. and station.xy. (x,y) read in from (build points) station.xy, and z values
!       read in from the depths in station.xy (not used for 2D variables), and then compute
!       along-channel angles (max. variance) at each point and along-channel vel.
!       Works with mixed tri/quads.
!       Inputs: station.xy; screen inputs
!       Outputs: fort.20 (along channel vel.); fort.18 (channel angle)

!       ifort -Bstatic -assume byterecl -O3 -o compute_alongchannel_vel.exe ../UtilLib/compute_zcor.f90 ../UtilLib/pt_in_poly_test.f90 compute_alongchannel_vel.f90

!											*
!****************************************************************************************
!
      program read_out
      use compute_zcor
      use pt_in_poly_test
      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:: i34(:),elnode(:,:)
      allocatable :: sigma(:),sigma_lcl(:,:),cs(:),ztot(:),x(:),y(:),dp(:),kbp00(:),kbp(:),ztmp2(:,:)
      allocatable :: out(:,:,:,:),out2(:,:,:),icum(:,:,:),eta2(:),node3(:,:),arco(:,:)
      allocatable :: ztmp(:),x00(:),y00(:),iep(:),out3(:,:,:),z00(:),theta_al(:),out4(:,:)
      dimension swild(3)
      integer :: nodel(3)
      
      pi=3.1415926
!      print*, 'Input file to read from (without *_):'
!      read(*,'(a30)')file63
      file63='hvel.64'
      
      print*, 'Input start and end file # to read:'
      read(*,*) iday1,iday2

      print*, 'Input start CORIE day (at T=0):'
      read(*,*) istart_day

      open(10,file='station.xy',status='old')
      read(10,*) 
      read(10,*) nxy
      allocate(x00(nxy),y00(nxy),z00(nxy),theta_al(nxy),stat=istat)
      if(istat/=0) stop 'Falied to allocate (1)'
      do i=1,nxy
        read(10,*)j,x00(i),y00(i),z00(i) !z00 ==0 is MSL; z00<0 is below MSL
      enddo !i
      close(10)

!      open(65,file='extract.out')
!      write(65,*)'(x,y)= ',x00,y00
      
!...  Header
!...
      open(63,file='1_'//file63,status='old',access='direct',recl=nbyte)
      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*, 'i23d=',i23d,' nrec= ',nrec

!     Vertical grid (obsolete)
      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
!
!      allocate(sigma(nvrt),cs(nvrt),ztot(nvrt),ztmp(nvrt),stat=istat)
!      if(istat/=0) stop 'Falied to allocate (2)'
!      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

!     Horizontal grid
      irec=irec+7+nvrt
      read(63,rec=irec+1) np
      read(63,rec=irec+2) ne
      irec=irec+2

      allocate(x(np),y(np),dp(np),kbp00(np),i34(ne),elnode(4,ne), &
     &out(nxy,3,nvrt,2),out2(nxy,nvrt,2),icum(np,nvrt,2),eta2(np),node3(nxy,3), &
     &arco(nxy,3),iep(nxy),stat=istat)
      if(istat/=0) stop 'Falied to allocate (3)'

      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)kbp00(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

      print*, 'last element',elnode(i34(ne),ne)

!     Read in vgrid.in
      allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp(np))
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp)
      allocate(ztmp(nvrt),ztmp2(nvrt,3))

!...  Find parent element for (x00,y00)
      iep=0
      do i=1,ne
        do l=1,nxy
          if(iep(l)/=0) cycle

          call pt_in_poly_single(i34(i),x(elnode(1:i34(i),i)),y(elnode(1:i34(i),i)),x00(l),y00(l),inside,arco(l,1:3),nodel)
          if(inside==1) then
            iep(l)=i
            !print*, 'Found:',l,arco(l,1:3),nodel
            node3(l,1:3)=elnode(nodel(1:3),i)
          endif !inside
        enddo !l; build pts
!          aa=0
!          ar=0 !area
!          do j=1,3
!            j1=j+1
!            j2=j+2
!            if(j1>3) j1=j1-3
!            if(j2>3) j2=j2-3
!            n0=elnode(j,i)
!            n1=elnode(j1,i)
!            n2=elnode(j2,i)
!            swild(j)=signa(x(n1),x(n2),x00(l),y(n1),y(n2),y00(l)) !temporary storage
!            aa=aa+abs(swild(j))
!            if(j==1) ar=signa(x(n1),x(n2),x(n0),y(n1),y(n2),y(n0))
!          enddo !j
!          if(ar<=0) then
!            print*, 'Negative area:',ar
!            stop
!          endif
!          ae=abs(aa-ar)/ar
!          if(ae<=1.e-5) then
!            iep(l)=i
!            node3(l,1:3)=elnode(1:3,i)
!            arco(l,1:3)=swild(1:3)/ar
!            arco(l,1)=max(0.,min(1.,arco(l,1)))
!            arco(l,2)=max(0.,min(1.,arco(l,2)))
!            if(arco(l,1)+arco(l,2)>1) then 
!              arco(l,3)=0
!              arco(l,2)=1-arco(l,1)
!            else
!              arco(l,3)=1-arco(l,1)-arco(l,2)
!            endif
!            cycle
!          endif

        ifl=0 !flag
        do l=1,nxy
          if(iep(l)==0) then
            ifl=1
            exit
          endif
        enddo !l
        if(ifl==0) exit
      enddo !i=1,ne

      do j=1,nxy
        if(iep(j)==0) then
          print*, 'Cannot find a parent for pt:',j,x00(j),y00(j)
          stop
        endif
      enddo !j

!...  Compute relative record # for a node and level for 3D outputs
!...
      icount=0
      do i=1,np
        do k=max0(1,kbp00(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_tot0=(iday1-1)*nrec
      it_tot=it_tot0
      nsteps=(iday2-iday1+1)*nrec
      allocate(out3(nsteps,nxy,2),out4(nsteps,2),stat=istat)
      if(istat/=0) stop 'Failed to allocate (5)'

      out3=0
      do iday=iday1,iday2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')iday
      open(63,file=it_char//'_'//file63,status='old',access='direct',recl=nbyte)

      irec=irec0
      do it1=1,nrec
!        read(63,rec=irec+1) time
        read(63,rec=irec+2) it
        irec=irec+2
        it_tot=it_tot+1
        !if(it_tot-it_tot0>mnit) stop 'Increase mnit'
        time=it_tot*dtout

!        print*, 'time=',time/86400

        do i=1,np
          read(63,rec=irec+i) eta2(i)
        enddo !i
        irec=irec+np

        out2=0
        if(i23d.eq.2) then
          do i=1,nxy
            do j=1,3 !nodes
              nd=node3(i,j)
              do m=1,ivs
                read(63,rec=irec+(nd-1)*ivs+m) tmp
                out2(i,1,m)=out2(i,1,m)+arco(i,j)*tmp
              enddo !m
            enddo !j
          enddo !i
          irec=irec+np*ivs
          write(18,'(e16.8,12(1x,f12.3))')time,(out2(i,1,1),i=1,nxy)
        else !i23d=3 
          do i=1,nxy
            do j=1,3 !nodes
              nd=node3(i,j)
              do k=max0(1,kbp00(nd)),nvrt
                do m=1,ivs
                  read(63,rec=irec+icum(nd,k,m)) out(i,j,k,m)
                enddo !m
              enddo !k
            enddo !j
          enddo !i
          irec=irec+icum(np,nvrt,ivs)

!         Do interpolation
          do i=1,nxy
            etal=0; dep=0; idry=0
            do j=1,3
              nd=node3(i,j)
              if(eta2(nd)+dp(nd)<h0) idry=1
              etal=etal+arco(i,j)*eta2(nd)
              dep=dep+arco(i,j)*dp(nd)
      
!             Debug
!              write(11,*)i,j,nd,dp(nd),arco(i,j)

            enddo !j
            if(idry==1) then
              if(file63(1:7).eq.'hvel.64') then
                out3(it_tot-it_tot0,i,1:2)=0
              else
                out3(it_tot-it_tot0,i,1:2)=-99
              endif
!              write(65,*)'Dry'
            else !element wet
!             Compute z-coordinates
              if(ivcor==1) then !localized
                do j=1,3
                  nd=node3(i,j)
                  do k=kbp(nd)+1,nvrt-1
                    ztmp2(k,j)=(eta2(nd)+dp(nd))*sigma_lcl(k,nd)+eta2(nd)
                  enddo !k
                  ztmp2(kbp(nd),j)=-dp(nd) !to avoid underflow
                  ztmp2(nvrt,j)=eta2(nd) !to avoid underflow
                enddo !j

                ztmp=0
                kbpl=minval(kbp(node3(i,1:3)))
                do k=kbpl,nvrt
                  do j=1,3
                    nd=node3(i,j)
                    ztmp(k)=ztmp(k)+arco(i,j)*ztmp2(max(k,kbp(nd)),j)
                  enddo !j
                enddo !k

              else if(ivcor==2) then !SZ
                call zcor_SZ_single(dep,etal,h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot, &
     &sigma,ztmp(:),idry2,kbpl)
              endif

              if(1==2) then
!-------------------------------------------------------------
              do k=kz,nvrt
                kin=k-kz+1
                hmod2=min(dep,h_s)
                if(hmod2<=h_c) then
                  ztmp(k)=sigma(kin)*(hmod2+etal)+etal
                else if(etal<=-h_c-(hmod2-h_c)*theta_f/sinh(theta_f)) then
                  write(*,*)'Pls choose a larger h_c (2):',etal,h_c
                  stop
                else
                  ztmp(k)=etal*(1+sigma(kin))+h_c*sigma(kin)+(hmod2-h_c)*cs(kin)
                endif

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

              if(dep<=h_s) then
                kbpl=kz
              else !z levels
!               Find bottom index
                kbpl=0
                do k=1,kz-1
                  if(-dep>=ztot(k).and.-dep<ztot(k+1)) then
                    kbpl=k
                    exit
                  endif
                enddo !k
                if(kbpl==0) then
                  write(*,*)'Cannot find a bottom level:',dep,i
                  stop
                endif
                ztmp(kbpl)=-dep
                do k=kbpl+1,kz-1
                  ztmp(k)=ztot(k)
                enddo !k
              endif

              do k=kbpl+1,nvrt
                if(ztmp(k)-ztmp(k-1)<=0) then
                  write(*,*)'Inverted z-level:',etal,dep,ztmp(k)-ztmp(k-1)
                  stop
                endif
              enddo !k
!-------------------------------------------------------------
              endif !1==2
       
              do k=kbpl,nvrt
                do m=1,ivs
                  do j=1,3
                    nd=node3(i,j)
                    kin=max(k,kbp00(nd))
                    out2(i,k,m)=out2(i,k,m)+arco(i,j)*out(i,j,kin,m)
                  enddo !j
                enddo !m
!                write(65,*)i,k,ztmp(k),(out2(i,k,m),m=1,ivs)     
              enddo !k

!             Interplate in vertical
              k0=0
              do k=kbpl,nvrt-1
                if(z00(i)>=ztmp(k).and.z00(i)<=ztmp(k+1)) then
                  k0=k
                  rat=(z00(i)-ztmp(k))/(ztmp(k+1)-ztmp(k))
                  exit
                endif
              enddo !k
              if(k0==0) then
                write(*,*)'Warning: failed to find a vertical level:',it,i,ztmp(kbpl),ztmp(nvrt)
                stop
!                Default for out3 is 0.
              else
                do m=1,ivs
                  out3(it_tot-it_tot0,i,m)=out2(i,k0,m)*(1-rat)+out2(i,k0+1,m)*rat !time index starts from 1
                enddo !m
              endif
            endif !dry/wet
          enddo !i=1,nxy
         
        endif !i23d
      enddo !it1=1,nrec

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

      it_out=it_tot-it_tot0
      if(nsteps/=it_out) stop 'Mismatch in # of records'

!     Along-channel angles
      do i=1,nxy
        av_u=sum(out3(1:it_out,i,1))/it_out
        av_v=sum(out3(1:it_out,i,2))/it_out
        out4(1:it_out,1)=out3(1:it_out,i,1)-av_u
        out4(1:it_out,2)=out3(1:it_out,i,2)-av_v
        var_u=sum(out4(1:it_out,1)*out4(1:it_out,1))/it_out
        var_v=sum(out4(1:it_out,2)*out4(1:it_out,2))/it_out
        var_uv=sum(out4(1:it_out,1)*out4(1:it_out,2))/it_out
        theta_al(i)=0.5*atan2(2*var_uv,var_u-var_v)
        write(18,*)theta_al(i)/pi*180,z00(i)

!       Debug
!        write(99,*)i,z00(i),av_u,av_v,var_u,var_v,var_uv
      enddo !i
      print*, 'done computing angles...'

      do it=1,it_out
        time=(it+it_tot0)*dtout/86400+istart_day
        write(20,'(e16.8,6000(1x,f12.3))')time, &
     &(out3(it,i,1)*cos(theta_al(i))+out3(it,i,2)*sin(theta_al(i)),i=1,nxy)
      enddo !it

      print*, 'Finished!'

      stop
      end

!      function signa(x1,x2,x3,y1,y2,y3)
!!...  Compute signed area formed by pts 1,2,3
!
!      signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2
!
!      return
!      end