!   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 nc outputs from scribe I/O versions for multiple files at all nodes 
!       Works for mixed tri/quad outputs on NODE/ELEMENT based vars.
!       Inputs: screen; nc file and out2d*.nc; vgrid.in (in this dir or ../)
!       Outputs: extract.out (ascii); optional: max/min/avg 2D output (*_[max|min|avg].gr3 and *_violators.bp)
!****************************************************************************************
!     ifort -O2 -assume byterecl -o read_output10_allnodes.exe ../UtilLib/extract_mod2.f90 ../UtilLib/schism_geometry.f90 ../UtilLib/compute_zcor.f90 read_output10_allnodes.f90 -I$NETCDF/include  -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff

      program read_out
      use netcdf
      use extract_mod2
      use schism_geometry_mod
      use compute_zcor

!      parameter(nbyte=4)
      character(len=30) :: file63,varname,outname(3),file62
      character(len=12) :: it_char
!      character*48 data_format
      logical::lexist
      integer :: nx(4,4,3),dimids(100),idims(100)
      allocatable :: sigma(:),cs(:),ztot(:)
      allocatable :: outvar(:,:),out(:,:,:),icum(:,:,:),eta2(:),ztmp(:,:)
      allocatable :: xcj(:),ycj(:),dp(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:)
      allocatable :: zs(:,:),ze(:,:),idry(:),outs(:,:,:),oute(:,:,:),rstat2d(:,:),rviolator(:)
      allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),kbp(:),sigma_lcl(:,:),kbp00(:)
      integer, allocatable :: elside(:,:),idry_e(:),nne(:),indel(:,:),iviolator(:)
      integer,allocatable :: i34(:),elnode(:,:)
      integer :: char_len,start_2d(2),start_3d(3),start_4d(4), &
     &count_2d(2),count_3d(3),count_4d(4)
      real*8,allocatable :: timeout(:),xnd(:),ynd(:)

      integer :: icomp_stats(3)

      print*, 'Input NODE-based variable name to read from nc (e.g. elevation):'
      read(*,'(a30)')varname
!!'
      print*, '<<<<<var name: ',varname

      print*, 'Is the var node-based or ele-based? 0: node based; 1: element based'
      read(*,*)inode_ele
!!'
      print*, '<<<<<inode_ele: ',inode_ele

      varname=adjustl(varname); len_var=len_trim(varname)
      
      print*, 'Input start and end stack # to read:'
      read(*,*) iday1,iday2
      print*, '<<<<<start and end stack #: ',iday1,iday2

      print*, 'Do you want to compute stats for 2D var? (0/1:min; 0/1:max; 0/1:time avg)'
!!'
      read(*,*) icomp_stats(1),icomp_stats(2),icomp_stats(3)
      print*, '<<<<<icomp_stats: ',icomp_stats
      outname(1)='min';outname(2)='max';outname(3)='avg'

      print*, 'Input a threshold: values below the threshold will be output as a bp file.'
!!'
      read(*,*) thres
      print*, '<<<<<threshold: ',thres

      print*, 'How do you want to initialize the variable values for min/max?'
!!'
      print*, '0: do nothing;  1: intialized to -dp (useful for maxelev)'
!!'
      read(*,*) iinitial
!!
      print*, '<<<<<initialization flag: ',iinitial
      if (iinitial.ne.0 .and. iinitial.ne.1) stop 'wrong initialization flag'

!!'
      if(varname(1:len_var).eq.'elev' .and. iinitial==1) then
        is_elev=1 
        print*, '<<<<<special treatment will be implemented for elev'
      else
        is_elev=0 
      endif

      open(65,file='extract.out')
      
