text:
module a2c
   IMPLICIT NONE

   ! cell: x(3), varh(5~7)
   real(kind=8),    dimension(:,:),   allocatable::real_table
   ! cell: level, cpu
   integer(kind=4), dimension(:,:),   allocatable::integer_table

   ! part: long id
   integer(kind=8), dimension(:,:),   allocatable::long_table
   ! part: family, tag
   integer(kind=1), dimension(:,:),   allocatable::byte_table

contains

subroutine a2c_count(ngridtot, nvarh, repository, xmin, xmax, ymin, &
                  & ymax, zmin, zmax, lmax, cpu_list)
  !--------------------------------------------------------------------------
  ! Ce programme calcule le cube cartesien pour les
  ! variables hydro d'une simulation RAMSES.
  ! Version F90 par R. Teyssier le 01/04/01.
  !--------------------------------------------------------------------------
  character(LEN=*),INTENT(IN)::repository
  character(LEN=128)::nomfich
  real(KIND=8),INTENT(IN)::xmin,xmax,ymin,ymax,zmin,zmax
  integer     ,INTENT(IN)::lmax
  integer,dimension(:),INTENT(IN)::cpu_list
  integer     ,INTENT(OUT)::ngridtot
  ! default value from python side is 0.

  integer::ndim,i,j,k,twotondim
  integer::ivar,ncpu,nboundary
  integer::nx,ny,nz,ilevel,idim
  integer::nlevelmax,ngrida
  integer::ind,ipos
  integer::ngridmax,icpu,ncpu_read
  real::boxlen

  integer::imin,imax,jmin,jmax,kmin,kmax
  integer,INTENT(OUT)::nvarh
  integer::ix,iy,iz
  integer::nx_full,ny_full,nz_full
  real(KIND=8)::xxmin,xxmax,yymin,yymax,zzmin,zzmax,dx,dx2
  real(KIND=8),dimension(:,:),allocatable::x,xg

  logical,dimension(:)  ,allocatable::ref
  integer,dimension(:,:),allocatable::son,ngridfile,ngridlevel,ngridbound
  real(KIND=8),dimension(1:8,1:3)::xc
  real(KIND=8),dimension(1:3)::xbound=(/0d0,0d0,0d0/)
  character(LEN=5)::nchar,ncharcpu
  character(LEN=80)::ordering
  logical::ok,ok_cell

  type level
     integer::ilevel
     integer::ngrid
     integer::imin
     integer::imax
     integer::jmin
     integer::jmax
     integer::kmin
     integer::kmax
  end type level

  type(level),dimension(1:50)::grid

  ncpu_read = size(cpu_list)

  !-----------------------------------------------
  ! Lecture du fichier hydro au format RAMSES
  !-----------------------------------------------
  ipos=INDEX(repository,'output_')
  nchar=repository(ipos+7:ipos+13)
  nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif
  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif

  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  open(unit=10,file=nomfich,status='old',form='unformatted')
  read(10)ncpu
  read(10)ndim
  read(10)nx,ny,nz
  read(10)nlevelmax
  read(10)ngridmax
  read(10)nboundary
  read(10)
  read(10)boxlen
  close(10)
  twotondim=2**ndim
  xbound=(/dble(nx/2),dble(ny/2),dble(nz/2)/)

  allocate(ngridfile(1:ncpu+nboundary,1:nlevelmax))
  allocate(ngridlevel(1:ncpu,1:nlevelmax))
  if(nboundary>0)allocate(ngridbound(1:nboundary,1:nlevelmax))

  xxmin=xmin ; xxmax=xmax
  yymin=ymin ; yymax=ymax
  zzmin=zmin ; zzmax=zmax

  !-----------------------------
  ! Compute hierarchy
  !-----------------------------
  do ilevel=1,lmax
     nx_full=2**ilevel
     ny_full=2**ilevel
     nz_full=2**ilevel
     imin=int(xxmin*dble(nx_full))+1
     imax=int(xxmax*dble(nx_full))+1
     jmin=int(yymin*dble(ny_full))+1
     jmax=int(yymax*dble(ny_full))+1
     kmin=int(zzmin*dble(nz_full))+1
     kmax=int(zzmax*dble(nz_full))+1
     grid(ilevel)%imin=imin
     grid(ilevel)%imax=imax
     grid(ilevel)%jmin=jmin
     grid(ilevel)%jmax=jmax
     grid(ilevel)%kmin=kmin
     grid(ilevel)%kmax=kmax
  end do

  !-----------------------------------------------
  ! Compute array size
  !----------------------------------------------
  ! Loop over processor files
  ngridtot = 0
  do k=1,ncpu_read
     icpu=cpu_list(k)
     call title(icpu,ncharcpu)

     ! Open AMR file and skip header
     nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
     open(unit=10,file=nomfich,status='old',form='unformatted')
     do i=1,21
        read(10)
     end do
     ! Read grid numbers
     read(10)ngridlevel
     ngridfile(1:ncpu,1:nlevelmax)=ngridlevel
     read(10)
     if(nboundary>0)then
        do i=1,2
           read(10)
        end do
        read(10)ngridbound
        ngridfile(ncpu+1:ncpu+nboundary,1:nlevelmax)=ngridbound
     endif
     read(10)
  ! ROM: comment the single follwing line for old stuff
     read(10)
     if(TRIM(ordering).eq.'bisection')then
        do i=1,5
           read(10)
        end do
     else
        read(10)
     endif
     read(10)
     read(10)
     read(10)

     ! Open HYDRO file and skip header
     nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
     open(unit=11,file=nomfich,status='old',form='unformatted')
     read(11)
     read(11)nvarh
     read(11)
     read(11)
     read(11)
     read(11)

     ! Loop over levels
     do ilevel=1,lmax

        ! Geometry
        dx=0.5**ilevel
        dx2=0.5*dx
        nx_full=2**ilevel
        ny_full=2**ilevel
        nz_full=2**ilevel
        do ind=1,twotondim
           iz=(ind-1)/4
           iy=(ind-1-4*iz)/2
           ix=(ind-1-2*iy-4*iz)
           xc(ind,1)=(dble(ix)-0.5D0)*dx
           xc(ind,2)=(dble(iy)-0.5D0)*dx
           xc(ind,3)=(dble(iz)-0.5D0)*dx
        end do

        ! Allocate work arrays
        ngrida=ngridfile(icpu,ilevel)
        grid(ilevel)%ngrid=ngrida
        if(ngrida>0)then
           allocate(xg(1:ngrida,1:ndim))
           allocate(son(1:ngrida,1:twotondim))
           allocate(x  (1:ngrida,1:ndim))
           allocate(ref(1:ngrida))
        endif

        ! Loop over domains
        do j=1,nboundary+ncpu

           ! Read AMR data
           if(ngridfile(j,ilevel)>0)then
              read(10) ! Skip grid index
              read(10) ! Skip next index
              read(10) ! Skip prev index
              ! Read grid center
              do idim=1,ndim
                 if(j.eq.icpu)then
                    read(10)xg(:,idim)
                 else
                    read(10)
                 endif
              end do
              read(10) ! Skip father index
              do ind=1,2*ndim
                 read(10) ! Skip nbor index
              end do
              ! Read son index
              do ind=1,twotondim
                 if(j.eq.icpu)then
                    read(10)son(:,ind)
                 else
                    read(10)
                 end if
              end do
              ! Skip cpu map
              do ind=1,twotondim
                 read(10)
              end do
              ! Skip refinement map
              do ind=1,twotondim
                 read(10)
              end do
           endif

           ! Read HYDRO data
           read(11)
           read(11)
           if(ngridfile(j,ilevel)>0)then
              ! Read hydro variables
              do ind=1,twotondim
                 do ivar=1,nvarh
                    read(11)
                 end do
              end do
           end if
        end do

        ! Compute map
        if(ngrida>0)then

           ! Loop over cells
           do ind=1,twotondim

              ! Compute cell center
              do i=1,ngrida
                 x(i,1)=(xg(i,1)+xc(ind,1)-xbound(1))
                 x(i,2)=(xg(i,2)+xc(ind,2)-xbound(2))
                 x(i,3)=(xg(i,3)+xc(ind,3)-xbound(3))
              end do
              ! Check if cell is refined
              do i=1,ngrida
                 ref(i)=son(i,ind)>0.and.ilevel<lmax
              end do
              ! Store data cube
              do i=1,ngrida
                 ok_cell= .not.ref(i).and. &
                      & (x(i,1)+dx2)>=xmin.and.&
                      & (x(i,2)+dx2)>=ymin.and.&
                      & (x(i,3)+dx2)>=zmin.and.&
                      & (x(i,1)-dx2)<=xmax.and.&
                      & (x(i,2)-dx2)<=ymax.and.&
                      & (x(i,3)-dx2)<=zmax
                 if(ok_cell)then
                   ngridtot = ngridtot + 1
                 end if
              end do

           end do
           ! End loop over cell

           deallocate(ref,x,son,xg)
        endif

     end do
     ! End loop over levels

     close(10)
     close(11)
  end do
  deallocate(ngridfile,ngridlevel)
  if(nboundary>0) deallocate(ngridbound)
