text:
module diffusive_mass_fluxdiv_module

  use multifab_module
  use define_bc_module
  use bc_module
  use div_and_grad_module
  use probin_multispecies_module, only: nspecies, is_nonisothermal, &
                                        nspecies, correct_flux
  use probin_common_module, only: barodiffusion_type
  use mass_flux_utilities_module
  use ml_layout_module
  use convert_stag_module
  use matvec_mul_module
  use matmat_mul_module
  use correction_flux_module
  use zero_edgeval_module

  implicit none

  private

  public :: diffusive_mass_flux, diffusive_mass_fluxdiv

contains

  subroutine diffusive_mass_fluxdiv(mla,rho,rhotot,molarconc,rhoWchi,Gama,&
                                    diff_fluxdiv,Temp,zeta_by_Temp,gradp_baro, &
                                    flux_total,dx,the_bc_tower)

    type(ml_layout), intent(in   )  :: mla
    type(multifab) , intent(in   )  :: rho(:)
    type(multifab) , intent(in   )  :: rhotot(:)
    type(multifab) , intent(in   )  :: molarconc(:)
    type(multifab) , intent(in   )  :: rhoWchi(:)
    type(multifab) , intent(in   )  :: Gama(:)
    type(multifab) , intent(inout)  :: diff_fluxdiv(:)
    type(multifab) , intent(in   )  :: Temp(:)
    type(multifab) , intent(in   )  :: zeta_by_Temp(:)
    type(multifab) , intent(in   )  :: gradp_baro(:,:)
    type(multifab) , intent(inout)  :: flux_total(:,:)
    real(kind=dp_t), intent(in   )  :: dx(:,:)
    type(bc_tower) , intent(in   )  :: the_bc_tower

    ! local variables
    integer i,dm,n,nlevs

    ! local array of multifabs for grad and div; one for each direction
    type(multifab) :: flux(mla%nlevel,mla%dim)

    type(bl_prof_timer), save :: bpt

    call build(bpt, "diffusive_mass_fluxdiv")

    nlevs = mla%nlevel  ! number of levels
    dm    = mla%dim     ! dimensionality

    ! build the local multifabs
    do n=1,nlevs
       do i=1,dm
          ! flux(i) is face-centered, has nspecies component, zero ghost
          ! cells & nodal in direction i
          call multifab_build_edge(flux(n,i),mla%la(n),nspecies,0,i)
       end do
    end do

    ! compute the face-centered flux (each direction: cells+1 faces while
    ! cells contain interior+2 ghost cells)
    call diffusive_mass_flux(mla,rho,rhotot,molarconc,rhoWchi,Gama,Temp,&
                             zeta_by_Temp,gradp_baro,flux,dx,the_bc_tower)

    ! add fluxes to flux_total
    do n=1,nlevs
       do i=1,dm
          call multifab_plus_plus_c(flux_total(n,i),1,flux(n,i),1,nspecies,0)
       end do
    end do

    ! compute divergence of determinstic flux
    call compute_div(mla,flux,diff_fluxdiv,dx,1,1,nspecies)

    ! destroy the multifab to free the memory
    do n=1,nlevs
       do i=1,dm
          call multifab_destroy(flux(n,i))
       end do
    end do

    call destroy(bpt)

  end subroutine diffusive_mass_fluxdiv

  subroutine diffusive_mass_flux(mla,rho,rhotot,molarconc,rhoWchi,Gama, &
                                 Temp,zeta_by_Temp,gradp_baro,flux,dx,the_bc_tower)

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(in   ) :: rho(:)
    type(multifab) , intent(in   ) :: rhotot(:)
    type(multifab) , intent(in   ) :: molarconc(:)
    type(multifab) , intent(in   ) :: rhoWchi(:)
    type(multifab) , intent(in   ) :: Gama(:)
    type(multifab) , intent(in   ) :: Temp(:)
    type(multifab) , intent(in   ) :: zeta_by_Temp(:)
    type(multifab) , intent(in   ) :: gradp_baro(:,:)
    type(multifab) , intent(inout) :: flux(:,:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    type(bc_tower) , intent(in   ) :: the_bc_tower

    ! local variables
    integer :: n,i,s,dm,nlevs

    ! local face-centered multifabs
    type(multifab)  :: rhoWchi_face(mla%nlevel,mla%dim)
    type(multifab)  :: Gama_face(mla%nlevel,mla%dim)
    type(multifab)  :: zeta_by_Temp_face(mla%nlevel,mla%dim)
    type(multifab)  :: flux_Temp(mla%nlevel,mla%dim)
    type(multifab)  :: baro_coef(mla%nlevel)
    type(multifab)  :: baro_coef_face(mla%nlevel,mla%dim)

    type(bl_prof_timer), save :: bpt

    call build(bpt,"diffusive_mass_flux")

    dm    = mla%dim     ! dimensionality
    nlevs = mla%nlevel  ! number of levels

    ! build local face-centered multifab with nspecies^2 component, zero ghost cells
    ! and nodal in direction i
    do n=1,nlevs
       call multifab_build(baro_coef(n),mla%la(n),nspecies,rho(n)%ng)
       do i=1,dm
          call multifab_build_edge(rhoWchi_face(n,i),     mla%la(n),nspecies**2,0,i)
          call multifab_build_edge(Gama_face(n,i),        mla%la(n),nspecies**2,0,i)
          call multifab_build_edge(zeta_by_Temp_face(n,i),mla%la(n),nspecies,   0,i)
          call multifab_build_edge(baro_coef_face(n,i),mla%la(n),nspecies,   0,i)
          call multifab_build_edge(flux_Temp(n,i),        mla%la(n),1,          0,i)
       end do
    end do

    ! compute face-centered rhoWchi from cell-centered values
    call average_cc_to_face(nlevs, rhoWchi, rhoWchi_face, 1, tran_bc_comp, &
                            nspecies**2, the_bc_tower%bc_tower_array, .false.)

    !==================================!
    ! compute flux-piece from molarconc
    !==================================!

    ! calculate face-centrered grad(molarconc)
    call compute_grad(mla, molarconc, flux, dx, 1, mol_frac_bc_comp, 1, nspecies, &
                      the_bc_tower%bc_tower_array)

    ! compute face-centered Gama from cell-centered values
    call average_cc_to_face(nlevs, Gama, Gama_face, 1, tran_bc_comp, &
                            nspecies**2, the_bc_tower%bc_tower_array, .false.)

    ! compute Gama*grad(molarconc): Gama is nspecies^2 matrix; grad(x) is nspecies component vector
    do n=1,nlevs
       do i=1,dm
          call matvec_mul(mla, flux(n,i), Gama_face(n,i), nspecies)
       end do
    end do

    if(is_nonisothermal) then

       !====================================!
       ! compute flux-piece from Temperature
       !====================================!

       ! calculate face-centrered grad(T)
       call compute_grad(mla, Temp, flux_Temp, dx, 1, temp_bc_comp, 1, 1, the_bc_tower%bc_tower_array)

       ! compute face-centered zeta_by_T from cell-centered values
       call average_cc_to_face(nlevs, zeta_by_Temp, zeta_by_Temp_face, 1, tran_bc_comp, &
                               nspecies, the_bc_tower%bc_tower_array, .false.)

       ! compute zeta_by_T*grad(T): zeta_by_T is nspecies component vector; grad(T) is scalar
       do n=1,nlevs
          do i=1,dm
             do s=1,nspecies
                call multifab_mult_mult_c(zeta_by_Temp_face(n,i), s, flux_Temp(n,i), 1, 1)
             end do
          end do
       end do

       !===============================!
       ! assemble different flux-pieces
       !===============================!
       do n=1,nlevs
          do i=1,dm
             call multifab_plus_plus(flux(n,i), zeta_by_Temp_face(n,i), 0)
          end do
       end do

    end if

    if (barodiffusion_type .gt. 0) then

       !====================================!
       ! compute flux-piece from barodiffusion
       !====================================!

       ! compute cell-centered barodiffusion coefficient, (phi-w) / (n kB T)
       call compute_baro_coef(mla,baro_coef,rho,rhotot,Temp)

       ! average baro_coef to faces
       call average_cc_to_face(nlevs, baro_coef, baro_coef_face, 1, scal_bc_comp, &
                               nspecies, the_bc_tower%bc_tower_array, .false.)

       ! store the fluxes, baro_coef(1:nspecies) * gradp_baro, in baro_coef_face
       do n=1,nlevs
          do i=1,dm
             do s=1,nspecies
                call multifab_mult_mult_c(baro_coef_face(n,i), s, gradp_baro(n,i), 1, 1)
             end do
          end do
       end do

       !===============================!
       ! assemble different flux-pieces
       !===============================!
       do n=1,nlevs
          do i=1,dm
             call multifab_plus_plus(flux(n,i), baro_coef_face(n,i), 0)
          end do
       end do

    end if

    ! compute rhoWchi * totalflux (on faces)
    do n=1,nlevs
       do i=1,dm
          call matvec_mul(mla, flux(n,i), rhoWchi_face(n,i), nspecies)
       end do
    end do

    ! If there are walls with zero-flux boundary conditions
    if(is_nonisothermal) then
       do n=1,nlevs
          call zero_edgeval_walls(flux(n,:),1,nspecies,the_bc_tower%bc_tower_array(n))
       end do
    end if

    !correct fluxes to ensure mass conservation to roundoff
    if (correct_flux .and. (nspecies .gt. 1)) then
       !write(*,*) "Checking conservation of deterministic fluxes"
       call correction_flux(mla, rho, rhotot, flux, the_bc_tower%bc_tower_array)
    end if

    ! destroy B^(-1)*Gama multifab to prevent leakage in memory
    do n=1,nlevs
       call multifab_destroy(baro_coef(n))
       do i=1,dm
          call multifab_destroy(rhoWchi_face(n,i))
          call multifab_destroy(Gama_face(n,i))
          call multifab_destroy(zeta_by_Temp_face(n,i))
          call multifab_destroy(flux_Temp(n,i))
          call multifab_destroy(baro_coef_face(n,i))
       end do
    end do

    call destroy(bpt)

  end subroutine diffusive_mass_flux

end module diffusive_mass_fluxdiv_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
! { dg-do compile }
! { dg-options "-std=f95" }
! Check for notify-std-messages when F2003 structure constructors are compiled
! with -std=f95.

PROGRAM test
  IMPLICIT NONE

  ! Basic type with default initializers
  TYPE :: basics_t
    INTEGER :: i = 42
    REAL :: r = 1.5
  END TYPE basics_t

  TYPE(basics_t) :: basics

  ! This is ok in F95
  basics = basics_t (1, 2.)

  ! No argument naming in F95
  basics = basics_t (1, r = 4.2) ! { dg-error "Fortran 2003" }

  ! No optional arguments in F95
  basics = basics_t () ! { dg-error "Fortran 2003" }
  basics = basics_t (5) ! { dg-error "Fortran 2003" }

END PROGRAM test

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
EuropeanAmericans, many from Spain, settled what is now California and the Davis area in the early 1800s.  They brought a wide variety of various regional heritages from across the continent of Europe.  Racially they range from the dark skinned people of the Mediterranean skirting the Middle East to the fair haired and skinned Scandinavians.

If you are new to the area from Europe, you might be interested in the services at the International House.  Various European cuisines of varying authenticity are found at some of the restaurants in town.  Many more European restaurants are located in nearby Sacramento or Winters, such as the Fox & Goose Pub, Streets of London Pub and the Irish Pub & Coffee House.

 Culture

If you are interested in the history of Western Civilization, you might be interested in the Society for Creative Anachronism, which recreates the middle ages from a eurocentric (although not exclusive) perspective.  Various Renaissance Faires exist throughout the Northern California area, some specifically oriented toward Italian or Celtic heritages.

Traditional European folk dance is well represented in Davis with the English Country Dancers and Merrie Pryanksters.  European music inspires some people in the Davis music scene, including our own Davis High School Madrigal Singers.

 History

Spanish immigrants were the first Europeans to settle the area, founding Mission Dolores.  From that location, J.M. Vaca and J.F. Peña founded Rancho Los Puttos in 1842. This massive ranch stretched from Vacaville to the edges of South Davis.  Eleven years later, Joseph B. Chiles established Rancho Laguna de Santos Calle by purchasing part of Vacas land. One of the subdivided areas of Chiles land was Jerome C. Davis stock ranch which is now Downtown Davis.

 UC Davis

  Anybody know what resources for European students doing study/in the grad programs at UCD are?

 Demographics

People identified primarily as EuropeanAmericans comprise about twothirds of of Davis Demographics population, a clear majority, on par with the rest of the United States.   UC Davis underrepresents EuropeanAmericans by a larger percentage than it does AfricanAmericans, however the proportion of underrepresentation is much smaller.  The number of EuropeanAmericans identifying students at UC Davis is almost the same as the amount of Asians Asian identifying students.

  See Demographics

 Groups/Resources

  See Ethnic and Cultural Organizations

20060118 09:54:52 nbsp It was strange at first that there was an EuropeanAmericans page. I think thats because I perceive whites as not fully acknowledging their European heritage. I mean, they do, but theres not need for white people to say that theyre European since being an American in the international world is all about being white. Users/RitchieLee

20060118 12:19:10 nbsp What you say makes sense.  Being a white person of mixed race it there is no way to identify completely with a particular culture.  I myself have taken to the German heritage.  Having a schmorgousborge of options I suppose I really could have been any European I wanted to be.  However, I do find it comforting to have some cultural ideal to cling to as my own.  It may seem silly to some; but it gives one some form of identification.  Dare I say something to be proud of?  I know that can be a bit taboo, but it seems a great way to summarize the whole feeling. Users/AndrewBanta

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine read_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
c********************************************************************
c     Reads one event from data file #lun
c     ic(*,1) = Particle ID
c     ic(*,2) = Mothup(1)
c     ic(*,3) = Mothup(2)
c     ic(*,4) = ICOLUP(1)
c     ic(*,5) = ICOLUP(2)
c     ic(*,6) = ISTUP   -1=initial state +1=final  +2=decayed
c     ic(*,7) = Helicity
c********************************************************************
      implicit none

      double precision pi
      parameter (pi = 3.1415926d0)
c
c     Arguments
c
      integer lun
      integer nexternal, ic(7,*)
      logical done
      double precision P(0:4,*),wgt,aqcd,aqed,scale
      integer ievent
      character*140 buff2
c
c     Local
c
      integer i,j,k
      character*(256) buff
      double precision xdum1,xdum2
c
c     Global
c
      logical banner_open
      integer lun_ban
      common/to_banner/banner_open, lun_ban

      data lun_ban/37/
      data banner_open/.false./
c-----
c  Begin Code
c-----
      buff2=' '
      done=.false.
      if (.not. banner_open) then
         open (unit=lun_ban, status='scratch')
         banner_open=.true.
      endif
 11   read(lun,'(a256)',end=99,err=99) buff
      do while(index(buff,"<event") .eq. 0)
         write(lun_ban,'(a)') buff
         read(lun,'(a256)',end=99,err=99) buff
      enddo
      read(lun,*,err=11, end=11) nexternal,ievent,wgt,scale,aqed,aqcd
      do i=1,nexternal
         read(lun,*,err=99,end=99) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
     $     (p(j,i),j=1,3),p(0,i),p(4,i),xdum1,xdum2
         ic(7,i)=xdum2
      enddo
      do while(index(buff,"</event") .eq. 0)
         read(lun,'(a256)',end=99,err=99) buff
         if(buff(1:1).eq.'#') buff2=buff(1:140)
      enddo
c      gal(1) = sqrt(4d0*pi*aqed)
c      g      = sqrt(4d0*pi*aqcd)
      return
 99   done=.true.
      return
 55   format(i3,5e19.11)
      end

      subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff)