!...  Get basic info from out2d*.nc
      !Returned vars: ne,np,ns,nrec,[xnd ynd dp](np),
      !elnode,i34,nvrt,h0,dtout,kbp
      call get_dims(iday1,np,ne,ns,nvrt,h0)
      allocate(xnd(np),ynd(np),dp(np),kbp(np),i34(ne),elnode(4,ne),stat=istat)
      if(istat/=0) stop 'alloc (1)'
      call readheader(iday1,np,ne,ns,kbp,i34,elnode,nrec,xnd,ynd,dp,dtout)

      print*
      print*, 'After header:',ne,np,nrec,i34(ne),elnode(1:i34(ne),ne),nvrt,h0,xnd(np),ynd(np),dp(np) !,start_time

      ! Compute geometry
      allocate(nne(np))
      do k=3,4
        do i=1,k
          do j=1,k-1
            nx(k,i,j)=i+j
            if(nx(k,i,j)>k) nx(k,i,j)=nx(k,i,j)-k
            if(nx(k,i,j)<1.or.nx(k,i,j)>k) then
              write(*,*)'nx wrong',i,j,k,nx(k,i,j)
              stop
            endif
          enddo !j
        enddo !i
      enddo !k

      nne=0
      do i=1,ne
        do j=1,i34(i)
          nd=elnode(j,i)
          nne(nd)=nne(nd)+1
        enddo
      enddo
      mnei=maxval(nne)

      allocate(indel(mnei,np),stat=istat)
      if(istat/=0) stop 'Failed to alloc. indel'
      nne=0
      do i=1,ne
        do j=1,i34(i)
          nd=elnode(j,i)
          nne(nd)=nne(nd)+1
          if(nne(nd)>mnei) then
            write(*,*)'Too many neighbors',nd
            stop
          endif
          indel(nne(nd),nd)=i
        enddo
      enddo !i

      call compute_nside(np,ne,i34,elnode(1:4,1:ne),ns2)
      if(ns2/=ns) then
        write(*,*)'Mismatch in side:',ns2,ns
        stop
      endif
      allocate(ic3(4,ne),elside(4,ne),isdel(2,ns2),isidenode(2,ns2),xcj(ns2),ycj(ns2),stat=istat)
      if(istat/=0) stop 'Allocation error: side(0)'
      call schism_geometry_single(np,ne,ns2,real(xnd),real(ynd),i34,elnode(1:4,1:ne),ic3(1:4,1:ne), &
     &elside(1:4,1:ne),isdel,isidenode,xcj,ycj)

      !For dimensioning purpose
      if(np>ns.or.ne>ns) stop 'ns is not largest'

      if (inode_ele==0) then
        last_dim=np
      elseif (inode_ele==1) then
        last_dim=ne
      endif

      print*, '<<<<<last_dim: ',last_dim

      allocate(idry_e(ne),rstat2d(3,ne),ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np), &
     &outvar(nvrt,last_dim),eta2(np),idry(np),ztmp(nvrt,np),timeout(nrec), &
     &rviolator(ne),iviolator(ne),kbp00(np))
      outvar=-huge(1.0) !test mem

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

!     Calculate kbp00
      if(ivcor==1) then
        kbp00=kbp
      else
        do i=1,np
          call zcor_SZ_single(dp(i),1.e8,h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma, &
     &ztmp(:,i),idry2,kbp00(i))
        enddo !i
      endif !ivcor
      