end subroutine a2c_count

subroutine a2c_load(repository &
                    &, xmin, xmax, ymin, ymax, zmin, zmax, lmax &
                    &, ngridtot, nvarh, cpu_list)
  implicit none

  integer,dimension(:),INTENT(IN)::cpu_list
  character(len=*),intent(IN)::repository
  character(len=128)::nomfich
  real(kind=8),intent(IN)::xmin,xmax,ymin,ymax,zmin,zmax
  integer     ,intent(IN)::lmax,ngridtot, nvarh
  !default value from python side is 0.

  !real(kind=8),intent(inout),dimension(ngridtot,3)::xarr
  !real(kind=8),intent(inout),dimension(ngridtot,nvarh)::varr
  !real(kind=8),intent(inout),dimension(ngridtot)::dxarr
  !integer,intent(inout),dimension(ngridtot)::cpuarr,refarr
  !logical,intent(inout),dimension(ngridtot)::refarr
  !f2py depend(ngridtot) xarr,dxarr,cpuarr,refarr
  !f2py depend(ngridtot, nvarh) varr

  integer::ndim,i,j,k,twotondim
  integer::ivar,ncpu,nboundary
  integer::nx,ny,nz,ilevel,idim
  integer::nlevelmax
  integer::ind,ipos,ngrida,icnt
  integer::ngridmax,icpu,ncpu_read
  !real::boxlen,t

  integer::imin,imax,jmin,jmax,kmin,kmax
  integer::nvarh_org
  integer::ix,iy,iz
  integer::nx_full,ny_full,nz_full
  real(kind=8)::xxmin,xxmax,yymin,yymax,zzmin,zzmax,dx,dx2
  real(kind=8),dimension(:,:),allocatable::x,xg
  real(kind=8),dimension(:,:,:),allocatable::var

  !real(kind=8),dimension(:)  ,allocatable::rho
  logical,dimension(:)  ,allocatable::ref
  integer,dimension(:,:),allocatable::son,ngridfile,ngridlevel,ngridbound
  real(kind=8),dimension(1:8,1:3)::xc
  real(kind=8),dimension(1:3)::xbound=(/0d0,0d0,0d0/)
  character(len=5)::nchar,ncharcpu
  character(len=80)::ordering
  logical::ok,ok_cell

  type level
     integer::ilevel
     integer::ngrid
     real(KIND=4),dimension(:,:,:),pointer::cube
     integer::imin
     integer::imax
     integer::jmin
     integer::jmax
     integer::kmin
     integer::kmax
  end type level

  type(level),dimension(1:50)::grid

  !-----------------------------------------------
  ! Lecture du fichier hydro au format RAMSES
  !-----------------------------------------------
  ipos=INDEX(repository,'output_')
  nchar=repository(ipos+7:ipos+13)
  nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif
  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif

  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  open(unit=10,file=nomfich,status='old',form='unformatted')
  read(10)ncpu
  read(10)ndim
  read(10)nx,ny,nz
  read(10)nlevelmax
  read(10)ngridmax
  read(10)nboundary

  close(10)
  twotondim=2**ndim
  xbound=(/dble(nx/2),dble(ny/2),dble(nz/2)/)

  allocate(ngridfile(1:ncpu+nboundary,1:nlevelmax))
  allocate(ngridlevel(1:ncpu,1:nlevelmax))
  if(nboundary>0)allocate(ngridbound(1:nboundary,1:nlevelmax))

  xxmin=xmin ; xxmax=xmax
  yymin=ymin ; yymax=ymax
  zzmin=zmin ; zzmax=zmax

  !-----------------------------
  ! Compute hierarchy
  !-----------------------------
  do ilevel=1,lmax
     nx_full=2**ilevel
     ny_full=2**ilevel
     nz_full=2**ilevel
     imin=int(xxmin*dble(nx_full))+1
     imax=int(xxmax*dble(nx_full))+1
     jmin=int(yymin*dble(ny_full))+1
     jmax=int(yymax*dble(ny_full))+1
     kmin=int(zzmin*dble(nz_full))+1
     kmax=int(zzmax*dble(nz_full))+1
     grid(ilevel)%imin=imin
     grid(ilevel)%imax=imax
     grid(ilevel)%jmin=jmin
     grid(ilevel)%jmax=jmax
     grid(ilevel)%kmin=kmin
     grid(ilevel)%kmax=kmax
  end do

  !------------------------------------------------
  ! Allocate arrays
  !------------------------------------------------
  icnt=0

  allocate(real_table(1:ngridtot, 1:ndim+nvarh))
  allocate(integer_table(1:ngridtot, 1:2))
  print *, "allocated"
  ! Loop over processor files
  ncpu_read = size(cpu_list)
  do k=1,ncpu_read
     icpu=cpu_list(k)
     call title(icpu,ncharcpu)
     print *, "CPU...", icpu
     ! Open AMR file and skip header
     nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
     open(unit=10,file=nomfich,status='old',form='unformatted')
     ! write(*,*)'Processing file '//TRIM(nomfich)
     do i=1,21
        read(10)
     end do
     ! Read grid numbers
     read(10)ngridlevel
     ngridfile(1:ncpu,1:nlevelmax)=ngridlevel
     read(10)
     if(nboundary>0)then
        do i=1,2
           read(10)
        end do
        read(10)ngridbound
        ngridfile(ncpu+1:ncpu+nboundary,1:nlevelmax)=ngridbound
     endif
     read(10)
  ! ROM: comment the single follwing line for old stuff
     read(10)
     if(TRIM(ordering).eq.'bisection')then
        do i=1,5
           read(10)
        end do
     else
        read(10)
     endif
     read(10)
     read(10)
     read(10)

     ! Open HYDRO file and skip header
     nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
     open(unit=11,file=nomfich,status='old',form='unformatted')
     read(11)
     read(11)nvarh_org
     read(11)
     read(11)
     read(11)
     read(11)

     ! Loop over levels
     do ilevel=1,lmax

      ! Geometry
        dx=0.5**ilevel
        dx2=0.5*dx
        nx_full=2**ilevel
        ny_full=2**ilevel
        nz_full=2**ilevel
        do ind=1,twotondim
           iz=(ind-1)/4
           iy=(ind-1-4*iz)/2
           ix=(ind-1-2*iy-4*iz)
           xc(ind,1)=(dble(ix)-0.5D0)*dx
           xc(ind,2)=(dble(iy)-0.5D0)*dx
           xc(ind,3)=(dble(iz)-0.5D0)*dx
        end do
      print *, "second allocation", ilevel, lmax, icnt
      ! Allocate work arrays
        ngrida=ngridfile(icpu,ilevel)
        grid(ilevel)%ngrid=ngrida
        if(ngrida>0)then
           allocate(xg(1:ngrida,1:ndim))
           allocate(son(1:ngrida,1:twotondim))
           allocate(var(1:ngrida,1:twotondim,1:nvarh))
           allocate(x  (1:ngrida,1:ndim))
           allocate(ref(1:ngrida))
        endif

        ! Loop over domains
        do j=1,nboundary+ncpu

           ! Read AMR data
           if(ngridfile(j,ilevel)>0)then
              read(10) ! Skip grid index
              read(10) ! Skip next index
              read(10) ! Skip prev index
              ! Read grid center
              do idim=1,ndim
                 if(j.eq.icpu)then
                    read(10)xg(:,idim)
                 else
                    read(10)
                 endif
              end do
              read(10) ! Skip father index
              do ind=1,2*ndim
                 read(10) ! Skip nbor index
              end do
              ! Read son index
              do ind=1,twotondim
                 if(j.eq.icpu)then
                    read(10)son(:,ind)
                 else
                    read(10)
                 end if
              end do
              ! Skip cpu map
              do ind=1,twotondim
                 read(10)
              end do
              ! Skip refinement map
              do ind=1,twotondim
                 read(10)
              end do
           endif

           ! Read HYDRO data
           read(11)
           read(11)
           if(ngridfile(j,ilevel)>0)then
              ! Read hydro variables
              do ind=1,twotondim
                 do ivar=1,nvarh_org
                    if(j.eq.icpu .and. ivar .le. nvarh)then
                       read(11)var(:,ind,ivar)
                    else
                       read(11)
                    end if
                 end do
              end do
           end if
        end do

        ! Compute map
        if(ngrida>0)then

           ! Loop over cells
           do ind=1,twotondim

              ! Compute cell center
              do i=1,ngrida
                 x(i,1)=(xg(i,1)+xc(ind,1)-xbound(1))
                 x(i,2)=(xg(i,2)+xc(ind,2)-xbound(2))
                 x(i,3)=(xg(i,3)+xc(ind,3)-xbound(3))
              end do
              ! Check if cell is refined
              do i=1,ngrida
                 ref(i)=son(i,ind)>0.and.ilevel<lmax
              end do
              ! Store data cube
              do i=1,ngrida
                 ok_cell= .not.ref(i).and. &
                      & (x(i,1)+dx2)>=xmin.and.&
                      & (x(i,2)+dx2)>=ymin.and.&
                      & (x(i,3)+dx2)>=zmin.and.&
                      & (x(i,1)-dx2)<=xmax.and.&
                      & (x(i,2)-dx2)<=ymax.and.&
                      & (x(i,3)-dx2)<=zmax
                 if(ok_cell)then
                    icnt = icnt + 1
                    real_table(icnt,1:3)=x(i,1:3)
                    real_table(icnt,4)=dx
                    integer_table(icnt,1) = icpu
                    print *, "3 ok", i, ind, nvarh

                    real_table(icnt,5:5+nvarh-1)=var(i,ind,1:nvarh)
                    print *, real_table(icnt,1:5+nvarh)
                    print *, "4 ok", i, ind, nvarh
                    !refarr(icnt) = ref(i) Do I ever use this??
                 end if
              end do

           end do
           ! End loop over cell

           deallocate(xg,son,var,ref,x)
        endif

     end do
     ! End loop over levels
     print *, "loop done", ilevel

     close(10)
     close(11)

  end do
  ! End loop over cpus

  deallocate(ngridfile,ngridlevel)
  if(nboundary>0) deallocate(ngridbound)