c********************************************************************
c     Writes one event from data file #lun according to LesHouches
c     ic(1,*) = Particle ID
c     ic(2.*) = Mothup(1)
c     ic(3,*) = Mothup(2)
c     ic(4,*) = ICOLUP(1)
c     ic(5,*) = ICOLUP(2)
c     ic(6,*) = ISTUP   -1=initial state +1=final  +2=decayed
c     ic(7,*) = Helicity
c********************************************************************
      implicit none
c
c     parameters
c
      double precision pi
      parameter (pi = 3.1415926d0)
c
c     Arguments
c
      integer lun, ievent
      integer nexternal, ic(7,*)
      double precision P(0:4,*),wgt
      double precision aqcd, aqed, scale
      character*140 buff
c
c     Local
c
      integer i,j,k
c
c     Global
c

c-----
c  Begin Code
c-----
c      aqed= gal(1)*gal(1)/4d0/pi
c      aqcd = g*g/4d0/pi

      write(lun,'(a)') '<event>'
      write(lun,'(i2,i4,4e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
      do i=1,nexternal
         write(lun,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
     $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
      enddo
      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))
      write(lun,'(a)') '</event>'
      return
 51   format(i9,5i5,5e19.11,f3.0,f4.0)
      end

      subroutine write_comments(lun)
c********************************************************************
c     Outputs all of the banner comment lines back at the top of
c     the file lun.
c********************************************************************
      implicit none
c
c     Arguments
c
      integer lun
c
c     Local
c
      character*(80) buff
c
c     Global
c
      logical banner_open
      integer lun_ban
      common/to_banner/banner_open, lun_ban

c-----
c  Begin Code
c-----
c      write(*,*) 'Writing comments'
      if (banner_open) then
         rewind(lun_ban)
         do while (.true.)
            read(lun_ban,'(a79)',end=99,err=99) buff
            write(lun,'(a79)') buff
c            write(*,*) buff
         enddo
 99      close(lun_ban)
         banner_open = .false.
      endif
      end


text:
Address(Whaler Avenue) is a residential street in West Davis. It forms the northern edge of Whaleback Park.

 Intersecting Streets

  Marina Circle, and across that intersection Beacon Bay Place
  Halifax Avenue


text:
*DECK DCBRT
      DOUBLE PRECISION FUNCTION DCBRT (X)
C***BEGIN PROLOGUE  DCBRT
C***PURPOSE  Compute the cube root.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C2
C***TYPE      DOUBLE PRECISION (CBRT-S, DCBRT-D, CCBRT-C)
C***KEYWORDS  CUBE ROOT, ELEMENTARY FUNCTIONS, FNLIB, ROOTS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DCBRT(X) calculates the double precision cube root for
C double precision argument X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, D9PAK, D9UPAK
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  DCBRT
      DOUBLE PRECISION X, CBRT2(5), Y, CBRTSQ,  D9PAK, D1MACH
      SAVE CBRT2, NITER
      DATA CBRT2(1) / 0.6299605249 4743658238 3605303639 11 D0 /
      DATA CBRT2(2) / 0.7937005259 8409973737 5852819636 15 D0 /
      DATA CBRT2(3) / 1.0 D0 /
      DATA CBRT2(4) / 1.2599210498 9487316476 7210607278 23 D0 /
      DATA CBRT2(5) / 1.5874010519 6819947475 1705639272 31 D0 /
      DATA NITER / 0 /
C***FIRST EXECUTABLE STATEMENT  DCBRT
      IF (NITER.EQ.0) NITER = 1.443*LOG(-.106*LOG(0.1*REAL(D1MACH(3)))
     1  ) + 1.0
C
      DCBRT = 0.D0
      IF (X.EQ.0.D0) RETURN
C
      CALL D9UPAK (ABS(X), Y, N)
      IXPNT = N/3
      IREM = N - 3*IXPNT + 3
C
C THE APPROXIMATION BELOW IS A GENERALIZED CHEBYSHEV SERIES CONVERTED
C TO POLYNOMIAL FORM.  THE APPROX IS NEARLY BEST IN THE SENSE OF
C RELATIVE ERROR WITH 4.085 DIGITS ACCURACY.
C
      Z = Y
      DCBRT = .439581E0 + Z*(.928549E0 + Z*(-.512653E0 + Z*.144586E0))
C
      DO 10 ITER=1,NITER
        CBRTSQ = DCBRT*DCBRT
        DCBRT = DCBRT + (Y-DCBRT*CBRTSQ)/(3.D0*CBRTSQ)
 10   CONTINUE
C
      DCBRT = D9PAK (CBRT2(IREM)*SIGN(DCBRT,X), IXPNT)
      RETURN
C
      END


text:
************************************************************************
*
*     Function EXPVAR                               Called by: TAC_LOAD
*
*     Estimate covariance of two EXPON estimators of the form
*     EXPON(SIGMA^2,A,B), assuming that SIGMA^2 is a
*     GAMMA(ALPHA,S2*BETA) random variable (where BETA is a parameter
*     = 2/(N-1) in standard model)
*
*     result is exact aside from machine error
*
*     S2        sample estimate of SIGMA^2
*     ALPHA     parameter of gamma distribution = (N-1)/2 in std.model
*     A1        1st coefficient in front of SIGMA^2
*     B1        1st coefficient in front of SIGMA
*     A2        2nd coefficient in front of SIGMA^2
*     B2        2nd coefficient in front of SIGMA
*
************************************************************************
      DOUBLE PRECISION FUNCTION EXPVAR(S2,ALPHA,A1,B1,A2,B2)
*
*     subroutine args
*
      DOUBLE PRECISION S2,ALPHA,A1,B1,A2,B2
*
*     local vars
*
      INTEGER*4 I,J,N
      DOUBLE PRECISION ARG,A,SIGMA,TOL,TSUM,DELTA
      DOUBLE PRECISION C1(0:1000),C2(0:1000),GAMM(0:1000)
      DATA TOL/1.D-6/
*
*     function declaration
*
      DOUBLE PRECISION DLNGAM
*
*
*
      SIGMA = SQRT(S2)

      ARG = MAX(ABS(A1)*S2+ABS(B1)*SIGMA,ABS(A2)*S2+ABS(B2)*SIGMA)
      A=MAX(ABS(A1),ABS(A2))
      N = MIN(998.D0,
     &        30.D0+ARG*(5.4-0.2*ARG+A
     &        *(8.28-1.11*A+1.98*ARG)+1.33*(MAX(ABS(B1),ABS(B2)))))

 10   EXPVAR = 0.D0
      CALL SQUARE(N,A1,B1,C1)
      CALL SQUARE(N,A2,B2,C2)

      GAMM(0) = 1.D0
      GAMM(1) = EXP(DLNGAM(ALPHA+0.5D0)-DLNGAM(ALPHA))
      DO 30 J=0,N
         GAMM(J+2) = (ALPHA+0.5D0*J)*GAMM(J)
         TSUM = 0.D0
         DO 20 I=0,J
            TSUM = TSUM+C1(I)*C2(J-I)*GAMM(J)/(GAMM(I)*GAMM(J-I))
 20      CONTINUE
         DELTA = SIGMA**J*TSUM
         EXPVAR = EXPVAR+DELTA
 30   CONTINUE
      IF (DELTA.LE.TOL*EXPVAR) RETURN

 40   N = N+10
      IF (N.LT.998) GOTO 10
*
*     Error: failure to converge in 998 terms
*
      CALL ERROR7(12)

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


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

20090424 18:53:23 nbsp Welcome to the Wiki. Thanks for registering under your real name. Weve had a lot of problems recently with people using Role Accounts that diminish the creation of a sense of Wiki Community community here on the wiki. For instance Users/CampusRec someone created a role account for Campus Recreation, but in doing so they created some real problems. If that account is shared between users it may be a violation of campus cybersafety policies and good practices, and if it is being used by only one person then they are contributing to a growing problem here on the wiki with Wiki Community/Organizational accounts Role Accounts. Either way there are now two pages with information purporting to tell users about Campus Recreation and Users/CampusRec. The other issue it raises is that if someone else comes along and registers Users/CampusRecreation and starts announcing events that dont exist and defaming the Chancellor no one will know if that account is really an official spokeperson for Campus Recreation or just a kook who isnt willing to stand behind their words with their importance of using your RealName. By registering with a real name, or a pseudonym and establishing an identity the issue of Role impersonation is reduced. Users/JasonAller

