text:
*DECK DVNLSD
      SUBROUTINE DVNLSD (Y, YH, LDYH, VSAV, SAVF, EWT, ACOR, IWM, WM,
     1                 F, JAC, PDUM, NFLAG, RPAR, IPAR)
      EXTERNAL F, JAC, PDUM
      DOUBLE PRECISION Y, YH, VSAV, SAVF, EWT, ACOR, WM, RPAR
      INTEGER LDYH, IWM, NFLAG, IPAR
      DIMENSION Y(*), YH(LDYH,*), VSAV(*), SAVF(*), EWT(*), ACOR(*),
     1          IWM(*), WM(*), RPAR(*), IPAR(*)
C-----------------------------------------------------------------------
C Call sequence input -- Y, YH, LDYH, SAVF, EWT, ACOR, IWM, WM,
C                        F, JAC, NFLAG, RPAR, IPAR
C Call sequence output -- YH, ACOR, WM, IWM, NFLAG
C COMMON block variables accessed:
C     /DVOD01/ ACNRM, CRATE, DRC, H, RC, RL1, TQ(5), TN, ICF,
C                JCUR, METH, MITER, N, NSLP
C     /DVOD02/ HU, NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
C
C Subroutines called by DVNLSD: F, DAXPY, DCOPY, DSCAL, DVJAC, DVSOL
C Function routines called by DVNLSD: DVNORM
C-----------------------------------------------------------------------
C Subroutine DVNLSD is a nonlinear system solver, which uses functional
C iteration or a chord (modified Newton) method.  For the chord method
C direct linear algebraic system solvers are used.  Subroutine DVNLSD
C then handles the corrector phase of this integration package.
C
C Communication with DVNLSD is done with the following variables. (For
C more details, please see the comments in the driver subroutine.)
C
C Y          = The dependent variable, a vector of length N, input.
C YH         = The Nordsieck (Taylor) array, LDYH by LMAX, input
C              and output.  On input, it contains predicted values.
C LDYH       = A constant .ge. N, the first dimension of YH, input.
C VSAV       = Unused work array.
C SAVF       = A work array of length N.
C EWT        = An error weight vector of length N, input.
C ACOR       = A work array of length N, used for the accumulated
C              corrections to the predicted y vector.
C WM,IWM     = Real and integer work arrays associated with matrix
C              operations in chord iteration (MITER .ne. 0).
C F          = Dummy name for user supplied routine for f.
C JAC        = Dummy name for user supplied Jacobian routine.
C PDUM       = Unused dummy subroutine name.  Included for uniformity
C              over collection of integrators.
C NFLAG      = Input/output flag, with values and meanings as follows:
C              INPUT
C                  0 first call for this time step.
C                 -1 convergence failure in previous call to DVNLSD.
C                 -2 error test failure in DVSTEP.
C              OUTPUT
C                  0 successful completion of nonlinear solver.
C                 -1 convergence failure or singular matrix.
C                 -2 unrecoverable error in matrix preprocessing
C                    (cannot occur here).
C                 -3 unrecoverable error in solution (cannot occur
C                    here).
C RPAR, IPAR = Dummy names for user's real and integer work arrays.
C
C IPUP       = Own variable flag with values and meanings as follows:
C              0,            do not update the Newton matrix.
C              MITER .ne. 0, update Newton matrix, because it is the
C                            initial step, order was changed, the error
C                            test failed, or an update is indicated by
C                            the scalar RC or step counter NST.
C
C For more details, see comments in driver subroutine.
C-----------------------------------------------------------------------
C
      include "vode.H"
C
C Type declarations for local variables --------------------------------
C
      DOUBLE PRECISION CCMAX, CRDOWN, CSCALE, DCON, DEL, DELP, ONE,
     1     RDIV, TWO, ZERO
      INTEGER I, IERPJ, IERSL, M, MAXCOR, MSBP
C
C Type declaration for function subroutines called ---------------------
C
      DOUBLE PRECISION DVNORM
      PARAMETER(CCMAX = 0.3D0, CRDOWN = 0.3D0, MAXCOR = 3, MSBP = 20)
      PARAMETER(RDIV  = 2.0D0, ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0)

C-----------------------------------------------------------------------
C On the first step, on a change of method order, or after a
C nonlinear convergence failure with NFLAG = -2, set IPUP = MITER
C to force a Jacobian update when MITER .ne. 0.
C-----------------------------------------------------------------------
      IF (JSTART .EQ. 0) NSLP = 0
      IF (NFLAG .EQ. 0) ICF = 0
      IF (NFLAG .EQ. -2) IPUP = MITER
      IF ( (JSTART .EQ. 0) .OR. (JSTART .EQ. -1) ) IPUP = MITER
C If this is functional iteration, set CRATE .eq. 1 and drop to 220
      IF (MITER .EQ. 0) THEN
        CRATE = ONE
        GO TO 220
      ENDIF