end subroutine

!=======================================================================
subroutine title(n,nchar)
!=======================================================================
  implicit none
  integer::n
  character*5::nchar

  character*1::nchar1
  character*2::nchar2
  character*3::nchar3
  character*4::nchar4
  character*5::nchar5

  if(n.ge.10000)then
     write(nchar5,'(i5)') n
     nchar = nchar5
  elseif(n.ge.1000)then
     write(nchar4,'(i4)') n
     nchar = '0'//nchar4
  elseif(n.ge.100)then
     write(nchar3,'(i3)') n
     nchar = '00'//nchar3
  elseif(n.ge.10)then
     write(nchar2,'(i2)') n
     nchar = '000'//nchar2
  else
     write(nchar1,'(i1)') n
     nchar = '0000'//nchar1
  endif

end subroutine title
!================================================================
!================================================================
!================================================================
!================================================================
subroutine hilbert3d(x,y,z,order,bit_length,npoint)
  implicit none

  integer     ,INTENT(IN)                     ::bit_length,npoint
  integer     ,INTENT(IN) ,dimension(1:npoint)::x,y,z
  real(kind=8),INTENT(OUT),dimension(1:npoint)::order
  !  real(kind=8),INTENT(OUT)::order

  logical,dimension(0:3*bit_length-1)::i_bit_mask
  logical,dimension(0:1*bit_length-1)::x_bit_mask,y_bit_mask,z_bit_mask
  integer,dimension(0:7,0:1,0:11)::state_diagram
  integer::i,ip,cstate,nstate,b0,b1,b2,sdigit,hdigit

  if(bit_length>bit_size(bit_length))then
     write(*,*)'Maximum bit length=',bit_size(bit_length)
     write(*,*)'stop in hilbert3d'
     stop
  endif

  state_diagram = RESHAPE( (/   1, 2, 3, 2, 4, 5, 3, 5,&
                            &   0, 1, 3, 2, 7, 6, 4, 5,&
                            &   2, 6, 0, 7, 8, 8, 0, 7,&
                            &   0, 7, 1, 6, 3, 4, 2, 5,&
                            &   0, 9,10, 9, 1, 1,11,11,&
                            &   0, 3, 7, 4, 1, 2, 6, 5,&
                            &   6, 0, 6,11, 9, 0, 9, 8,&
                            &   2, 3, 1, 0, 5, 4, 6, 7,&
                            &  11,11, 0, 7, 5, 9, 0, 7,&
                            &   4, 3, 5, 2, 7, 0, 6, 1,&
                            &   4, 4, 8, 8, 0, 6,10, 6,&
                            &   6, 5, 1, 2, 7, 4, 0, 3,&
                            &   5, 7, 5, 3, 1, 1,11,11,&
                            &   4, 7, 3, 0, 5, 6, 2, 1,&
                            &   6, 1, 6,10, 9, 4, 9,10,&
                            &   6, 7, 5, 4, 1, 0, 2, 3,&
                            &  10, 3, 1, 1,10, 3, 5, 9,&
                            &   2, 5, 3, 4, 1, 6, 0, 7,&
                            &   4, 4, 8, 8, 2, 7, 2, 3,&
                            &   2, 1, 5, 6, 3, 0, 4, 7,&
                            &   7, 2,11, 2, 7, 5, 8, 5,&
                            &   4, 5, 7, 6, 3, 2, 0, 1,&
                            &  10, 3, 2, 6,10, 3, 4, 4,&
                            &   6, 1, 7, 0, 5, 2, 4, 3 /), &
                            & (/8 ,2, 12 /) )

  do ip=1,npoint

     ! convert to binary
     do i=0,bit_length-1
        x_bit_mask(i)=btest(x(ip),i)
        y_bit_mask(i)=btest(y(ip),i)
        z_bit_mask(i)=btest(z(ip),i)
     enddo

     ! interleave bits
     do i=0,bit_length-1
        i_bit_mask(3*i+2)=x_bit_mask(i)
        i_bit_mask(3*i+1)=y_bit_mask(i)
        i_bit_mask(3*i  )=z_bit_mask(i)
     end do

     ! build Hilbert ordering using state diagram
     cstate=0
     do i=bit_length-1,0,-1
        b2=0 ; if(i_bit_mask(3*i+2))b2=1
        b1=0 ; if(i_bit_mask(3*i+1))b1=1
        b0=0 ; if(i_bit_mask(3*i  ))b0=1
        sdigit=b2*4+b1*2+b0
        nstate=state_diagram(sdigit,0,cstate)
        hdigit=state_diagram(sdigit,1,cstate)
        i_bit_mask(3*i+2)=btest(hdigit,2)
        i_bit_mask(3*i+1)=btest(hdigit,1)
        i_bit_mask(3*i  )=btest(hdigit,0)
        cstate=nstate
     enddo

     ! save Hilbert key as double precision real
     order(ip)=0.
     do i=0,3*bit_length-1
        b0=0 ; if(i_bit_mask(i))b0=1
        order(ip)=order(ip)+dble(b0)*dble(2)**i
     end do

  end do

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

!#####################################################################
    subroutine clear_all()
!#####################################################################
        ! Clean old data table is used more then once
        implicit none
        if(allocated(real_table)) deallocate(real_table)
        if(allocated(integer_table)) deallocate(integer_table)
        if(allocated(byte_table)) deallocate(byte_table)
        if(allocated(long_table)) deallocate(long_table)
    end subroutine clear_all