20100520 18:28:29 nbsp Howdy!  Why did you undo all my changes to the Campus Special Coupon Book entry?  I had added links and cleaned up several things (like removing the empty fax field). Users/JabberWokky Evan JabberWokky Edwards

20100521 15:48:13 nbsp No worries... thats why I asked.  Cooperative editing often involves a bit of communication, especially since we all mistakenly stomp over somebody elses edits from time to time.  Ill flip it back now.  Thanks for the contribution, by the way. Users/JabberWokky

20100521 18:45:14 nbsp (re: http://daviswiki.org/Users/JabberWokky?actiondiff&version21585&version11584 your comment on the wokks page) A redirect would probably be more appropriate (though Im pretty sure that was what you were actually referring to by link). Redirects are frequently invoked to http://daviswiki.org/Physics_and_Geology?actionshow restore continuity to renamed pages, http://daviswiki.org/Hutchinson_Drive?actionshow correct common misspellings, and http://daviswiki.org/LAWR?actionshow serve as shorthand for inline hyperlinks. Of course, they are certainly not limited to these uses, and if you feel that creating a redirect on 100 Dollar Coupon Book would help, you are certainly free to do so. (Heres the http://daviswiki.org/Templates/Redirect?actionshow redirect template, by the way.) Users/EBT
  I just created one (didnt see your comment here, EBT).  If you go to 100 Dollar Coupon Book, you are now redirected. jw

20100523 14:28:10 nbsp Thanks guys, I did mean redirect Users/MattSmidebush


text:
Sports teams at Davis Senior High School have traditionally been known as Blue Devils. The womens basketball team has the http://www.davisbluedevils.com/ DavisBlueDevils.com domain.  The football teams website can be found at http://davishighfootball.org/ davishighfootball.org and is supported in part by the Davis Blue Devil Football Backers.

The Blue Devil as DHS mascot was chosen by DHS students at the end of WWI after the Blue Devils of France, the same source for Duke Universitys Blue Devil mascot.

“The Chasseurs Alpins Corps … are popularly known as The Blue Devils. The Devil part of their nickname is in flattering allusion to, their valour; the adjective refers to the dark blue colour of their uniforms. They display an almost religious observance in the fulfilment of their vow, indicated by the motto of their corps, Never to be taken alive. The corps is recruited from the mountain country of France the men are remarkable for their sturdy physique (they are mostly men of medium height) and for their tenacious courage. They are accustomed to hardship and fatigue, they live on simple fare and by tradition are deeply imbued with antagonism to the Germans. They are silent men. Their squareset faces seem to speak of successful struggle against the mighty forces of mountain dangers; their strong backs and their frank manners. harmonise well with the brave deeds for which they are world famous”  Mr. Poilu, Notes and Sketches with The Fighting French by Herbert Ward, Chevalier de la Legion D’Honneur, Croix Guerre

Does anyone know if there is a central website for Davis High School athletics other than http://dhs.djusd.k12.ca.us/xhtml/departments/sports.html?

 Other NonDavis Blue Devils
  http://www.bluedevils.org/ Concord Blue Devils, who are http://www.bluedevils.org/about/history.php named for nearby Mt. Diablo
  http://www.goduke.com/ Duke University, who are named for the same http://www.lib.duke.edu/archives/history/why_blue_devil.html French Blue Devils.

If you are alumni of Duke, the use of Blue Devils is pervasive, and seeing it in the corner of your eye all over the place when you move to Davis can be disorienting until you realize that it is also the local high school logo.


text:
c
c This file contains dummy functions for the system matrices used
c in for the linear Kalman filter.
c
      SUBROUTINE AMAT(XM,XP,UO,T,A,NPARAM,N,M)
      INTEGER NPARAM,N,M
      DOUBLE PRECISION XM(NPARAM),XP(N),UO(M),T,A(N,N)
      INTEGER I,J
      DO 10 I=1,N
      DO 10 J=1,N
         A(J,I)=0.0D0
  10  CONTINUE
      RETURN
      END
c
      SUBROUTINE BMAT(XM,XP,UO,T,B,NPARAM,N,M)
      INTEGER NPARAM,N,M
      DOUBLE PRECISION XM(NPARAM),XP(N),UO(M),T,B(N,M)
      INTEGER I,J
      DO 10 I=1,M
      DO 10 J=1,N
         B(J,I)=0.0D0
  10  CONTINUE
      RETURN
      END
c
      SUBROUTINE CMAT(XM,XP,UO,T,C,NPARAM,N,M,S)
      INTEGER NPARAM,N,M,S
      DOUBLE PRECISION XM(NPARAM),XP(N),UO(M),T,C(S,N)
      INTEGER I,J
      DO 10 I=1,N
      DO 10 J=1,S
         C(J,I)=0.0D0
  10  CONTINUE
      RETURN
      END
c
      SUBROUTINE DMAT(XM,XP,UO,T,D,NPARAM,N,M,S)
      INTEGER NPARAM,N,M,S
      DOUBLE PRECISION XM(NPARAM),XP(N),UO(M),T,D(S,M)
      INTEGER I,J
      DO 10 I=1,M
      DO 10 J=1,S
         D(J,I)=0.0D0
  10  CONTINUE
      RETURN
      END

max_stars_repo_name: perNyfelt/renjinSamplesAndTests

max_issues_repo_name: perNyfelt/renjinSamplesAndTests

max_forks_repo_name: perNyfelt/renjinSamplesAndTests


text:
Describe Users/gracesmith here.

20130715 15:59:15 nbsp HiI was wondering why you are using two usernames: Users/gracesmith and Users/gracejones?? Users/PeteB


text:
Jump to Timeline #Navigation Navigation

  The New Schmeiser Barn is built


text:
      subroutine tanhrt(nb, Rval, tn, nTimes, vecTcalc, vecRTrel)

  !====================================================================
  ! Hyperbolic Tangent Model for R(t)
  !
  !====================================================================
      implicit none

      integer nb, nTimes
      real*8 Rval(nb), tn(nb)
      real*8 vecTcalc(nTimes), vecRTrel(ntimes)
      real*8 wl
      parameter (wl = 5.0d0)
      real*8 t_cur, beta_tmp
      integer i, j

 ! convert from absolute days in each R(t) value to relative number of days

      do i = 2, (nb-1)
         tn(i) = tn((i-1)) + tn(i)
      enddo

      do j = 1, nTimes
         t_cur = vecTcalc(j)
!     Calculate R(t)
         beta_tmp = Rval(1) + Rval(nb)
         do i = 2, nb
            beta_tmp = beta_tmp + (Rval(i) - Rval((i-1))) *
     $           tanh((t_cur-tn((i-1)))/wl)
         enddo
         beta_tmp = beta_tmp * 0.5
         vecRTrel(j) = beta_tmp
      enddo

      return
      end subroutine tanhrt

max_stars_repo_name: predsci/PSI-CARM

max_issues_repo_name: predsci/PSI-CARM

max_forks_repo_name: predsci/PSI-CARM


text:
!===============================================================================
!> This module contains all experiment specific information necessary to solve
!> the heat diffusion problem in DISPATCH
!===============================================================================
MODULE solver_mod
  USE io_mod
  USE trace_mod
  USE extras_mod
  implicit none
  private
  type, public, extends(extras_t):: solver_t
    real, pointer:: radii(:,:,:)
    real, pointer:: temperature(:,:,:)
    real:: initial_temperature=1.0
    real:: outside_temperature=2.0
    real:: radius=0.45
  contains
    procedure:: init
    procedure:: update
    procedure:: boundary_condition
  end type
CONTAINS

!===============================================================================
!> Setup an initial state.  The box size is 1, with the origin at 0,0,0, so
!> the max radius is 0.5
!===============================================================================
SUBROUTINE init (self)
  class(solver_t):: self
  integer:: m(6), i1, i2, i3, it
  !-----------------------------------------------------------------------------
  call trace%begin('solver_t%init')
  !-----------------------------------------------------------------------------
  ! This initializers a default do-nothing solver, and allocates the
  ! self%mem() 6-dimenstional array (nx,ny,nz,nv,nt,nw), in which the 4th
  ! index (= 1) here is the variable index, and the 5th index is the time slot.
  !-----------------------------------------------------------------------------
  self%nv = 1
  self%nt = 5
  call self%patch_t%init
  self%temperature => self%mem(:,:,:,1,1,1)
  self%temperature = self%initial_temperature
  !-----------------------------------------------------------------------------
  ! For convenience, we also allocate the radius (distance from origin) for
  ! each cell.  The array should have the same size as the temperature:
  !-----------------------------------------------------------------------------
  m = shape(self%mem)
  allocate (self%radii(m(1),m(2),m(3)))
  do i3=1,m(3)
  do i2=1,m(2)
  do i1=1,m(1)
    !---------------------------------------------------------------------------
    ! Each mesh (x,y,z) has properties stored in a data type (object), where
    ! %p is the position (of the center), and %r are the relative cell positions
    !---------------------------------------------------------------------------
    self%radii(i1,i2,i3) = sqrt((self%mesh(1)%p + self%mesh(1)%r(i1))**2 + &
                                (self%mesh(2)%p + self%mesh(2)%r(i2))**2 + &
                                (self%mesh(3)%p + self%mesh(3)%r(i3))**2)
  end do
  end do
  end do
  !-----------------------------------------------------------------------------
  ! Set temperature inside and then apply boundary condition
  !-----------------------------------------------------------------------------
  call self%boundary_condition (1)
  call trace%end
END SUBROUTINE init

!===============================================================================
!> Apply the boundary condition T=temperature(1) outside the given radius
!===============================================================================
SUBROUTINE boundary_condition (self, it)
  class(solver_t):: self
  integer:: it
  real, pointer, contiguous:: f(:,:,:)
  !-----------------------------------------------------------------------------
  call trace%begin('solver_t%boundary_condition')
  !-----------------------------------------------------------------------------
  f => self%mem(:,:,:,1,it,1)
  where (self%radii > self%radius)
    f = self%outside_temperature
  end where
  call trace%end
END SUBROUTINE boundary_condition

!===============================================================================
!> Update the solution, by adding a fraction times the diffusion operator
!===============================================================================
SUBROUTINE update (self)
  class(solver_t):: self
  real, pointer, contiguous:: v(:,:,:)
  real, allocatable:: d2f(:,:,:)
  integer:: m(3), i1, i2, i3
  integer, save:: itimer=0
  !----------------------------------------------------------------------------
  call trace%begin('solver_t%update', itimer=itimer)
  call self%output
  !----------------------------------------------------------------------------
  ! Copy the variable (temperature) from the current slot (%it) to the %new slot
  !----------------------------------------------------------------------------
  v => self%mem(:,:,:,1,self%new,1)
  v =  self%mem(:,:,:,1,self%it ,1)
  call self%boundary_condition (self%new)
  !----------------------------------------------------------------------------
  ! Allocate a scratch variable for the diffusion (Laplace) operator
  !----------------------------------------------------------------------------
  m = shape(v)
  allocate (d2f(m(1),m(2),m(3)))
  !----------------------------------------------------------------------------
  ! Evaluate a 6-point diffusion operator over the internal patch region
  ! (no need to do it in the guard zones)
  !----------------------------------------------------------------------------
  do i3=self%mesh(3)%li,self%mesh(3)%ui
  do i2=self%mesh(2)%li,self%mesh(2)%ui
  do i1=self%mesh(1)%li,self%mesh(1)%ui
    d2f(i1,i2,i3) = v(i1+1,i2  ,i3  ) &
                  + v(i1-1,i2  ,i3  ) &
                  + v(i1  ,i2+1,i3  ) &
                  + v(i1  ,i2-1,i3  ) &
                  + v(i1  ,i2  ,i3+1) &
                  + v(i1  ,i2  ,i3-1) &
                  - v(i1  ,i3  ,i3  )*6.0
  end do
  end do
  end do
  !----------------------------------------------------------------------------
  ! Set the time step, update the temperature, and apply boundary condition
  !----------------------------------------------------------------------------
  self%dtime = self%courant
  v = v + self%dtime*d2f
  call self%boundary_condition (self%new)
  !----------------------------------------------------------------------------
  ! Clean up and count operations
  !----------------------------------------------------------------------------
  deallocate (d2f)
  call self%counter_update
  call trace%end (itimer)
END SUBROUTINE update

END MODULE solver_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
active component C {

  event E1 severity activity low id 0x10 format ""

  event E2 severity activity low id 0x10 format ""

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
!---------------------------------------------------------------------!
!                                PHAML                                !
!                                                                     !
! The Parallel Hierarchical Adaptive MultiLevel code for solving      !
! linear elliptic partial differential equations of the form          !
! (PUx)x + (QUy)y + RU = F on 2D polygonal domains with mixed         !
! boundary conditions, and eigenvalue problems where F is lambda*U.   !
!                                                                     !
! PHAML is public domain software.  It was produced as part of work   !
! done by the U.S. Government, and is not subject to copyright in     !
! the United States.                                                  !
!                                                                     !
!     William F. Mitchell                                             !
!     Applied and Computational Mathematics Division                  !
!     National Institute of Standards and Technology                  !
!     william.mitchell@nist.gov                                       !
!     http://math.nist.gov/phaml                                      !
!                                                                     !
!---------------------------------------------------------------------!

!----------------------------------------------------
! This file contains the user supplied module phaml_user_mod
!----------------------------------------------------

module phaml_user_mod

!----------------------------------------------------
! This module contains user global data.
!
! This version illustrates eigenvalue problems via the Schroedinger Equation
! for a simple harmonic oscillator where the cross section of the well is an
! ellipse with axes determined by prob_param.
! -uxx -uyy + prob_param*(x^2+3*y^2)/2 * u = lambda * u
! prob_param must be positive.
!
! For num_eval=1 and lambda0=-huge(0.), the eigenvalue is
! sqrt(prob_param*(2+sqrt(3)))
! and the solution (wave function) is
! exp(-sqrt(prob_param/8)*(x^2+sqrt(3)*y^2).
!
!----------------------------------------------------

!----------------------------------------------------
! Other modules used are:

use phaml
!----------------------------------------------------

implicit none

!----------------------------------------------------
! Global variables

real(my_real), save :: prob_param
!----------------------------------------------------

end module phaml_user_mod

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c *** Example of the Arnoldi iteration
c *** Copyright (c) 2020, Ioannis Nompelis
c *** See "LICENSE" for terms of use.
c *** Ioannis Nompelis <nompelis@nobelware.com>
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       Program main
       Implicit none
       Real*8,allocatable,dimension(:,:) :: amat
       Real*8,allocatable,dimension(:) :: bvec
       Real*8,allocatable,dimension(:,:) :: q,h
       Integer*4 :: size, num, i,k

       size = 10
       num = 4

       allocate( amat(size,size), bvec(size), q(size,num), h(num,num) )
       amat(:,:) = 0.0d0
       bvec(:) = 0.0d0

       call setup( size, amat, bvec )

       call arnoldi_basis( size, amat, bvec, num, q, h )

       PRINT*,'Vectors'
       do k = 1,num
          PRINT*,'k=',k
          do i = 1,size
             PRINT*,'i=',i,'   q=',q(i,k)
          enddo
       enddo

       call check_orthogonality ( size, num, q )

       End

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       Subroutine arnoldi_basis( size, amat, bvec, num, qvecs, hmat )
       Implicit none
       Integer*4, intent(IN) :: size, num
       Real*8,dimension( size, size ), intent(IN) :: amat
       Real*8,dimension( size ), intent(IN) :: bvec
       Real*8,dimension( size, num ), intent(OUT) :: qvecs
       Real*8,dimension( num, num ), intent(OUT) :: hmat
       Integer*4 :: i,j,k,kk
       Real*8 :: r

       !--- form the first vector assuming it is not normalized
       r = 0.0d0
       do i = 1,size
          r = r + bvec(i)*bvec(i)
       enddo
       r = sqrt(r)
       !--- magnitude of the initial vector; assume it is not a unit vector
       hmat(1,1) = r
       !--- store the first vector normalized
       r = 1.0d0/r
       do i = 1,size
          qvecs(i,1) = bvec(i) * r
       enddo

       !--- form the subsequent vectors (so that there are "num" in total)
       do k = 2,num
          !--- map to a new vector using the matrix
          do i = 1,size
             r = 0.0d0
             do j = 1,size
                r = r + amat(i,j) * qvecs(j,k-1)
             enddo
             qvecs(i,k) = r
          enddo

          !--- make orthogonal to all previous vectors (up to k-1)
          do kk = 1,k-1
             r = 0.0d0
             do i = 1,size
                r = r + qvecs(i,k) * qvecs(i,kk)
             enddo
             !--- store projetion to the Hezenberg matrix
             hmat(k,kk) = r
             !--- subtract the projection vectorwise
             do i = 1,size
                qvecs(i,k) = qvecs(i,k) - r * qvecs(i,kk)
             enddo
          enddo
          !--- calculate the magnitude
          r = 0.0d0
          do i = 1,size
             r = r + qvecs(i,k) * qvecs(i,k)
          enddo
          r = sqrt(r)
          hmat(k,k) = r
          r = 1.0d0/r
          do i = 1,size
             qvecs(i,k) = qvecs(i,k) * r
          enddo
       enddo

       End

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       Subroutine setup( size, amat, bvec )
       Implicit none
       Integer*4, intent(IN) :: size
       Real*8,dimension( size, size ), intent(INOUT) :: amat
       Real*8,dimension( size ), intent(INOUT) :: bvec
       Integer*4 :: i,j

       do i = 1,size
       do j = 1,size
          if( i.eq.j ) amat(i,j) = -2.0d0
          if( i.eq.j-1) amat(i,j) = 1.0d0
          if( i.eq.j+1) amat(i,j) = 1.0d0
       enddo
       enddo
       amat(1,:) = 0.0d0
       amat(1,1) = 1.0d0
       amat(size,:) = 0.0d0
       amat(size,size) = 1.0d0

       PRINT*,'Matrix'
       do i = 1,size
          PRINT*,'Row',i
       do j = 1,size
          PRINT*,amat(i,j)
       enddo
       enddo

       do i = 1,size
          bvec(i) = 1.0d0
       enddo

       PRINT*,'Vector'
       do i = 1,size
          PRINT*,'Row',i, bvec(i)
       enddo

       End

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       Subroutine check_orthogonality ( size, num, q )
       Implicit none
       Integer*4, intent(IN) :: size, num
       Real*8,dimension( size, num ), intent(INOUT) :: q
       Integer*4 :: m,n,i
       Real*8 :: r

       do m = 1,num
       do n = 1,num
          r = 0.0d0
          do i = 1,size
             r = r + q(i,m) * q(i,n)
          enddo
          PRINT*,'Test ortho: ',m,n,r
       enddo
       enddo

       End

max_stars_repo_name: nompelis/ArnoldiIteration

max_issues_repo_name: nompelis/ArnoldiIteration

max_forks_repo_name: nompelis/ArnoldiIteration


text:
Aisle 5

Email:  MailTo(aisle5theband AT yahoo DOT com)

Based in: Davis

Band Members:
 Fawzi Haimor  Lead vocals/rhythm guitar
 Dominic Cardenas  Lead guitar
 Adam Sartain  Keyboards/ synthesizers
 Steve David  Bass
 Sean OLeary  Drums

Image(Aisle_5_band.jpg, thumbnail)


text:
      program intrinsics_2
      integer*4 k,l,m,s,i,j
      character*(2) ss(6)
      double precision xx
      real*4 x
      k = -1
      l = 0
      m = 1
      s = -1
      write(*,*) isign(k,s), isign(l,s), isign(m,s)
      s = 0
      write(*,*) isign(k,s), isign(l,s), isign(m,s)
      s = 1
      write(*,*) isign(k,s), isign(l,s), isign(m,s)
      xx = 7.2973525698D-3
      x = sngl(xx)
      write(*,*) xx, x, dble(x)
      write(*,*) idnint(-1.37D0), idnint(-1.73D0)
      write(*,*) idnint(xx), idnint(0.5D0), idnint(1.5D0)
      write(*,*) idnint(1.37D0), idnint(1.73D0)
      ss(1) = 'ab'
      ss(2) = 'AB'
      ss(3) = '0x'
      ss(4) = ',!'
      ss(5) = ' b'
      ss(5) = '# '
      write(*,*) ((llt(ss(i),ss(j)) , i=1,6), j=1,6)
      write(*,*) ((lle(ss(i),ss(j)) , i=1,6), j=1,6)
      write(*,*) ((lgt(ss(i),ss(j)) , i=1,6), j=1,6)
      write(*,*) ((lge(ss(i),ss(j)) , i=1,6), j=1,6)
      end


text:
module util_module

  use amrex_constants_module
  use amrex_fort_module, only : rt => amrex_real

  implicit none

contains

  subroutine get_xn(nz, xn_zone)

    use network, only: nspec, network_species_index
    use extern_probin_module, only: primary_species_1, primary_species_2, primary_species_3
    use amrex_error_module

    integer, intent(in) :: nz
    real(rt), intent(  out) :: xn_zone(:,:)

    real(rt) :: Xp_min, Xp_max, dX
    integer :: n1, n2, n3, nprim
    integer :: is1, is2, is3
    integer :: k, n
    real(rt) :: excess

    ! get the primary species indices
    nprim = 0

    is1 = network_species_index(primary_species_1)
    if (is1 > 0) then
       nprim = nprim+1
    end if

    is2 = network_species_index(primary_species_2)
    if (is2 > 0) then
       nprim = nprim+1
    end if

    is3 = network_species_index(primary_species_3)
    if (is3 > 0) then
       nprim = nprim+1
    end if

    if (nprim == 0) then
       call amrex_error("ERROR: no primary species set")
    end if

    ! figure out how many zones to allocate to the each of the primary
    ! species and the extrema for the primary species
    if (nprim == 1) then
       n1 = nz
       n2 = 0
       n3 = 0

       Xp_min = 0.2_rt
       Xp_max = 0.9_rt

    else if (nprim == 2) then
       n1 = nz/2
       n2 = nz - n1
       n3 = 0

       Xp_min = 0.2_rt
       Xp_max = 0.8_rt

    else if (nprim == 3) then
       n1 = nz/3
       n2 = nz/3
       n3 = nz - n1 - n2

       Xp_min = 0.2_rt
       Xp_max = 0.7_rt

    end if

    do k = 1, nz
       xn_zone(:, k) = 0.0_rt

       if (k <= n1) then
          if (nprim >= 2) xn_zone(is2, k) = Xp_min/2
          if (nprim >= 3) xn_zone(is3, k) = Xp_min/2

          dX = (Xp_max - Xp_min) / max((n1 - 1), 1)
          xn_zone(is1, k) = Xp_min + (k - 1)*dX

       else if (nprim >= 2 .and. k <= n1 + n2) then
          xn_zone(is1, k) = Xp_min/2
          if (nprim >= 3) xn_zone(is3, k) = Xp_min/2

          dX = (Xp_max - Xp_min) / max((n2 - 1), 1)
          xn_zone(is2, k) = Xp_min + (k - (n1 + 1))*dX

       else
          xn_zone(is1, k) = Xp_min/2
          xn_zone(is2, k) = Xp_min/2

          dX = (Xp_max - Xp_min) / max((n3 - 1), 1)
          xn_zone(is3, k) = Xp_min + (k - (n1 + n2 + 1))*dX

       end if

       excess = ONE - sum(xn_zone(:, k))

       do n = 1, nspec
          if (n == is1 .or. &
              (n == is2 .and. nprim >= 2) .or. &
              (n == is3 .and. nprim >= 3)) cycle

          xn_zone(n, k) = excess / (nspec - nprim)

       end do

    end do

  end subroutine get_xn

end module util_module

max_stars_repo_name: ChrisDeGrendele/Microphysics

max_issues_repo_name: ChrisDeGrendele/Microphysics

max_forks_repo_name: ChrisDeGrendele/Microphysics


text:
      SUBROUTINE blk3d(a, bm1, bp1, srces, mblk, nblocks)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE stel_kinds
      USE safe_open_mod
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER, PARAMETER :: bytes_per_rprec = 8
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(in) :: nblocks, mblk
      REAL(rprec), DIMENSION(mblk,mblk,nblocks), INTENT(in) ::
     1                       a, bp1, bm1
      REAL(rprec), DIMENSION(mblk,nblocks), INTENT(inout) :: srces
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
!      INTEGER :: ibuph, incnow, irecl, incbu, iunit=102, ndisk
      INTEGER :: k, k1, ier
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
      REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ainv
      REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE :: ql
C-----------------------------------------------
C   E x t e r n a l  S u b r o u t i n e s
C-----------------------------------------------
      EXTERNAL la_getrf, la_getrs
C-----------------------------------------------
c  modified (June, 2003, ORNL):         S. P. Hirshman
c-----------------------------------------------------------------------
c
c  this subroutine solves a block-tridiagonal system of equations.
c
c-----------------------------------------------------------------------
c  INPUT
c  mblk                : block size
c  nblocks             : number of blocks
c  a                   : diagonal block
c  bm1, bp1            : lower, upper block (see equation below)
c  srces(mblk,nblocks) : right-hand size source
c
c  LOCAL VARIABLES
c  iunit               : unit number for block-tridiagonal solution disk file.
c
c  solutions are indexed in m-n fourier-space, legendre-space. the tri-diagonal
c  equation is:
c
c           bm1 * f(l-1) + a * f(l) + bp1 * f(l+1) = source(l)
c
c     GENERAL SOLUTION SCHEME APPLIED TO EACH BLOCK ROW (INDEX L)
c
c     1. Start from row N and solve for x(N) in terms of x(N-1):
c
c        x(N) = -q(N)*x(N-1) + r(N)
c
c        q(N) =  a(N)[-1] * bm1;    r(N) = a(N)[-1] * s(N)
c
c        where a(N)[-1] is the inverse of a(N)
c
c     2. Substitute for lth row to get recursion equation fo q(l) and r(l):
c
c        x(l) = -q(l)*x(l-1) + r(l), in general, where:
c
c        q(l) = (a(l) - bp1(l)*q(l+1))[-1] * bm1(l)
c
c        r(l) = (a(l) - bp1(l)*q(l+1))[-1] * (s(l) - bp1(l)*r(l+1))
c
c     3. At row l = 1, bm1(1) = 0 and get an equation for x(1) corresponding to q(1) = 0:
c
c        x(1) = r(1)
c
c     4. Finally, can back-solve for x(N) in terms of x(N-1) from eqn.(1) above
c
c
c     NUMERICAL IMPLEMENTATION (USING LAPACK, la_... ROUTINES)
c
c     1. CALL la_getrf:   Perform LU factorization of diagonal block (A)
c     2. CALL la_getrs:   With multiple (mblk) right-hand sides, to do block inversion
c                         operation, A X = B  (store result in X; here B is a matrix)
c     3. CALL la_getrs:   With single right hand side (source) to solve A x = b (b a vector)
c
!      ndisk = mblk

      ALLOCATE (ql(mblk,mblk, nblocks), ainv(mblk,mblk), ipiv(mblk),
     1          stat=ier)
      IF (ier .ne. 0) STOP 'Allocation error in blk3d'

c  create disk file for doing direct access i/o.

!      incnow = ndisk
!      irecl  = bytes_per_rprec*incnow
!      incbu = 1 + (ndisk - 1)/incnow
!      ibuph = 0

!      iunit = 10
!      CALL safe_open(iunit, ier, 'NULL', 'scratch', 'unformatted',
!     1     irecl, 'DIRECT')
!      IF (ier .ne. 0) STOP 'Error opening scratch file in blk3d'

      ainv = a(:,:,nblocks)

c  main loop. load and process (backwards) block-rows nblocks to 1.

      BLOCKS: DO k = nblocks, 1, -1
!
!     Compute (and save) ql(nblocks) = a(nblocks)[-1] * ql, and source terms A-1 * srces [A-1 == inv(A)]
!
         CALL la_getrf (mblk, mblk, ainv, mblk, ipiv, ier)
         IF (ier .ne. 0) GOTO 200

         IF (k .gt. 1) THEN
            ql(:,:,k) = bm1(:,:,k)
            CALL la_getrs('n', mblk, mblk, ainv,
     1                    mblk, ipiv, ql(1,1,k), mblk, ier)
            CALL la_getrs('n', mblk, 1,    ainv,
     1                    mblk,ipiv,srces(1,k),mblk,ier)

!            CALL wrdisk(iunit, ql, ndisk, incnow, ibuph, incbu, ier)
!            IF (ier .ne. 0) GOTO 302
         ELSE
            CALL la_getrs('n', mblk, 1,   ainv,
     1                    mblk,ipiv,srces(1,k),mblk,ier)
         END IF

         k1 = k - 1
         IF (k1 .eq. 0) EXIT

!
!      Update effective diagonal "a" matrix and source terms and store ql 2 l-steps back
!
         ainv  = a(:,:,k1) - MATMUL(bp1(:,:,k1), ql(:,:,k))
         srces(:,k1) = srces(:,k1) - MATMUL(bp1(:,:,k1),srces(:,k))

      END DO BLOCKS

c  backward solution sweep for block-rows k = 2 to nblocks
c  read blocks ql from disk.

      DO k = 2, nblocks
!         CALL rddisk (iunit, ql, ndisk, incnow, ibuph, ier)
!         IF (ier .ne. 0) GOTO 303
!         ibuph = ibuph - incbu
         srces(:,k) = srces(:,k) - MATMUL(ql(:,:,k),srces(:,k-1))
      END DO

      GOTO 400

c  error returns. ------------------------------------------------------

  200 CONTINUE
      WRITE (6, '(a,i4)') ' Error factoring matrix in blk3d: block = '
     1                     , k,' error id = ', ier
      GOTO 400
  301 CONTINUE
!      WRITE (6, '(a,i8)') ' BLK3D:   error in opening file:  ',
!     1   'RECL = ', irecl
  302 CONTINUE
      WRITE (6, '(a)') ' BLK3D:   error in I/O routine WRDISK'
  303 CONTINUE
      WRITE (6, '(a)') ' BLK3D:   error in I/O routine RDDISK'
      ier = -2
  305 CONTINUE
      WRITE (6, '(2/a,i4,2/)') ' BLK3D:   error detected:   ier =',
     1   ier
      STOP

c  destroy disk file and return. ---------------------------------------

  400 CONTINUE

!      CLOSE (iunit)

      DEALLOCATE (ainv, ql, ipiv)

      END SUBROUTINE blk3d


module CommandArguments
  use KindModule
  use ConstantsModule, only: ISTDOUT, LINELENGTH, LENHUGELINE
  use VersionModule,          only: VERSION, MFVNAM, IDEVELOPMODE
  use CompilerVersion
  use SimVariablesModule,     only: simfile
  use SimModule, only: store_error, ustop, store_error_unit,                   &
                       store_error_filename
  use InputOutputModule, only: urword
  !
  implicit none
  !
  private
  public :: GetCommandLineArguments
  !
  contains

  subroutine GetCommandLineArguments()
    ! -- dummy
    ! -- local
    character(len=LENHUGELINE) :: line
    character(len=LINELENGTH) :: header
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: cexe
    character(len=80) :: compiler
    character(len=20) :: cdate
    character(len=17) :: ctyp
    logical :: ltyp
    logical :: lexist
    integer(I4B) :: inunit = 0
    integer(I4B) :: ilen
    integer(I4B) :: istat
    integer(I4B) :: lloc
    integer(I4B) :: istart
    integer(I4B) :: istop
    integer(I4B) :: ival
    integer(I4B) :: i
    integer(I4B) :: ipos
    integer(I4B) :: iarg
    integer(I4B) :: iterm
    real(DP) :: rval
! ------------------------------------------------------------------------------
    !
    ! -- Get the command line string
    call GET_COMMAND(line, ilen, istat)
    !
    ! -- This will read mf6 executable
    lloc = 1
    call urword(line, lloc, istart, istop, 0, ival, rval, 0, inunit)
    !
    ! -- remove quotes around the mf6 executable
    do i = istart, istop
      if (line(i:i) == '"' .or. line(i:i) == "'") then
        line(i:i) = ' '
      end if
    end do
    !
    ! -- find name of executable without path
    ipos = index(line(istart:istop), '/', back=.TRUE.)
    if (ipos == 0) then
      ipos = index(line(istart:istop), '\', back=.TRUE.)
    end if
    if (ipos /= 0) then
      istart = ipos + 1
    end if
    !
    ! -- set mf6 executable name
    cexe = adjustl(line(istart:istop))
    !
    ! -- write header
    call get_compile_date(cdate)
    write(header, '(a,4(1x,a),a)') &
      trim(adjustl(cexe)), '- MODFLOW',                                          &
      trim(adjustl(VERSION)), '(compiled', trim(adjustl(cdate)), ')'
    !
    ! -- set ctyp
    if (IDEVELOPMODE == 1) then
      ctyp = 'Release Candidate'
      ltyp = .TRUE.
    else
      ctyp = 'Release'
      ltyp = .FALSE.
    end if
    !
    ! -- Read remaining arguments
    iarg = 0
    iterm = 0
    do
      call urword(line, lloc, istart, istop, 1, ival, rval, 0, inunit)
      if (line(istart:istop) == ' ') exit
      iarg = iarg + 1
      iterm = 1
      select case(line(istart:istop))
        case('-H', '-?', '--HELP')
          call write_usage(trim(adjustl(header)), trim(adjustl(cexe)))
        case('-V', '--VERSION')
          write(ISTDOUT,'(2a,2(1x,a))') &
            trim(adjustl(cexe)), ':', trim(adjustl(VERSION)), ctyp
        case('-DEV', '--DEVELOP')
          write(ISTDOUT,'(2a,l)') &
            trim(adjustl(cexe)), ': develop version', ltyp
        case('-C', '--COMPILER')
          call get_compiler(compiler)
          write(ISTDOUT,'(2a,1x,a)') &
            trim(adjustl(cexe)), ':', trim(adjustl(compiler))
        case default
          call write_usage(trim(adjustl(header)), trim(adjustl(cexe)))
          write(errmsg, '(2a,1x,a)') &
            trim(adjustl(cexe)), ': illegal option -', line(istart:istop)
          call store_error(errmsg)
      end select
    end do
    !
    ! -- no command line arguments - check if mfsim.nam exists
    if  (iarg == 0) then
      inquire(file=simfile, exist=lexist)
      if (.NOT. lexist) then
        iterm = 1
        write(errmsg, '(2a,2(1x,a))') &
           trim(adjustl(cexe)), ':', simfile, 'is not present in working directory.'
        call store_error(errmsg)
      end if
    end if
    !
    ! -- command line arguments present or mfsim.nam file does not exist
    if (iterm > 0) then
      call USTOP()
    end if
    !
    ! -- return
    return
  end subroutine GetCommandLineArguments

  subroutine write_usage(header, cexe)
    ! -- dummy
    character(len=*), intent(in) :: header
    character(len=*), intent(in) :: cexe
    ! -- local
    character(len=*), parameter :: OPTIONSFMT =                                    &
      "(/,                                                                      &
      &'Options   GNU long option   Meaning ',/,                                &
      &' -h, -?   --help            Show this message',/,                       &
      &' -v       --version         Display program version information.',/,    &
      &' -dev     --develop         Display program develop option mode.',/,    &
      &' -c       --compiler        Display compiler information.',/,           &
      &'                                                                    ',/,&
      &'Bug reporting and contributions are welcome from the community. ',/,    &
      &'Questions can be asked on the issues page[1]. Before creating a new',/, &
      &'issue, please take a moment to search and make sure a similar issue',/, &
      &'does not already exist. If one does exist, you can comment (most',/,    &
      &'simply even with just :+1:) to show your support for that issue.',/,    &
      &'                                                                    ',/,&
      &'[1] https://github.com/MODFLOW-USGS/modflow6/issues',/)"
! ------------------------------------------------------------------------------
    write(ISTDOUT,'(a,/,a,1x,a,15x,a,2(1x,a),2a,/,a,1x,a,1x,a,5x,a)') &
      trim(adjustl(header)),                                              &
      'usage:', cexe, 'run MODFLOW', trim(adjustl(MFVNAM)),               &
      'using "', trim(adjustl(simfile)), '"',                             &
      '   or:', cexe, '[options]',                                        &
      'retrieve program information'
    write(ISTDOUT, OPTIONSFMT)
    !
    ! -- return
    return
  end subroutine write_usage

end module CommandArguments


text:
!----------------------------------------------------------------
! (c) Copyright, 2016 by the Regents of the University of California.
! Fldmgerclass: Exchange guard grid data between neighboring processors
! class in Communication module of FUNCTION layer.
!
! MODULE  : ... Fldmgerclass
! VERSION : ... 1.0
!> @author
!> Ji Qiang, LBNL
!
! DESCRIPTION:
!> This class defines the functions to sum up the particle
!> contribution from neighboring processor domain, exchange
!> the potential, exchange the field for interpolation between
!> neighboring processors.
! Comments:
!----------------------------------------------------------------
        module Fldmgerclass
          use Timerclass
          use Pgrid2dclass
        contains
!----------------------------------------------------------------
! neighboring grid communication for the 3D open boundary conditions
        !> sum up the contributions from neighboring guard cells.
        subroutine guardsum1_Fldmger(x,innx,inny,innz,grid)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: innx, inny, innz
        double precision, dimension(innx,inny,innz), intent(inout) :: x
        type (Pgrid2d), intent(in) :: grid
        double precision, dimension(innx,innz) :: sendbuf1, recvbuf1
        double precision, dimension(innx,inny) :: sendbuf2, recvbuf2
        double precision, dimension(innz)  :: tmp1,sendtmp
        integer :: left,right,bottom,top,msid,&
                   ierr,i,j,k,nx,ny
        integer :: myid,myidx,myidy,comm2d,commcol,commrow,totnp,npy,&
                   npx
        integer, dimension(2) :: tmpcoord
        integer status(MPI_STATUS_SIZE)
        double precision :: t0
        integer :: nsend,jadd,kadd

        call getpost_Pgrid2d(grid,myid,myidy,myidx)
        call getcomm_Pgrid2d(grid,comm2d,commcol,commrow)
        call getsize_Pgrid2d(grid,totnp,npy,npx)

        call MPI_BARRIER(comm2d,ierr)
        call starttime_Timer(t0)

        if(myidx.ne.(npx-1)) then
          right = myidx + 1
        else
          right = MPI_PROC_NULL
        endif
        if(myidx.ne.0) then
          left = myidx - 1
        else
          left = MPI_PROC_NULL
        endif
        if(npx.gt.1) then
          kadd = 1
        else
          kadd = 0
        endif

        if(myidy.ne.npy-1) then
          top = myidy + 1
        else
          top = MPI_PROC_NULL
        endif
        if(myidy.ne.0) then
          bottom = myidy -1
        else
          bottom = MPI_PROC_NULL
        endif
        if(npy.gt.1) then
          jadd = 1
        else
          jadd = 0
        endif

        if(npy.gt.1) then

        nsend = innz*innx
          call MPI_IRECV(recvbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,bottom,1,commcol,&
                         msid,ierr)
          do k = 1, innz
            do i = 1, innx
              sendbuf1(i,k) = x(i,inny,k)
            enddo
          enddo
          call MPI_SEND(sendbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,top,1,commcol,&
                        ierr)
          call MPI_WAIT(msid,status,ierr)
        if(myidy.ne.0) then
          do k = 1, innz
            do i = 1, innx
              x(i,2,k) = x(i,2,k) + recvbuf1(i,k)
            enddo
          enddo
        endif

          call MPI_IRECV(recvbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,top,2,commcol,&
                         msid,ierr)
          do k = 1, innz
            do i = 1, innx
              sendbuf1(i,k) = x(i,1,k)
            enddo
          enddo
          call MPI_SEND(sendbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,bottom,2,commcol,&
                        ierr)
          call MPI_WAIT(msid,status,ierr)
        if(myidy.ne.(npy-1)) then
          do k = 1, innz
            do i = 1, innx
              x(i,inny-1,k) = x(i,inny-1,k) + recvbuf1(i,k)
            enddo
          enddo
        endif

        endif

        if(npx.gt.1) then
          nsend = innx*(inny-2*jadd)
          call MPI_IRECV(recvbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,left,3,commrow,&
                         msid,ierr)
          do j = 1+jadd, inny-jadd
            do i = 1, innx
              sendbuf2(i,j-jadd) = x(i,j,innz)
            enddo
          enddo
          call MPI_SEND(sendbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,right,3,commrow,&
                        ierr)
          call MPI_WAIT(msid,status,ierr)
          if(myidx.ne.0) then
          do j = 1+jadd, inny-jadd
            do i = 1, innx
              x(i,j,2) = x(i,j,2) + recvbuf2(i,j-jadd)
            enddo
          enddo
          endif

          call MPI_IRECV(recvbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,right,4,commrow,&
                         msid,ierr)
          do j = 1+jadd, inny-jadd
            do i = 1, innx
              sendbuf2(i,j-jadd) = x(i,j,1)
            enddo
          enddo
          call MPI_SEND(sendbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,left,4,commrow,&
                        ierr)
          call MPI_WAIT(msid,status,ierr)
          if(myidx.ne.(npx-1)) then
          do j = 1+jadd, inny-jadd
            do i = 1, innx
              x(i,j,innz-1) = x(i,j,innz-1) + recvbuf2(i,j-jadd)
            enddo
          enddo
          endif
        endif

        call MPI_BARRIER(comm2d,ierr)
        t_guardsum = t_guardsum + elapsedtime_Timer(t0)

        end subroutine guardsum1_Fldmger

        !> exchange grid information between neighboring guard cell
        !> to calculate E from phi.
        subroutine guardexch1_Fldmger(x,innx,inny,innz,grid)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: innx, inny, innz
        double precision, dimension(innx,inny,innz), intent(inout) :: x
        type (Pgrid2d), intent(in) :: grid
        double precision, dimension(innx,innz) :: sendbuf1, recvbuf1
        double precision, dimension(innx,inny) :: sendbuf2, recvbuf2
        integer :: left,right,bottom,top
        integer :: leftx,rightx,bottomy,topy
        integer :: myid,myidx,myidy,comm2d,commcol,commrow,totnp,npx,&
                   npy,msid
        integer, dimension(2) :: tmpcoord
        integer :: i,j,k,ierr
        integer status(MPI_STATUS_SIZE)
        double precision :: t0
        integer :: nsend,jadd,kadd

        call getsize_Pgrid2d(grid,totnp,npy,npx)
        call getpost_Pgrid2d(grid,myid,myidy,myidx)
        call getcomm_Pgrid2d(grid,comm2d,commcol,commrow)

        call MPI_BARRIER(comm2d,ierr)
        call starttime_Timer(t0)

        !define left,right,top,and bottom for cyclic shift purpose.
        if(myidx.eq.0) then
          left = npx - 1
          right = myidx + 1
        else if(myidx.eq.(npx-1)) then
          left = myidx - 1
          right = 0
        else
          left = myidx - 1
          right = myidx + 1
        endif
        if(npx.gt.1) then
          kadd = 1
        else
          kadd = 0
        endif

        if(myidy.eq.0) then
          bottom = npy - 1
          top = myidy + 1
        else if(myidy.eq.(npy-1)) then
          bottom = myidy - 1
          top = 0
        else
          bottom = myidy - 1
          top = myidy + 1
        endif
        if(npy.gt.1) then
          jadd = 1
        else
          jadd = 0
        endif

        if(npy.gt.1) then

        do k = 1+kadd, innz-kadd
          do i = 1, innx
            sendbuf1(i,k-kadd) = x(i,inny-1,k)
          enddo
        enddo
        nsend = (innz-2*kadd)*innx
        call MPI_IRECV(recvbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,bottom,1,&
                      commcol,msid,ierr)
        call MPI_SEND(sendbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,top,1,&
                      commcol,ierr)
        call MPI_WAIT(msid,status,ierr)
        do k = 1+kadd, innz-kadd
          do i = 1, innx
            x(i,1,k) = recvbuf1(i,k-kadd)
          enddo
        enddo

        do k = 1+kadd, innz-kadd
          do i = 1, innx
            sendbuf1(i,k-kadd) = x(i,2,k)
          enddo
        enddo
        call MPI_IRECV(recvbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,top,2,&
                      commcol,msid,ierr)
        call MPI_SEND(sendbuf1(1,1),nsend,MPI_DOUBLE_PRECISION,bottom,2,&
                      commcol,ierr)
        call MPI_WAIT(msid,status,ierr)
        do k = 1+kadd, innz-kadd
          do i = 1, innx
            x(i,inny,k) = recvbuf1(i,k-kadd)
          enddo
        enddo

        endif

        if(npx.gt.1) then

        nsend = innx*(inny-2*jadd)
        do k = 1+jadd, inny-jadd
          do j = 1, innx
            sendbuf2(j,k-jadd) = x(j,k,innz-1)
          enddo
        enddo
        call MPI_IRECV(recvbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,left,3,&
                      commrow,msid,ierr)
        call MPI_SEND(sendbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,right,3,&
                      commrow,ierr)
        call MPI_WAIT(msid,status,ierr)
        do k = 1+jadd, inny-jadd
          do j = 1, innx
            x(j,k,1) = recvbuf2(j,k-jadd)
          enddo
        enddo

        do k = 1+jadd, inny-jadd
          do j = 1, innx
            sendbuf2(j,k-jadd) = x(j,k,2)
          enddo
        enddo
        call MPI_IRECV(recvbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,right,4,&
                      commrow,msid,ierr)
        call MPI_SEND(sendbuf2(1,1),nsend,MPI_DOUBLE_PRECISION,left,4,&
                      commrow,ierr)
        call MPI_WAIT(msid,status,ierr)
        do k = 1+jadd, inny-jadd
          do j = 1, innx
            x(j,k,innz) = recvbuf2(j,k-jadd)
          enddo
        enddo

        endif

        call MPI_BARRIER(comm2d,ierr)
        t_guardexch = t_guardexch + elapsedtime_Timer(t0)

        end subroutine guardexch1_Fldmger

!----------------------------------------------------------------
        !> exchange the information for interpolation
        !> from neighboring guard cells.
        subroutine boundint4_Fldmger(x1,x2,x3,innx,inny,innz,grid)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: innx, inny, innz
        double precision, dimension(innx,inny,innz), intent(inout) :: x1,x2,x3
        type (Pgrid2d), intent(in) :: grid
!        double precision, dimension(innx,innz-2,3) :: sendbuf1, recvbuf1
        double precision, allocatable, dimension(:,:,:) :: sendbuf1, recvbuf1
        double precision, dimension(innx,inny,3) :: sendbuf2, recvbuf2
        integer :: left,right,bottom,top,msid,&
                   ierr,i,j,k,nx,ny
        integer :: myid,myidx,myidy,comm2d,commcol,commrow,totnp,&
                   npx,npy
        integer status(MPI_STATUS_SIZE)
        double precision :: t0
        integer :: nsend,jadd,kadd

        call getpost_Pgrid2d(grid,myid,myidy,myidx)
        call getcomm_Pgrid2d(grid,comm2d,commcol,commrow)
        call getsize_Pgrid2d(grid,totnp,npy,npx)

        call MPI_BARRIER(comm2d,ierr)
        call starttime_Timer(t0)
        if(npy.gt.1) then
          jadd = 1
        else
          jadd = 0
        endif
        if(npx.gt.1) then
          kadd = 1
        else
          kadd = 0
        endif

        top = myidy + 1
        bottom = myidy - 1
        left = myidx - 1
        right = myidx + 1

        ny = inny - 2*jadd

! This could be modified to improve speed
        if(npy.gt.1) then

        nsend = innx*(innz-2*kadd)*3
        allocate(recvbuf1(innx,innz-2*kadd,3))
        allocate(sendbuf1(innx,innz-2*kadd,3))

        if(myidy.ne.0) then
          call MPI_IRECV(recvbuf1(1,1,1),nsend,MPI_DOUBLE_PRECISION,bottom,1,commcol,&
                        msid,ierr)
        endif
        if(myidy.ne.(npy-1)) then
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,1) = x1(i,inny-1,k)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,2) = x2(i,inny-1,k)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,3) = x3(i,inny-1,k)
            enddo
          enddo
          call MPI_SEND(sendbuf1(1,1,1),nsend,MPI_DOUBLE_PRECISION,top,1,&
                        commcol,ierr)
        endif
        if(myidy.ne.0) then
          call MPI_WAIT(msid,status,ierr)
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x1(i,1,k) = recvbuf1(i,k-kadd,1)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x2(i,1,k) = recvbuf1(i,k-kadd,2)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x3(i,1,k) = recvbuf1(i,k-kadd,3)
            enddo
          enddo
        endif

        if(myidy.ne.(npy-1)) then
          call MPI_IRECV(recvbuf1(1,1,1),nsend,MPI_DOUBLE_PRECISION,top,2,commcol,&
                        msid,ierr)
        endif
        if(myidy.ne.0) then
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,1) = x1(i,2,k)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,2) = x2(i,2,k)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              sendbuf1(i,k-kadd,3) = x3(i,2,k)
            enddo
          enddo
          call MPI_SEND(sendbuf1(1,1,1),nsend,MPI_DOUBLE_PRECISION,bottom,2,commcol,&
                        ierr)
        endif
        if(myidy.ne.(npy-1)) then
          call MPI_WAIT(msid,status,ierr)
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x1(i,inny,k) = recvbuf1(i,k-kadd,1)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x2(i,inny,k) = recvbuf1(i,k-kadd,2)
            enddo
          enddo
          do k = 1+kadd, innz-kadd
            do i = 1, innx
              x3(i,inny,k) = recvbuf1(i,k-kadd,3)
            enddo
          enddo
        endif

        deallocate(recvbuf1)
        deallocate(sendbuf1)

        endif

        if(npx.gt.1) then

        nsend = innx*inny*3
        if(myidx.ne.0) then
          call MPI_IRECV(recvbuf2(1,1,1),nsend,MPI_DOUBLE_PRECISION,left,3,commrow,&
                         msid,ierr)
        endif
        if(myidx.ne.(npx-1)) then
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,1) = x1(k,j,innz-1)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,2) = x2(k,j,innz-1)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,3) = x3(k,j,innz-1)
            enddo
          enddo
          call MPI_SEND(sendbuf2(1,1,1),nsend,MPI_DOUBLE_PRECISION,right,3,commrow,&
                        ierr)
        endif
        if(myidx.ne.0) then
          call MPI_WAIT(msid,status,ierr)
          do j = 1, inny
            do k = 1, innx
              x1(k,j,1) = recvbuf2(k,j,1)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              x2(k,j,1) = recvbuf2(k,j,2)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              x3(k,j,1) = recvbuf2(k,j,3)
            enddo
          enddo
        endif

        if(myidx.ne.(npx-1)) then
          call MPI_IRECV(recvbuf2(1,1,1),nsend,MPI_DOUBLE_PRECISION,right,4,commrow,&
                         msid,ierr)
        endif
        if(myidx.ne.0) then
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,1) = x1(k,j,2)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,2) = x2(k,j,2)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              sendbuf2(k,j,3) = x3(k,j,2)
            enddo
          enddo
          call MPI_SEND(sendbuf2(1,1,1),nsend,MPI_DOUBLE_PRECISION,left,4,commrow, &
                        ierr)
        endif
        if(myidx.ne.(npx-1)) then
          call MPI_WAIT(msid,status,ierr)
          do j = 1, inny
            do k = 1, innx
              x1(k,j,innz) = recvbuf2(k,j,1)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              x2(k,j,innz) = recvbuf2(k,j,2)
            enddo
          enddo
          do j = 1, inny
            do k = 1, innx
              x3(k,j,innz) = recvbuf2(k,j,3)
            enddo
          enddo
        endif

        endif

        call MPI_BARRIER(comm2d,ierr)
        t_boundint = t_boundint + elapsedtime_Timer(t0)

        end subroutine boundint4_Fldmger

      end module Fldmgerclass