C-----------------------------------------------------------------------
C RC is the ratio of new to old values of the coefficient H/EL(2)=h/l1.
C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
C to force DVJAC to be called, if a Jacobian is involved.
C In any case, DVJAC is called at least every MSBP steps.
C-----------------------------------------------------------------------
      DRC = ABS(RC-ONE)
      IF (DRC .GT. CCMAX .OR. NST .GE. NSLP+MSBP) IPUP = MITER
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken.  A convergence test is
C made on the r.m.s. norm of each correction, weighted by the error
C weight vector EWT.  The sum of the corrections is accumulated in the
C vector ACOR(i).  The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
 220  M = 0
      DELP = ZERO
      CALL DCOPY (N, YH(1,1), 1, Y, 1 )
      CALL F (N, TN, Y, SAVF, RPAR, IPAR)
      NFE = NFE + 1
      IF (IPUP .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, the matrix P = I - h*rl1*J is reevaluated and
C preprocessed before starting the corrector iteration.  IPUP is set
C to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
      CALL DVJAC (Y, YH, LDYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, IERPJ,
     1           RPAR, IPAR)
      IPUP = 0
      RC = ONE
      DRC = ZERO
      CRATE = ONE
      NSLP = NST
C If matrix is singular, take error return to force cut in step size. --
      IF (IERPJ .NE. 0) GO TO 430
 250  DO 260 I = 1,N
 260    ACOR(I) = ZERO
C This is a looping point for the corrector iteration. -----------------
 270  IF (MITER .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation.
C-----------------------------------------------------------------------
      DO 280 I = 1,N
 280    SAVF(I) = RL1*(H*SAVF(I) - YH(I,2))
      DO 290 I = 1,N
 290    Y(I) = SAVF(I) - ACOR(I)
      DEL = DVNORM (N, Y, EWT)
      DO 300 I = 1,N
 300    Y(I) = YH(I,1) + SAVF(I)
      CALL DCOPY (N, SAVF, 1, ACOR, 1)
      GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix.  The correction is scaled by the factor
C 2/(1+RC) to account for changes in h*rl1 since the last DVJAC call.
C-----------------------------------------------------------------------
 350  DO 360 I = 1,N
 360    Y(I) = (RL1*H)*SAVF(I) - (RL1*YH(I,2) + ACOR(I))
      CALL DVSOL (WM, IWM, Y, IERSL)
      NNI = NNI + 1
      IF (IERSL .GT. 0) GO TO 410
      IF (METH .EQ. 2 .AND. RC .NE. ONE) THEN
        CSCALE = TWO/(ONE + RC)
        CALL DSCAL (N, CSCALE, Y, 1)
      ENDIF
      DEL = DVNORM (N, Y, EWT)
      CALL DAXPY (N, ONE, Y, 1, ACOR, 1)
      DO 380 I = 1,N
 380    Y(I) = YH(I,1) + ACOR(I)
C-----------------------------------------------------------------------
C Test for convergence.  If M .gt. 0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C-----------------------------------------------------------------------
 400  IF (M .NE. 0) CRATE = MAX(CRDOWN*CRATE,DEL/DELP)
      DCON = DEL*MIN(ONE,CRATE)/TQ(4)
      IF (DCON .LE. ONE) GO TO 450
      M = M + 1
      IF (M .EQ. MAXCOR) GO TO 410
      IF (M .GE. 2 .AND. DEL .GT. RDIV*DELP) GO TO 410
      DELP = DEL
      CALL F (N, TN, Y, SAVF, RPAR, IPAR)
      NFE = NFE + 1
      GO TO 270
C
 410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
      ICF = 1
      IPUP = MITER
      GO TO 220
C
 430  CONTINUE
      NFLAG = -1
      ICF = 2
      IPUP = MITER
      RETURN
C
C Return for successful step. ------------------------------------------
 450  NFLAG = 0
      JCUR = 0
      ICF = 0
      IF (M .EQ. 0) ACNRM = DEL
      IF (M .GT. 0) ACNRM = DVNORM (N, ACOR, EWT)
      RETURN
C----------------------- End of Subroutine DVNLSD ----------------------
      END


text:
Address(Rutgers Drive) is a residential street in Central Davis.

 Intersecting Streets

  Anderson Road  at a crossing for Cesar Chavez Elementary School
  Oeste Drive
  Reed Drive
  Redwood Lane

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
  !! Local Data - Offsets
  INTEGER, DIMENSION(:), ALLOCATABLE :: send_per_process
  INTEGER, DIMENSION(:), ALLOCATABLE :: send_offsets
  INTEGER, DIMENSION(:), ALLOCATABLE :: recv_per_process
  INTEGER, DIMENSION(:), ALLOCATABLE :: recv_offsets
  !! Local Data - Send/Recv Buffers
  INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_row
  INTEGER, DIMENSION(:), ALLOCATABLE :: send_buffer_col
  INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_row
  INTEGER, DIMENSION(:), ALLOCATABLE :: recv_buffer_col
  !! ETC
  INTEGER :: num_processes
  INTEGER :: counter, inner_counter, insert_pt
  INTEGER :: mpi_error

  !! Allocate Size Buffers
  CALL MPI_COMM_SIZE(comm, num_processes, mpi_error)
  ALLOCATE(send_per_process(num_processes))
  ALLOCATE(send_offsets(num_processes))
  ALLOCATE(recv_per_process(num_processes))
  ALLOCATE(recv_offsets(num_processes))

  !! Figure Out How Much Data Gets Sent
  DO counter = 1, num_processes
     send_per_process(counter) = triplet_lists(counter)%CurrentSize
  END DO
  send_offsets(1) = 0
  DO counter = 2, num_processes
     send_offsets(counter) = send_offsets(counter-1) + &
          & send_per_process(counter-1)
  END DO

  !! Figure Out How Much Data Gets Received
  CALL MPI_ALLTOALL(send_per_process, 1, MPINTINTEGER, recv_per_process, 1, &
       & MPINTINTEGER, comm, mpi_error)
  recv_offsets(1) = 0
  DO counter = 2, num_processes
     recv_offsets(counter) = recv_offsets(counter-1) + &
          & recv_per_process(counter-1)
  END DO

  !! Allocate And Fill Send Buffers
  ALLOCATE(send_buffer_row(SUM(send_per_process)))
  ALLOCATE(send_buffer_col(SUM(send_per_process)))
  ALLOCATE(send_buffer_val(SUM(send_per_process)))
  ALLOCATE(recv_buffer_row(SUM(recv_per_process)))
  ALLOCATE(recv_buffer_col(SUM(recv_per_process)))
  ALLOCATE(recv_buffer_val(SUM(recv_per_process)))

  !! Fill Send Buffer
  insert_pt = 1
  DO counter = 1, num_processes
     DO inner_counter = 1, triplet_lists(counter)%CurrentSize
        CALL GetTripletAt(triplet_lists(counter), inner_counter, temp_triplet)
        send_buffer_row(insert_pt) = temp_triplet%index_row
        send_buffer_col(insert_pt) = temp_triplet%index_column
        send_buffer_val(insert_pt) = temp_triplet%point_value
        insert_pt = insert_pt + 1
     END DO
  END DO

  !! Do Actual Send
  CALL MPI_Alltoallv(send_buffer_col, send_per_process, send_offsets, &
       & MPINTINTEGER, recv_buffer_col, recv_per_process, recv_offsets, &
       & MPINTINTEGER, comm, mpi_error)
  CALL MPI_Alltoallv(send_buffer_row, send_per_process, send_offsets, &
       & MPINTINTEGER, recv_buffer_row, recv_per_process, recv_offsets, &
       & MPINTINTEGER, comm, mpi_error)
  CALL MPI_Alltoallv(send_buffer_val, send_per_process, send_offsets, &
       & MPIDATATYPE, recv_buffer_val, recv_per_process, recv_offsets, &
       & MPIDATATYPE, comm, mpi_error)

  !! Unpack Into The Output Triplet List
  CALL ConstructTripletList(local_data_out, size_in=SUM(recv_per_process))
  DO counter = 1, SUM(recv_per_process)
     local_data_out%data(counter)%index_column = recv_buffer_col(counter)
     local_data_out%data(counter)%index_row = recv_buffer_row(counter)
     local_data_out%data(counter)%point_value = recv_buffer_val(counter)
  END DO

  !! Cleanup
  DEALLOCATE(send_per_process)
  DEALLOCATE(send_offsets)
  DEALLOCATE(recv_per_process)
  DEALLOCATE(recv_offsets)
  DEALLOCATE(send_buffer_row)
  DEALLOCATE(send_buffer_col)
  DEALLOCATE(send_buffer_val)
  DEALLOCATE(recv_buffer_row)
  DEALLOCATE(recv_buffer_col)
  DEALLOCATE(recv_buffer_val)

max_stars_repo_name: shoshijak/NTPoly

max_issues_repo_name: omni-compiler/NTPoly

max_forks_repo_name: omni-compiler/NTPoly


text:
	SUBROUTINE DSSPCL ( szspcl, ispwid, size, jspwid, iret )
C************************************************************************
C* DSSPCL								*
C*									*
C* This subroutine sets the special symbol parameters.  If these	*
C* parameters are not positive, no changes are made.			*
C*									*
C* DSSPCL ( SZSPCL, ISPWID, SIZE, JSPWID, IRET )			*
C*									*
C* Input parameters:							*
C*	SZSPCL		REAL		Special symbol size		*
C*	ISPWID		INTEGER		Special symbol line width	*
C*									*
C* Output parameters:							*
C*	SIZE		REAL		Special symbol size		*
C*	JSPWID		INTEGER		Special symbol line width	*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. Linda/GSC		 8/96	Based on DSWTHR				*
C* S. Jacobs/NCEP	 3/97	Added check for VG driver		*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVACT.CMN'
	INCLUDE		'DEVCHR.CMN'
C-----------------------------------------------------------------------
	iret = NORMAL
C
	IF ( szspcl .gt. 0. ) tsprsz = szspcl
	IF ( ispwid .gt. 0  ) mspwid = ispwid
	size   = tsprsz
	jspwid = mspwid
C
C*	If the driver is VG, send the attributes to the device
C*	for output.
C
	IF  ( ddev .eq. 'VG' )  THEN
	    CALL HSSYMB ( 7, 0, size, jspwid, iret )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module m_Read
  use m_Useful
  implicit none
  private
!
  public :: InitRead
  public :: ReadFrames
!
contains
!
  subroutine InitRead (file, mxz, unitout)
    character (len=*), parameter :: myname = "InitRead"
    character (len=k_mw), intent (inout) :: file
    type (mxzType), intent (inout) :: mxz
    integer, intent (inout) :: unitout
!
    integer :: errno, n, k, frames
    character (len=k_ml) :: line
    integer :: unitinp
!
    unitinp = GetUnit ()
    open (unit=unitinp, file=trim(file), status="old", action="read", iostat=errno)
!
    if (errno /= 0) then
      call error ("Error reading file "//trim(file), myname, .true., unitout)
    end if
! determine no of atoms and no of frames
    read (unitinp,*, iostat=errno) n
    if (errno /= 0) then
      call error ("Error reading no of atoms from file "//trim(file), myname, .true., unitout)
    end if
    mxz%natoms = n
    k = 1
    do while (GetLine(unitinp, line))
      k = k + 1
    end do
    if (Mod(k, n+2) == 0) then
      mxz%nframes = k / (n+2)
    else
      call error ("Non integer number of frames in file "//trim(file), myname, .true., unitout)
    end if
    close (unitinp)
  end subroutine InitRead
!
!
  subroutine ReadFrames (file, mxz, unitout)
    character (len=*), parameter :: myname = "ReadFrames"
    character (len=k_mw), intent (inout) :: file
    type (mxzType), intent (inout) :: mxz
    integer, intent (inout) :: unitout
!
    integer :: errno, err, i, n, j
    character (len=k_ml) :: line
    character (len=k_mw) :: saux, a, b, c
    real (k_pr) :: x, y, z, dx, dy, dz, q, tdx, tdy, tdz, charge
    integer :: unitinp
!
    unitinp = GetUnit ()
    open (unit=unitinp, file=trim(file), status="old", action="read", iostat=errno)
    do i = 1, mxz%nframes
      read (unitinp,*, iostat=errno) n
      write (saux, '(i0)') i
      if (errno /= 0) then
        call error ("Error reading first line of frame "//trim(saux), myname, .true., unitout)
      end if
      read (unitinp, fmt='(a)', iostat=errno) a
      if (errno /= 0) then
        call error ("Error reading label of frame "//trim(saux), myname, .true., unitout)
      end if
      mxz%frames(i)%frameLabel = trim (a)
      read (a,*, iostat=errno) b, c, q
      if (errno /= 0) then
        call error ("could non determine the timestamp of frame "//trim(saux), myname, .false., unitout)
      end if
      mxz%frames(i)%timestamp = q
      err = 0
      do j = 1, mxz%natoms
        read (unitinp,*, iostat=errno) mxz%frames(i)%element(j), mxz%frames(i)%x(j), mxz%frames(i)%y(j), mxz%frames(i)%z(j), &
       & mxz%frames(i)%q(j), mxz%frames(i)%dx(j), mxz%frames(i)%dy(j), mxz%frames(i)%dz(j)
        err = err + Abs (errno)
      end do
      if (err /= 0) then
        call error ("error readind atoms data from frame "//trim(saux), myname, .true., unitout)
      end if
!! now compute the dipole it may seem like a waste of loops to do it here and not in previous loop but thing what mwy happen it
!! reading fails
      tdx = 0.0_k_pr
      tdy = 0.0_k_pr
      tdz = 0.0_k_pr
      charge = 0.0_k_pr
!$omp  parallel  private(j)
!$omp do schedule(static) reduction(+ : tdx, tdy,tdz,charge)
      do j = 1, mxz%natoms
        tdx = mxz%frames(i)%dx(j) + tdx
        tdy = mxz%frames(i)%dy(j) + tdy
        tdz = mxz%frames(i)%dz(j) + tdz
        charge = charge + mxz%frames(i)%q(j)
      end do
!$omp end do
!$omp end parallel
      mxz%frames(i)%tdx = tdx
      mxz%frames(i)%tdy = tdy
      mxz%frames(i)%tdz = tdz
      mxz%frames(i)%charge = charge
      mxz%frames(i)%dipole = Sqrt (mxz%frames(i)%tdx**2+mxz%frames(i)%tdy**2+mxz%frames(i)%tdz**2)
      if (Abs(mxz%frames(i)%dipole) > tiny(1.0_k_pr)) then
        mxz%frames(i)%tdx = mxz%frames(i)%tdx / mxz%frames(i)%dipole
        mxz%frames(i)%tdy = mxz%frames(i)%tdy / mxz%frames(i)%dipole
        mxz%frames(i)%tdz = mxz%frames(i)%tdz / mxz%frames(i)%dipole
      end if
      mxz%frames(i)%projDipole = mxz%frames(i)%dipole * &
     & (mxz%frames(1)%tdx*mxz%frames(i)%tdx+mxz%frames(1)%tdy*mxz%frames(i)%tdy+mxz%frames(1)%tdz*mxz%frames(i)%tdz)
    end do
    close (unitinp)
  end subroutine ReadFrames
!
!> \brief logical function reads a line from a unit
!> \details
!> \author Alin M Elena
!> \date 14th of January 2006
!> \param uno integer, unit number from where to read the line
!> \param line character(len=*), the line that was read
!> \return .true. if successfull in reading the line, .false. otherwise
!
  logical function GetLine (uno, line)
    character (len=*), parameter :: myname = 'GetLine'
    character (len=k_ml), intent (out) :: line
    integer, intent (in) :: uno
    integer :: errno
!
    inquire (unit=uno, iostat=errno)
    GetLine = .false.
    if (errno /= 0) then
      write (*,*) "Unexpected error opening the input file(s)"
      stop
    end if
    read (uno, fmt='(a)', iostat=errno) line
    if (errno == 0) GetLine = .true.
  end function GetLine
!
end module m_Read

max_stars_repo_name: naumruso/tdtbuj

max_issues_repo_name: naumruso/tdtbuj

max_forks_repo_name: naumruso/tdtbuj


!NORMAL, MATSUBARA GREEN'S FUNCTIONS
subroutine ed_get_gimp_matsubara_1(Gmats)
  complex(8),dimension(Nspin,Nspin,Norb,Norb,Lmats),intent(inout) :: Gmats
  Gmats(:,:,:,:,:) = impGmats(:,:,:,:,:)
end subroutine ed_get_gimp_matsubara_1

subroutine ed_get_gimp_matsubara_2(Gmats)
  complex(8),dimension(Nspin*Norb,Nspin*Norb,Lmats),intent(inout) :: Gmats
  integer  :: io,jo,iorb,jorb,ispin,jspin
  do ispin=1,Nspin
     do jspin=1,Nspin
        do iorb=1,Norb
           do jorb=1,Norb
              io = iorb + (ispin-1)*Norb
              jo = jorb + (jspin-1)*Norb
              Gmats(io,jo,:) = impGmats(ispin,jspin,iorb,jorb,:)
           enddo
        enddo
     enddo
  enddo
end subroutine ed_get_gimp_matsubara_2

subroutine ed_get_gimp_matsubara_3(Gmats,ispin,jspin,iorb,jorb)
  complex(8),dimension(Lmats),intent(inout) :: Gmats
  integer                                   :: iorb,jorb,ispin,jspin
  Gmats(:) = impGmats(ispin,jspin,iorb,jorb,:)
end subroutine ed_get_gimp_matsubara_3

subroutine ed_get_gimp_matsubara_lattice_1(Gmats,Nsites)
  integer                                                                :: Nsites
  complex(8),dimension(Nsites,Nspin,Nspin,Norb,Norb,Lmats),intent(inout) :: Gmats
  Gmats(1:Nsites,:,:,:,:,:) = Gmatsii(1:Nsites,:,:,:,:,:)
end subroutine ed_get_gimp_matsubara_lattice_1

subroutine ed_get_gimp_matsubara_lattice_2(Gmats,Nsites)
  integer                                                                :: Nsites
  complex(8),dimension(Nsites,Nspin*Norb,Nspin*Norb,Lmats),intent(inout) :: Gmats
  integer                                                                :: io,jo,iorb,jorb,ispin,jspin,ilat
  do ilat=1,Nsites
     do ispin=1,Nspin
        do jspin=1,Nspin
           do iorb=1,Norb
              do jorb=1,Norb
                 io = iorb + (ispin-1)*Norb
                 jo = jorb + (jspin-1)*Norb
                 Gmats(ilat,io,jo,:) = Gmatsii(ilat,ispin,jspin,iorb,jorb,:)
              enddo
           enddo
        enddo
     enddo
  enddo
end subroutine ed_get_gimp_matsubara_lattice_2

subroutine ed_get_gimp_matsubara_lattice_3(Gmats,Nsites,ispin,jspin,iorb,jorb)
  integer                                          :: Nsites
  complex(8),dimension(Nsites,Lmats),intent(inout) :: Gmats
  integer                                          :: iorb,jorb,ispin,jspin
  Gmats(1:Nsites,:) = Gmatsii(1:Nsites,ispin,jspin,iorb,jorb,:)
end subroutine ed_get_gimp_matsubara_lattice_3

subroutine ed_get_gimp_matsubara_lattice_11(Gmats,ilat)
  integer                                                         :: ilat
  complex(8),dimension(Nspin,Nspin,Norb,Norb,Lmats),intent(inout) :: Gmats
  Gmats(:,:,:,:,:) = Gmatsii(ilat,:,:,:,:,:)
end subroutine ed_get_gimp_matsubara_lattice_11

subroutine ed_get_gimp_matsubara_lattice_21(Gmats,ilat)
  integer                                                         :: ilat
  complex(8),dimension(Nspin*Norb,Nspin*Norb,Lmats),intent(inout) :: Gmats
  integer                                                         :: io,jo,iorb,jorb,ispin,jspin
  do ispin=1,Nspin
     do jspin=1,Nspin
        do iorb=1,Norb
           do jorb=1,Norb
              io = iorb + (ispin-1)*Norb
              jo = jorb + (jspin-1)*Norb
              Gmats(io,jo,:) = Gmatsii(ilat,ispin,jspin,iorb,jorb,:)
           enddo
        enddo
     enddo
  enddo
end subroutine ed_get_gimp_matsubara_lattice_21

subroutine ed_get_gimp_matsubara_lattice_31(Gmats,ilat,ispin,jspin,iorb,jorb)
  integer                                   :: ilat
  complex(8),dimension(Lmats),intent(inout) :: Gmats
  integer                                   :: iorb,jorb,ispin,jspin
  Gmats(:) = Gmatsii(ilat,ispin,jspin,iorb,jorb,:)
end subroutine ed_get_gimp_matsubara_lattice_31


module IOClass
    use iso_fortran_env
    use MathClass
    use StringClass
    implicit none

    integer(int32),parameter :: PF_JSON=1
    integer(int32),parameter :: PF_CSV=2

    type :: IO_
        integer :: fh=100
        logical :: active=.false.
        logical :: EOF=.true.
        character(1) :: state
        character(200)::path,name,extention
        character(:),allocatable:: title
        character(:),allocatable:: xlabel,ylabel,zlabel
        character(:),allocatable :: filename
    contains
        procedure,public :: unit => unitIO

        procedure,public :: numLine => numLineIO

        procedure,public :: flush => flushIO

        procedure,pass   :: openIOchar
        procedure,pass   :: openIOstring

        procedure, pass :: parseIOChar200
        procedure, pass :: parseIO2keysChar200
        generic,public :: parse =>parseIOChar200,parseIO2keysChar200

        generic,public :: open => openIOchar, openIOstring
        !procedure,public :: open => openIO
        procedure,pass :: writeIOchar,writeIOcharchar,writeIOcharcharchar
        procedure,pass :: writeIOstring,writeIOstringstring,writeIOstringstringstring

        ! int-char-int
        procedure,pass :: writeIOint32re64

        procedure,pass :: writeIOint32re64vector
        procedure,pass :: writeIOint32int32vector

        procedure,pass :: writeIOint32
        procedure,pass :: writeIOint32int32
        procedure,pass :: writeIOint32int32int32
        procedure,pass :: writeIOint32int32int32int32
        procedure,pass :: writeIOint32int32int32int32int32
        procedure,pass :: writeIOint32int32int32int32int32int32

        procedure,pass :: writeIOint32Vector
        procedure,pass :: writeIOint32Array

        procedure,pass :: writeIOre64
        procedure,pass :: writeIOre64re64
        procedure,pass :: writeIOre64re64re64
        procedure,pass :: writeIOre64re64re64re64
        procedure,pass :: writeIOre64re64re64re64re64
        procedure,pass :: writeIOre64re64re64re64re64re64

        procedure,pass :: writeIOre64Vector
        procedure,pass :: writeIOre64Array
        procedure,pass :: writeIOcomplex64
        procedure,pass :: writeIOcomplex64Vector
        procedure,pass :: writeIOcomplex64Array
        generic,public :: write => writeIOchar,writeIOstring,writeIOre64,writeIOre64Vector,writeIOre64Array,&
            writeIOint32,writeIOint32Vector,writeIOint32Array,&
            writeIOre64re64,writeIOre64re64re64,writeIOre64re64re64re64,&
            writeIOre64re64re64re64re64,writeIOre64re64re64re64re64re64,&
            writeIOint32int32,writeIOint32int32int32,writeIOint32int32int32int32,&
            writeIOint32int32int32int32int32,writeIOint32int32int32int32int32int32,&
            writeIOstringstring,writeIOstringstringstring,&
            writeIOcharchar,writeIOcharcharchar,&
            writeIOint32re64vector,writeIOint32int32vector,&
            writeIOint32re64,&
            writeIOcomplex64,writeIOcomplex64Vector,writeIOcomplex64Array
        !procedure,public :: write => writeIO
        procedure,pass :: readIOchar
        procedure,pass :: readIOInt
        procedure,pass :: readIOIntVector
        procedure,pass :: readIOIntArray
        procedure,pass :: readIOReal64
        procedure,pass :: readIOReal64Vector
        procedure,pass :: readIOReal64Array
        generic,public :: read => readIOchar,readIOInt,readIOIntVector,readIOIntArray&
            ,readIOReal64,readIOReal64Vector,readIOReal64Array

        procedure,public :: plot => plotIO
        procedure,public :: splot => splotIO

        procedure,public :: readline => readlineIO
        procedure,public :: close => closeIO
    end type

    interface print
        module procedure printChar, printReal64, printReal32, printInt64, printInt32
    end interface print

    interface disp
        module procedure printChar, printReal64, printReal32, printInt64, printInt32
    end interface disp

    interface plot
        module procedure plotRealArray
    end interface

    interface spy
        module procedure spyRealArray
    end interface
contains

! ===========================================
subroutine readIOInt(obj,val)
    class(IO_),intent(in) :: obj
    integer(int32),intent(inout) :: val

    read(obj%fh,*) val

end subroutine
! ===========================================

! ===========================================
subroutine readIOIntVector(obj,val)
    class(IO_),intent(in) :: obj
    integer(int32),intent(inout) :: val(:)

    read(obj%fh,*) val(:)

end subroutine
! ===========================================

! ===========================================
subroutine readIOIntArray(obj,val)
    class(IO_),intent(in) :: obj
    integer(int32),intent(inout) :: val(:,:)

    read(obj%fh,*) val(:,:)

end subroutine
! ===========================================

! ===========================================
subroutine readIOReal64(obj,val)
    class(IO_),intent(in) :: obj
    real(real64),intent(inout) :: val

    read(obj%fh,*) val

end subroutine
! ===========================================

! ===========================================
subroutine readIOReal64Vector(obj,val)
    class(IO_),intent(in) :: obj
    real(real64),intent(inout) :: val(:)

    read(obj%fh,*) val(:)

end subroutine
! ===========================================

! ===========================================
subroutine readIOReal64Array(obj,val)
    class(IO_),intent(in) :: obj
    real(real64),intent(inout) :: val(:,:)

    read(obj%fh,*) val(:,:)

end subroutine
! ===========================================

function numLineIO(obj,name) result(line)
    class(IO_),intent(inout) :: obj
    type(IO_) :: f
    character(*),intent(in) :: name
    integer(int32) :: line
    character(len=1) :: content

    call f%open(name)

    line=1
    do
        content = f%readline()
        if(f%EOF .eqv. .true.) then
            line=line-1
            exit
        endif
        line = line+1
    enddo

    call f%close()

end function

! #############################################
function readlineIO(obj) result(ret)
    class(IO_),intent(inout) :: obj
    character(len=:),allocatable :: ret

    if(obj%EOF .eqv. .true.)then
        print *, "ERROR :: file is not opened or EOF"
        allocate(character(len=30000) :: ret )
        ret = " "
        return
    endif

    allocate(character(len=30000) :: ret )
    read(obj%fh,'(A)',end=100) ret
    ret = trim(adjustl(ret) )
    return

100 ret = " "
    obj%EOF =.true.
end function
! #############################################

! #############################################
function unitIO(obj) result(unit)
    class(IO_),intent(inout) :: obj
    integer(int32) :: unit
    unit=obj%fh
end function
! #############################################

! #############################################
subroutine openIOchar(obj,path,state,name,extention,fh)
    class(IO_),intent(inout) :: obj
    character(1),optional,intent(in) :: state ! w or r
    character(*),optional,intent(in)::path,name,extention
    integer(int32),optional,intent(in) :: fh
    logical :: yml=.False.

    if(present(state) )then
        obj%state = state
    else
        obj%state="w"
    endif

!    if(present(extention) )then
!        if( trim(extention) == "yml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".yml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".YML" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".yaml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".YAML" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == "yaml" )then
!            yml=.True.
!        endif
!    endif
!    if( index(path,"yml") /= 0 )then
!        yml=.True.
!    endif
!    if( index(path,"yaml") /= 0 )then
!        yml=.True.
!    endif
!
!    if(yml .eqv. .true.)then
!
!        return
!    endif

    if(obj%active .eqv. .true.)then

        print *, "ERROR :: "//trim(obj%path)//trim(obj%name)//trim(obj%extention)//" is already opened."
        stop
    endif

    obj%active=.true.
    if(present(fh) )then
        obj%fh=fh
    else
        obj%fh=10
    endif

    obj%path="./"
    obj%name="untitled"
    obj%name=".txt"

    if(present(path) )then
        obj%path=trim(path)
        if(present(name) )then
            obj%name=trim(name)
            if(present(extention) )then
                obj%extention=trim(extention)
                if(present(state) )then
                    if(state=="r")then
                        open(newunit=obj%fh,file=trim(path)//trim(name)//trim(extention),status='old')
                        obj%filename=trim(path)//trim(name)//trim(extention)
                    elseif(state=="w")then
                        open(newunit=obj%fh,file=trim(path)//trim(name)//trim(extention),status='replace')
                        obj%filename = trim(path)//trim(name)//trim(extention)
                    else
                        call print("Error :: IOClass % open >> argument <state> should be w or r ")
                        stop
                    endif
                else
                    open(newunit=obj%fh,file=trim(path)//trim(name)//trim(extention) )
                    obj%filename = trim(path)//trim(name)//trim(extention)
                endif
            else
                if(present(state) )then
                    if(state=="r")then
                        open(newunit=obj%fh,file=trim(path)//trim(name),status='old' )
                        obj%filename = trim(path)//trim(name)
                    elseif(state=="w")then
                        open(newunit=obj%fh,file=trim(path)//trim(name),status='replace' )
                        obj%filename=trim(path)//trim(name)
                    else
                        call print("Error :: IOClass % open >> argument <state> should be w or r ")
                        stop
                    endif
                else
                    open(newunit=obj%fh,file=trim(path)//trim(name) )
                    obj%filename = trim(path)//trim(name)
                endif
            endif
        else
            if(present(state) )then
                if(state=="r")then
                    open(newunit=obj%fh,file=trim(path),status='old' )
                    obj%filename = trim(path)
                elseif(state=="w")then
                    open(newunit=obj%fh,file=trim(path),status='replace' )
                    obj%filename = trim(path)
                else
                    call print("Error :: IOClass % open >> argument <state> should be w or r ")
                    stop
                endif
            else
                open(newunit=obj%fh,file=trim(path) )
                obj%filename = trim(path)
            endif
        endif
    else
        open(newunit=obj%fh,file="./untitled.txt",status="replace" )
        obj%filename = "./untitled.txt"
    endif

    obj%EOF = .false.

end subroutine openIOchar
! #############################################

! #############################################
subroutine openIOstring(obj,path_s,state,name_s,extention_s,fh)
    class(IO_),intent(inout) :: obj
    character(1),optional,intent(in) :: state ! w or r
    type(String_),intent(in) ::path_s
    type(String_),optional,intent(in)::name_s,extention_s
    character(len=:),allocatable::path,name,extention
    integer(int32),optional,intent(in) :: fh
    logical :: yml=.False.

    if(present(state) )then
        obj%state = state
    else
        obj%state="w"
    endif
!    if(present(extention) )then
!        if( trim(extention) == "yml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".yml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".YML" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".yaml" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == ".YAML" )then
!            yml=.True.
!        endif
!    endif
!    if(present(extention) )then
!        if( trim(extention) == "yaml" )then
!            yml=.True.
!        endif
!    endif
!    if( index(path,"yml") /= 0 )then
!        yml=.True.
!    endif
!    if( index(path,"yaml") /= 0 )then
!        yml=.True.
!    endif
!
!    if(yml .eqv. .true.)then
!
!        return
!    endif

    if(obj%active .eqv. .true.)then

        print *, "ERROR :: "//trim(obj%path)//trim(obj%name)//trim(obj%extention)//" is already opened."
        stop
    endif

    obj%active=.true.
    if(present(fh) )then
        obj%fh=fh
    else
        obj%fh=10
    endif

    obj%path="./"
    obj%name="untitled"
    obj%name=".txt"

    if(present(state) )then
        if(state == "w")then
            open(newunit=obj%fh,file=trim(path_s%str()),status="replace" )
            obj%filename = trim(path_s%str())
        elseif(state == "r")then
            open(newunit=obj%fh,file=trim(path_s%str()),status="old" )
            obj%filename = trim(path_s%str())
        else
            call print("Error :: IOClass % open >> argument <state> should be w or r ")
        endif
    endif

    obj%path=trim(path_s%str() )
    if(present(name_s) )then
        obj%name=trim(name_s%str())
        if(present(extention_s) )then
            obj%extention=trim(extention_s%str())
            open(newunit=obj%fh,file=trim(path_s%str())//trim(name_s%str())//trim(extention_s%str()) )
            obj%filename = trim(path_s%str())//trim(name_s%str())//trim(extention_s%str())
        else
            open(newunit=obj%fh,file=trim(path_s%str())//trim(name_s%str()) )
            obj%filename = trim(path_s%str())//trim(name_s%str())
        endif
    else
        open(newunit=obj%fh,file=trim(path_s%str()) )
        obj%filename = trim(path_s%str())
    endif

    obj%EOF = .false.

end subroutine openIOstring
! #############################################

! #############################################
subroutine writeIOchar(obj,char)
    class(IO_),intent(inout) :: obj
    character(*),intent(in) :: char

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') char

end subroutine writeIOchar
! #############################################

! #############################################
subroutine writeIOcharchar(obj,char1,char2)
    class(IO_),intent(inout) :: obj
    character(*),intent(in) :: char1,char2

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') char1//" "//char2

end subroutine writeIOcharchar
! #############################################

! #############################################
subroutine writeIOcharcharchar(obj,char1,char2,char3)
    class(IO_),intent(inout) :: obj
    character(*),intent(in) :: char1,char2,char3

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') char1//" "//char2//" "//char3

end subroutine writeIOcharcharchar
! #############################################

! #############################################
subroutine writeIOint32(obj,in32)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32))

end subroutine writeIOint32
! #############################################

! #############################################
subroutine writeIOint32re64(obj,in32,re64)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32
    real(real64),intent(in) :: re64

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32))//" "//trim(str(re64))

end subroutine writeIOint32re64
! #############################################

! #############################################

subroutine writeIOint32int32(obj,in32_1,in32_2)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32_2

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32_1))//" "//trim(str(in32_2))

end subroutine

! ####################################################

! #############################################

subroutine writeIOint32int32int32(obj,in32_1,in32_2,in32_3)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32_2,in32_3

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32_1))//" "//trim(str(in32_2))//" "//trim(str(in32_3))

end subroutine

! ####################################################

! #############################################

subroutine writeIOint32int32int32int32(obj,in32_1,in32_2,in32_3,in32_4)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32_2,in32_3,in32_4

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32_1))//" "//trim(str(in32_2))//" "//trim(str(in32_3))//" "//trim(str(in32_4))

end subroutine

! ####################################################

! #############################################

subroutine writeIOint32int32int32int32int32(obj,in32_1,in32_2,in32_3,in32_4,in32_5)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32_2,in32_3,in32_4,in32_5

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32_1))//" "//trim(str(in32_2))//" "&
        //trim(str(in32_3))//" "//trim(str(in32_4))//" "//trim(str(in32_5))

end subroutine

! ####################################################

! #############################################

subroutine writeIOint32int32int32int32int32int32(obj,in32_1,in32_2,in32_3,in32_4,in32_5,in32_6)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32_2,in32_3,in32_4,in32_5,in32_6

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(in32_1))//" "//trim(str(in32_2))//" "&
        //trim(str(in32_3))//" "//trim(str(in32_4))//" "//trim(str(in32_5))&
        //" "//trim(str(in32_6))

end subroutine

! ####################################################

! #############################################
subroutine writeIOint32Vector(obj,in32)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32(:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    do i=1,size(in32)
        write(obj%fh, '(A)') trim(str(in32(i) ))
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOint32int32vector(obj,in32_1,in32)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1,in32(:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)',advance="no") trim(str(in32_1))
    do i=1,size(in32)
        write(obj%fh, '(A)',advance="no") " "//trim(str(in32(i) ))
    enddo
    write(obj%fh, '(A)',advance="yes") " "
end subroutine
! #############################################

! #############################################
subroutine writeIOint32re64Vector(obj,in32_1,re64)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32_1
    real(real64),intent(in) :: re64(:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)',advance="no") trim(str(in32_1))
    do i=1,size(re64)
        write(obj%fh, '(A)',advance="no") " "//trim(str(re64(i) ))
    enddo
    write(obj%fh, '(A)',advance="yes") " "
end subroutine
! #############################################

! #############################################
subroutine writeIOint32Array(obj,in32)
    class(IO_),intent(inout) :: obj
    integer(int32),intent(in) :: in32(:,:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    do i=1,size(in32,1)
        write(obj%fh, *) in32(i,:)
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOre64(obj,re64)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    write(obj%fh, '(A)') trim(str(re64))

end subroutine writeIOre64
! #############################################

subroutine writeIOre64re64(obj,re64_1,re64_2)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64_1,re64_2

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    if( isnan(re64_1) .or. abs(re64_1) > HUGE(real64)  )then
        write(obj%fh, '(A)') "NaN "//trim(str(re64_2))
    elseif( isnan(re64_2) .or. abs(re64_2) > HUGE(real64) )then
        write(obj%fh, '(A)') trim(str(re64_1))//" NaN"
    else
        write(obj%fh, '(A)') trim(str(re64_1))//" "//trim(str(re64_2))
    endif
end subroutine

! ####################################################

! #############################################

subroutine writeIOre64re64re64(obj,re64_1,re64_2,re64_3)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64_1,re64_2,re64_3

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(re64_1))//" "//trim(str(re64_2))//" "//trim(str(re64_3))

end subroutine

! ####################################################

! #############################################

subroutine writeIOre64re64re64re64(obj,re64_1,re64_2,re64_3,re64_4)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64_1,re64_2,re64_3,re64_4

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(re64_1))//" "//trim(str(re64_2))//" "//trim(str(re64_3))//" "//trim(str(re64_4))

end subroutine

! ####################################################

! #############################################

subroutine writeIOre64re64re64re64re64(obj,re64_1,re64_2,re64_3,re64_4,re64_5)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64_1,re64_2,re64_3,re64_4,re64_5

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(re64_1))//" "//trim(str(re64_2))//" "&
        //trim(str(re64_3))//" "//trim(str(re64_4))//" "//trim(str(re64_5))

end subroutine

! ####################################################

! #############################################

subroutine writeIOre64re64re64re64re64re64(obj,re64_1,re64_2,re64_3,re64_4,re64_5,re64_6)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64_1,re64_2,re64_3,re64_4,re64_5,re64_6

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    write(obj%fh, '(A)') trim(str(re64_1))//" "//trim(str(re64_2))//" "&
        //trim(str(re64_3))//" "//trim(str(re64_4))//" "//trim(str(re64_5))&
        //" "//trim(str(re64_6))

end subroutine

! ####################################################

! #############################################
subroutine writeIOre64Vector(obj,re64,sparse)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64(:)
    logical,optional,intent(in) :: sparse
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif

    if(present(sparse) )then
        if(sparse)then
            do i=1,size(re64,1)
                if(re64(i)==0.0d0)then
                    write(obj%fh, '(A)', advance='no') '0 '
                else
                    write(obj%fh, '(A)', advance='no') '* '
                endif
                write(obj%fh, '(A)', advance='yes') ' '
            enddo
            return
        endif
    endif

    do i=1,size(re64)
        write(obj%fh, '(A)') trim(str(re64(i) ))
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOre64Array(obj,re64,sparse)
    class(IO_),intent(inout) :: obj
    real(real64),intent(in) :: re64(:,:)
    logical,optional,intent(in) :: sparse
    integer(int32) :: i,j

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    if(present(sparse) )then
        if(sparse)then

            do i=1,size(re64,1)
                do j=1,size(re64,2)
                    if(re64(i,j)==0.0d0)then
                        write(obj%fh, '(A)', advance='no') '0 '
                    else
                        write(obj%fh, '(A)', advance='no') '* '
                    endif
                enddo
                write(obj%fh, '(A)', advance='yes') ' '
            enddo
            return
        endif
    endif
    do i=1,size(re64,1)
        write(obj%fh, *) re64(i,:)
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOcomplex64(obj,complex64)
    class(IO_),intent(inout) :: obj
    complex(kind(0d0) ),intent(in) :: complex64

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    write(obj%fh, '(A)') trim(str(complex64))

end subroutine writeIOcomplex64
! #############################################

! #############################################
subroutine writeIOcomplex64Vector(obj,complex64)
    class(IO_),intent(inout) :: obj
    complex(kind(0d0) ),intent(in) :: complex64(:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    do i=1,size(complex64,1)
        write(obj%fh, '(A)') trim(str(complex64(i) ))
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOcomplex64Array(obj,complex64)
    class(IO_),intent(inout) :: obj
    complex(kind(0d0) ),intent(in) :: complex64(:,:)
    integer(int32) :: i

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    do i=1,size(complex64,1)
        write(obj%fh, *) complex64(i,:)
    enddo
end subroutine
! #############################################

! #############################################
subroutine writeIOstring(obj,string)
    class(IO_),intent(inout) :: obj
    type(String_),intent(in) :: string

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    write(obj%fh, '(A)') str(string)

end subroutine writeIOstring
! #############################################

! #############################################
subroutine writeIOstringstring(obj,string1,string2)
    class(IO_),intent(inout) :: obj
    type(String_),intent(in) :: string1,string2

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    write(obj%fh, '(A)') str(string1)//str(string2)

end subroutine
! #############################################

! #############################################
subroutine writeIOstringstringstring(obj,string1,string2,string3)
    class(IO_),intent(inout) :: obj
    type(String_),intent(in) :: string1,string2,string3

    if(obj%state=="r")then
        call print("IOClass >> Error >> This file is readonly. ")
        call print("Nothing is written.")
        return
    endif
    write(obj%fh, '(A)') str(string1)//" "//str(string2)//" "//str(string3)

end subroutine
! #############################################

! #############################################
subroutine readIOchar(obj,char)
    class(IO_),intent(inout) :: obj
    character(200) :: char

    read(obj%fh,'(A)' ) char

end subroutine readIOchar
! #############################################

! #############################################
!function readIOstring(obj) result(char)
!    class(IO_),intent(inout) :: obj
!    type(String_) :: char
!
!    read(obj%fh,'(A)' ) char%all
!
!end function readIOstring
! #############################################

! #############################################
subroutine closeIO(obj)
    class(IO_),intent(inout) :: obj

    if(obj%active .eqv. .false.)then
        print *, "ERROR :: "//"file is already closed. filename = "//obj%filename
        return
    endif
    close(obj%fh)
    obj%fh=0
    obj%active=.false.

end subroutine closeIO
! #############################################

! #############################################
subroutine printChar(char)
    character(*),intent(in) :: char

    write(*,'(A)' ) trim(char)

end subroutine
! #############################################

! #############################################
subroutine printReal64(re64)
    real(real64),intent(in) :: re64
    character(20) :: char

    write(char, '(f20.10)') re64
    write(*,'(A)' ) trim(adjustl(char))

end subroutine
! #############################################

! #############################################
subroutine printReal32(re32)
    real(real32),intent(in) :: re32

    character(20) :: char

    write(char, '(f20.5)') re32
    write(*,'(A)' ) trim(adjustl(char))

end subroutine
! #############################################

! #############################################
subroutine printint64(in64)
    integer(int64),intent(in) :: in64

    character(20) :: char

    write(char, '(i15.7)') in64
    write(*,'(A)' ) trim(adjustl(char))

end subroutine
! #############################################
! #############################################
subroutine printint32(in32)
    integer(int32),intent(in) :: in32

    character(20) :: char

    write(char, '(i10)') in32
    write(*,'(A)' ) trim(adjustl(char))
end subroutine
! #############################################

subroutine plotRealArray(x,y,z,xr,yr,zr)
    real(real64),intent(in) :: x(:),y(:)
    real(real64),optional,intent(in) :: z(:)
    character(*),optional,intent(in) :: xr, yr, zr
    type(IO_) :: f
    integer(int32) :: i

    if(present(z) )then

        call f%open("__plotRealArray__buf_.txt")
        do i=1,size(x)
            call f%write(str(x(i))//" "//str(y(i)) //" "//str(z(i)) )
        enddo
        call f%close()

        call print("")
        call print("Press Ctrl+c to close window.")
        call f%open("__plotRealArray__buf_.gp")
        call f%write('unset key' )
        do i=1,size(x)
            if(present(xr) )then
                call f%write('set xr'//trim(xr) )
            endif
            if(present(yr) )then
                call f%write('set yr'//trim(yr) )
            endif
            if(present(zr) )then
                call f%write('set zr'//trim(zr) )
            endif
            call f%write('splot "__plotRealArray__buf_.txt" w l')
            call f%write("pause -1")
        enddo
        call f%close()

        call system("gnuplot __plotRealArray__buf_.gp")

        call system("rm __plotRealArray__buf_.txt")
        call system("rm __plotRealArray__buf_.gp")

        return
    endif

    call f%open("__plotRealArray__buf_.txt")
    do i=1,size(x)
        call f%write(str(x(i))//" "//str(y(i)) )
    enddo
    call f%close()

    call print("")
    call print("Press Ctrl+c to close window.")
    call f%open("__plotRealArray__buf_.gp")
    call f%write('unset key' )

    if(present(xr) )then
        call f%write('set xr'//trim(xr) )
    endif
    if(present(yr) )then
        call f%write('set yr'//trim(yr) )
    endif
    call f%write('plot "__plotRealArray__buf_.txt" w l')
    call f%write("pause 10")
    call f%close()

    call system("gnuplot __plotRealArray__buf_.gp")

    call system("rm __plotRealArray__buf_.txt")
    call system("rm __plotRealArray__buf_.gp")

end subroutine

subroutine spyRealArray(array)
    real(real64),intent(in) :: array(:,:)
    real(real64),allocatable :: x(:),y(:)
    integer(int32) :: i,j,non_zero

    non_zero=0
    do i=1,size(array,1)
        do j=1,size(array,2)
            if(array(i,j)/=0.0d0 )then
                non_zero=non_zero+1
            endif
        enddo
    enddo
    allocate(x(non_zero),y(non_zero) )

    non_zero=0
    do i=1,size(array,1)
        do j=1,size(array,2)
            if(array(i,j)/=0.0d0 )then
                non_zero=non_zero+1
                x(non_zero) = dble(i)
                y(non_zero) = dble(j)
            endif
        enddo
    enddo

    call plot(x=x, y=y)

end subroutine

subroutine plotIO(obj,name,option)
    class(IO_),intent(inout) ::  obj
    character(*),optional,intent(in) :: name
    character(*),optional,intent(in) :: option
    type(IO_) :: gp_script

    if(present(name) )then
        obj%filename = name
    endif
    call obj%open(obj%filename,"r")
    call gp_script%open(trim(obj%filename)//"_gp_script.gp","w")
    call gp_script%write("set xlabel '"//obj%xlabel//"'")
    call gp_script%write("set ylabel '"//obj%ylabel//"'")
    call gp_script%write("set title '"//obj%title//"'")
    if(present(option) )then
        call gp_script%write("plot '"//obj%filename//"' "//option)
    else
        call gp_script%write("plot '"//obj%filename//"' ")
    endif

    call gp_script%close()
    call system("gnuplot "//trim(obj%filename)//"_gp_script.gp -pause")
    call obj%close()
end subroutine

subroutine splotIO(obj,name,option)
    class(IO_),intent(inout) ::  obj
    character(*),optional,intent(in) :: name
    character(*),optional,intent(in) :: option
    type(IO_) :: gp_script

    if(present(name) )then
        obj%filename = name
    endif
    call obj%open(obj%filename,"r")
    call gp_script%open(trim(obj%filename)//"_gp_script.gp","w")
    call gp_script%write("set xlabel '"//obj%xlabel//"'")
    call gp_script%write("set ylabel '"//obj%ylabel//"'")
    call gp_script%write("set zlabel '"//obj%zlabel//"'")
    call gp_script%write("set title '"//obj%title//"'")
    if(present(option) )then
        call gp_script%write("splot '"//obj%filename//"' "//option)
    else
        call gp_script%write("splot '"//obj%filename//"' ")
    endif

    call gp_script%close()
    call system("gnuplot "//trim(obj%filename)//"_gp_script.gp -pause")

end subroutine

subroutine flushIO(obj)
    class(IO_),intent(in) :: obj

    flush(obj%fh)

end subroutine

function parseIOChar200(obj,filename,fileformat,key1,debug) result(ret)
    class(IO_),intent(inout) :: obj
    integer(int32),optional,intent(in) :: fileformat
    character(*),intent(in) :: key1
    character(*),optional,intent(in) :: filename
    character(:),allocatable :: line
    integer(int32)::blcount=0
    integer(int32)::rmc,id,fformat
    character(200) :: ret
    logical,optional,intent(in) :: debug

    ret = " "
    fformat = input(default=1,option=fileformat)

    if(present(filename) )then
        call obj%open(filename,"r")
        if(index(filename,".json")/=0 )then
            fformat = 1
        endif
    endif

    if(fformat==PF_JSON)then
        do
            line = obj%readline()
            if(present(debug) )then
                if(debug)then
                    print *, trim(line)
                endif
            endif
            if( adjustl(trim(line))=="{" )then
                blcount=1
                cycle
            endif
            if( adjustl(trim(line))=="}" )then
                exit
            endif

            if(blcount==1)then
                if(index(line,trim(key1) )/=0)then
                    rmc=index(line,",")
                    ! カンマがあれば除く
                    if(rmc /= 0)then
                        line(rmc:rmc)=" "
                    endif
                    id = index(line,":")
                    read(line(id+1:),*) ret
                    ret = adjustl(ret)
                    exit
                else
                    cycle
                endif
            endif

            if(obj%EOF) exit

        enddo
    endif

    if(present(filename) )then
        call obj%close()
    endif

end function

function parseIO2keysChar200(obj,filename,fileformat,key1,key2,debug) result(ret)
    class(IO_),intent(inout) :: obj
    integer(int32),optional,intent(in) :: fileformat
    character(*),intent(in) :: key1,key2
    character(*),optional,intent(in) :: filename
    character(:),allocatable :: line
    integer(int32)::blcount=0
    integer(int32)::rmc,id,fformat
    character(200) :: ret
    logical,optional,intent(in) :: debug
    ret = " "

    fformat = input(default=1,option=fileformat)
    if(present(filename) )then
        call obj%open(filename,"r")
        if(index(filename,".json")/=0 )then
            fformat = 1
        endif
    endif

    if(fformat==PF_JSON)then
        do
            line = obj%readline()
            if(present(debug) )then
                if(debug)then
                    print *, trim(line)
                endif
            endif
            if( adjustl(trim(line))=="{" )then
                blcount=1
                cycle
            endif
            if( adjustl(trim(line))=="}" )then
                exit
            endif

            if(blcount==1)then
                if(index(line,trim(key1) )/=0)then
                    ! search for the second key
                    do
                        line = obj%readline()
                        if(present(debug) )then
                            if(debug)then
                                print *, trim(line)
                            endif
                        endif

                        if( index(line,"}")/=0 )then
                            exit
                        endif

                        if(index(line,trim(key2))/=0 )then
                            rmc=index(line,",")
                            if(rmc /= 0)then
                                line(rmc:rmc)=" "
                            endif
                            id = index(line,":")
                            read(line(id+1:),*) ret
                            ret = adjustl(ret)
                            exit
                        endif

                    enddo
                else
                    cycle
                endif
            endif

            if(obj%EOF) exit

        enddo
    endif

    if(present(filename) )then
        call obj%close()
    endif

end function

end module IOClass


module time_info_M

    USE global_parameters_M

    implicit none

!---Time information
    logical :: time_solve     ! decide if we are doing time solve or not
    real(dp)   alpha          ! time solve
    real(dp)   t0             ! starting time
    real(dp)   delta_t        ! time step
    real(dp)   tmax           ! max time
    real(dp)   t_initial      ! starting time
    real(dp)   ramp_start_time
    real(dp)   ramp_end_time
    real(dp)   step_start_time
    real(dp)  step_end_time
!---Time constant for flow slowdown
    real(dp)  time_constant
    real(dp) flow_reduction_percent

end module time_info_M


text:
C
C	$Id: gcelcd.f,v 1.5 2008-07-27 00:20:56 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
        SUBROUTINE GCELCD(WKID,IOS,STATUS)
C
C  Process the control elements.
C
C  At present the working control elements include:
C     VDC INTEGER PRECISION
C     CLIP RECTANGLE
C     CLIP INDICATOR
C  All other control elements will be ignored.
C
      include 'trinst.h'
      include 'trpars.h'
      include 'trstat.h'
      include 'trcode.h'
      include 'gkscom.h'
C
      INTEGER WKID, IOS, STATUS
      INTEGER TMP, XMIN, XMAX, YMIN, YMAX, CLPDAT(4)
C
      STATUS = 0
C
      IF (OPID .EQ. CELMVI) THEN
C
C  VDC integer precision.
C
        CALL GOPDEC(MOPLEN,MWHCPR,1,IOS,STATUS)
      ELSE IF (OPID .EQ. CELMCR) THEN
C
C  Set the clip rectangle for normalization transformation 1 that
C  we are currently using.  Transform the clipping rectangle if
C  IGSGCP is non-zero.
C
        CALL GOPDEC(CLPDAT,MOPLEN,4,IOS,STATUS)
        XMIN = CLPDAT(1)
        XMIN = MIN(XMIN,CLPDAT(3))
        XMAX = CLPDAT(3)
        XMAX = MAX(XMAX,CLPDAT(1))
        YMIN = CLPDAT(2)
        YMIN = MIN(YMIN,CLPDAT(4))
        YMAX = CLPDAT(4)
        YMAX = MAX(YMAX,CLPDAT(2))
        IF (IGSGCP .EQ. 0) THEN
          CALL GSWN(1,REAL(XMIN)/32767.,REAL(XMAX)/32767.,
     +                   REAL(YMIN)/32767.,REAL(YMAX)/32767.)
          CALL GSVP(1,REAL(XMIN)/32767.,REAL(XMAX)/32767.,
     +                   REAL(YMIN)/32767.,REAL(YMAX)/32767.)
        ELSE
          XMN = REAL(XMIN)/32767.
          XMX = REAL(XMAX)/32767.
          YMN = REAL(YMIN)/32767.
          YMX = REAL(YMAX)/32767.
          XM1 = CURTM(1,1)*XMN + CURTM(1,2)*YMN + CURTM(1,3)
          YM1 = CURTM(2,1)*XMN + CURTM(2,2)*YMN + CURTM(2,3)
          XM2 = CURTM(1,1)*XMX + CURTM(1,2)*YMX + CURTM(1,3)
          YM2 = CURTM(2,1)*XMX + CURTM(2,2)*YMX + CURTM(2,3)
          XM1 = MAX(0.,XM1)
          XM2 = MIN(1.,XM2)
          YM1 = MAX(0.,YM1)
          YM2 = MIN(1.,YM2)
          CALL GSWN(1,XM1,XM2,YM1,YM2)
          CALL GSVP(1,XM1,XM2,YM1,YM2)
        ENDIF
C
      ELSE IF (OPID.EQ.CELMCI) THEN
C
C  Set the clipping indicator.
C
        CALL GOPDEC(TMP,MENCPR,1,IOS,STATUS)
        CALL GSCLIP(TMP)
      END IF
C
C  Check the instruction status.
C
      IF (STATUS .NE. 0) RETURN
C
C  If the instruction is continued, then loop to get the
C  full operand list.
C
  10  CONTINUE
      IF (CNTINU) THEN
        CALL GNPART(IOS,STATUS)
        IF (STATUS .NE. 0) RETURN
        GO TO 10
      END IF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
Quail Ridge Wilderness Conservancy lies on a peninsula at the southwest end of Lake Berryessa, in Napa County.  It is home to a wide variety of flora and fauna native to Californias Inner Coast Range.  The reserves website indicates it is home to black oak, blue oak, interior live oak, oracle oak, scrub oak, and valley oak, as well as about 15 different species of bunchgrasses that are quite rare todayincluding Californias state grass (as of January 1, 2005), nassella pulchra, or Purple Needlegrass.  As a relatively pristine wilderness reserve, it is studied by biologists and ecologists.  It is one of 36 http://nrs.ucdavis.edu/ natural reserves operated by the University of California.  Three webcams look out over different parts of the reserve and can be viewed from the link above.  They also have several informational videos, listed on their main website.

While the Ridge is not freely open to the public, they offer a free guided walk each month, generally taking place from 2 to 6pm on a Saturday.  Those interested in reserving a spot on the walk should call QRWCs Executive Director, Frank Maurer, at (530) 2194477 (cell) or (530) 7581387 at least one day in advance.  Carpools for the walks leave from the Research Farm between 1 and 2pm, depending on the season.

The 2011 guided walk schedule is as follows:
  8 January
  6 February
  19 March
  16 April
  14 May (2:30  6:30)
  18 June
  9 July
  6 August
  17 September
  5 November
  10 December

You can also arrange for private guided walks.  They ask for a minimum donation of $200 for a walk for up to 10 people, which will include complimentary wine and finger foods.