end module a2c

max_stars_repo_name: Hoseung/pyRamAn

max_issues_repo_name: Hoseung/pyRamAn

max_forks_repo_name: Hoseung/pyRamAn


text:
LOGICAL FUNCTION KFINDX( FNAME, COL, ROW, ECOUNT, SDATES, STIMES, KFLENS, EVENTS )  &
        RESULT( KFFLAG )

    !!***********************************************************************
    !! Version "$Id: kfindx.f90 1 2017-06-10 18:05:20Z coats $"
    !! EDSS/Models-3 I/O API.
    !! Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr.,
    !! (C) 2003-2013 Baron Advanced Meteorological Systems,
    !! (C) 2007-2013 Carlie J. Coats, Jr., and
    !! (C) 2014-2015 UNC Institute for the Environment.
    !! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
    !! See file "LGPL.txt" for conditions of use.
    !!.........................................................................
    !!  function body starts at line  82
    !!
    !!  FUNCTION:  reads the event descriptions for the indicated cell
    !! 	from the KF-Cloud-Event file with logical name FNAME.
    !!
    !!  RETURN VALUE:  TRUE iff the operation succeeds (and the data is available)
    !!
    !!  PRECONDITIONS REQUIRED:  FNAME is a  KF-Cloud-Event file already
    !!       opened by kfindx()
    !!
    !!  REVISION  HISTORY:
    !!      Adapted  4/1996 by CJC from READ3().
    !!
    !!      Modified  5/1998 by CJC for OpenMP thread-safety
    !!
    !!      Modified  1/2002 by CJC:  check TRIMLEN() of FNAME
    !!
    !!      Modified 7/2003 by CJC:  bugfix -- clean up critical sections
    !!      associated with INIT3()
    !!
    !!      Modified 03/2010 by CJC: F9x changes for I/O API v3.1
    !!
    !!      Modified 02/2015 by CJC for I/O API 3.2: USE M3UTILIO
    !!
    !!      Modified 10/2015 by CJC for I/O API 3.2: USE MODNCFIO,
    !!      use NF_*() interfaces.  F90 "free" source-format
    !!***********************************************************************

    USE M3UTILIO
    USE MODNCFIO

    IMPLICIT NONE

    !!...........   INCLUDES:

    INCLUDE 'STATE3.EXT'

    !!...........   ARGUMENTS and their descriptions:

    CHARACTER*(*), INTENT(IN   ) :: FNAME      !  logical file name
    INTEGER      , INTENT(IN   ) :: COL        !  column number for this event
    INTEGER      , INTENT(IN   ) :: ROW        !  row    number for this event
    INTEGER      , INTENT(  OUT) :: ECOUNT     !  # of events for this col-row
    INTEGER      , INTENT(  OUT) :: SDATES(*)  !  starting date,  formatted YYYYDDD
    INTEGER      , INTENT(  OUT) :: STIMES(*)  !  starting time,  formatted HHMMSS
    INTEGER      , INTENT(  OUT) :: KFLENS(*)  !  event duration, formatted HHMMSS
    INTEGER      , INTENT(  OUT) :: EVENTS(*)  !  event numbers

    !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

    INTEGER         FID             !  subscript  for STATE3 arrays
    INTEGER         FNUM            !  CDFID3( FID )
    INTEGER         IERR            !  netCDF error status return
    INTEGER         DIMT( 5 )       !  corner   for NF_GET_VARA_*()
    INTEGER         DELT( 5 )       !  diagonal for NF_GET_VARA_*()
    CHARACTER*16    FIL16           !  scratch file-name buffer
    CHARACTER*256   MESG            !  scratch message buffer
    LOGICAL         EFLAG

    !!***********************************************************************
    !!   begin body of function  KFINDX

    !!.......   Check that Models-3 I/O has been initialized:

    EFLAG = .FALSE.
!$OMP   CRITICAL( S_INIT )
    IF ( .NOT. FINIT3 ) THEN
        LOGDEV = INIT3()
        EFLAG   = .TRUE.
    END IF          !  if not FINIT3
!$OMP   END CRITICAL( S_INIT )
    IF ( EFLAG ) THEN
        CALL M3MSG2(  'KFINDX:  I/O API not yet initialized.' )
        KFFLAG = .FALSE.
        RETURN
    END IF

    IF ( LEN_TRIM( FNAME ) .GT. 16 ) THEN
        MESG = 'File "'// FNAME //'"'
        CALL M3MSG2( MESG )
        WRITE( MESG, '( A, I10 )' ) 'Max file name length 16; actual:', LEN_TRIM( FNAME )
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        RETURN
    END IF          !  if len( fname ) > 16, or if len( vname ) > 16

    !!.......   Find netCDF index for the file, and check time step availability:

    FIL16 = FNAME   !  fixed-length-16 scratch copy of name
    FID   = INDEX1( FIL16, COUNT3, FLIST3 )

    IF ( FID .EQ. 0 ) THEN  !  file not available
        MESG = 'File:  '//FIL16// ' not yet opened.'
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        RETURN
    END IF          !  if file not available

    !!.......   Check col/row for this cell:

    IF      ( COL .LT. 1 .OR. COL .GT. NCOLS3( FID ) ) THEN
        WRITE( MESG,91010 )                             &
                'Column requested:    ', COL,           &
                'out of range; max is ', NCOLS3( FID )
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        RETURN
    ELSE IF ( ROW .LT. 1 .OR. ROW .GT. NROWS3( FID ) ) THEN
        WRITE( MESG,91010 )                             &
                'Row requested:       ', ROW,           &
                'out of range; max is ', NROWS3( FID )
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        RETURN
    END IF

    FNUM = CDFID3( FID )

!$OMP   CRITICAL( S_NC )

    IF ( VOLAT3( FID ) ) THEN     !  volatile file:  synch with disk

        IERR = NF_SYNC( CDFID3( FID ) )
        IF ( IERR .NE. 0 ) THEN
            WRITE( MESG,91010 ) 'netCDF error number', IERR,    &
                'Error with disk synchronization for file:  ' // FIL16
            CALL M3WARN( 'KFINDX', 0,0, MESG )
            KFFLAG = .FALSE.
            GO TO  999        !  return from kfindx()
        END IF      !  if ncsnc() failed

    END IF          !  if file not available, or if file is volatile

    !!.......   Read number of events for this cell:

    DIMT( 1 ) = COL
    DIMT( 2 ) = ROW
    IERR = NF_GET_VAR1_INT( FNUM,  NINDX3( FID ), DIMT, ECOUNT )
    IF ( IERR .NE. 0 ) THEN

        WRITE( MESG,91010 ) 'Requested col:  ', COL, 'Requested row:  ', ROW
        CALL M3MESG( MESG )
        WRITE( MESG,91010 ) 'netCDF error number', IERR, 'netCDF ID:  ', FNUM
        CALL M3MESG( MESG )
        MESG = 'Error reading KFCOUNT from ' // FIL16
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        GO TO  999        !  return from kfindx()

    END IF

    !!.......   Read the event-specification parameters for ths cell:

    DIMT( 1 ) = 1
    DELT( 1 ) = ECOUNT
    DIMT( 2 ) = COL
    DELT( 2 ) = 1
    DIMT( 3 ) = ROW
    DELT( 3 ) = 1

    IERR = NF_GET_VARA_INT( FNUM, SINDX3( FID ), DIMT, DELT, EVENTS )
    IF ( IERR .NE. 0 ) THEN     !  timestep flag not yet written

        WRITE( MESG,91010 ) 'Requested col:  ', COL, 'Requested row:  ', ROW
        CALL M3MESG( MESG )
        WRITE( MESG,91010 ) 'netCDF error number', IERR, 'netCDF ID:  ', FNUM
        CALL M3MESG( MESG )
        MESG = 'Error reading KFEVENT from ' // FIL16
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        GO TO  999        !  return from kfindx()

    END IF          !  if ierr bad or if timestep flags bad

    IERR = NF_GET_VARA_INT( FNUM, LINDX3( FID ), DIMT, DELT, SDATES )
    IF ( IERR .NE. 0 ) THEN     !  timestep flag not yet written

        WRITE( MESG,91010 ) 'Requested col:  ', COL, 'Requested row:  ', ROW
        CALL M3MESG( MESG )
        WRITE( MESG,91010 ) 'netCDF error number', IERR, 'netCDF ID:  ', FNUM
        CALL M3MESG( MESG )
        MESG = 'Error reading KFSDATE from ' // FIL16
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        GO TO  999        !  return from kfindx()

    END IF          !  if ierr bad or if timestep flags bad

    IERR = NF_GET_VARA_INT( FNUM, XINDX3( FID ), DIMT, DELT, STIMES )
    IF ( IERR .NE. 0 ) THEN     !  timestep flag not yet written

        WRITE( MESG,91010 ) 'Requested col:  ', COL, 'Requested row:  ', ROW
        CALL M3MESG( MESG )
        WRITE( MESG,91010 ) 'netCDF error number', IERR, 'netCDF ID:  ', FNUM
        CALL M3MESG( MESG )
        MESG = 'Error reading KFSTIME from ' // FIL16
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        GO TO  999        !  return from kfindx()

    END IF          !  if ierr bad or if timestep flags bad

    IERR = NF_GET_VARA_INT( FNUM, YINDX3( FID ), DIMT, DELT, KFLENS )
    IF ( IERR .NE. 0 ) THEN     !  timestep flag not yet written

        WRITE( MESG,91010 ) 'Requested col:  ', COL, 'Requested row:  ', ROW
        CALL M3MESG( MESG )
        WRITE( MESG,91010 ) 'netCDF error number', IERR, 'netCDF ID:  ', FNUM
        CALL M3MESG( MESG )
        MESG = 'Error reading KFLNGTH from ' // FIL16
        CALL M3WARN( 'KFINDX', 0,0, MESG )
        KFFLAG = .FALSE.
        GO TO  999        !  return from kfindx()

    END IF          !  if ierr bad or if timestep flags bad

    KFFLAG = .TRUE.