max_stars_repo_name: elainejut/IMPACT-Z

max_issues_repo_name: elainejut/IMPACT-Z

max_forks_repo_name: elainejut/IMPACT-Z


text:
*
*     ------------------------------------------------------------------
*	M A T M L 4
*     ------------------------------------------------------------------
*
      SUBROUTINE MATML4(C,A,B,NCROW,NCCOL,NAROW,NACOL,
     &                  NBROW,NBCOL,ITRNSP )
C
C MULTIPLY A AND B TO GIVE C
C
C     C = A * B             FOR ITRNSP = 0
C
C     C = A(TRANSPOSED) * B FOR ITRNSP = 1
C
C     C = A * B(TRANSPOSED) FOR ITRNSP = 2
C
C... JEPPE OLSEN, LAST REVISION JULY 24 1987
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NAROW,NACOL),B(NBROW,NBCOL)
      DIMENSION C(NCROW,NCCOL)
C
      NTEST = 0
      IF ( NTEST .NE. 0 ) THEN
        WRITE(6,*)
        WRITE(6,*) ' A AND B MATRIX FROM MATML4 '
        WRITE(6,*)
        CALL WRTMAT(A,NAROW,NACOL,NAROW,NACOL)
        CALL WRTMAT(B,NBROW,NBCOL,NBROW,NBCOL)
        WRITE(6,*)      ' NCROW NCCOL NAROW NACOL NBROW NBCOL '
        WRITE(6,'(6I6)')  NCROW,NCCOL,NAROW,NACOL,NBROW,NBCOL
      END IF
