PROGRAM pyramid_test
    USE Precision
    USE vtk_datasets,   ONLY : polygonal_data
    USE vtk_attributes, ONLY : scalar, attributes
    USE vtk_cells,      ONLY : polygon
    USE vtk,            ONLY : vtk_legacy_write
    IMPLICIT NONE
    !! author: Ian Porter
    !! date: 12/31/2017
    !!
    !! This is a test of a pyramid geometry using polygonal data
    !!
    INTEGER(i4k), PARAMETER     :: n_params_to_write = 1
    TYPE (polygonal_data)       :: pyramid
    TYPE (attributes), DIMENSION(n_params_to_write) :: point_vals_to_write, cell_vals_to_write
    INTEGER(i4k)                :: i = 0_i4k
    INTEGER(i4k),     PARAMETER :: n_x = 3, n_y = 3, n_z = 3, n_faces = 20, unit = 20
    CHARACTER(LEN=*), PARAMETER :: filename = 'pyramid.vtk'
    CHARACTER(LEN=*), PARAMETER :: title    = 'Testing of pyramid geometry'
    INTEGER(i4k), DIMENSION(3)  :: dims
    REAL(r8k), DIMENSION(n_faces,1:n_params_to_write)     :: cell_vals
    REAL(r8k), DIMENSION(n_x*n_y*n_z,1:n_params_to_write) :: point_vals
    REAL(r8k), DIMENSION(3,n_x*n_y*n_z), PARAMETER :: points = RESHAPE ( &
      & [ 0.0_r8k,   0.0_r8k,   0.0_r8k,  &
      &   0.5_r8k,   0.0_r8k,   0.0_r8k,  &
      &   1.0_r8k,   0.0_r8k,   0.0_r8k,  &
      &   0.0_r8k,   0.5_r8k,   0.0_r8k,  &
      &   0.5_r8k,   0.5_r8k,   0.0_r8k,  &
      &   1.0_r8k,   0.5_r8k,   0.0_r8k,  &
      &   0.0_r8k,   1.0_r8k,   0.0_r8k,  &
      &   0.5_r8k,   1.0_r8k,   0.0_r8k,  &
      &   1.0_r8k,   1.0_r8k,   0.0_r8k,  &
      &   0.25_r8k,  0.25_r8k,  0.5_r8k,  &
      &   0.50_r8k,  0.25_r8k,  0.5_r8k,  &
      &   0.75_r8k,  0.25_r8k,  0.5_r8k,  &
      &   0.25_r8k,  0.50_r8k,  0.5_r8k,  &
      &   0.50_r8k,  0.50_r8k,  0.5_r8k,  &
      &   0.75_r8k,  0.50_r8k,  0.5_r8k,  &
      &   0.25_r8k,  0.75_r8k,  0.5_r8k,  &
      &   0.50_r8k,  0.75_r8k,  0.5_r8k,  &
      &   0.75_r8k,  0.75_r8k,  0.5_r8k,  &
      &   0.375_r8k, 0.375_r8k, 0.75_r8k, &
      &   0.50_r8k,  0.375_r8k, 0.75_r8k, &
      &   0.625_r8k, 0.375_r8k, 0.75_r8k, &
      &   0.375_r8k, 0.50_r8k,  0.75_r8k, &
      &   0.50_r8k,  0.50_r8k,  0.75_r8k, &
      &   0.625_r8k, 0.50_r8k,  0.75_r8k, &
      &   0.375_r8k, 0.625_r8k, 0.75_r8k, &
      &   0.50_r8k,  0.625_r8k, 0.75_r8k, &
      &   0.625_r8k, 0.625_r8k, 0.75_r8k ], [3,n_x*n_y*n_z] )
    REAL(r8k), DIMENSION(n_x*n_y*n_z), PARAMETER :: power = &
      & [ 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, &
      &   100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 200.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, &
      &   100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k, 100.0_r8k ]
    INTEGER(i4k), DIMENSION(n_faces), PARAMETER :: cellID = &
      & [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 ]
    TYPE(polygon), DIMENSION(n_faces) :: polygon_faces
    CHARACTER(LEN=10), DIMENSION(n_params_to_write), PARAMETER :: cell_dataname = &
      & [ 'cellIDs   ' ]
    CHARACTER(LEN=10), DIMENSION(n_params_to_write), PARAMETER :: point_dataname = &
      & [ 'Power_(W) ' ]

    CALL polygon_faces(1)%setup ( (/ 0, 1, 4, 3 /) )
    CALL polygon_faces(2)%setup ( (/ 1, 2, 5, 4 /) )
    CALL polygon_faces(3)%setup ( (/ 3, 4, 7, 6 /) )
    CALL polygon_faces(4)%setup ( (/ 4, 5, 8, 7 /) )
    CALL polygon_faces(5)%setup ( (/ 9, 10, 13, 12 /) )
    CALL polygon_faces(6)%setup ( (/ 10, 11, 14, 13 /) )
    CALL polygon_faces(7)%setup ( (/ 12, 13, 16, 15 /) )
    CALL polygon_faces(8)%setup ( (/ 13, 14, 17, 16 /) )
    CALL polygon_faces(9)%setup ( (/ 18, 19, 22, 21 /) )
    CALL polygon_faces(10)%setup ( (/ 19, 20, 23, 22 /) )
    CALL polygon_faces(11)%setup ( (/ 21, 22, 25, 24 /) )
    CALL polygon_faces(12)%setup ( (/ 22, 23, 26, 25 /) )
    CALL polygon_faces(13)%setup ( (/ 0, 1, 10, 19, 18, 9 /) )
    CALL polygon_faces(14)%setup ( (/ 1, 2, 11, 20, 19, 10 /) )
    CALL polygon_faces(15)%setup ( (/ 0, 3, 12, 21, 18, 9 /) )
    CALL polygon_faces(16)%setup ( (/ 3, 6, 15, 24, 21, 12 /) )
    CALL polygon_faces(17)%setup ( (/ 6, 7, 16, 25, 24, 15 /) )
    CALL polygon_faces(18)%setup ( (/ 7, 8, 17, 26, 25, 16 /) )
    CALL polygon_faces(19)%setup ( (/ 8, 5, 14, 23, 26, 17 /) )
    CALL polygon_faces(20)%setup ( (/ 5, 2, 11, 20, 23, 14 /) )

    cell_vals(:,1)  = REAL(cellID(:))
    point_vals(:,1) = power(:)

    dims = (/ n_x, n_y, n_z /)
    CALL pyramid%init (points=points, polygons=polygon_faces)
    DO i = 1, n_params_to_write
        ! Cell values
        IF (.NOT. ALLOCATED(cell_vals_to_write(i)%attribute))THEN
            ALLOCATE(scalar::cell_vals_to_write(i)%attribute)
            cell_vals_to_write(1)%n = SIZE(cell_vals(:,1))
        END IF
        CALL cell_vals_to_write(i)%attribute%init (TRIM(cell_dataname(i)), numcomp=1, values1d=cell_vals(:,i))
        ! Point values
        IF (.NOT. ALLOCATED(point_vals_to_write(i)%attribute))THEN
            ALLOCATE(scalar::point_vals_to_write(i)%attribute)
            point_vals_to_write(1)%n = SIZE(point_vals(:,1))
        END IF
        CALL point_vals_to_write(i)%attribute%init (TRIM(point_dataname(i)), numcomp=1, values1d=point_vals(:,i))
    END DO

    CALL vtk_legacy_write (unit, pyramid, celldatasets=cell_vals_to_write, pointdatasets=point_vals_to_write, &
      &                    filename=filename, title=title)

    WRITE(*,*) 'Finished'

END PROGRAM pyramid_test