999     CONTINUE        !  target of "exit from routine"

!$OMP   END CRITICAL( S_NC )

    RETURN

    !!******************  FORMAT  STATEMENTS   ******************************

    !!...........   Error and warning message formats..... 91xxx

91010   FORMAT ( 5 ( A , :, I5, :, 2X ) )

END FUNCTION KFINDX

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
Subroutine Cell2node(Fg,Fc,scalar)
use Aero2DCOM
implicit none
integer i,j
real(8) Fg(Ip,Jp),Fc(Ic,Jc)
character(*) scalar
DO j=1,Jc
  DO i=2,Ic
  if(j==1.and.i>=Ib1.and.i<=Ib2+1) then
   if((scalar=='U'.or.scalar=='V').and.Turmod/='inv'.or.scalar=='Tn'.or.scalar=='Tk') then
    Fg(i,j)=0
   else if(scalar=='T') then
    Fg(i,j)=Tf
   else if(scalar=='rho'.and.Proctrl=='incom') then
    Fg(i,j)=rhoi
   else
    Fg(i,j)=0.5*(Fc(i,j)+Fc(i-1,j))
   end if
  else if(j==1) then
   Fg(i,j)=(Fc(i,j)+Fc(i-1,j)+Fc(Ic+1-i,j)+Fc(Ic+2-i,j))/4
  else
   Fg(i,j)=(Fc(i,j)+Fc(i-1,j)+Fc(i,j-1)+Fc(i-1,j-1))/4
  end if
  end DO
end DO
if(scalar=='U') then
  Fg(:,Jp)=Ui
  Fg(1,:)=Ui
  Fg(Ip,:)=Ui
else if(scalar=='V') then
  Fg(:,Jp)=Vi
  Fg(1,:)=Vi
  Fg(Ip,:)=Vi
else if(scalar=='P') then
  Fg(:,Jp)=0
  Fg(1,:)=0
  Fg(Ip,:)=0
else if(scalar=='T') then
  Fg(:,Jp)=Ta
  Fg(1,:)=Ta
  Fg(Ip,:)=Ta
else if(scalar=='rho') then
  Fg(:,Jp)=rhoi
  Fg(1,:)=rhoi
  Fg(Ip,:)=rhoi
else if(scalar=='Tn') then
  Fg(:,Jp)=Tni
  Fg(1,:)=Tni
  Fg(Ip,:)=Tni
else if(scalar=='Tk') then
  Fg(:,Jp)=Tki
  Fg(1,:)=Tki
  Fg(Ip,:)=Tki
else if(scalar=='Te') then
  Fg(:,Jp)=Tei
  Fg(1,:)=Tei
  Fg(Ip,:)=Tei
else if(scalar=='Tw') then
  Fg(:,Jp)=Twi
  Fg(1,:)=Twi
  Fg(Ip,:)=Twi
end if
end Subroutine Cell2node

max_stars_repo_name: sdhzhs/ransfoil

max_issues_repo_name: sdhzhs/ransfoil

max_forks_repo_name: sdhzhs/ransfoil


program threePoint
implicit none
include 'mpif.h'

real, allocatable :: x(:), y(:)
integer :: i, j, nLocal, nGlobal
character(len=32) :: arg
integer :: ierr, myRank, numRanks
integer :: leftNeigh, rghtNeigh
integer :: status(MPI_STATUS_SIZE)

call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, myRank, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, numRanks, ierr)

call GET_COMMAND_ARGUMENT(1,arg)
read(arg,*) nGlobal

! Now setup the _local_ array
nLocal = int(nGlobal/numRanks)
if (myRank == 0) nLocal = nLocal + mod(nGlobal,numRanks)

allocate(x(0:nLocal+1))
allocate(y(nLocal))

! intialize vector
do i=1,nLocal
    x(i) = i + myRank*(nLocal + mod(nGlobal,numRanks))
end do

! now deal with the "ghost" zones
leftNeigh = myRank-1
if (myRank==0) leftNeigh = numRanks-1
rghtNeigh = myRank+1
if (myRank==numRanks-1) rghtNeigh = 0

if (mod(myRank,2)==0) then
    ! exchange with left neighbor rank
    call MPI_SENDRECV(x(1),1,MPI_REAL,leftNeigh,0,x(0),1,MPI_REAL,leftNeigh,0,MPI_COMM_WORLD,status,ierr)
else
    ! exchnage with right neighbor rank
    call MPI_SENDRECV(x(nLocal),1,MPI_REAL,rghtNeigh,0,x(nLocal+1),1,MPI_REAL,rghtNeigh,0,MPI_COMM_WORLD,status,ierr)
end if

! Now do the other set of exchanges
if (mod(myRank,2)==1) then
    ! exchange with left neighbor rank
    call MPI_SENDRECV(x(1),1,MPI_REAL,leftNeigh,1,x(0),1,MPI_REAL,leftNeigh,1,MPI_COMM_WORLD,status,ierr)
else
    ! exchnage with right neighbor rank
    call MPI_SENDRECV(x(nLocal),1,MPI_REAL,rghtNeigh,1,x(nLocal+1),1,MPI_REAL,rghtNeigh,1,MPI_COMM_WORLD,status,ierr)
end if

! perform averaging
do i=1,nLocal
    y(i) = (x(i-1) + x(i) + x(i+1))/3.
end do

print *, myRank, y(1:nLocal)

call MPI_FINALIZE(ierr)

end program threePoint


! written by jxzou at 20210213: replace Cartesian coordinates in a input file
! written by jxzou at 20210217: support Molpro

program main
 implicit none
 integer :: i, itype
 integer, parameter :: iout = 6
 character(len=11) :: str
 character(len=240) :: xyzname, inpname

 i = iargc()
 if(i /= 3) then
  write(iout,'(/,A)')' ERROR in subroutine replace_xyz_in_inp: wrong command&
                     & line arguments!'
  write(iout,'(A)')  ' Example 1 (OpenMolcas): replace_xyz_in_inp a.xyz a.input -molcas'
  write(iout,'(A)')  ' Example 2 (OpenMolcas): replace_xyz_in_inp a.out a.input -molcas'
  write(iout,'(A)')  ' Example 3 (Molpro)    : replace_xyz_in_inp a.xyz a.com -molpro'
  write(iout,'(A,/)')' Example 4 (Molpro)    : replace_xyz_in_inp a.out a.com -molpro'
  stop
 end if

 str = ' '
 call getarg(1, xyzname)
 call getarg(2, inpname)
 call getarg(3, str)
 call require_file_exist(xyzname)
 call require_file_exist(inpname)

 select case(TRIM(str))
 case('-molcas','-openmolcas')
  itype = 1
 case('-molpro')
  itype = 2
 case('-gamess','-gms')
  itype = 3
 case default
  itype = 0
  write(iout,'(A)') 'ERROR in subroutine replace_xyz_in_inp: wrong command&
                   & line argument str='//TRIM(str)
  stop
 end select

 call replace_xyz_in_inp(xyzname, inpname, itype)
 stop