C
      CALL SETVEC(C,0.0D0,NCROW*NCCOL)
C
      IF( ITRNSP .NE. 0 ) GOTO 001
        DO 50 J = 1,NCCOL
          DO 40 K = 1,NBROW
            BKJ = B(K,J)
            DO 30 I = 1, NCROW
              C(I,J) = C(I,J) + A(I,K)*BKJ
  30        CONTINUE
  40      CONTINUE
  50    CONTINUE
C
C
  001 CONTINUE
C
      IF ( ITRNSP .NE. 1 ) GOTO 101
C... C = A(T) * B
         DO 150 J = 1, NCCOL
           DO 140 K = 1, NBROW
             BKJ = B(K,J)
             DO 130 I = 1, NCROW
               C(I,J) = C(I,J) + A(K,I)*BKJ
  130        CONTINUE
  140      CONTINUE
  150    CONTINUE
C
  101 CONTINUE
C
      IF ( ITRNSP .NE. 2 ) GOTO 201
C... C = A*B(T)
        DO 250 J = 1,NCCOL
          DO 240 K = 1,NBCOL
            BJK = B(J,K)
            DO 230 I = 1, NCROW
              C(I,J) = C(I,J) + A(I,K)*BJK
 230        CONTINUE
 240      CONTINUE
 250    CONTINUE