!...  Time iteration
!...
      outvar=-99 !init.
      ztmp=-99
      if(iinitial==1) then
        print*,'setting initial elev = -dp for all (dry and wet) nodes'
        rstat2d(1,1:np)=-dp !for elev, min is -h
        rstat2d(2,1:np)=-dp !for elev, min is -h
        rstat2d(3,1:np)=-dp !for elev, min is -h
      else
        if (icomp_stats(1)==1) then  !min
          rstat2d(1,:)=huge(1.0)
        elseif (icomp_stats(2)==1) then !max
          rstat2d(2,:)=-huge(1.0)
        elseif (icomp_stats(3)==1) then !avg
          rstat2d(3,:)=0.0
        endif
      endif

      do iday=iday1,iday2 !1,ndays
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')iday
      it_char=adjustl(it_char)
      leng=len_trim(it_char)
      file62='out2d_'//it_char(1:leng)//'.nc'
      iret=nf90_open(trim(adjustl(file62)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid4)
      !time is double
      iret=nf90_inq_varid(ncid4,'time',itime_id)
      iret=nf90_get_var(ncid4,itime_id,timeout,(/1/),(/nrec/))
      print*, 'time=',timeout !,trim(adjustl(file63))
 
      !Find nc file
      file63=varname(1:len_var)//'_'//it_char(1:leng)//'.nc'
      inquire(file=file63,exist=lexist)
      if(lexist) then
        i23d=2 !3D var
      else
        i23d=1 !2D
        file63=file62
      endif

      iret=nf90_open(trim(adjustl(file63)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid)
      iret=nf90_inq_varid(ncid,varname(1:len_var),ivarid1)
      if(iret/=nf90_NoErr) stop 'Var not found'
      iret=nf90_Inquire_Variable(ncid,ivarid1,ndims=ndims,dimids=dimids)
      if(ndims>100) stop 'increase dimension of dimids & idims'
      do i=1,ndims
       iret=nf90_Inquire_Dimension(ncid,dimids(i),len=idims(i))
      enddo !i
      npes=idims(ndims-1) !np|ne|ns
      if(npes/=np.and.npes/=ne) stop 'can only handle node- or elem-based'
!'
      if(idims(ndims)/=nrec) stop 'last dim is not time'

!        iret=nf90_get_att(ncid,ivarid1,'i23d',i23d)
!        if(i23d<=0.or.i23d>6) stop 'wrong i23d'
!        if(i23d>3.and.ics==2) stop 'ics=2 with elem-based var'
!        iret=nf90_get_att(ncid,ivarid1,'ivs',ivs)
!        !print*, 'i23d:',i23d,ivs,idims(1:ndims)
!
      do irec=1,nrec
        !Get elev
        iret=nf90_inq_varid(ncid4,'elevation',itmp)
        start_2d(1)=1; start_2d(2)=irec
        count_2d(1)=np; count_2d(2)=1
        iret=nf90_get_var(ncid4,itmp,eta2,start_2d,count_2d)

        if(i23d==1) then !2D
          start_2d(1)=1; start_2d(2)=irec
          count_2d(1)=npes; count_2d(2)=1
          iret=nf90_get_var(ncid,ivarid1,outvar(1,1:npes),start_2d,count_2d)
        else !3D
          start_3d(1:2)=1; start_3d(3)=irec
          count_3d(1)=nvrt; count_3d(2)=npes; count_3d(3)=1
          iret=nf90_get_var(ncid,ivarid1,outvar(:,1:npes),start_3d,count_3d)
        endif 

        !Available now: outvar(nvrt,np|ne), eta2(np)
          if(i23d==1) then !2D
!           Output: time, 2D variable at all nodes
            if(maxval(icomp_stats(1:3))/=0) &
     &write(65,'(e14.6,1000000(1x,e14.4))')timeout(irec)/86400,(outvar(1,i),i=1,np)

            !Compute stats (magnitude for vectors)
            if(sum(icomp_stats(:))/=0) then
              if(is_elev==1 .and. icomp_stats(2)==1) then !maxelev

                !Mark wet elem
                do i=1,ne
                  if(minval(outvar(1,elnode(1:i34(i),i))+dp(elnode(1:i34(i),i)))>0) then
                    idry_e(i)=0
                  else
                    idry_e(i)=1
                  endif 
                enddo !i
                do i=1,ne
                  if(idry_e(i)==0) then !wet
                    ifl=1 !isolated wet
                    loop2: do j=1,i34(i)
                      nd=elnode(j,i)
                      do m=1,nne(nd)
                        ie=indel(m,nd)
                        if(i/=ne.and.idry_e(ie)==0) then
                          ifl=0
                          exit loop2
                        endif
                      enddo !m
                    enddo loop2 !j

                    if(ifl==0) then ! not isolated wet
                      do j=1,i34(i)
                        nd=elnode(j,i)
                        tmp=outvar(1,nd)
                        if(tmp+dp(nd)>0) rstat2d(2,nd)=max(rstat2d(2,nd),tmp)
                      enddo
                    endif 
                  endif !idry_e
                enddo !i=1,ne 
              else !other variables, or min_ele or avg_ele
                do i=1,last_dim
                  tmp=outvar(1,i)
                  if (icomp_stats(1)==1) then !min
                    rstat2d(1,i)=min(rstat2d(1,i),tmp)
                  endif
                  if(icomp_stats(2)==1) then !max
                    rstat2d(2,i)=max(rstat2d(2,i),tmp)
                  endif
                  if(icomp_stats(3)==1) then !min
                    rstat2d(3,i)=tmp+rstat2d(3,i)
                  endif
                enddo !i 
              endif !is_elev
            endif ! icomp_stats
          else !if(i23d==3) then !3D 
            !Compute z coordinates
            do i=1,last_dim
              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)
              endif

              do k=max0(1,kbp00(i)),nvrt
!               Output: time, node #, level #, z-coordinate, 3D variable 
                if(idry(i)==1) then
                  write(65,*)timeout(irec)/86400,i,k,-1.e6,outvar(k,i)
                else
                  write(65,*)timeout(irec)/86400,i,k,ztmp(k,i),outvar(k,i)
                endif
              enddo !k
 
              !Debug
              !write(65,*)x(i),y(i),out(i,1,1:ivs) 

            enddo !i
          endif !i23d
      enddo !irec

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

      
      nrec_total=nrec*(iday2-iday1+1)
      if(icomp_stats(3)==1) then
        rstat2d(3,:)=rstat2d(3,:)/real(nrec_total)
      endif
      print*, '<<<<<total rec: ', nrec_total

      do k=1,3
        if(icomp_stats(k)==0) cycle

        if (inode_ele==0) then !node-based, write gr3
          open(12,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'.gr3',status='replace')
          write(12,*)iday1,iday2
          write(12,*)ne,np
          do i=1,np
            write(12,*)i,xnd(i),ynd(i),rstat2d(k,i)
          enddo !i
          do i=1,ne
            write(12,*)i,i34(i),elnode(1:i34(i),i)
          enddo !i
          close(12)
        elseif(inode_ele==1) then !ele based, write prop
          open(12,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'.prop',status='replace')
          do i=1,ne
            write(12,*)i,rstat2d(k,i)
          enddo
        endif

        close(12)
      enddo !k=1,3

      !record violators
      do k=1,3
        nrec=0; rviolator=0.0; iviolator=0
        if (icomp_stats(k)==0) cycle
        do i=1,last_dim
          if (rstat2d(k,i)<thres) then
            nrec=nrec+1
            rviolator(nrec)=rstat2d(k,i) !val
            iviolator(nrec)=i !id
          endif
        enddo
        
        open(13,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'_violators.bp',status='replace')
        write(13,*) 'violators <', thres
        write(13,*) nrec

        if (inode_ele==0) then !node-based
          do i=1,nrec
            ip=iviolator(i)
            write(13,*) i,xnd(ip),ynd(ip),rviolator(i),ip
          enddo
        elseif (inode_ele==1) then !ele based
          do i=1,nrec
            ie=iviolator(i)
            write(13,*) i,sum(xnd(elnode(1:i34(ie),ie)))/i34(ie),&
              & sum(ynd(elnode(1:i34(ie),ie)))/i34(ie), rviolator(i),ie
          enddo
        endif
        close(13)
      enddo !k=1,3

      print*, 'Finished!'

      stop
      end