end program main

! replace Cartesian coordinates in a input file
subroutine replace_xyz_in_inp(xyzname, inpname, itype)
 implicit none
 integer :: i, natom
 integer, intent(in) :: itype ! 1/2 for Molcas/Molpro/GAMESS
 integer, parameter :: iout = 6
 real(kind=8), allocatable :: coor(:,:)
 character(len=240), intent(in) :: xyzname, inpname

 i = LEN_TRIM(xyzname)

 if(xyzname(i-3:i) == '.xyz') then
  call read_natom_from_xyz(xyzname, natom)
  allocate(coor(3,natom))
  call read_coor_from_xyz(xyzname, natom, coor)
 else
  select case(itype)
  case(1)
   call read_natom_from_molcas_out(xyzname, natom)
   allocate(coor(3,natom))
   call read_coor_from_molcas_out(xyzname, natom, coor)
  case(2)
   call read_natom_from_molpro_out(xyzname, natom)
   allocate(coor(3,natom))
   call read_coor_from_molpro_out(xyzname, natom, coor)
  case default
   write(iout,'(A)') 'ERROR in subroutine replace_xyz_in_inp: file format not&
                     & supported.'
   write(iout,'(A)') 'Filename='//TRIM(xyzname)
   stop
  end select
 end if

 select case(itype)
 case(1)
  call replace_coor_in_molcas_inp(inpname, natom, coor)
 case(2)
  call replace_coor_in_molpro_inp(inpname, natom, coor)
 case default
  write(iout,'(A)') 'ERROR in subroutine replace_xyz_in_inp: file format not&
                   & supported.'
  write(iout,'(A)') 'Filename='//TRIM(inpname)
  stop
 end select

 deallocate(coor)
 return
end subroutine replace_xyz_in_inp

! replace Cartesian coordinates in a (Open)Molcas input file
subroutine replace_coor_in_molcas_inp(inpname, natom, coor)
 implicit none
 integer :: i, k, fid, fid1
 integer, intent(in) :: natom
 integer, parameter :: iout = 6
 real(kind=8), intent(in) :: coor(3,natom)
 character(len=6) :: elem
 character(len=16) :: key
 character(len=240) :: buf0, buf, inpname1
 character(len=240), intent(in) :: inpname

 i = index(inpname, '.inp', back=.true.)
 inpname1 = inpname(1:i-1)//'_new.input'
 open(newunit=fid,file=TRIM(inpname),status='old',position='rewind')
 open(newunit=fid1,file=TRIM(inpname1),status='replace')

 i = 0
 do while(.true.)
  read(fid,'(A)',iostat=k) buf
  if(k /= 0) exit
  write(fid1,'(A)') TRIM(buf)

  key = buf(1:16)
  call lower(key)
  if(key == 'end of basis set') then
   read(buf0,*) elem
   BACKSPACE(fid1)
   BACKSPACE(fid1)
   i = i + 1
   write(fid1,'(A,3(2X,F15.8),A)') TRIM(elem),coor(1:3,i),'   Angstrom'
   write(fid1,'(A)') 'End of basis set'
   if(i == natom) exit
  end if

  buf0 = buf
 end do ! for while

 if(k /= 0) then
  write(iout,'(A)') 'ERROR in subroutine replace_coor_in_molcas_inp: insuffic&
                    &ient number of atoms found in file '//TRIM(inpname)
  close(fid1,status='delete')
  close(fid)
  stop
 end if

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  write(fid1,'(A)') TRIM(buf)
 end do ! for while

 close(fid)
 close(fid1)
 return
end subroutine replace_coor_in_molcas_inp

! replace Cartesian coordinates in a Molpro input file
subroutine replace_coor_in_molpro_inp(inpname, natom, coor)
 implicit none
 integer :: i, fid, fid1
 integer, intent(in) :: natom
 integer, parameter :: iout = 6
 real(kind=8), intent(in) :: coor(3,natom)
 character(len=6) :: elem
 character(len=8) :: key
 character(len=240) :: buf, inpname1
 character(len=240), intent(in) :: inpname

 i = index(inpname, '.com', back=.true.)
 inpname1 = inpname(1:i-1)//'_new.com'
 open(newunit=fid,file=TRIM(inpname),status='old',position='rewind')
 open(newunit=fid1,file=TRIM(inpname1),status='replace')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  write(fid1,'(A)') TRIM(buf)

  key = buf(1:8)
  call lower(key)
  if(key == 'geometry') exit
 end do ! for while

 if(i /= 0) then
  write(iout,'(A)') 'ERROR in subroutine replace_coor_in_molcas_inp: insuffic&
                    &ient number of atoms found'
  write(iout,'(A)') 'in file '//TRIM(inpname)
  close(fid1,status='delete')
  close(fid)
  stop
 end if

 do i = 1, natom, 1
  elem = ' '
  read(fid,*) elem
  write(fid1,'(A,3(2X,F15.8))') TRIM(elem),coor(1:3,i)
 end do ! for i

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  write(fid1,'(A)') TRIM(buf)
 end do ! for while

 close(fid)
 close(fid1)
 return
end subroutine replace_coor_in_molpro_inp


text:
	SUBROUTINE IN_BDTA  ( iret )
C************************************************************************
C* IN_BDTA								*
C*									*
C* This subroutine serves as a BLOCKDATA statement, initializing	*
C* variables in GEMPAK common blocks.  This subroutine is called	*
C* by IP_INIT.  If a GEMPAK program does not call IP_INIT, it must	*
C* call IN_BDTA directly.						*
C*									*
C* IN_BDTA  ( IRET )							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 2/87						*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* M. desJardins/GSFC	 7/90	Initialize LUNGEM for lun flags		*
C* M. desJardins/NMC	 1/92	Initiaize mbchan			*
C* K. Brill/NMC		01/92	Remove CALL IN_BOLD			*
C* M. desJardins/NMC	 8/94	Initialize IESTAT			*
C* C. Lin/EAI		 8/95	Added IMFTYP				*
C* S. Jacobs/NCEP	 7/96	Added IFINVD				*
C* K. Tyle/GSC		12/96	Replace IESTAT with call to ER_STAT;	*
C*				added ER_INIT call			*
C* R. Tian/SAIC		 1/04	Added GD_INIT call			*
C* T. Piper/SAIC	12/07	Removed GD_INIT				*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE		'ADBUFF.CMN'
	INCLUDE		'IMGDEF.CMN'
C------------------------------------------------------------------------
	iret= 0
C
C*	Initialize IM variables.
C
	imftyp = IFINVD
C
C*	Initialize DP variables.
C
	DO  i = 1, MMFLDP
	    iflgdp (i) = 0
	END DO
C
C*	Initialize IP variables.
C
	lunp = 0
	lunf = 0
C
C*	Initialize PC variables.
C
	tablrd = .false.
	dsflg  = .false.
	stnflg = .false.
	pptabl = .false.
C
C*	Initialize SF, SN, GD DM file numbers and LUNs.
C
	DO  i = 1, MMFILE
	    isfcfn (i) = -1
	    isndfn (i) = -1
	    igrdfn (i) = -1
	    lundm (i)  = -1
	END DO
C
C*	Initialize DM variables
C
	kclast = 0
	DO i = 1, MCACHE
	    kwrite (i) = .false.
	    kcflno (i) = 0
	    kcrecn (i) = 0
	END DO
C
C*	Initialize logical unit number flags.
C
	DO  i = 1, 10
	    lungem (i) = 0
	END DO