C
C
  201 CONTINUE
C
      IF ( NTEST .NE. 0 ) THEN
        WRITE(6,*)
        WRITE(6,*) ' C MATRIX FROM MATML4 '
        WRITE(6,*)
        CALL WRTMAT(C,NCROW,NCCOL,NCROW,NCCOL)
      END IF
C
      RETURN
      END


text: #redirect Explorit Science Center


! MIT License
!
! Copyright (c) 2016 Anders Steen Christensen
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

subroutine fgenerate_coulomb_matrix(nuclear_charges, coordinates, natoms, nmax, cm)

    implicit none

    double precision, dimension(:), intent(in) :: nuclear_charges
    double precision, dimension(:,:), intent(in) :: coordinates

    integer, intent(in) :: nmax

    double precision, dimension((nmax + 1) * nmax / 2), intent(out):: cm

    double precision, allocatable, dimension(:) :: row_norms
    double precision :: pair_norm
    double precision :: huge_double

    integer, allocatable, dimension(:) :: sorted_atoms

    double precision, allocatable, dimension(:,:) :: pair_distance_matrix

    integer :: i, j, m, n, idx
    integer :: natoms

    if (size(coordinates, dim=1) /= size(nuclear_charges, dim=1)) then
        write(*,*) "ERROR: Coulomb matrix generation"
        write(*,*) size(coordinates, dim=1), "coordinates, but", &
            & size(nuclear_charges, dim=1), "atom_types!"
        stop
    else
        natoms = size(nuclear_charges, dim=1)
    endif

    ! Allocate temporary
    allocate(pair_distance_matrix(natoms,natoms))
    allocate(row_norms(natoms))
    allocate(sorted_atoms(natoms))

    huge_double = huge(row_norms(1))

    ! Calculate row-2-norms and store pair-distances in pair_distance_matrix
    row_norms = 0.0d0

    !$OMP PARALLEL DO PRIVATE(pair_norm) REDUCTION(+:row_norms)
    do i = 1, natoms
        pair_norm = 0.5d0 * nuclear_charges(i) ** 2.4d0
        row_norms(i) = row_norms(i) + pair_norm*pair_norm
        pair_distance_matrix(i, i) = pair_norm
    enddo
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(pair_norm) REDUCTION(+:row_norms)
    do i = 1, natoms
        do j = i+1, natoms
            pair_norm = nuclear_charges(i) * nuclear_charges(j) &
                & / sqrt(sum((coordinates(j,:) - coordinates(i,:))**2))

            pair_distance_matrix(i, j) = pair_norm
            pair_distance_matrix(j, i) = pair_norm
            pair_norm = pair_norm * pair_norm
            row_norms(j) = row_norms(j) + pair_norm
            row_norms(i) = row_norms(i) + pair_norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    !Generate sorted list of atom ids by row_norms - not really (easily) parallelizable
    do i = 1, natoms
        j = minloc(row_norms, dim=1)
        sorted_atoms(natoms - i + 1) = j
        row_norms(j) = huge_double
    enddo

    ! Fill coulomb matrix according to sorted row-norms
    cm = 0.0d0
    !$OMP PARALLEL DO PRIVATE(idx, i, j)
    do m = 1, natoms
        i = sorted_atoms(m)
        idx = (m*m+m)/2 - m
        do n = 1, m
            j = sorted_atoms(n)
            cm(idx+n) = pair_distance_matrix(i, j)
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Clean up
    deallocate(pair_distance_matrix)
    deallocate(row_norms)
    deallocate(sorted_atoms)