text: #redirect Users/DanielQuincy

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module mod_output
  use mpi
  use decomp_2d_io
  use mod_param     , only: dims,dx,dy,dz
  use mod_common_mpi, only:ierr,myid,coord
  use mod_types
  implicit none
  private
  public out0d,out1d,out1d_2,out2d,out3d,write_log_output,write_visu_2d,write_visu_3d
  contains
  subroutine out0d(fname,n,var)
    !
    ! appends the first n entries of an array
    ! var to a file
    ! fname -> name of the file
    ! n     -> number of entries
    ! var   -> input array of real values
    !
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in) :: n
    real(rp), intent(in), dimension(:) :: var
    integer :: iunit
    character(len=30) :: cfmt
    integer :: i
    !
    write(cfmt,'(A,I3,A)') '(',n,'E15.7)'
    iunit = 10
    if (myid .eq. 0) then
      open(iunit,file=fname,position='append')
      write(iunit,trim(cfmt)) (var(i),i=1,n)
      close(iunit)
    endif
    return
  end subroutine out0d
  !
  subroutine out1d(fname,n,idir,z,dzlzi,p)
    !
    ! writes the profile of a variable averaged
    ! over two domain directions
    !
    ! fname -> name of the file
    ! n     -> size of the input array
    ! idir  -> direction of the profile
    ! z     -> z coordinate (grid is non-uniform in z)
    ! dzlzi -> dz/lz weight of a grid cell for averaging over z
    ! p     -> 3D input scalar field
    !
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in), dimension(3) :: n
    integer , intent(in) :: idir
    real(rp), intent(in), dimension(0:) :: z,dzlzi
    real(rp), intent(in), dimension(0:,0:,0:) :: p
    real(rp), allocatable, dimension(:) :: p1d
    integer :: i,j,k,ii,jj
    integer :: iunit
    integer, dimension(3) :: ng
    !
    ng(:) = n(:)
    ng(1:2) = n(1:2)*dims(1:2)
    iunit = 10
    select case(idir)
    case(3)
      allocate(p1d(n(3)))
      do k=1,ng(3)
        p1d(k) = 0.
        do j=1,n(2)
          do i=1,n(1)
            p1d(k) = p1d(k) + p(i,j,k)
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,p1d(1),ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      p1d(:) = p1d(:)/(1.*ng(1)*ng(2))
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do k=1,n(3)
          write(iunit,'(2E15.7)') z(k),p1d(k)
        enddo
        close(iunit)
      endif
    case(2)
      allocate(p1d(ng(2)))
      p1d(:) = 0.
      do j=1,n(2)
        jj = coord(2)*n(2)+j
        p1d(jj) = 0.
        do k=1,n(3)
          do i=1,n(1)
            p1d(jj) = p1d(jj) + p(i,j,k)*dzlzi(k)
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,p1d(1),ng(2),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      p1d(:) = p1d(:)/(1.*ng(1))
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do j=1,ng(2)
          write(iunit,'(2E15.7)') (1.*j-.5)/(1.*ng(2)),p1d(j)
        enddo
        close(iunit)
      endif
    case(1)
      allocate(p1d(ng(1)))
      p1d(:) = 0.
      do i=1,n(1)
        ii = coord(1)*n(1)+i
        p1d(i) = 0.
        do k=1,n(3)
          do j=1,n(2)
            p1d(ii) = p1d(ii) + p(i,j,k)*dzlzi(k)
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,p1d(1),ng(1),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      p1d(:) = p1d(:)/(1.*ng(2))
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do i=1,ng(1)
          write(iunit,'(2E15.7)') (1.*i-.5)/(1.*n(1)),p1d(j)
        enddo
        close(iunit)
      endif
    end select
    deallocate(p1d)
  end subroutine out1d
  !
  !
  subroutine out2d(fname,inorm,islice,p)
    !
    ! saves a planar slice of a scalar field into a binary file
    !
    ! fname  -> name of the output file
    ! inorm  -> plane is perpendicular to direction
    !           inorm (1,2,3)
    ! islice -> plane is of constant index islice
    !           in direction inorm
    ! p      -> 3D input scalar field
    !
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in) :: inorm,islice
    real(rp), intent(in), dimension(:,:,:) :: p
    !
    select case(inorm)
    case(1) !normal to x --> yz plane
       call decomp_2d_write_plane(3,p,inorm,islice,fname)
    case(2) !normal to y --> zx plane
       call decomp_2d_write_plane(3,p,inorm,islice,fname)
    case(3) !normal to z --> xy plane
       call decomp_2d_write_plane(3,p,inorm,islice,fname)
    end select
    return
  end subroutine out2d
  !
  subroutine out3d(fname,nskip,p)
    !
    ! saves a 3D scalar field into a binary file
    !
    ! fname  -> name of the output file
    ! nskip  -> array with the step size for which the
    !           field is written; i.e.: (/1,1,1/)
    !           writes the full field
    ! p      -> 3D input scalar field
    !
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in), dimension(3) :: nskip
    real(rp), intent(in), dimension(:,:,:) :: p
    integer :: fh
    integer(kind=MPI_OFFSET_KIND) :: filesize,disp
    !
    call MPI_FILE_OPEN(MPI_COMM_WORLD, fname, &
         MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL,fh, ierr)
    filesize = 0_MPI_OFFSET_KIND
    call MPI_FILE_SET_SIZE(fh,filesize,ierr)
    disp = 0_MPI_OFFSET_KIND
    call decomp_2d_write_every(3,p,nskip(1),nskip(2),nskip(3),fname,.true.)
    call MPI_FILE_CLOSE(fh,ierr)
    return
  end subroutine out3d
  !
  subroutine write_log_output(fname,fname_fld,varname,nmin,nmax,nskip,time,istep)
    !
    ! appends information about a saved binary file to a log file
    ! this file is used to generate a xdmf file for visualization of field data
    !
    ! fname     -> name of the output log file
    ! fname_fld -> name of the saved binary file (excluding the directory)
    ! varname   -> name of the variable that is saved
    ! nmin      -> first element of the field that is saved in each direction, e.g. (/1,1,1/)
    ! nmax      -> last  element of the field that is saved in each direction, e.g. (/ng(1),ng(2),ng(3)/)
    ! nskip     -> step size between nmin and nmax, e.g. (/1,1,1/) if the whole array is saved
    ! time      -> physical time
    ! istep     -> time step number
    !
    implicit none
    character(len=*), intent(in) :: fname,fname_fld,varname
    integer , intent(in), dimension(3) :: nmin,nmax,nskip
    real(rp), intent(in)               :: time
    integer , intent(in)               :: istep
    character(len=100) :: cfmt
    integer :: iunit
    !
    iunit = 10
    write(cfmt, '(A)') '(A,A,A,9i5,E15.7,i7)'
    if (myid .eq. 0) then
      open(iunit,file=fname,position='append')
      write(iunit,trim(cfmt)) trim(fname_fld),' ',trim(varname),nmin,nmax,nskip,time,istep
      close(iunit)
    endif
  end subroutine write_log_output
  !
  subroutine write_visu_3d(datadir,fname_bin,fname_log,varname,nmin,nmax,nskip,time,istep,p)
    !
    ! wraps the calls of out3d and write_log_output into the same subroutine
    !
    implicit none
    character(len=*), intent(in)          :: datadir,fname_bin,fname_log,varname
    integer , intent(in), dimension(3)    :: nmin,nmax,nskip
    real(rp), intent(in)                  :: time
    integer , intent(in)                  :: istep
    real(rp), intent(in), dimension(:,:,:) :: p
    !
    call out3d(trim(datadir)//trim(fname_bin),nskip,p)
    call write_log_output(trim(datadir)//trim(fname_log),trim(fname_bin),trim(varname),nmin,nmax,nskip,time,istep)
    return
  end subroutine write_visu_3d
  !
  subroutine write_visu_2d(datadir,fname_bin,fname_log,varname,inorm,nslice,ng,time,istep,p)
    !
    ! wraps the calls of out2d and write-log_output into the same subroutine
    !
    implicit none
    character(len=*), intent(in)          :: datadir,fname_bin,fname_log,varname
    integer , intent(in)                  :: inorm,nslice
    integer , intent(in), dimension(3)    :: ng
    real(rp), intent(in)                  :: time
    integer , intent(in)                  :: istep
    real(rp), intent(in), dimension(:,:,:) :: p
    integer , dimension(3) :: nmin_2d,nmax_2d
    !
    call out2d(trim(datadir)//trim(fname_bin),inorm,nslice,p)
    select case(inorm)
    case(1)
      nmin_2d(:) = (/nslice,1    ,1    /)
      nmax_2d(:) = (/nslice,ng(2),ng(3)/)
    case(2)
      nmin_2d(:) = (/1    ,nslice,1    /)
      nmax_2d(:) = (/ng(1),nslice,ng(3)/)
    case(3)
      nmin_2d(:) = (/1    ,1    ,nslice/)
      nmax_2d(:) = (/ng(1),ng(2),nslice/)
    end select
    call write_log_output(trim(datadir)//trim(fname_log),trim(fname_bin),trim(varname),nmin_2d,nmax_2d,(/1,1,1/),time,istep)
    return
  end subroutine write_visu_2d
  !
  subroutine out1d_2(fname,n,idir,z,u,v,w) ! e.g. for a channel with streamwise dir in x
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in), dimension(3) :: n
    integer , intent(in) :: idir
    real(rp), intent(in), dimension(0:) :: z
    real(rp), intent(in), dimension(0:,0:,0:) :: u,v,w
    real(rp), allocatable, dimension(:) :: um,vm,wm,u2,v2,w2,uw
    integer :: i,j,k
    integer :: iunit
    integer, dimension(3) :: ng
    integer :: q
    !
    ng(:) = n(:)
    ng(1:2) = n(1:2)*dims(1:2)
    iunit = 10
    select case(idir)
    case(3)
      q = n(3)
      allocate(um(0:q+1),vm(0:q+1),wm(0:q+1),u2(0:q+1),v2(0:q+1),w2(0:q+1),uw(0:q+1))
      do k=1,n(3)
        um(k) = 0.
        vm(k) = 0.
        wm(k) = 0.
        u2(k) = 0.
        v2(k) = 0.
        w2(k) = 0.
        uw(k) = 0.
        do j=1,n(2)
          do i=1,n(1)
            um(k) = um(k) + u(i,j,k)
            vm(k) = vm(k) + v(i,j,k)
            wm(k) = wm(k) + 0.50*(w(i,j,k-1) + w(i,j,k))
            u2(k) = u2(k) + u(i,j,k)**2
            v2(k) = v2(k) + v(i,j,k)**2
            w2(k) = w2(k) + 0.50*(w(i,j,k)**2+w(i,j,k-1)**2)
            uw(k) = uw(k) + 0.25*(u(i-1,j,k) + u(i,j,k))* &
                                   (w(i,j,k-1) + w(i,j,k))
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,um(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,vm(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,wm(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,u2(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,v2(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,w2(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,uw(1),n(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      um(:) = um(:)/(1.*ng(1)*ng(2))
      vm(:) = vm(:)/(1.*ng(1)*ng(2))
      wm(:) = wm(:)/(1.*ng(1)*ng(2))
      u2(:) = sqrt(u2(:)/(1.*ng(1)*ng(2)) - um(:)**2)
      v2(:) = sqrt(v2(:)/(1.*ng(1)*ng(2)) - vm(:)**2)
      w2(:) = sqrt(w2(:)/(1.*ng(1)*ng(2)) - wm(:)**2)
      uw(:) = uw(:)/(1.*ng(1)*ng(2)) - um(:)*wm(:)
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do k=1,n(3)
          write(iunit,'(8E15.7)') z(k),um(k),vm(k),wm(k), &
                                       u2(k),v2(k),w2(k), &
                                       uw(k)
        enddo
        close(iunit)
      endif
      deallocate(um,vm,wm,u2,v2,w2,uw)
    case(2)
    case(1)
    end select
  end subroutine out1d_2
  !
  subroutine out2d_2(fname,n,idir,z,u,v,w) ! e.g. for a duct
    implicit none
    character(len=*), intent(in) :: fname
    integer , intent(in), dimension(3) :: n
    integer , intent(in) :: idir
    real(rp), intent(in), dimension(0:) :: z
    real(rp), intent(in), dimension(0:,0:,0:) :: u,v,w
    real(rp), allocatable, dimension(:,:) :: um,vm,wm,u2,v2,w2,uv,uw,vw
    integer :: i,j,k,ii,jj,kk
    integer :: iunit
    integer, dimension(3) :: ng
    integer :: p,q
    real(rp) :: x,y
    !
    ng(:) = n(:)
    ng(1:2) = n(1:2)*dims(1:2)
    iunit = 10
    select case(idir)
    case(3)
    case(2)
      p = ng(1)
      q = ng(3)
      allocate(um(p,q),vm(p,q),wm(p,q),u2(p,q),v2(p,q),w2(p,q),uv(p,q),vw(p,q))
      !
      um(:,:) = 0.
      vm(:,:) = 0.
      wm(:,:) = 0.
      u2(:,:) = 0.
      v2(:,:) = 0.
      w2(:,:) = 0.
      uv(:,:) = 0.
      vw(:,:) = 0.
      do k=1,n(3)
        kk = k
        do i=1,n(1)
          ii = i+coord(1)*n(1)
          um(ii,kk) = 0.
          vm(ii,kk) = 0.
          wm(ii,kk) = 0.
          u2(ii,kk) = 0.
          v2(ii,kk) = 0.
          w2(ii,kk) = 0.
          vw(ii,kk) = 0.
          uv(ii,kk) = 0.
          do j=1,n(2)
            jj = j+coord(2)*n(2)
            um(ii,kk) = um(ii,kk) + 0.5*(u(i-1,j,k)+u(i,j,k))
            vm(ii,kk) = vm(ii,kk) + v(i,j,k)
            wm(ii,kk) = wm(ii,kk) + 0.5*(w(i,j,k-1)+w(i,j,k))
            u2(ii,kk) = u2(ii,kk) + 0.5*(u(i-1,j,k)**2+u(i,j,k)**2)
            v2(ii,kk) = v2(ii,kk) + v(i,j,k)**2
            w2(ii,kk) = w2(ii,kk) + 0.5*(w(i,j,k-1)**2+w(i,j,k)**2)
            vw(ii,kk) = vw(ii,kk) + 0.25*(v(i,j-1,k) + v(i,j,k))* &
                                           (w(i,j,k-1) + w(i,j,k))
            uv(ii,kk) = uv(ii,kk) + 0.25*(u(i-1,j,k) + u(i,j,k))* &
                                           (v(i,j-1,k) + v(i,j,k))
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,um(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,vm(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,wm(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,u2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,v2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,w2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,vw(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,uv(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      um(:,:) =      um(:,:)/(1.*ng(2))
      vm(:,:) =      vm(:,:)/(1.*ng(2))
      wm(:,:) =      wm(:,:)/(1.*ng(2))
      u2(:,:) = sqrt(u2(:,:)/(1.*ng(2)) - um(:,:)**2)
      v2(:,:) = sqrt(v2(:,:)/(1.*ng(2)) - vm(:,:)**2)
      w2(:,:) = sqrt(w2(:,:)/(1.*ng(2)) - wm(:,:)**2)
      vw(:,:) =      vw(:,:)/(1.*ng(2)) - vm(:,:)*wm(:,:)
      uv(:,:) =      uv(:,:)/(1.*ng(2)) - um(:,:)*vm(:,:)
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do k=1,ng(3)
          do i=1,ng(1)
            y = (i-.5)*dx
            write(iunit,'(10E15.7)') y,z(k),um(i,k),vm(i,k),wm(i,k), &
                                            u2(i,k),v2(i,k),w2(i,k), &
                                            vw(i,k),uv(i,k)
          enddo
        enddo
        close(iunit)
      endif
      deallocate(um,vm,wm,u2,v2,w2,vw,uv)
    case(1)
      p = ng(2)
      q = ng(3)
      allocate(um(p,q),vm(p,q),wm(p,q),u2(p,q),v2(p,q),w2(p,q),uv(p,q),uw(p,q))
      !
      um(:,:) = 0.
      vm(:,:) = 0.
      wm(:,:) = 0.
      u2(:,:) = 0.
      v2(:,:) = 0.
      w2(:,:) = 0.
      uv(:,:) = 0.
      uw(:,:) = 0.
      do k=1,n(3)
        kk =  k
        do j=1,n(2)
          jj = j+coord(2)*n(2)
          um(jj,kk) = 0.
          vm(jj,kk) = 0.
          wm(jj,kk) = 0.
          u2(jj,kk) = 0.
          v2(jj,kk) = 0.
          w2(jj,kk) = 0.
          uv(jj,kk) = 0.
          uw(jj,kk) = 0.
          do i=1,n(1)
            ii = i+coord(1)*n(1)
            um(jj,kk) = um(jj,kk) + u(i,j,k)
            vm(jj,kk) = vm(jj,kk) + 0.5*(v(i,j-1,k)+v(i,j,k))
            wm(jj,kk) = wm(jj,kk) + 0.5*(w(i,j,k-1)+w(i,j,k))
            u2(jj,kk) = u2(jj,kk) + u(i,j,k)**2
            v2(jj,kk) = v2(jj,kk) + 0.5*(v(i,j-1,k)**2+v(i,j,k)**2)
            w2(jj,kk) = w2(jj,kk) + 0.5*(w(i,j,k-1)**2+w(i,j,k)**2)
            uv(jj,kk) = uv(jj,kk) + 0.25*(u(i-1,j,k) + u(i,j,k))* &
                                         (v(i,j-1,k) + v(i,j,k))
            uw(jj,kk) = uw(jj,kk) + 0.25*(u(i-1,j,k) + u(i,j,k))* &
                                         (w(i,j,k-1) + w(i,j,k))
          enddo
        enddo
      enddo
      call mpi_allreduce(MPI_IN_PLACE,um(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,vm(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,wm(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,u2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,v2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,w2(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,uv(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      call mpi_allreduce(MPI_IN_PLACE,uw(1,1),ng(1)*ng(3),MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
      um(:,:) =      um(:,:)/(1.*ng(2))
      vm(:,:) =      vm(:,:)/(1.*ng(2))
      wm(:,:) =      wm(:,:)/(1.*ng(2))
      u2(:,:) = sqrt(u2(:,:)/(1.*ng(2)) - um(:,:)**2)
      v2(:,:) = sqrt(v2(:,:)/(1.*ng(2)) - vm(:,:)**2)
      w2(:,:) = sqrt(w2(:,:)/(1.*ng(2)) - wm(:,:)**2)
      uv(:,:) =      uv(:,:)/(1.*ng(2)) - um(:,:)*vm(:,:)
      uw(:,:) =      uw(:,:)/(1.*ng(2)) - um(:,:)*wm(:,:)
      if(myid.eq.0) then
        open(unit=iunit,file=fname)
        do kk=1,ng(3)
          do jj=1,ng(2)
            y = (jj-.5)*dy
            write(iunit,'(10E15.7)') y,z(kk),um(jj,kk),vm(jj,kk),wm(jj,kk), &
                                             u2(jj,kk),v2(jj,kk),w2(jj,kk), &
                                             uv(jj,kk),uw(jj,kk)
          enddo
        enddo
        close(iunit)
      endif
      deallocate(um,vm,wm,u2,v2,w2,vw,uv)
    end select
    return
  end subroutine out2d_2
end module mod_output

max_stars_repo_name: maxcuda/CaNS

max_issues_repo_name: M-Abbaszadeh/CaNS

max_forks_repo_name: M-Abbaszadeh/CaNS


text:
! ------------------------------------------------------------ !
! Finding selected elements of the inverse of sparce matrix.
! Implementation based on the Ref. "Knitting algorithm for calculating
! Green functions in quantum systems".
! * K. Kazymyrenko and X. Waintal, Phys. Rev. B 77, 115119
! How to use:
! 1. Call selinv%init_memory(no_nodes,no_vals,vals,rc)
! 2. and find inverses for selected elements: call selinv%invert(elements,evals)
! Bugs:
! - does not work for some matrices
! - for now cannot be used to find first element of the matrix
! ------------------------------------------------------------ !
module modskminv
use modcommons
implicit none
private
integer :: SKMINV_DEBUG_LEVEL = 0 ! show 1, hide 0 debug information

! -----------------------------------------------------------
! Structure which holds the value of the n-th element in the
! matrix and the coupling values to neightbours. This type
! is generated from the rows(),cols(),values() information of
! the matrix.
! -----------------------------------------------------------
type node
    complex*16,allocatable :: Tij(:) ! hoping from current i-th element to neightboors T(i->j)
    complex*16,allocatable :: Tji(:) ! inverse of above T(j->i)
    integer   ,allocatable :: Iij(:) ! global indices of neightbours
    integer   ,allocatable :: Iij_less(:) ! global indices of neightbours with id smaller than current node
    complex*16 :: H0                 ! on site value
    integer    :: no_bonds           ! number of hopings, size of the Tij,Tji, and Iij vectors
    integer    :: no_used            ! number of connections in dynamic interface (used in algorithm)
    integer    :: fl_id              ! dynamically contains mapping from global id to local id in the interface (here called frontline)
    integer    :: no_bonds_less      ! number of neightbours with id smaller than current node id. Size of Iij_less
    contains
    procedure,pass(this) :: reset_node ! reset variables and deallocate arrays
    procedure,pass(this) :: alloc_mem  ! allocate all arrays
end type node

! -----------------------------------------------------------
! Describes the moving interface of not fully connected nodes
! Here we use one structure to define the moving frontline and
! Greens function of selected elements.
! -----------------------------------------------------------
type frontline
    complex*16,allocatable :: Gij(:,:) ! Greens function of the moving interface/frontline
    integer   ,allocatable :: Iij(:)   ! global id-s of nodes in the interface/frontline
    integer    :: no_sites             ! number of sites in the frontline
!    complex*16 :: Gp,Gxa,Gax           ! Gp - holds greens function at specific point (i,j), Gxa,Gax - are temporal variables
!    integer    :: i,j                  ! defines the element to be inversed
    contains
    procedure,pass(this) :: free_frontline!free memory
    procedure,pass(this) :: init_frontline!allocate memory
!    procedure,pass(this) :: init_point!allocate memory for element (i,j)
    procedure,pass(this) :: copy_frontline!(this,scr_ft)
end type frontline

type element
    complex*16,allocatable :: Gij(:,:) ! Greens function of the moving interface/frontline
    integer    :: no_sites             ! number of sites in the frontline
    complex*16 :: Gp,Gxa,Gax           ! Gp - holds greens function at specific point (i,j), Gxa,Gax - are temporal variables
    integer    :: i,j                  ! defines the element to be inversed
    contains
    procedure,pass(this) :: free_element!free memory
    procedure,pass(this) :: init_element!allocate memory for element (i,j)
    procedure,pass(this) :: copy_element!(this,scr_ft)
end type element

! -----------------------------------------------------------
! Main type responsible for the selected inversion of the
! sparse matrix. Contains structure of the matrix represented
! by nodes and seleted elements to be inversed.
! -----------------------------------------------------------
type skminv
    type(node) , allocatable    :: nodes(:)            ! node representation of the inversed matrix
!    type(frontline)             :: flNew,flOld         ! dynamic frontline
    type(element),allocatable :: elOld(:),elNew(:)   ! dynamic Greens functions of selected elements
    integer                     :: no_nodes,no_elements! size of the matrix, and number of points to be inversed
    contains
    procedure,pass(this) :: init_memory
    procedure,pass(this) :: invert
    procedure,pass(this) :: free_memory
end type skminv

! -----------------------------------------------------------
! Make public some of the variables
! -----------------------------------------------------------
public :: SKMINV_DEBUG_LEVEL
public :: skminv
contains

! -----------------------------------------------------------
!                           SKMINV
! Free all allocated memory: matrix and selected points.
! -----------------------------------------------------------
subroutine free_memory(this)
    class(skminv) :: this
    integer :: i
    ! deallocate matrix
    do i = 1 , this%no_nodes
        call this%nodes(i)%reset_node()
    enddo
    if(allocated(this%nodes))deallocate(this%nodes)
    ! just in case deallocate selected points, they
    ! should be deallocated in invert function
    do i = 1 , this%no_elements
        call this%elNew(i)%free_element()
        call this%elOld(i)%free_element()
    enddo
    if(allocated(this%elNew))deallocate(this%elNew)
    if(allocated(this%elOld))deallocate(this%elOld)

    ! free moving frontlines
!    call this%flNew%free_frontline()
!    call this%flOld%free_frontline()
    this%no_nodes    = 0
    this%no_elements = 0
end subroutine free_memory

subroutine realloc_fl(new_size,gij,Iij,old_size)
    integer     :: new_size,old_size
    complex(16),allocatable :: gij(:,:)
    integer,allocatable     :: Iij(:)

    if(new_size /= old_size) then
        deallocate(gij,Iij)
        allocate(Gij(new_size,new_size))
        allocate(Iij(new_size))
        old_size = new_size
    endif
end subroutine realloc_fl
! -----------------------------------------------------------
!                           SKMINV
! Find the selected elements of the inverse of the input matrix.
! elements(no_elements,2) - array which contains pairs (row,col)
!               to elements which are suppose to be inverted
! values(no_elements) - contains values of selected elements.
!
! TODO:
! call init_memory before this function to allocate memory.
! -----------------------------------------------------------
subroutine invert(this,elements,values)
    class(skminv)           :: this
    integer,dimension(:,:)  :: elements
    complex*16,dimension(:) :: values
    ! internal variables:
    integer    :: no_elems,i
    complex*16 :: GAA,GXA,GAX,tmpC
    integer    :: c_node,c_elem,no_fl,k,i1,i2,s1,s2,p,q
    type(element) :: tmp_element
    ! profiling variables
    doubleprecision :: progress,progress_it
    doubleprecision :: t_step0,t_step1,t_step1e,t_step2,t_step2e,t_tmp,t_step3,t_step4,t_total
    ! elements
    complex(16),allocatable :: oGij1(:,:),nGij1(:,:)
    complex(16),allocatable :: oGij2(:,:),nGij2(:,:)
    complex(16),allocatable :: e_Gax(:),e_Gxa(:),e_Gp(:)

    complex(16),allocatable ::fl_new_Gij(:,:),fl_old_Gij(:,:)
    integer,allocatable     ::fl_new_Iij(:),fl_old_Iij(:),fl_id(:),fl_id_less(:)
    integer :: fl_new_no_sites,fl_old_no_sites
    integer :: no_bonds_less
    complex(16),allocatable  :: Tij(:) , Tji(:)
    if(SKMINV_DEBUG_LEVEL > 0) then
        print*,"skminv::starting calculation for:"
        print*,"        Num. elements:",size(values,1)
        print*,"        Matrix. size :",this%no_nodes
    endif
    t_step0  = 0
    t_step1  = 0
    t_step1e = 0
    t_step2  = 0
    t_step2e = 0
    t_step3  = 0
    t_step4  = 0
    t_total  = 0

    t_tmp   = get_clock()
    t_total = t_tmp
    ! allocate elements
    no_elems         = size(values,1)
    this%no_elements = no_elems
    allocate(this%elNew(no_elems))
    allocate(this%elOld(no_elems))

    allocate(oGij1(1,no_elems))
    allocate(oGij2(1,no_elems))
    allocate(nGij1(1,no_elems))
    allocate(nGij2(1,no_elems))
    allocate(e_Gax(no_elems))
    allocate(e_Gxa(no_elems))
    allocate(e_Gp(no_elems ))
    oGij1  = 0
    oGij2  = 0
    nGij1  = 0
    nGij2  = 0
    e_Gp  = 0
    e_Gxa = 0
    e_Gax = 0

    ! initialize elements to invert
    do i = 1,no_elems
        call this%elOld(i)%init_element(1)
        this%elOld(i)%i      = elements(i,1)
        this%elOld(i)%j      = elements(i,2)
        this%elOld(i)%Gp     = 0
!        this%elOld(i)%Iij(1) = 1 ! set starting element
        call this%elNew(i)%init_element(1)
        this%elNew(i)%i  = this%elOld(i)%i
        this%elNew(i)%j  = this%elOld(i)%j
        this%elNew(i)%Gp = 0
!        this%elNew(i)%Iij(1) = 1 ! set starting element
    enddo

    ! initialize first element to be inverted
!    call this%flNew%free_frontline()
!    call this%flOld%free_frontline()
!    call this%flNew%init_frontline(1)
!    call this%flOld%init_frontline(1)
!
!    this%flNew%Gij(1,1)  = 1/this%nodes(1)%H0
!    this%flNew%Iij(1)    = 1
    this%nodes(1)%fl_id  = 1
!    call this%flOld%copy_frontline(this%flNew)

    allocate(fl_new_Gij(1,1))
    allocate(fl_old_Gij(1,1))
    allocate(fl_new_Iij(1))
    allocate(fl_old_Iij(1))
    allocate(fl_id(1))
    allocate(fl_id_less(1))

    fl_new_Gij = 1/this%nodes(1)%H0
    fl_new_Iij = 1
    fl_old_Gij = fl_new_Gij
    fl_old_Iij = fl_new_Iij

    fl_new_no_sites = 1
    fl_old_no_sites = 1

    fl_id = 1
    allocate(Tij(1),Tji(1))

    t_step0 = get_clock() - t_tmp
    ! Perform recursive iteration through all the matrix
    progress_it = 0.1
    do c_node = 2 , this%no_nodes

        no_bonds_less = this%nodes(c_node)%no_bonds_less
        deallocate(Tij,Tji,fl_id_less)
        allocate(Tij(no_bonds_less),Tji(no_bonds_less),fl_id_less(no_bonds_less))
        do s1 = 1 , this%nodes(c_node)%no_bonds_less
            Tij(s1) = this%nodes(c_node)%Tij(s1)
            Tji(s1) = this%nodes(c_node)%Tji(s1)
            fl_id_less(s1) = this%nodes(this%nodes(c_node)%Iij_less(s1))%fl_id
        enddo

        progress = c_node/(this%no_nodes+0.0)
        if(progress > progress_it) then
            progress_it = progress_it + 0.05
            if(SKMINV_DEBUG_LEVEL == 2) then
                print"(A,I3,A)"," skminv::solving stage:",nint(progress*100),"%"
            endif
        endif
        t_tmp = get_clock()
        ! Calculate the Gaa from Eq. (0)
        tmpC = 0
        do s1 = 1 , no_bonds_less
            i1 = fl_id_less(s1)
            do s2 = 1 , no_bonds_less
                i2   = fl_id_less(s2)
!                tmpC = tmpC + this%nodes(c_node)%Tij(s1)*this%nodes(c_node)%Tji(s2)*this%flOld%Gij(i2,i1)
                tmpC = tmpC + Tij(s1)*Tji(s2)*fl_old_Gij(i2,i1)
            enddo ! end of s2
        enddo ! end of s1

        ! Updade first step of algorithm
        GAA = 1.0/( this%nodes(c_node)%H0 - tmpC )

        ! Actualize the number of coupling on the frontline after addition
        ! of new site.
        do s1 = 1 , no_bonds_less
            k = this%nodes(c_node)%Iij_less(s1)
            this%nodes(k)%no_used      = this%nodes(k)%no_used      + 1
            this%nodes(c_node)%no_used = this%nodes(c_node)%no_used + 1
        enddo

        if(SKMINV_DEBUG_LEVEL == 1) then
            ! For the last step this should be equal to last element
            print"(A,i,A,2f12.6,A)","Current node:",c_node," Gaa value=(",GAA,")"
        endif
        t_step1  = t_step1 + get_clock() - t_tmp
        t_tmp    = get_clock()
        ! -----------------------------------------------
        ! Updating the Gaa values of the selected elements
        do c_elem = 1 , no_elems
            GXA   = 0
            GAX   = 0
            ! special treatment is done if current element is equal c_node
            ! if yes Eq. (9) cannot be use and we use (8) instead
            if( elements(c_elem,2) == c_node  )then
                GXA = GAA
            else ! normal treatement
                do s1 = 1 , no_bonds_less
                    i2  = fl_id_less(s1)
!                    GXA = GXA - this%nodes(c_node)%Tij(s1)*this%elOld(c_elem)%Gij(i2,1)*GAA
                    GXA = GXA - Tij(s1)*oGij1(i2,c_elem)*GAA
                enddo
            endif
            ! same for i-th component
            if( elements(c_elem,1) == c_node  )then
                GAX = GAA
            else
                do s1 = 1 , no_bonds_less
                    i2 = fl_id_less(s1)
!                    GAX = GAX - this%nodes(c_node)%Tji(s1)*this%elOld(c_elem)%Gij(i2,2)*GAA
                    GAX = GAX - Tji(s1)*oGij2(i2,c_elem)*GAA
                enddo
            endif
            ! Update element value and temporal Gax,Gxa for later use.
!            this%elOld(c_elem)%Gax = GAX
!            this%elOld(c_elem)%Gxa = GXA
!            this%elOld(c_elem)%Gp  = this%elOld(c_elem)%Gp + GXA*GAX/GAA

            e_Gax(c_elem) = GAX
            e_Gxa(c_elem) = GXA
            e_Gp (c_elem) = e_Gp (c_elem) + GXA*GAX/GAA
        enddo ! end of loop over elements

        ! -----------------------------------------------
        ! Updating the number of sites on the frontline
        k = 0
!        do p = 1 , this%flOld%no_sites
        do p = 1 , fl_old_no_sites
            ! count only those sites which number of used connection is smaller than number of bounds
!            if(this%nodes(this%flOld%Iij(p))%no_used < this%nodes(this%flOld%Iij(p))%no_bonds )then
            if(this%nodes(fl_old_Iij(p))%no_used < this%nodes(fl_old_Iij(p))%no_bonds )then
                k = k + 1
            endif
        enddo
        if(this%nodes(c_node)%no_used < this%nodes(c_node)%no_bonds) k = k + 1
        no_fl = k ! new number of sites in the fronline

        if(no_fl == 0) exit ! this will exit for the last element where frontline is closed
!        call this%flNew%init_frontline(no_fl) ! allocate memory for new fronline
        call realloc_fl(no_fl,fl_new_Gij,fl_new_Iij,fl_new_no_sites)

        k = 0
!        do p = 1, this%flOld%no_sites
        do p = 1, fl_old_no_sites
            if(this%nodes(fl_old_Iij(p))%no_used < this%nodes(fl_old_Iij(p))%no_bonds )then
                k = k + 1
                fl_new_Iij(k) = fl_old_Iij(p)
            endif
        enddo
        ! Add new site to the frontline
        if(this%nodes(c_node)%no_used < this%nodes(c_node)%no_bonds) then
            k = k + 1
            fl_new_Iij(k) = c_node
!            this%flNew%Iij(k)         = c_node
        endif

        if(no_fl /= size(fl_id))then
            deallocate(fl_id)
            allocate(fl_id(no_fl))
        endif
        fl_id = 0
        do p = 1 , fl_new_no_sites
            fl_id(p) = this%nodes(fl_new_Iij(p))%fl_id
        enddo

        if(SKMINV_DEBUG_LEVEL == 1) then
            ! For the last step this should be equal to last element
            print"(A,i,A,10000i)","New fron line size:",no_fl," fl. ids:",fl_new_Iij
        endif

        t_step1e = t_step1e + get_clock() - t_tmp
        t_tmp    = get_clock()
        ! Update of the rest elements of the frontline Greens function:
        ! See Eq. (9) and (10), note that we have minus sign. Which
        ! was working for us for any matrix.
        do k = 1 , no_fl-1
!            p   = this%flNew%Iij(k) ! atom do ktorego nalezy nowy gj
!            p   = fl_new_Iij(k) ! atom do ktorego nalezy nowy gj
            GXA = 0
            GAX = 0
!            i1  = this%nodes(p)%fl_id
            i1  = fl_id(k)
!            print*,k, i1 , size(fl_id)
            do s1 = 1 , no_bonds_less
!                q  = this%nodes(c_node)%Iij_less(s1)
                i2 = fl_id_less(s1)
!                i2 = this%nodes(q)%fl_id
                GXA = GXA + Tij(s1)*fl_old_Gij(i1,i2)
                GAX = GAX + Tji(s1)*fl_old_Gij(i2,i1)
            enddo
            fl_new_Gij(k,no_fl) = -GXA*GAA
            fl_new_Gij(no_fl,k) = -GAX*GAA
        enddo
        ! Update rest of the matrix: See Eq. (11)
        do k = 1 , no_fl-1
!            s1 = fl_new_Iij(k)
!            i1 = this%nodes(s1)%fl_id
            i1 = fl_id(k)
        do p = 1 , no_fl-1
!            s2 = fl_new_Iij(p)
!            i2 = this%nodes(s2)%fl_id
            i2 = fl_id(p)
            fl_new_Gij(k,p) = fl_old_Gij(i1,i2)  &
                                + fl_new_Gij(k,no_fl)*fl_new_Gij(no_fl,p)/GAA
        enddo
        enddo
        ! Fill last element of the matrix
        fl_new_Gij(no_fl,no_fl) = GAA

        t_step2 = t_step2 + get_clock() - t_tmp
        t_tmp   = get_clock()

        ! Same calculation as above but for the selected elements

        if(no_fl /= size(nGij1,2))then
            deallocate(nGij1)
            deallocate(nGij2)
            allocate(nGij1(no_fl,c_elem))
            allocate(nGij2(no_fl,c_elem))
        endif

        do c_elem = 1 , no_elems
            nGij1(no_fl,c_elem) = e_Gxa(c_elem)
            nGij2(no_fl,c_elem) = e_Gax(c_elem)
            GAX = e_Gxa(c_elem)/GAA
            GXA = e_Gax(c_elem)/GAA
            do k = 1 , no_fl-1
!                s1 = fl_new_Iij(k)
                i1 = fl_id(k)
                nGij1(k,c_elem) = oGij1(i1,c_elem) + GAX*fl_new_Gij(no_fl,k)
                nGij2(k,c_elem) = oGij2(i1,c_elem) + fl_new_Gij(k,no_fl)*GXA
            enddo
        enddo

        if(no_fl /= size(oGij1,2))then
            deallocate(oGij1)
            deallocate(oGij2)
            allocate(oGij1(no_fl,c_elem))
            allocate(oGij2(no_fl,c_elem))
        endif
        oGij1 = nGij1
        oGij2 = nGij2

        t_step2e = t_step2e + get_clock() - t_tmp
        t_tmp   = get_clock()
        ! recalculate mapping from global to local in the new frontline

        k = 0
!        do p = 1, this%flOld%no_sites
        do p = 1, fl_old_no_sites
            if(this%nodes(fl_old_Iij(p))%no_used < this%nodes(fl_old_Iij(p))%no_bonds )then
                k = k + 1
                this%nodes(fl_old_Iij(p))%fl_id = k
!                fl_id(p) = k
            endif
        enddo
        if(this%nodes(c_node)%no_used < this%nodes(c_node)%no_bonds) then
            k = k + 1
            this%nodes(c_node)%fl_id = k
!            fl_id(no_fl) = k
        endif

        ! copy new to old
!        call this%flOld%copy_frontline(this%flNew)
        call realloc_fl(no_fl,fl_old_Gij,fl_old_Iij,fl_old_no_sites)
        fl_old_Gij = fl_new_Gij
        fl_old_Iij = fl_new_Iij

        t_step3 = t_step3 + get_clock() - t_tmp
    enddo ! end of loop c_node

    ! Copy calulated values to output array
    do c_elem = 1 , no_elems
        values(c_elem) = e_Gp(c_elem)
    enddo

    ! Freeing memory
    do i = 1 , this%no_elements
        call this%elNew(i)%free_element()
        call this%elOld(i)%free_element()
    enddo
    if(allocated(this%elNew))deallocate(this%elNew)
    if(allocated(this%elOld))deallocate(this%elOld)
    this%no_elements = 0
    t_total = get_clock() - t_total
    if(SKMINV_DEBUG_LEVEL == 2) then
        print*,"skminv::Profilling info calculation times for selected stages:"
        print*,"    Initialization time  :",t_step0,"[s]"
        print*,"    Step 1. for interface:",t_step1,"[s]"
        print*,"    Step 1. for elements :",t_step1e,"[s]"
        print*,"    Step 2. for interface:",t_step2,"[s]"
        print*,"    Step 2. for elements :",t_step2e,"[s]"
        print*,"    Renumbering nodes    :",t_step3,"[s]"
        print*,"    Total time           :",t_total,"[s]"
    endif
    if(SKMINV_DEBUG_LEVEL > 0) then
        print*,"skminv::finished"
    endif
!    if(allocated(Gij)) deallocate(Gij)
end subroutine invert

! ------------------------------------------------------------------------
!                           SKMINV
! Initialize memory, and convert matrix from row-col representation
! to arrays of nodes.
! no_nodes - number of nodes, i.e. size of the matrix.
! no_vals  - lenght of the array vals(:) and rc(:,1:2)
! vals(:)  - array which contains non-zero elements of the matrix
! rc(:,row/col)- contains row=1 and col=2 id to non zero values of matrix.
! ------------------------------------------------------------------------
subroutine init_memory(this,no_nodes,no_vals,vals,rc)
    class(skminv) :: this
    integer       :: no_nodes,no_vals
    complex*16    :: vals(:)
    integer       :: rc(:,:)
    ! local variables
    integer    :: i,r,c,k

    ! Free memory before initializating another one
    call this%free_memory()
    this%no_nodes = no_nodes
    allocate(this%nodes(no_nodes))

    ! clear memory of nodes
    do i = 1 , no_nodes
        call this%nodes(i)%reset_node()
    enddo
    ! caclulate number of hopings for each row/node
    do i = 1 , no_vals
        r = rc(i,1)
        c = rc(i,2)
        if(r /= c) this%nodes(r)%no_bonds = this%nodes(r)%no_bonds + 1
    enddo
    ! reallocate memory to calculated number of bonds in previous step
    ! alloc_mem uses this%nodes(i)%no_bonds to allocate memory
    do i = 1 , no_nodes
        call this%nodes(i)%alloc_mem()
    enddo

    ! Fill matrices using provided arrays
    do i = 1 , no_vals
        r = rc(i,1)
        c = rc(i,2)
        if( r == c ) then
            this%nodes(r)%H0 = vals(i)
        else
            this%nodes(r)%no_used = this%nodes(r)%no_used + 1
            k = this%nodes(r)%no_used
            this%nodes(r)%Tij(k) = vals(i)
            this%nodes(r)%Iij(k) = c
        endif
    enddo
    ! Fill the Tji vectors, this increases the number of memory
    ! usage but in the main loop one does not have to search for
    ! the Tij coupling.
    do i = 1 , no_nodes
        this%nodes(i)%no_used = 0
        do k = 1 , this%nodes(i)%no_bonds
            r = this%nodes(i)%Iij(k)
            do c = 1 , this%nodes(r)%no_bonds
                if( this%nodes(r)%Iij(c) == i ) exit
            enddo
            this%nodes(i)%Tji(k) = this%nodes(r)%Tij(c)
        enddo
    enddo
    ! Calculate number of neightboors with id smaller than current node id.
    do i = 1 , no_nodes
        this%nodes(i)%no_bonds_less = 0
        do k = 1 , this%nodes(i)%no_bonds
            if(this%nodes(i)%Iij(k) < i) this%nodes(i)%no_bonds_less = this%nodes(i)%no_bonds_less + 1
        enddo
        deallocate(this%nodes(i)%Iij_less)
        if(this%nodes(i)%no_bonds_less > 0) allocate(this%nodes(i)%Iij_less(this%nodes(i)%no_bonds_less))
    enddo
    ! Fill Iij_less array.
    do i = 1 , no_nodes
        r = 0
        if(this%nodes(i)%no_bonds_less /= 0) then
        do k = 1 , this%nodes(i)%no_bonds
            if(this%nodes(i)%Iij(k) < i) then
                r = r + 1
                this%nodes(i)%Iij_less(r) = this%nodes(i)%Iij(k)
            endif
        enddo
        endif
    enddo
    if(SKMINV_DEBUG_LEVEL == 1) then
        print*,"skminv::printing couplings values:"
        do i = 1 , no_nodes
            print"(A,i,A,1000i)","row=",i," Iij:", this%nodes(i)%Iij
            do k = 1 , this%nodes(i)%no_bonds
                print"(A,i,A,i,A,2f12.5)","from ",i," to:",this%nodes(i)%Iij(k)," t(i,j)=",(this%nodes(i)%Tij(k))
            enddo
            do k = 1 , this%nodes(i)%no_bonds
                print"(A,i,A,i,A,2f12.5)","from ",this%nodes(i)%Iij(k)," to:",i," t(j,i)=",(this%nodes(i)%Tji(k))
            enddo
        enddo
    endif
!    call this%flNew%init_frontline(1)
!    call this%flOld%init_frontline(1)

end subroutine init_memory

! -------------------------------------------------------------- !
!
! -------------------------------------------------------------- !
subroutine reset_node(this)
    class(node) :: this
    this%no_bonds      = 0
    this%no_used       = 0
    this%no_bonds_less = 0
    this%fl_id         = 0
    this%H0            = 0
    if(allocated(this%Iij))      deallocate(this%Iij)
    if(allocated(this%Tij))      deallocate(this%Tij)
    if(allocated(this%Tji))      deallocate(this%Tji)
    if(allocated(this%Iij_less)) deallocate(this%Iij_less)
endsubroutine reset_node

! -------------------------------------------------------------- !
!
! -------------------------------------------------------------- !
subroutine alloc_mem(this)
    class(node) :: this
    allocate(this%Iij(this%no_bonds))
    allocate(this%Iij_less(this%no_bonds_less))
    allocate(this%Tij(this%no_bonds))
    allocate(this%Tji(this%no_bonds))
    this%Iij      = 0
    if(this%no_bonds_less > 0) this%Iij_less = 0

    this%Tij = 0
    this%Tji = 0
endsubroutine alloc_mem

! -------------------------------------------------------------- !
!
! -------------------------------------------------------------- !
subroutine init_frontline(this,ns)
    class(frontline) :: this
    integer :: ns
    if(ns /= this%no_sites) then
        this%no_sites = ns
        if(allocated(this%Iij)) deallocate(this%Iij)
        if(allocated(this%Gij)) deallocate(this%Gij)
        allocate(this%Iij(this%no_sites))
        allocate(this%Gij(this%no_sites,this%no_sites))
    endif
    this%Iij = 0
    this%Gij = 0
endsubroutine init_frontline

subroutine free_frontline(this)
    class(frontline) :: this
    this%no_sites = 0
    if(allocated(this%Iij)) deallocate(this%Iij)
    if(allocated(this%Gij)) deallocate(this%Gij)
endsubroutine free_frontline
! -------------------------------------------------------------- !
!
! -------------------------------------------------------------- !
subroutine copy_frontline(this,scr_ft)
    class(frontline) :: this
    type(frontline) :: scr_ft
    call this%init_frontline(scr_ft%no_sites)
    this%Gij      = scr_ft%Gij
    this%Iij      = scr_ft%Iij
    this%no_sites = scr_ft%no_sites
endsubroutine copy_frontline

! -------------------------------------------------------------- !
!
! -------------------------------------------------------------- !
subroutine init_element(this,ns)
    class(element) :: this
    integer :: ns
    if(ns /= this%no_sites) then
        this%no_sites = ns
        if(allocated(this%Gij)) deallocate(this%Gij)
        allocate(this%Gij(this%no_sites,2))
    endif
!    this%Gij = 0
endsubroutine init_element

subroutine free_element(this)
    class(element) :: this
    this%no_sites = 0
    this%i        = 0
    this%j        = 0
    this%Gp       = 0
    this%Gax      = 0
    this%Gxa      = 0
    if(allocated(this%Gij)) deallocate(this%Gij)
endsubroutine free_element

subroutine copy_element(this,scr_ft)
    class(element) :: this
    type(element) :: scr_ft

    call this%init_element(scr_ft%no_sites)
    this%Gij      = scr_ft%Gij
    this%no_sites = scr_ft%no_sites
endsubroutine copy_element

end module modskminv

max_stars_repo_name: kmkolasinski/Bubel

max_issues_repo_name: kmkolasinski/Bubel

max_forks_repo_name: kmkolasinski/Bubel


text:
use plantfem
implicit none

type(IO_) :: f
character(:),allocatable :: file_path
real(real64),allocatable :: accel_wave(:,:)

file_path = "../Data/Earthquake/20070716/Knet/NIG0190707161013.EW"

accel_wave = f%import(format=KNT_ASCII,name=file_path)

call f%plot(x=accel_wave(:,1),Fx=accel_wave(:,2),option="with lines" )

end

max_stars_repo_name: kazulagi/faba

max_issues_repo_name: kazulagi/faba

max_forks_repo_name: kazulagi/faba


!***********************************/
!	Name:ユーザー指定の時間刻みがCFL条件に抵触する場合，自動的に時間刻みを変更するプログラム
!	Alias:CheckCFL4FixedTime
!	Description:変更された時間刻みは，元の時間刻みの整数分の1になるようにして，なるべく統一時間刻みに戻す
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.02.01
!	Update:
!	Other:
!***********************************/
subroutine CheckCFL4FixedTime(Geom,CC,iSplit)
use StructVar_Mod
use ConstantVar_Mod
implicit none

    type(Geometry),intent(in) :: Geom
    type(CellCenter), intent(inout) :: CC
    integer, intent(inout) :: iSplit

    call CalcTimeWidth(Geom,CC)
    iSplit = ceiling(DefaultTimeStep / minval(CC%TimeWidth(:,1,1),1)) !ceiling:天井= 切り上げ

return
end subroutine CheckCFL4FixedTime


text:
c--- NB: phi-dagger amplitudes related directly to phi amplitudes
c---     using parity and charge conjugation
      double complex function A0Hqarbmppm(j1,j2,j3,j4,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4
      double complex A0phiqarbmppm

      A0Hqarbmppm=A0phiqarbmppm(j1,j2,j3,j4,za,zb)
     .           +A0phiqarbmppm(j2,j1,j4,j3,zb,za)

      return
      end

      double complex function A0Hqarbmpmp(j1,j2,j3,j4,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4
      double complex A0phiqarbmpmp

      A0Hqarbmpmp=A0phiqarbmpmp(j1,j2,j3,j4,za,zb)
     .           +A0phiqarbmpmp(j2,j1,j4,j3,zb,za)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text: #redirect Campus Unions and Recreation Board


text:
      SUBROUTINE CHRGE(P,Q)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION P(*),Q(*)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
C***********************************************************************
C
C      CHRGE STORES IN Q THE TOTAL ELECTRON DENSITIES ON THE ATOMS
C
C      ON INPUT P      = DENSITY MATRIX
C
C      ON OUTPUT Q     = ATOM ELECTRON DENSITIES
C
C***********************************************************************
      K=0
      DO 10 I=1,NUMAT
         IA=NFIRST(I)
         IB=NLAST(I)
         Q(I)=0.D0
         DO 10 J=IA,IB
            K=K+J
   10 Q(I)=Q(I)+P(K)
      RETURN
      END

ext: for

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


program test_initialize

  use bmif_1_2, only: BMI_FAILURE
  use bmigiplf
  use fixtures, only: status, config_file

  implicit none

  type (bmi_gipl) :: m
  integer :: status1, status2, s
  real*8 :: end_time

  status1 = m%initialize(config_file)
  status = m%finalize()

  if (status1 .ne. 0) then
     stop BMI_FAILURE
  end if

end program test_initialize


  !
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 1998 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                          Spherepack                           *
!     *                                                               *
!     *       A Package of Fortran Subroutines and Programs           *
!     *                                                               *
!     *              for Modeling Geophysical Processes               *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *                  John Adams and Paul Swarztrauber             *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!
!
!
!     1/976
!
!     a program for testing all vector laplacian and its inverse
!
!     (1) first set the vector function (rotation about equator) using
!         v = cos(phi) and w = -cos(theta)*sin(phi)
!
!     (2) set vector laplacian ananlytically
!         vlap = -2.*cos(phi)=-2.*v, wlap = -2.*w
!         (i.e., L(v, w) = -2.*(v, w) so (v, w) is an eigenfunction for the
!         vector Laplacian with eigenvalue -2.
!
!     (3) compute the coefficients br, bi, cr, ci of (v, w) using vector analysis
!
!     (3) compute the vector laplacian of (v, w) using vlapec, vlapes, vlapgc, vlapgs
!
!     (4) compare (3) with (2)
!
!     (5) invert (4) and compare with (v, w)
!
program tvlap

    use, intrinsic :: ISO_Fortran_env, only: &
        stdout => OUTPUT_UNIT

    use spherepack, only: &
        ip, & ! Integer precision
        wp, & ! Working precision
        Sphere, &
        Regularsphere, &
        GaussianSphere

    ! Explicit typing only
    implicit none

    ! Dictionary
    class(Sphere), allocatable :: solver

    ! Test gaussian grid
    allocate (GaussianSphere :: solver)
    call test_vector_laplacian_routines(solver)
    deallocate (solver)

    ! Test regular grid
    allocate (RegularSphere :: solver)
    call test_vector_laplacian_routines(solver)
    deallocate (solver)

contains

    subroutine test_vector_laplacian_routines( sphere_type)

        ! Dummy arguments
        class(Sphere), intent(inout)  :: sphere_type

        ! Local variables
        integer(ip), parameter        :: NLONS = 16
        integer(ip), parameter        :: NLATS = 29
        integer(ip)                   :: i, j, k ! Counters
        real(wp)                      :: original_polar_component(NLATS, NLONS)
        real(wp)                      :: original_azimuthal_component(NLATS, NLONS)
        real(wp)                      :: exact_polar_laplacian(NLATS, NLONS)
        real(wp)                      :: exact_azimuthal_laplacian(NLATS, NLONS)
        real(wp)                      :: polar_component(NLATS, NLONS)
        real(wp)                      :: azimuthal_component(NLATS, NLONS)
        real(wp)                      :: approximate_polar_laplacian(NLATS, NLONS)
        real(wp)                      :: approximate_azimuthal_laplacian(NLATS, NLONS)
        character(len=:), allocatable :: previous_polar_laplacian_error
        character(len=:), allocatable :: previous_polar_inversion_error
        character(len=:), allocatable :: previous_azimuthal_laplacian_error
        character(len=:), allocatable :: previous_azimuthal_inversion_error

        !  Set up workspace arrays
        select type(sphere_type)
            type is (GaussianSphere)

            !  Initialize gaussian sphere object
            sphere_type = GaussianSphere(NLATS, NLONS)

            ! Allocate known error from previous platform
            allocate (previous_polar_laplacian_error, source='     polar laplacian error     = 1.761258e-12')
            allocate (previous_azimuthal_laplacian_error, source='     azimuthal laplacian error = 8.710810e-13')
            allocate (previous_polar_inversion_error, source='     polar inversion error     = 5.551115e-16')
            allocate (previous_azimuthal_inversion_error, source='     azimuthal inversion error = 7.216450e-16')

            type is (RegularSphere)

            ! Initialize regular sphere
            sphere_type = RegularSphere(NLATS, NLONS)

            ! Allocate known error from previous platform
            allocate (previous_polar_laplacian_error, source='     polar laplacian error     = 3.113065e-13')
            allocate (previous_azimuthal_laplacian_error, source='     azimuthal laplacian error = 4.702905e-13')
            allocate (previous_polar_inversion_error, source='     polar inversion error     = 5.551115e-16')
            allocate (previous_azimuthal_inversion_error, source='     azimuthal inversion error = 5.551115e-16')
        end select

        !
        !  test all vector laplacian and inverse vector laplacian subroutines
        !
        associate (&
            r => sphere_type%unit_vectors%radial, &
            phi => sphere_type%unit_vectors%azimuthal, &
            ve => original_polar_component, &
            we => original_azimuthal_component, &
            velap => exact_polar_laplacian, &
            welap => exact_azimuthal_laplacian &
           )
            do j=1, NLONS
                do i=1, NLATS
                    associate (&
                        cost => r(i, j)%z, &
                        cosp => phi(i, j)%y, &
                        sinp => -phi(i, j)%x &
                       )
                        !
                        !  set vector field v, w
                        !
                        ve(i, j) = cosp
                        we(i, j) = -cost*sinp
                        !
                        !  set vector laplacian vlap, wlap
                        !
                        velap(i, j) = -2.0_wp * ve(i, j)
                        welap(i, j) = -2.0_wp * we(i, j)
                    end associate
                end do
            end do
        end associate

        associate (&
            ve => original_polar_component, &
            we => original_azimuthal_component, &
            vlap => approximate_polar_laplacian, &
            wlap => approximate_azimuthal_laplacian &
           )
            !
            !  Compute vector laplacian
            !
            call sphere_type%get_laplacian(ve, we, vlap, wlap)

        end associate

        !
        !  Compute laplacian error
        !
        associate (&
            velap => exact_polar_laplacian, &
            welap => exact_azimuthal_laplacian, &
            vlap => approximate_polar_laplacian, &
            wlap => approximate_azimuthal_laplacian &
           )
            associate (&
                err2v => maxval(abs(vlap - velap)), &
                err2w => maxval(abs(wlap - welap)) &
               )

                !
                !  Print earlier output from platform with 64-bit floating point
                !    arithmetic followed by the output from this computer
                !
                write (stdout, '(/a/)') '     tvlap *** TEST RUN *** '
                write (stdout, '(a)') '     grid type = '//sphere_type%grid%grid_type
                write (stdout, '(a)') '     Testing vector laplacian'
                write (stdout, '(2(a, i3))') '     nlat = ', NLATS, ' nlon = ', NLONS
                write (stdout, '(a)') '     Previous 64 bit floating point arithmetic result '
                write (stdout, '(a)') previous_polar_laplacian_error
                write (stdout, '(a)') previous_azimuthal_laplacian_error
                write (stdout, '(a)') '     The output from your computer is: '
                write (stdout, '(a, 1pe15.6)') '     polar laplacian error     = ', err2v
                write (stdout, '(a, 1pe15.6/)') '     azimuthal laplacian error = ', err2w
            end associate
        end associate
        !
        !  Now recompute (v, w) inverting (velap, welap)
        !
        associate (&
            v => polar_component, &
            w => azimuthal_component, &
            velap => exact_polar_laplacian, &
            welap => exact_azimuthal_laplacian &
           )
            call sphere_type%invert_laplacian( velap, welap, v, w)
        end associate

        !
        !  compare this v, w with original
        !
        associate (&
            ve => original_polar_component, &
            we => original_azimuthal_component, &
            v => polar_component, &
            w => azimuthal_component &
           )
            associate (&
                err2v => maxval(abs(v - ve)), &
                err2w => maxval(abs(w - we)) &
               )
                !
                !  Print earlier output from platform with 64-bit floating point
                !    arithmetic followed by the output from this computer
                !
                write (stdout, '(/a/)') '     tvlap *** TEST RUN *** '
                write (stdout, '(a)') '     grid type = '//sphere_type%grid%grid_type
                write (stdout, '(a)') '     Testing vector laplacian inversion'
                write (stdout, '(2(a, i3))') '     nlat = ', NLATS, ' nlon = ', NLONS
                write (stdout, '(a)') '     Previous 64 bit floating point arithmetic result '
                write (stdout, '(a)') previous_polar_inversion_error
                write (stdout, '(a)') previous_azimuthal_inversion_error
                write (stdout, '(a)') '     The output from your computer is: '
                write (stdout, '(a, 1pe15.6)') '     polar inversion error     = ', err2v
                write (stdout, '(a, 1pe15.6/)') '     azimuthal inversion error = ', err2w
            end associate
        end associate
        !
        !  Release memory
        deallocate (previous_polar_laplacian_error)
        deallocate (previous_polar_inversion_error)
        deallocate (previous_azimuthal_laplacian_error)
        deallocate (previous_azimuthal_inversion_error)
        call sphere_type%destroy()

    end subroutine test_vector_laplacian_routines

end program tvlap


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

Welcome to the Wiki Howdy, Mr. or Miss Ultomega, and Welcome to the Wiki!   You might want to check out the importance of using your RealName, just so we can get to know you (or not: its your choice, but people were pretty friendly here).  My names Evan, pleased to meet you, and congratulations on your engagement, and I hope the best for your wedding and marriage!  Users/JabberWokky Evan JabberWokky Edwards

20080213 20:04:42 nbsp I am sorry you went some where els. but custom design jewelry is more expensive i also give a life time warranty with all my diamond jewelry is that something you received at the other jewelers. ok well best of luck and in the future if there is some thing i can do come see us. Thank you for your feed back.  Users/3rdstreetjeweler