C
C*	Initialize channel number for GPLT mailbox, so that an error
C*	can be generated for missing call to GINITP.
C
	mbchan = -999
C
C*	Initialize ERCMN.CMN and ERCMN.H files.
C
	CALL ER_STAT ( 0, 0, .false., iret )
	CALL ER_INIT ( iret )
C*
 	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


module rzvv_grid
    use constants

    implicit none

    save

    real :: R_range, z_Range, vPer_range, vPar_range
    real :: R_binSize, z_binSize, vPer_binSize, vPar_binSize
    real, allocatable :: R_binEdges (:), &
            z_binEdges (:), &
            vPer_binEdges (:), &
            vPar_binEdges (: )
    real, allocatable :: R_binCenters (:), &
            z_binCenters (:), &
            vPer_binCenters (:), &
            vPar_binCenters (: )
    real :: z_min, z_max, R_min, R_max

    real, allocatable :: f_rzvv (:,:,:,:), f_rzvv_global(:,:,:,:)
    real, allocatable :: density (:,:)
    real :: normFac_
    real(kind=dbl) :: v_sigma, x_sigma

contains
    subroutine init_rzvv_grid ()
        use read_particle_list
        use eqdsk
        use constants
        use read_namelist
        use init_mpi
        implicit none

        integer :: i
        i = 0

        allocate ( R_binEdges ( R_nBins + 1 ), &
                   z_binEdges ( z_nBins + 1 ), &
                   vPer_binEdges ( vPer_nBins + 1 ), &
                   vPar_binEdges ( vPar_nBins + 1 ), &
                   R_binCenters ( R_nBins ), &
                   z_binCenters ( z_nBins ), &
                   vPer_binCenters ( vPer_nBins ), &
                   vPar_binCenters ( vPar_nBins ) )
        allocate ( f_rzvv ( R_nBins, z_nBins, vPer_nBins, vPar_nBins ) )
        allocate ( density ( R_nBins, z_nBins ) )

        R_min   = minVal ( rbbbs )
        R_max   = maxVal ( rbbbs )
        R_range = R_max - R_min
        R_binSize   = R_range / R_nBins
        R_binEdges  = (/ (i*R_binSize,i=0,R_nBins) /) + R_min
        R_binCenters    = R_binEdges(1:R_nBins) + R_binSize / 2.0

        z_min   = minVal ( zbbbs )
        z_max = maxVal ( zbbbs )
        z_range = ( z_max - z_min )
        z_binSize   = z_range / z_nBins
        z_binEdges  = (/ (i*z_binSize,i=0,z_nBins) /) + z_min
        z_binCenters    = z_binEdges(1:z_nBins) + z_binSize / 2.0

        if(mpi_pId==1) write(*,*) 'R_min: ', R_min
        if(mpi_pId==1) write(*,*) 'R_max: ', R_max
        if(mpi_pId==1) write(*,*) 'z_min: ', z_min
        if(mpi_pId==1) write(*,*) 'z_max: ', z_max

        vPer_range = vPerInRange / 100.0 * c
        vPer_binSize   = vPer_range / vPer_nBins
        vPer_binEdges   = (/ (i*vPer_binSize,i=0,vPer_nBins) /)
        vPer_binCenters    = vPer_binEdges(1:vPer_nBins) &
            + vPer_binSize / 2.0

        vPar_range  = vParInRange / 100.0 * c
        vPar_binSize    = 2.0 * vPar_range / vPar_nBins
        vPar_binEdges   = (/ (i*vPar_binSize,i=0,vPar_nBins) /) - vPar_range
        vPar_binCenters = vPar_binEdges(1:vPar_nBins) + vPar_binSize / 2.0

        f_rzvv  = 0.0

    end subroutine init_rzvv_grid

end module rzvv_grid


text:
Hoffmann Automotive is a local automobile repair shop in business since 1979. Owned by the Hoffmann family, they have a good reputation for customer service. Service is a high priority at Hoffmann Automotive.  Customers who need a ride to or from home/work will be accommodated by most any available technician, or by Jeff Hoffmann himself.

Hoffmann isnt one of those cheapo budget places, but you definitely get what you pay for.  The service is fantastic and they stand behind their work.  The parts they set you up with generally come with long warranties and easy replacement (on the off chance its needed).  Its a pretty popular place  call ahead to make an appointment, since they can get booked up several days in advance.

While they do have a waiting area in the shop, they offer you a ride home when you drop off your vehicle, and pick you up when the service is done.  Its worth taking them up on the offerit typically takes a few hours for, e.g., oil change + headlight replacement.

History

They used to be operated out of the Hoffmann 76 Gas Stations gas station. Both facilities were owned by the Hoffmann family and garnered a good reputation for customer service. Hoffmann family owned Davis Muffler & Hitch Center, which now shares (shared?) the 2nd Street Second Street location.

Not to be confused with the former Hoffmans Garage.

   Having car problems?  There are many automobile repair automotive repair services in town!

Also see /Hoffmann additional comments Hoffmann additional comments

20100216 18:27:59 nbsp Im taking my vehicle in tomorrow, which will be the first time I have used Hoffmann for service.  I have a failed window regulator, with my window stuck in the down position.  They were polite on the phone, and asked me if my car is garaged, which it is.  I am assuming they were going to let me bring it in tonight instead of tomorrow morning so it would be garaged.  Anyway, I looked up the entry in the wiki to get the phone number, and I read through all of the comments.  I have two interesting observations (make your own inferences):

1)  Every single negative comment on this page (except Truman) was made by a person that created an account, left a negative comment, and then disappeared from the Wiki.

2)  Trumans contribution to the wiki comprises two comments on this page, and one edit to correct the pricing of a competitor... Users/DavidGrundler

  Just picked up my car.  Diagnosis was proper, and fix was fast and well within the expected price range.  They were polite and knowledgeable.  I will definitely use them again.  They also conducted a complimentary 30point inspection.   Users/DavidGrundler

20100309 13:26:41 nbsp I used to work across the street from Rick and Jeff when they owned the G & 5th street property and can tell you that they were never underhanded or shady, as some of the ghost posters here have claimed. I have seen these people work their hindquarters off trying to help even the most rude of people and have always seen them to be friendly, honest, fair and professional in every way. They have been here for over 30 years and do very little, if any, advertising which means they have been doing word of mouth business for over 30 years. You tel me what kind of business can do that if they are shady cheating or any of the other slanderous statements I have seen here.

I recommend this business and these people VERY HIGHLY.

Dont be fooled by the ghost poster that will pop n here, make a derogatory statement and never be seen or heard from again or post directly on a competitive businesses site.  Users/WesOne

20100412 15:13:30 nbsp I went in to Hoffman Automotive today for the first time.  Needed new tires, and one of my headlights went out last week, so that was on the list.  And it was just about time for an oil change, so... hey, why not.  When I called this morning to check on availability and prices of tires, I spoke with Jeff Hoffman, the owner.  He got my info and said hed do a little research and call me back shortly.  It took less than 10 minutes, and I made an appointment for shortly thereafter.  Went on in and Jeff was super friendly.  Got me set up with an estimate (which was about $10 off and didnt include tax... close enough, although certainly not exact), and then drove me home personally.  A few hours later I got a call from Jason there, who came on out to get me.  If you have a chance, ask Jason about New Orleans right after Katrinahes got crazy stories!

The prices arent cheap, but they also arent outrageous.  Extremely friendly, and as far as I can tell so far, the quality of the work is good.  The headlight seems to be nicely aligned and all that.  Ill definitely be coming back regularly. Users/TomGarberson

20100412 21:50:11 nbsp As a UCD Student years ago, I used do my own car repair.  For the past twentyfive years, I’ve relinquished auto maintenance to Rick and Jeff Hoffman, who’ve serviced many cars for me.  They’re knowledgeable, professional, practical and run a firstclass shop.   Users/pschrupp

20100707 10:23:13 nbsp I was working on my own brakes at 11:00 at night on June 30 and thought that I could change them out.  My buddy omar and I had replaced front disk breaks many times before so we thought that the rears would be the same.  Afer an hour into it and lots of reading repair sites on the interet we found that there were some special tools needed to push in the calipers.  The bad thing was that We were leaving for camping the next morning for the fourth of July weekend.  In the morning I took the car to Hoffmann Automotive and they expitided the repairs with in a few hours for me. We were able to leave by noon.  I want to thank Jeff and Rick for making this happen.  The service was excellent.  The bill was hefty because I asked them to put in OEM parts.  The labor charge was totally resonable. Users/Khanh

20100802 16:10:54 nbsp Just went to Hoffman for a tire repair.  Cheap, friendly, pretty quick (I didnt make an appointment).  Ill be back! Users/DanLawyer

20101117 09:05:05 nbsp Took my car here yesterday.  It was great, timely, informative service.  The quality of work and explanations made me very happy.  I would recommend this place to anyone. Over all excellent Users/KatilinAdams

20110511 11:27:09 nbsp I just had an very, very positive experience with Hoffman Automotive. A piece of my exhaust fell off my car on the freeway  that was exciting. They diagnosed my problem and did enough of a temporary fix that I can get home to the East Bay, where Ill do my own repairs at home. They saved me from a lot of trouble and expense  they didnt charge me at all! I hope I dont break down in Davis again, but if I do, Ill definitely take my car back to Hoffman. They were very polite, friendly, and kind. Users/KatieHarris

20110613 15:28:00 nbsp Ive now been here several times for various services.  In addition to a full set of new tires, which I mentioned in an earlier review, theyve done my 45,000 mile service which included flushing the power steering and transmission lines.  I got a new battery, another new headlight (the other one... Im not that hard on my headlights!) and some oil changes there.  100% of the service has been fantastic, and everyone Ive talked to there has been extremely friendly.  This last time my wife picked up the car for me while I was at work, and she was impressed with the folks she talked to.  In the past, at other service shops/dealerships, shes felt like she was talked down to because she was a woman, but there was no sense of that here.

It looks like there were a few negative comments about this place years ago from single edit users accusing it of shady business practices.  Its up to you to make your own judgment call, but personally, I suspect most of them are BS.  Every business can make mistakes, so I wont go so far as to say theyre all bogus, but Hoffmanns is an extremely high quality repair shop that, from my experience, is completely trustworthy.  If you want high quality work, high quality parts, and friendly service, give them a visit. Users/TomGarberson

20110808 13:48:10 nbsp Taking my car to Hoffmann Automotive was a relief. Im a fairly new car owner as well as a young woman, and what I fear from unfamiliar repair shops is having my ignorance of fair cost and car mechanics exploited. I dont believe it was, here  the first visit they performed a thorough inspection of my car, blew the dust from my brake drums in an attempt to silence a disturbing scraping sound, kept me informed during the process, and finally assured me that the brakes were safe. Following the inspection the owner reviewed the results with me, informed me of upcoming maintenance checkpoints for my vehicle, and charged for the hour+ labor what I have gathered from my history with less reliable cars is reasonable for a AAAbacked auto shop.

On my second visit, just today, I brought to them an oil leak I thought was severe. Despite being busy they got my car up, examined, and out of there in about an hour, and, especially because Im on crutches right now, I was thrilled with the effort made to square me away quickly. When explaining to me what was done (bolttightening, removal of the excess oil I added) and what I should do/expect, the owner has never insulted my intelligence, and I really value that. After all this, there was no charge for the hour they spent on my vehicle while I got to study and catch endearing little exchanges between the family members.

The business has shown great consideration to me during my tiny crises, and it has secured a return a customer. Users/Willikers

20110810 12:09:36 nbsp BAD BAD BAD EXPERIENCE! AND VERY EXPENSIVE PLACE!!!!

1) Just to ckech why the engine light was on they charge me: 150$!!!! Just to check!!! All the other mechanics do it for free!!! (operation that takes less then 2 minutes!)

2) In order to fix the car they ask me 650$

3) I run away from this place and I went to another mechanic (he check the engine light for free) and he discover that was a very simple work, with a cost of 5$20$ to fix it!!!

..............NO COMMENT...........

20120219 12:56:52 nbsp I took my car in to get a new starter. The staff was friendly and professional (Jason in particular made me feel like a valued customer, even though it was my first visit), and they kept in touch with me throughout the replacement. Service was quick and competitively priced. If my car needs anything in the future, Im definitely calling these guys back! Users/BombayCollins

20120317: I have taken my car here three times in the past year, and I have always been very happy with the results. They are friendly, knowledgeable, and professional. The first time I went there I asked them to do an overall check of my car and they came back with a list of things that could be done, from an oil change to a new fan belt. I told them I only had a certain amount of money to spend and asked them to recommend what was most important. They did so without any hesitation and without trying to convince me to spend more. The next time I came in, I wanted to buy windshield wipers from them and install them myself, but when they heard that I was going to do it myself, the guy behind the counter ran out and did it for me, for free! They have proven to be honest and to stand behind their work, and I am very glad that I have found an auto repair company in town that I can trust.Elizabeth C.

20120318 19:26:49 nbsp I always have a great experience bringing my car here. They are very helpful and always answer my questions. I just took my car here last week to get if fixed after another automotive place was unable to fix a noise it was making. They fixed it, as well as fixed a few other things, all in one day.  Users/KMcGuire

20121016 17:11:52 nbsp I have been using Hoffmanns to service my cars for more than 20 years. They are trustworthy, funny, and will go out of their way to help someone in need. I would highly recommend taking your car to them if it needs service.  Users/PeriSiepman

20121121 18:11:03 nbsp Its Thanksgiving eve and I will be heading home tomorrow, thanks to Hoffmann Automotive. I am a firsttime customer who was discouraged with the service and honesty of some other shops in Davis. After a free 30point inspection, their estimate for work was lower and more thorough than anything Id seen so far.  Moreover, in addition to getting me back on the road before the holiday, when there was a mix up about when to do the additional services, they did not charge me for them at all. My previous experience led me to expect, when authorizing service over the phone, I would be forced to pay (my word against theirs). Instead, Hoffmann Automotive treated me with integrity and respect, even though it cost them money. I am thankful to have found an auto shop in Davis that I feel good recommending to my friends.  Users/rlreeves

20121227 16:21:29 nbsp I feel weird leaving a third review for this place, but I had yet another good experience.  I recently took my car in for a substantial service and they accidentally doublebooked the courtesy car they loan out.  They apologized and immediately resolved the situation by giving me a ride over to Enterprise, where they covered a rental for me for the day.  The inconvenience was minimal, and they gave me 10% off the labor as an apology.

Friendly, knowledgeable folks with a lot of integrity.  Highly recommended. Users/TomGarberson

20130125 18:14:49 nbsp The Hoffman brothers have always been kind and considerate as well as professional and very experienced.  Ive known them for over 30 years and have NEVER heard anyone who has had a bad experience there.  Im so glad they are here to care for our cars. Users/LoisRichter

20130404 18:01:33 nbsp $191 to change ONE HID LIGHT BULB. ouch.

Dear college students, save yourself some money and stay away from these heartless people. Go somewhere else where they wont rip you off. My biggest mistake was trusting their 5 star reviews on yelp and not ask anyone else about the price I had to pay for beforehand. OH DID I FORGET TO MENTION HOW RUDE THEY WERE AFTER my friends busted them for charging $150 for a part that should cost less than half of that price?
 Users/honeybun

        Didnt they give you an estimate before the work was performed that showed you the proposed charge?  they should have I believe and you could have certainly decided at that point to take your car elsewhere.
           users/Loneshark
        This sounds like BS.  Ive had two headlights replaced there. Both times were with other service  oil changes, one with a larger maintenance. Both times they charged for the bulb only ($10ish I think) and didnt charge a dime for the labor. TomGarberson
        The cost of just the HID bulb for my car runs $160.00$220.00 at cheaper aftermarket places.  Much higher through dealers.  HID bulbs can be very expensive, especially the xenon variety.   Users/DavidGrundler

20130405 10:04:09 nbsp Honey Bun you are an idiot. Do you accuse 5 star restaurants of being a rip off for charging $50 for a plate of food that  should cost less than half of that price? Everything you pay for is marked up from Advil at the hospital to the cherries the guy at the fruit stand sells. If you didnt like the cost of parts and labor that BY LAW they are required to provide you with BEFORE they do any work, you should have gone someplace else.  Users/grumpyoldcodger

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search