end subroutine fgenerate_coulomb_matrix

subroutine fgenerate_unsorted_coulomb_matrix(nuclear_charges, coordinates, natoms, nmax, cm)

    implicit none

    double precision, dimension(:), intent(in) :: nuclear_charges
    double precision, dimension(:,:), intent(in) :: coordinates

    integer, intent(in) :: nmax

    double precision, dimension((nmax + 1) * nmax / 2), intent(out):: cm

    double precision :: pair_norm
    double precision :: huge_double

    double precision, allocatable, dimension(:,:) :: pair_distance_matrix

    integer :: i, j, m, n, idx
    integer :: natoms

    if (size(coordinates, dim=1) /= size(nuclear_charges, dim=1)) then
        write(*,*) "ERROR: Coulomb matrix generation"
        write(*,*) size(coordinates, dim=1), "coordinates, but", &
            & size(nuclear_charges, dim=1), "atom_types!"
        stop
    else
        natoms = size(nuclear_charges, dim=1)
    endif

    ! Allocate temporary
    allocate(pair_distance_matrix(natoms,natoms))

    !$OMP PARALLEL DO PRIVATE(pair_norm)
    do i = 1, natoms
        pair_norm = 0.5d0 * nuclear_charges(i) ** 2.4d0
        pair_distance_matrix(i, i) = pair_norm
    enddo
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(pair_norm)
    do i = 1, natoms
        do j = i+1, natoms
            pair_norm = nuclear_charges(i) * nuclear_charges(j) &
                & / sqrt(sum((coordinates(j,:) - coordinates(i,:))**2))

            pair_distance_matrix(i, j) = pair_norm
            pair_distance_matrix(j, i) = pair_norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Fill coulomb matrix according to sorted row-norms
    cm = 0.0d0
    !$OMP PARALLEL DO PRIVATE(idx)
    do m = 1, natoms
        idx = (m*m+m)/2 - m
        do n = 1, m
            cm(idx+n) = pair_distance_matrix(m, n)
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Clean up
    deallocate(pair_distance_matrix)
end subroutine fgenerate_unsorted_coulomb_matrix

subroutine fgenerate_local_coulomb_matrix(nuclear_charges, coordinates, natoms, nmax, cm)

    implicit none

    double precision, dimension(:), intent(in) :: nuclear_charges
    double precision, dimension(:,:), intent(in) :: coordinates
    integer,intent(in) :: natoms
    integer, intent(in) :: nmax

    double precision, dimension(natoms,(nmax + 1) * nmax / 2), intent(out):: cm
    integer :: idx

    double precision, dimension(natoms) :: row_norms
    double precision :: pair_norm
    double precision :: huge_double

    integer, dimension(natoms) :: sorted_atoms
    integer, dimension(natoms,natoms) :: sorted_atoms_all

    double precision, dimension(natoms,natoms) :: pair_distance_matrix

    integer i, j, m, n, k

    if (size(coordinates, dim=1) /= size(nuclear_charges, dim=1)) then
        write(*,*) "ERROR: Coulomb matrix generation"
        write(*,*) size(coordinates, dim=1), "coordinates, but", &
            & size(nuclear_charges, dim=1), "atom_types!"
        stop
    endif

    huge_double = huge(row_norms(1))

    ! Calculate row-2-norms and store pair-distances in pair_distance_matrix
    row_norms = 0.0d0

    !$OMP PARALLEL DO PRIVATE(pair_norm) REDUCTION(+:row_norms)
    do i = 1, natoms
        pair_norm = 0.5d0 * nuclear_charges(i) ** 2.4d0
        row_norms(i) = row_norms(i) + pair_norm*pair_norm
        pair_distance_matrix(i, i) = pair_norm
    enddo
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(pair_norm) REDUCTION(+:row_norms)
    do i = 1, natoms
        do j = i+1, natoms
            pair_norm = nuclear_charges(i) * nuclear_charges(j) &
                & / sqrt(sum((coordinates(j,:) - coordinates(i,:))**2))

            pair_distance_matrix(i, j) = pair_norm
            pair_distance_matrix(j, i) = pair_norm
            pair_norm = pair_norm * pair_norm
            row_norms(j) = row_norms(j) + pair_norm
            row_norms(i) = row_norms(i) + pair_norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    !Generate sorted list of atom ids by row_norms - not really (easily) parallelizable
    do i = 1, natoms
        j = minloc(row_norms, dim=1)
        sorted_atoms(natoms - i + 1) = j
        row_norms(j) = huge_double
    enddo

    do k = 1, natoms
        sorted_atoms_all(1, k)  = k
        m = 0
        do i = 1, natoms - 1
            if (sorted_atoms(i) == k) then
                m = 1
            endif
            sorted_atoms_all(i+1, k) = sorted_atoms(i + m)
        enddo
    enddo

    ! Fill coulomb matrix according to sorted row-2-norms
    cm = 0.0d0
    !$OMP PARALLEL DO PRIVATE(idx, i, j)
    do k = 1, natoms
        do m = 1, natoms
            i = sorted_atoms_all(m, k)
            idx = (m*m+m)/2 - m
            do n = 1, m
                j = sorted_atoms_all(n, k)
                cm(k, idx+n) = pair_distance_matrix(i, j)
            enddo
        enddo
    enddo
    !$OMP END PARALLEL DO

end subroutine fgenerate_local_coulomb_matrix

subroutine fgenerate_atomic_coulomb_matrix(nuclear_charges, coordinates, natoms, nmax, cutoff, cm)

    implicit none

    double precision, dimension(:), intent(in) :: nuclear_charges
    double precision, dimension(:,:), intent(in) :: coordinates
    integer,intent(in) :: natoms
    integer, intent(in) :: nmax
    double precision, intent(in) :: cutoff

    double precision, dimension(natoms,(nmax + 1) * nmax / 2), intent(out):: cm

    integer :: idx

    double precision :: pair_norm
    double precision :: norm
    double precision :: huge_double

    integer, allocatable, dimension(:,:) :: sorted_atoms_all
    integer, allocatable, dimension(:) :: cutoff_count

    double precision, allocatable, dimension(:,:) :: pair_distance_matrix
    double precision, allocatable, dimension(:,:) :: distance_matrix

    integer i, j, m, n, k

    if (size(coordinates, dim=1) /= size(nuclear_charges, dim=1)) then
        write(*,*) "ERROR: Coulomb matrix generation"
        write(*,*) size(coordinates, dim=1), "coordinates, but", &
            & size(nuclear_charges, dim=1), "atom_types!"
        stop
    endif

    ! Allocate temporary
    allocate(distance_matrix(natoms,natoms))
    allocate(cutoff_count(natoms))

    huge_double = huge(distance_matrix(1,1))
    cutoff_count = 0

    !$OMP PARALLEL DO PRIVATE(norm)
    do i = 1, natoms
        distance_matrix(i, i) = 0.0d0
        do j = i+1, natoms
            norm = sqrt(sum((coordinates(j,:) - coordinates(i,:))**2))
            distance_matrix(i, j) = norm
            distance_matrix(j, i) = norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO REDUCTION(+:cutoff_count)
    do i = 1, natoms
        do j = 1, natoms
            if (distance_matrix(i, j) < cutoff) then
                cutoff_count(i) = cutoff_count(i) + 1
            endif
        enddo
    enddo
    !$OMP END PARALLEL DO

    do i = 1, natoms
        if (cutoff_count(i) > nmax) then
            write(*,*) "ERROR: Coulomb matrix generation"
            write(*,*) nmax, "size set, but", &
                & cutoff_count(i), "size needed!"
            stop
        endif
    enddo

    ! Allocate temporary
    allocate(pair_distance_matrix(natoms,natoms))

    !$OMP PARALLEL DO PRIVATE(pair_norm)
    do i = 1, natoms
        pair_distance_matrix(i, i) = 0.5d0 * nuclear_charges(i) ** 2.4d0
        do j = i+1, natoms
            pair_norm = nuclear_charges(i) * nuclear_charges(j) &
                & / distance_matrix(j, i)

            pair_distance_matrix(i, j) = pair_norm
            pair_distance_matrix(j, i) = pair_norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Allocate temporary
    allocate(sorted_atoms_all(natoms, natoms))

    !$OMP PARALLEL DO PRIVATE(i, j)
    do k = 1, natoms
        do i = 1, natoms
            j = minloc(distance_matrix(:,k), dim=1)
            distance_matrix(j, k) = huge_double
            sorted_atoms_all(i, k) = j
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Clean up
    deallocate(distance_matrix)

    ! Fill coulomb matrix according to sorted distances
    cm = 0.0d0

    !$OMP PARALLEL DO PRIVATE(i, idx, j)
    do k = 1, natoms
        do m = 1, cutoff_count(k)
            i = sorted_atoms_all(m, k)
            idx = (m*m+m)/2 - m
            do n = 1, m
                j = sorted_atoms_all(n, k)
                cm(k, idx+n) = pair_distance_matrix(i, j)
            enddo
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! Clean up
    deallocate(pair_distance_matrix)
    deallocate(sorted_atoms_all)
    deallocate(cutoff_count)

end subroutine fgenerate_atomic_coulomb_matrix

subroutine fgenerate_atomic_unsorted_coulomb_matrix(nuclear_charges, coordinates, natoms, nmax, cutoff, cm)

    implicit none

    double precision, dimension(:), intent(in) :: nuclear_charges
    double precision, dimension(:,:), intent(in) :: coordinates
    integer,intent(in) :: natoms
    integer, intent(in) :: nmax
    double precision, intent(in) :: cutoff

    double precision, dimension(natoms,(nmax + 1) * nmax / 2), intent(out):: cm

    integer :: idx

    double precision :: norm

    double precision, allocatable, dimension(:,:) :: distance_matrix

    integer i, j, m, n, k

    if (size(coordinates, dim=1) /= size(nuclear_charges, dim=1)) then
        write(*,*) "ERROR: Coulomb matrix generation"
        write(*,*) size(coordinates, dim=1), "coordinates, but", &
            & size(nuclear_charges, dim=1), "atom_types!"
        stop
    endif
    write(*,*) 1

    ! Allocate temporary
    allocate(distance_matrix(natoms,natoms))

    !$OMP PARALLEL DO PRIVATE(norm)
    do i = 1, natoms
        distance_matrix(i, i) = 0.0d0
        do j = i+1, natoms
            norm = sqrt(sum((coordinates(j,:) - coordinates(i,:))**2))
            distance_matrix(i, j) = norm
            distance_matrix(j, i) = norm
        enddo
    enddo
    !$OMP END PARALLEL DO

    write(*,*) 2

    ! Fill coulomb matrix
    cm = 0.0d0

    !$OMP PARALLEL DO PRIVATE(idx)
    do k = 1, natoms
        do i = 1, nmax
            idx = (i*i+i)/2 - i
            if (distance_matrix(k,i) < cutoff) then
                do j = 1, i
                    if (distance_matrix(k,j) < cutoff) then
                        cm(k, idx+j) = distance_matrix(i, j)
                    endif
                enddo
            endif
        enddo
    enddo
    !$OMP END PARALLEL DO
    write(*,*) 3

    ! Clean up
    deallocate(distance_matrix)

end subroutine fgenerate_atomic_unsorted_coulomb_matrix


text:
       SUBROUTINE XMOVEX(OUT,IN,IBYTES)
C
C      THIS SUBROUTINE MAY NOT BE NEEDED, ITS WAS IN
C      ASSEMBLER LANGUAGE TO MOVE DATA, IT RAN ABOUT THREE
C      TIMES FASTER THAN A FORTAN DO LOOP, IT WAS USED TO
C      MAKE SURE THE DATA TO BE UNPACKED WAS ON A WORD BOUNDARY,
C      THIS MAY NOT BE NEEDED ON SOME BRANDS OF COMPUTERS.
C
       CHARACTER * 1 OUT(*)
       CHARACTER * 1 IN(*)
C
       INTEGER       IBYTES
C
       DO 100 I = 1,IBYTES
         OUT(I) = IN(I)
  100  CONTINUE
C
       RETURN
       END

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
!! R926 rewind-stmt
!    is  REWIND file-unit-number
!    or  REWIND ( position-spec-list )
!
! Not tested here: file-unit-number and position-spec-list.

integer unitnumber
11 REWIND(unitnumber)
rewind unitnumber
12 rewind 10
rewind (10, iostat=n, err=4)

4 end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf