text:
subroutine condinit

 ! ------------------------------------ !
 !           Initial Condition          !
 !    [ Kelvi-Helmholtz instability ]   !
 ! ------------------------------------ !

 !**************************************************!
 !  MPI PARALLEL VERSION: VALENTINI-FAGANELLO 2009  !
 !   3D PARALLEL VERSION: FAGANELLO 2010            !
 !**************************************************!

 !**************************************************!
 !  FLR corrected MHD equilibria:   CERRI 2011      !
 !                                                  !
 !  for details see:                                !
 !                                                  !
 !  Cerri et al., Phys. Plasmas 20, 112112 (2013)   !
 !                                                  !
 !**************************************************!

use parameter_mod
use dom_distr_mod
use box_mod
use fields_UJ_mod
use fields_DP_mod
use fields_EB_mod
use parallel_mod

IMPLICIT NONE

integer :: ix, iy, iz, ir, iky, ikz

REAL(dp) :: ky, kz, zkk, pph, xx, ab, zz
REAL(dp) :: n0, B0
REAL(dp) :: T0i_para, T0i_perp, T0e_para, T0e_perp
REAL(dp) :: gam_i_perp, gam_i_para, gam_e_perp, gam_e_para
REAL(dp) :: gaminv_i_perp, gam_tilde
REAL(dp) :: beta_i_perp, beta_e_perp
REAL(dp) :: beta_i_perp_tilde, beta_e_perp_tilde
REAL(dp) :: C_0, C_1, C_inf
REAL(dp) :: Den_inf, Dinv_inf!, DeltaDen
REAL(dp) :: pi_para_inf, pi_perp_inf, pe_para_inf, pe_perp_inf, anis_i_0, anis_e_0
REAL(dp) :: Ti_para_inf, Ti_perp_inf, Te_para_inf, Te_perp_inf

REAL(dp),ALLOCATABLE :: equil(:)
REAL(dp),ALLOCATABLE :: f_0(:), fnct(:), equilDen(:), equilB(:)
REAL(dp),ALLOCATABLE :: gx1(:,:,:)
REAL(dp),ALLOCATABLE :: gy1(:), phs(:)
REAL(dp),ALLOCATABLE :: zx(:)
REAL(dp),ALLOCATABLE :: g0(:,:,:), gt(:,:,:)

 CHARACTER*13  filen_DT
 CHARACTER*1 F_1
 CHARACTER*2 F_2
 CHARACTER*3 F_3
 CHARACTER*4 F_4

allocate(equil(nxl))
allocate(f_0(nxl))
allocate(fnct(nxl))
allocate(equilB(nxl))
allocate(equilDen(nxl))
allocate(g0(nxl,ny,nz))
allocate(gx1(nxl,ny,nz))
allocate(gy1(ny))
allocate(phs(ny))
allocate(gt(nx,nyt,nz))
allocate(zx(nx))

SELECT CASE (istart)

 CASE (0)   !*** istart=0 ==> first RUN ***!

  !**********************************************!
  !  MHD one-fluid velocity: shear flow along x  !
  !                                              !
  !  Uy(x) = Aeq * tanh[ (x - Lx/2) / Leq ]      !
  !                                              !
  !  (Leq in d_i units, Aeq in v_A units)        !
  !**********************************************!

   Ux = 0.
   Uz = 0.

   equil = Aeq * dtanh((x(ixlg:ixlg+nxl-1) - xl*0.5d0 ) / Leq)

   if (mpime==0) then
   write(6,*) ' ----------------------------- '
   write(6,*) ' U_0, L_U  = ', Aeq, Leq
   write(6,*) ' ----------------------------- '
   endif

   do iz=1, nz
    do iy=1, ny  !  nyl == ny
     Uy(:,iy,iz) = equil(:)
    enddo
   enddo

  !------------------------------------------------------------!
  ! 2D incompressible perturbations on U: dU = e_z x grad(phi) !
  !------------------------------------------------------------!

   ir = 1

   call random_seed(ir)

   call random_number(phs)

  ! Sum up modes on phi (2D --> no kz)

   g0 = 0.

   do iky = 1, MIN(Ny/8,16)
     ky  = float(iky)
     zkk = 1.0 / ky
     pph = (phs(ir) - 0.5) * 2.0 * pgreco
     ir = ir + 1
     do iz = 1, nz
       do iy = 1, ny
         do ix = 1, nxl
           xx = x(ixlg+ix-1) - xl/2.0
           ab = abs(xx)/Leq
           IF (ab .GE. 9.0)  then
             g0(ix,iy,iz)= 0.0
           ELSE
             g0(ix,iy,iz) = g0(ix,iy,iz) + cos(ky * y(iy) * hyl + pph) * zkk *&
             exp(-1.0 * xx ** 2.0 / Leq ** 2.0)
           ENDIF
         enddo
       enddo
     enddo
   enddo

  ! compute: dU = ez x grad(phi)
  ! (multiplied by ampl, so that max(|dU|) ~ ampl)

   CALL traspdist( g0, gt, 1 )

   do iz = 1, nz
     do iy = 1, nyt
       call  derx_1( gt(:,iy,iz), zx )
       gt(:,iy,iz) = zx
     enddo
   enddo

   CALL traspdist( gx1, gt, -1 )

   Uy = Uy + ampl * gx1

   do iz = 1, nz
     do ix = 1, nxl
       call  dery_1( g0(ix,:,iz), gy1 )
       Ux(ix,:,iz) = - ampl * gy1
     enddo
   enddo

   Uz = 0.0

deallocate(g0)
deallocate(gx1)
deallocate(gy1)
deallocate(equil)
deallocate(phs)
deallocate(gt)
deallocate(zx)

  !**************************************************!
  ! Compute FLR-corrected MHD equilibrium            !
  ! ------------------------------------------------ !
  !  for details see:                                !
  !                                                  !
  !  Cerri et al., Phys. Plasmas 20, 112112 (2013)   !
  !                                                  !
  !**************************************************!

   !-- politropic index of the equilibrium
   !   (default: isothermal initial condition)

   gam_i_perp = 1.0d0
   gam_i_para = 1.0d0

   gam_e_perp = 1.0d0
   gam_e_para = 1.0d0

   gaminv_i_perp = 1.0d0 / gam_i_perp

   !-- equilibrium quantities "at infinity"

   !density
   Den_inf = 1.0d0
   Dinv_inf = 1.0d0 / Den_inf
   !ion pressures and temperatures
   anis_i_0 = 0.0d0     ! ion anisotropy at infinity (A = Pperp/Ppara - 1)
   pi_para_inf = 0.5d0  ! (note: beta_inf = 0.5*p_inf/B0 => 0.5d0 -> beta = 1)
   pi_perp_inf = (1.0d0 + anis_i_0)*pi_para_inf
   Ti_para_inf = Dinv_inf * pi_para_inf
   Ti_perp_inf = Dinv_inf * pi_perp_inf
   !electron pressures and temperatures
   anis_e_0 = 0.0d0     ! electron anisotropy at infinity (A = Pperp/Ppara - 1)
   pe_para_inf = 0.5d0 ! (0.5d0 --> beta = 1)
   pe_perp_inf = (1.0d0 + anis_e_0)*pe_para_inf
   Te_para_inf = Dinv_inf * pe_para_inf
   Te_perp_inf = Dinv_inf * pe_perp_inf
   !perpendicular plasma betas
   beta_i_perp = (2.0d0 * pi_perp_inf) / B00    ! ions
   beta_e_perp = (2.0d0 * pe_perp_inf) / B00    ! elettrons

   ! useful quantities: "reduced" betas and gamma
   ! [see: Cerri et al., PoP 20, 112112 (2013)]
   gam_tilde = (gam_e_perp - gam_i_perp)*gaminv_i_perp
   beta_i_perp_tilde = beta_i_perp / (1.0d0 + beta_i_perp + beta_e_perp)
   beta_e_perp_tilde = beta_e_perp / (1.0d0 + beta_i_perp + beta_e_perp)

   if (mpime==0) then
   write(6,*) ' ------------------------------------------------------------------ '
   write(6,*) ' gam_i_perp,  gam_i_para,  gam_e_perp,  gam_e_para,  gam_tilde  '
   write(6,*) ' ',gam_i_perp, gam_i_para, gam_e_perp, gam_e_para, gam_tilde
   write(6,*) ' Den_inf,  pi_para_inf,  pi_perp_inf,  pe_para_inf,  pe_perp_inf  '
   write(6,*) ' ',Den_inf, pi_para_inf, pi_perp_inf, pe_para_inf, pe_perp_inf
   write(6,*) ' Ti_para_inf,  Ti_perp_inf,  Te_para_inf,  Te_perp_inf  '
   write(6,*) ' ',Ti_para_inf, Ti_perp_inf, Te_para_inf, Te_perp_inf
   write(6,*) ' beta_i_perp,  beta_e_perp,  beta_i_perp_tilde,  beta_e_perp_tilde  '
   write(6,*) ' ',beta_i_perp, beta_e_perp, beta_i_perp_tilde, beta_e_perp_tilde
   write(6,*) ' ------------------------------------------------------------------ '
   endif

   if (mpime==0) then
   write(6,*) ' ----------------------------- '
   write(6,*) ' Delta_N  = ', DeltaDen
   write(6,*) ' ----------------------------- '
   endif

  ! compute FLR-correction function for equivalent MHD profiles
  ! [see: Cerri et al., PoP 20, 112112 (2013)]

  if (flr_on .eq. 1) then
    C_0 = (0.5d0/dsqrt(B00))*(Aeq / Leq)*beta_i_perp_tilde
    equilDen(:) = 1.0d0 - 0.5d0*DeltaDen*Dinv_inf*(1.0d0 - dtanh((x(ixlg:ixlg+nxl-1) - xl*0.5d0 ) / Leq) )
    f_0(:) = 1.0d0 / (1.0d0 - C_0*equilDen(:)*(1.0d0 - (dtanh((x(ixlg:ixlg+nxl-1) - xl*0.5d0 ) / Leq) )**2.0)  )
  else
    f_0(:) = 1.0d0
  endif

  fnct(:) = equilDen(:) * f_0(:)

  !*** Compute actual equilibrium profiles ***!

   !density
   do iz=1, nz
     do iy=1, ny  !  nyl == ny
      Den(:,iy,iz) = Den_inf * (fnct(:)**gaminv_i_perp)
     enddo
   enddo

   Dinv = 1.0d0 / Den

   !gyrotropic ion pressures
   do iz = 1, nz
     do iy = 1, ny  ! nyl == ny
       pi_perp(:,iy,iz) = pi_perp_inf * fnct(:)
     enddo
   enddo
   do iz=1, nz
     do iy=1, ny  ! nyl == ny
       pi_para(:,iy,iz) = pi_para_inf * (fnct(:)**(gam_i_para*gaminv_i_perp))
     enddo
   enddo

   !gyrotropic electron pressures
   do iz = 1, nz
     do iy = 1, ny  ! nyl == ny
       pe_perp(:,iy,iz) = pe_perp_inf * (fnct(:)**(gam_e_perp*gaminv_i_perp))
     enddo
   enddo
   do iz=1, nz
     do iy=1, ny  ! nyl == ny
       pe_para(:,iy,iz) = pe_para_inf * (fnct(:)**(gam_e_para*gaminv_i_perp))
     enddo
   enddo

  !gyrotropic temperatures: simply T = p/n
  Ti_perp = Dinv * pi_perp
  Ti_para = Dinv * pi_para
  Te_perp = Dinv * pe_perp
  Te_para = Dinv * pe_para

 !*** equilibrium magnetic field

  !make sure that B is in (y,z) plane
  !(i.e., not perpendicular to the flow)
  Bx = 0.d0

  !FLR-corrected equilibrium profile
  equilB(:) = (1.0d0+beta_i_perp+beta_e_perp)*f_0(:) - (beta_i_perp+beta_e_perp)*fnct(:)
  do iz=1, nz
    do iy=1, ny  !  nyl == ny
      do ix=1, nxl
        Bz(ix,iy,iz) = dsqrt(B00*equilB(ix)) * dcos(angolo)
      enddo
    enddo
  enddo
  do iz=1, nz
    do iy=1, ny  !  nyl == ny
      do ix=1, nxl
        By(ix,iy,iz) = dsqrt(B00*equilB(ix)) * dsin(angolo)
      enddo
    enddo
  enddo

deallocate(f_0)
deallocate(fnct)
deallocate(equilB)
deallocate(equilDen)

 !*****************************************************************!
 ! Compute: passive tracer (defined on the initial velocity shear) !
 !   [ note: this is problem dependent, adjust at will ]           !
 !*****************************************************************!

  do iz = 1, nz
    do ix = 1, nxl
      Tracciante(ix,:,iz) = 0.6 + 0.4 * dtanh(( x(ix + ixlg - 1)-xl*0.5d0) / Leq)
    enddo
  enddo

 !***********************************!
 ! Compute: initial current density  !
 !***********************************!

  call corrente

 !********************************************************!
 ! Compute: initial fluid velocities (ions and electrons) !
 !********************************************************!

  call u_ei

 !*********************************!
 ! Compute: initial electric field !
 !*********************************!

  call ohm

 !**************************!
 ! Compute: inital momentum !
 !**************************!

  nU_x = Den * Ux
  nU_y = Den * Uy
  nU_z = Den * Uz

 CASE(+1)   !istart = 1 -> Restart from existing RUN using the same timestep dt

  !***********!
  !  RESTART  !
  ! (same dt) !
  !***********!

  l_rst_me=l_rst+mpime

 if ( mpime.lt.10 ) then
   write( F_1, '(1i1)' ) mpime
   filen_DT  = 'DATA_TF_000'//F_1
 endif

 if ( mpime.ge.10.and.mpime.le.99 ) then
   write( F_2, '(1i2)' ) mpime
   filen_DT  = 'DATA_TF_00'//F_2
 endif

 if ( mpime.gt.99 ) then
   write( F_3, '(1i3)' ) mpime
   filen_DT  = 'DATA_TF_0'//F_3
 endif

 if ( (mpime.gt.999) .and. (mpime.le.9999) ) then
   write( F_4, '(1i4)' ) mpime
   filen_DT  = 'DATA_TF_'//F_4
 endif

 OPEN(l_rst_me,status='unknown', form='unformatted',file=filen_DT)

 READ(l_rst_me) ioutt, ioutx
 READ(l_rst_me) tempo, tt_last, tx_last
 READ(l_rst_me) Ti_para, Ti_perp, Te_para, Te_perp
 READ(l_rst_me) Uex, Uey, Uez, Uix, Uiy, Uiz, Ux, Uy, Uz, Jx, Jy, Jz
 READ(l_rst_me) Den, Dinv, pe_para, pe_perp, pi_para, pi_perp, nU_x, nU_y, nU_z
 READ(l_rst_me) Ex, Ey, Ez, Bx, By, Bz, Tracciante!, Tracciante_e

 CLOSE(l_rst_me)

 CASE(+2)   !istart = 2 -> Restart from existing RUN, but using a different timestep dt

  !****************!
  !     RESTART    !
  ! (different dt) !
  !****************!

  l_rst_me=l_rst+mpime

 if ( mpime.lt.10 ) then
   write( F_1, '(1i1)' ) mpime
   filen_DT  = 'DATA_TF_000'//F_1
 endif

 if ( (mpime.ge.10) .and. (mpime.le.99) ) then
   write( F_2, '(1i2)' ) mpime
   filen_DT  = 'DATA_TF_00'//F_2
 endif

 if ( (mpime.gt.99) .and. (mpime.le.999) ) then
   write( F_3, '(1i3)' ) mpime
   filen_DT  = 'DATA_TF_0'//F_3
 endif

 if ( (mpime.gt.999) .and. (mpime.le.9999) ) then
   write( F_4, '(1i4)' ) mpime
   filen_DT  = 'DATA_TF_'//F_4
 endif

 OPEN(l_rst_me,status='unknown', form='unformatted',file=filen_DT)

 READ(l_rst_me) ioutt, ioutx
 READ(l_rst_me) tempo, tt_last, tx_last
 READ(l_rst_me) Ti_para, Ti_perp, Te_para, Te_perp
 READ(l_rst_me) Uex, Uey, Uez, Uix, Uiy, Uiz, Ux, Uy, Uz, Jx, Jy, Jz
 READ(l_rst_me) Den, Dinv, pe_para, pe_perp, pi_para, pi_perp, nU_x, nU_y, nU_z
 READ(l_rst_me) Ex, Ey, Ez, Bx, By, Bz, Tracciante

 CLOSE(l_rst_me)

END SELECT

 !*************************************!
 ! Compute: initial gyroviscous tensor !
 !*************************************!

  if (flr_on .eq. 1) then
    write(6,*) ''
    write(6,*) ' IONS FLR: ON! '
    call flr_i
  else !if (flr_on .eq. 0) then
    write(6,*) ''
    write(6,*) ' NO IONS FLR: G_i set to zero! '
    Gi_xx = 0.0d0
    Gi_xy = 0.0d0
    Gi_xz = 0.0d0
    Gi_yz = 0.0d0
  !else
  !  write(6,*) ''
  !  write(6,*) '*****************************************'
  !  write(6,*) '***          !!! WARNING !!!          ***'
  !  write(6,*) '*****************************************'
  !  write(6,*) '*  invalid value for FLR_ON parameter!  *'
  !  write(6,*) '* ------------------------------------- *'
  !  write(6,*) '*    FLR_ON must have 1 or 0 value...   *'
  !  write(6,*) '*****************************************'
  endif

 if (mpime==0) then
 write(6,*) ' --------------------------------------------- '
 write(6,*) 'SUBROUTINE CONDINIT, B00, angolo, Lx, Ly, Lz :'
 write(6,805) B00, angolo,  x(nx)-x(1), y(ny) - y(1), z(nz) - z(1)
 write(6,*) ' --------------------------------------------- '
 endif

622  format(1x, 2(2x, 1i4), 2(2x, 1e11.4))
822  format(1x, 2(2x, 1f5.1), 2(2x, 1e11.4))
802  format(1x, 2(2x, 1e11.4))
803  format(1x, 3(2x, 1e11.4))
805  format(1x, 5(2x, 1e11.4))

end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
! ###################################################################
! Copyright (c) 2013-2016, Marc De Graef/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:EMKAM.f90
!--------------------------------------------------------------------------
!
! PROGRAM: EMKAM
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief create a tiff file with an OSM map
!
!> @date 07/29/16 MDG 1.0 original
!--------------------------------------------------------------------------
program EMKAM

use local
use typedefs
use NameListTypedefs
use NameListHandlers
use files
use constants
use dictmod
use io
use error
use ECPmod
use EBSDiomod
use EBSDmod
use EBSDDImod
use HDF5
use HDFsupport
use TIFF_f90
use rotations
use stringconstants

IMPLICIT NONE

character(fnlen)                        :: nmldeffile, progname, progdesc
type(KAMNameListType)                   :: enl
type(EBSDIndexingNameListType)          :: ebsdnl

logical                                 :: stat, readonly, noindex
integer(kind=irg)                       :: hdferr, nlines, FZcnt, Nexp, nnm, nnk, Pmdims, i, j, k, olabel, Nd, Ne, ipar(10), &
                                           ipar2(6), pgnum, ipat, ipf_wd, ipf_ht, idims2(2)
character(fnlen)                        :: groupname, dataset, dpfile, energyfile, masterfile, efile, fname
integer(HSIZE_T)                        :: dims2(2)
!type(dicttype),pointer                  :: dict
type(dicttype)                          :: dict
real(kind=sgl)                          :: testeu(3), eu1(3), eu2(3), ro1(4), ro2(4), s, da

character(fnlen),allocatable            :: stringarray(:)
integer(kind=irg),allocatable           :: tmi(:,:), tmitmp(:,:), indexmain(:,:)
real(kind=sgl),allocatable              :: kam(:,:), eulers(:,:), Eulerstmp(:,:), Eulervals(:,:), avEuler(:,:), &
                                           dplist(:,:), dplisttmp(:,:)

type(HDFobjectStackType),pointer        :: HDF_head

nullify(HDF_head)

nmldeffile = 'EMKAM.nml'
progname = 'EMKAM.f90'
progdesc = 'Generate a nearest-neighbor Kernel Average Misorientation (KAM) map'

! print some information
call EMsoft(progname, progdesc)

! deal with the command line arguments, if any
call Interpret_Program_Arguments(nmldeffile,1,(/ 84 /), progname)

! deal with the namelist stuff
call GetKAMNameList(nmldeffile,enl)

!====================================
! read the relevant fields from the dot product HDF5 file

! open the fortran HDF interface
call h5open_EMsoft(hdferr)

dpfile = trim(EMsoft_getEMdatapathname())//trim(enl%dotproductfile)
dpfile = EMsoft_toNativePath(dpfile)

! is this a proper HDF5 file ?
call h5fis_hdf5_f(trim(dpfile), stat, hdferr)

if (stat) then
! open the dot product file
  readonly = .TRUE.
  hdferr =  HDF_openFile(dpfile, HDF_head, readonly)

! get the energyfile and masterfile parameters from NMLParameters
groupname = SC_NMLparameters
    hdferr = HDF_openGroup(groupname, HDF_head)
groupname = SC_EBSDIndexingNameListType
    hdferr = HDF_openGroup(groupname, HDF_head)

dataset = SC_energyfile
    call HDF_readDatasetStringArray(dataset, nlines, HDF_head, hdferr, stringarray)
    energyfile = trim(stringarray(1))
    deallocate(stringarray)

dataset = SC_nnk
    call HDF_readDatasetInteger(dataset, HDF_head, hdferr, nnk)

dataset = SC_ipfwd
    call HDF_readDatasetInteger(dataset, HDF_head, hdferr, ipf_wd)
dataset = SC_ipfht
    call HDF_readDatasetInteger(dataset, HDF_head, hdferr, ipf_ht)

! and close the NMLparameters group
    call HDF_pop(HDF_head)
    call HDF_pop(HDF_head)

! open the Scan 1/EBSD/Data group
    groupname = SC_Scan1
    hdferr = HDF_openGroup(groupname, HDF_head)
groupname = SC_EBSD
    hdferr = HDF_openGroup(groupname, HDF_head)
groupname = SC_Data
    hdferr = HDF_openGroup(groupname, HDF_head)

! integers
dataset = SC_FZcnt
    call HDF_readDatasetInteger(dataset, HDF_head, hdferr, FZcnt)

dataset = SC_NumExptPatterns
    call HDF_readDatasetInteger(dataset, HDF_head, hdferr, Nexp)

dataset = SC_TopDotProductList
    call HDF_readDatasetFloatArray2D(dataset, dims2, HDF_head, hdferr, dplisttmp)
    allocate(dplist(dims2(1),Nexp))
    do i=1,Nexp
      dplist(1:dims2(1),i) = dplisttmp(1:dims2(1),i)
    end do
    deallocate(dplisttmp)
    nnm = dims2(1)

dataset = SC_EulerAngles
    call HDF_readDatasetFloatArray2D(dataset, dims2, HDF_head, hdferr, Eulerstmp)
    allocate(Eulervals(3,FZcnt))
    do i=1,FZcnt
      Eulervals(1:3,i) = Eulerstmp(1:3,i)
    end do
    deallocate(Eulerstmp)
    Eulervals = Eulervals * sngl(cPi)/180.0

! arrays
dataset = SC_TopMatchIndices
    call HDF_readDatasetIntegerArray2D(dataset, dims2, HDF_head, hdferr, tmitmp)
    allocate(tmi(nnk,Nexp))
    do i=1,Nexp
      tmi(1:nnk,i) = tmitmp(1:nnk,i)
    end do
    idims2 = (/ nnk, Nexp /)
    deallocate(tmitmp)
  call HDF_pop(HDF_head,.TRUE.)

! close the fortran HDF interface
  call h5close_EMsoft(hdferr)

  call Message('dot product HDF5 file read')
end if

!===================================
! we will also need some of the crystallographic data, so that requires
! extracting the xtalname from the energyfile
  efile = trim(EMsoft_getEMdatapathname())//trim(energyfile)
  efile = EMsoft_toNativePath(efile)

! first, we need to check whether or not the input file is of the HDF5 format type; if
! it is, we read it accordingly, otherwise we use the old binary format.
!
  call h5fis_hdf5_f(efile, stat, hdferr)

  if (stat) then
! open the fortran HDF interface
    call h5open_EMsoft(hdferr)

    nullify(HDF_head)

! open the MC file using the default properties.
    readonly = .TRUE.
    hdferr =  HDF_openFile(efile, HDF_head, readonly)

! open the namelist group
groupname = SC_NMLparameters
    hdferr = HDF_openGroup(groupname, HDF_head)
groupname = SC_MCCLNameList
    hdferr = HDF_openGroup(groupname, HDF_head)

! read all the necessary variables from the namelist group
dataset = SC_xtalname
    call HDF_readDatasetStringArray(dataset, nlines, HDF_head, hdferr, stringarray)
    ebsdnl%MCxtalname = trim(stringarray(1))
    deallocate(stringarray)

    call HDF_pop(HDF_head,.TRUE.)

! close the fortran HDF interface
    call h5close_EMsoft(hdferr)

    call Message('energy file read')
  else
    efile = 'File '//trim(efile)//' is not an HDF5 file'
    call FatalError('EMAverageOrient',efile)
  end if
  pgnum = GetPointGroup(ebsdnl%MCxtalname)

! and next we compute the KAM map (kam)
allocate(kam(ipf_wd,ipf_ht),eulers(3,Nexp))

! do we need to do an orientation average first ?
if (enl%orav.ne.0) then
  ipar2(1) = pgnum
  ipar2(2) = FZcnt
  ipar2(3) = Nexp
  ipar2(4) = nnk
  ipar2(5) = Nexp*ceiling(float(ipf_wd*ipf_ht)/float(Nexp))
  ipar2(6) = enl%orav
  call Message('Computing orientation averages ... ')
  call EBSDgetAverageOrientations(ipar2, Eulervals, tmi, dplist, eulers)
  eulers = eulers*sngl(cPi)/180.0
else
  do i=1,Nexp
    eulers(1:3,i) = Eulervals(1:3,tmi(1,i))
  end do
end if

! compute the Kernel Average Misorientation map
dict%Num_of_init = 3
dict%Num_of_iterations = 30
dict%pgnum = pgnum
call DI_Init(dict,'nil')

call Message('Computing KAM map... ')
call EBSDgetKAMMap(Nexp, eulers, ipf_wd, ipf_ht, dict, kam)
kam = kam*180.0/sngl(cPi)

write (*,*) 'KAM range = ',minval(kam), maxval(kam)

where (kam.gt.enl%kamcutoff) kam = enl%kamcutoff
where (kam.lt.0.0) kam = 0.0
kam = kam/maxval(kam)
kam = kam*255.0

TIFF_filename = trim(EMsoft_getEMdatapathname())//trim(enl%kamtiff)
TIFF_filename = EMsoft_toNativePath(TIFF_filename)

TIFF_nx = ipf_wd
TIFF_ny = ipf_ht
! allocate memory for image
allocate(TIFF_image(0:TIFF_nx-1,0:TIFF_ny-1))
! fill the image with whatever data you have (between 0 and 255)
 do i=0,TIFF_nx-1
  do j=0,TIFF_ny-1
   TIFF_image(i,j) = kam(i+1,ipf_ht-j)
  end do
 end do
! create the file
 call TIFF_Write_File

call Message('KAM map written to '//trim(TIFF_filename))

end program EMKAM

max_stars_repo_name: imikejackson/EMsoft

max_issues_repo_name: imikejackson/EMsoft

max_forks_repo_name: imikejackson/EMsoft


!======================================================================
!======================================================================
! Utility Program for Listing, Deleting, Relabeling Labeled Solutions
!                        in AUTO97 Data Files
!======================================================================
!======================================================================
!
PROGRAM KEEPSP
  USE UTILITY
  LOGICAL, EXTERNAL :: DELETEFN
  CALL KEEPMAIN(DELETEFN)
END PROGRAM KEEPSP

LOGICAL FUNCTION DELETEFN(I,ITP)
  INTEGER I,ITP
  DELETEFN = (ITP==4 .OR. ITP==-4)
END FUNCTION DELETEFN


text:
Ron Pina is a former elder at Grace Valley Christian Center and has since left the church.  He has written an essay http://daviswiki.org/Grace_Valley_Christian_Center/Talk?actiondiff&version230&version129 Regarding Qualifications For Ministry.

20080514 19:14:50 nbsp Welcome to the Davis Wiki, Ron!  May the Lord bless you and your family as you live in HIS freedom from day to day! Users/RichLindvall

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


PROGRAM main

  USE mo_kind,        ONLY: dp, sp, i4
  USE mo_integrate,   ONLY: int_regular
  use mo_ansi_colors, only: color, c_red, c_green

  IMPLICIT NONE

  INTEGER(i4), PARAMETER  :: nn = 9
  REAL(dp), DIMENSION(nn) :: dat
  REAL(sp), DIMENSION(nn) :: sat
  INTEGER(i4)             :: ii
  real(sp)                :: stmp
  real(dp)                :: dtmp

  LOGICAL :: isgood

  Write(*,*) ''
  Write(*,*) 'Test mo_integrate.f90'

  ! Double precision
  forall(ii=1:nn) dat(ii) = ii
  isgood = .true.
  dtmp   = int_regular(dat)
  isgood = isgood .and. (nint(dtmp) == 40)

  if (isgood) then
     write(*,*) 'mo_integrate double precision ', color('o.k.', c_green)
  else
     write(*,*) 'mo_integrate double precision ', color('failed!', c_red)
  endif

  ! Double precision
  forall(ii=1:nn) sat(ii) = ii
  isgood = .true.
  stmp   = int_regular(sat)
  isgood = isgood .and. (nint(stmp) == 40)

  if (isgood) then
     write(*,*) 'mo_integrate single precision ', color('o.k.', c_green)
  else
     write(*,*) 'mo_integrate single precision ', color('failed!', c_red)
  endif

END PROGRAM main


!
!     -------------------------------------------------------------
!      D L S A 5
!     -------------------------------------------------------------
!                                                                  *
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville             October 1996   *
!
      SUBROUTINE DLSA5(K, JA1, KA, IRE, IAT, REC)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONSTS_C
      USE MEDEFN_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  11:15:35  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE ixjtik_I
      USE sixj_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: K
      INTEGER , INTENT(IN) :: JA1
      INTEGER  :: KA
      INTEGER , INTENT(IN) :: IRE
      INTEGER , INTENT(OUT) :: IAT
      REAL(DOUBLE) , INTENT(OUT) :: REC
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: K1, ITI1, ITI1S, ITI, ITIS, JI
      REAL(DOUBLE) :: A3
!-----------------------------------------------
      REC = ZERO
      K1 = IHSH + IHSH - 1
      ITI1 = J1QN1(K1,K) - 1
      ITI1S = J1QN2(K1,K) - 1
      K1 = K1 - 1
      IF (JA1 == IHSH) THEN
         ITI = J1QN1(IHSH,K) - 1
         ITIS = J1QN2(IHSH,K) - 1
         JI = J1QN1(K1,K) - 1
      ELSE
         JI = J1QN1(IHSH,K) - 1
         ITI = J1QN1(K1,K) - 1
         ITIS = J1QN2(K1,K) - 1
      ENDIF
      IF (IRE == 0) THEN
         IF (IXJTIK(KA,ITIS,ITI,JI,ITI1,ITI1S) /= 0) IAT = 1
      ELSE
         CALL SIXJ (KA, ITIS, ITI, JI, ITI1, ITI1S, 0, A3)
         REC = A3*DSQRT(DBLE((ITI + 1)*(ITI1S + 1)))
         IF (MOD(KA + JI + ITIS + ITI1,4) /= 0) REC = -REC
         IAT = 1
         IF (JA1 == IHSH) RETURN
         IF (MOD(ITI + ITIS - ITI1S - ITI1 + 2*JI,4) /= 0) REC = -REC
      ENDIF
      RETURN
      END SUBROUTINE DLSA5


text:
!*robodoc*m* CGPACK/cgca_m2mpiio
!  NAME
!    cgca_m2mpiio
!  SYNOPSIS

!$Id: cgca_m2mpiio.f90 380 2017-03-22 11:03:09Z mexas $

module cgca_m2mpiio

!  AUTHOR
!    Anton Shterenlikht
!  COPYRIGHT
!    See LICENSE
!  DESCRIPTION
!    Module with MPI/IO related routines
!  CONTAINS
!    Subroutines: cgca_pswci2
!  USES
!    cgca_m1co
!  USED BY
!    cgca
!  SOURCE

use cgca_m1co
use mpi
implicit none
private
public :: cgca_pswci2

contains

!*roboend*

!*robodoc*s* cgca_m2mpiio/cgca_pswci2
!  NAME
!    cgca_pswci2
!  SYNOPSIS

subroutine cgca_pswci2( coarray, stype, fname )

!  INPUTS

    integer( kind=iarr ), allocatable, intent( in ) ::                 &
      coarray(:,:,:,:)[:,:,:]
    integer( kind=idef ),intent( in ) :: stype
    character( len=* ), intent( in ) :: fname

!  OUTPUTS
!    None
!  SIDE EFFECTS
!    A single binary file is created using MPI/IO
!    with contents of coarray.
!  AUTHOR
!    Anton Shterenlikht, adapted from the code
!    written by David Henty, EPCC
!  COPYRIGHT
!   Note that this routine has special Copyright conditions.
!
!    !----------------------------------------------------------------------------!
!    !                                                                            !
!    !  MPI-IO routine for Fortran Coarrays                                       !
!    !                                                                            !
!    !  David Henty, EPCC; d.henty@epcc.ed.ac.uk                                  !
!    !                                                                            !
!    !  Copyright 2013 the University of Edinburgh                                !
!    !                                                                            !
!    !  Licensed under the Apache License, Version 2.0 (the "License");           !
!    !  you may not use this file except in compliance with the License.          !
!    !  You may obtain a copy of the License at                                   !
!    !                                                                            !
!    !      http://www.apache.org/licenses/LICENSE-2.0                            !
!    !                                                                            !
!    !  Unless required by applicable law or agreed to in writing, software       !
!    !  distributed under the License is distributed on an "AS IS" BASIS,         !
!    !  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  !
!    !  See the License for the specific language governing permissions and       !
!    !  limitations under the License.                                            !
!    !                                                                            !
!    !----------------------------------------------------------------------------!
!
!  DESCRIPTION
!    Parallel Stream Write Coarray of Integers:
!    - coarray - what array to dump
!    - stype - what cell state type to dump
!    - fname - what file name to use
!  NOTES
!    All images must call this routine!
!
!    MPI must be initialised prior to calling this routine,
!    most probably in the main program.
!    Likewise MPI must be terminated only when no further MPI
!    routines can be called. This will most likely be in the
!    main program.
!    There are some assumptions about the shape of the passed array.
!
!    The default integer is assumed for the array at present!
!  USES
!    cgca_m1co, MPI library
!  USED BY
!    none, end user.
!  SOURCE

integer, parameter :: totdim = 4, arrdim = totdim-1, coardim = 3

integer :: img, nimgs, comm, ierr=0, rank=0, mpisize=0, filetype,      &
           mpi_subarray, fh, funit
integer, dimension(totdim) :: asizehal
integer, dimension(arrdim) :: arrsize, arstart, artsize
integer, dimension(coardim) :: coarsize, copos
integer( kind=MPI_OFFSET_KIND ) :: disp = 0
integer, dimension(MPI_STATUS_SIZE) :: mpistat

character( len=80 ) :: iomsg

  img = this_image()
nimgs = num_images()

asizehal(:) = shape( coarray )
   copos(:) = this_image( coarray )

! Subtract halos
 arrsize(:) = asizehal(1:arrdim) - 2
coarsize(:) = ucobound(coarray) - lcobound(coarray) + 1

! Does the array fit exactly?
if ( product( coarsize ) .ne. nimgs) then
 write(*,*) 'ERROR: cgca_m2mpiio/cgca_pswci2: non-conforming coarray'
 error stop
end if

comm = MPI_COMM_WORLD
call MPI_Comm_size( comm, mpisize, ierr )
call MPI_Comm_rank( comm, rank, ierr )

! Sanity check
if ( mpisize .ne. nimgs .or. rank .ne. img-1 ) then
  write(*,*) 'ERROR: cgca_m2mpiio/cgca_pswci2: MPI/coarray mismatch'
  error stop
end if

! Define filetype for this process, ie what portion of the global array
! this process owns. Starting positions use C-indexing
! (ie counting from 0).
artsize(:) = arrsize(:) * coarsize(:)
arstart(:) = arrsize(:) * (copos(:)-1)

! debug
!write (*,*) "image",img, "asizehal", asizehal, "copos", copos,         &
!  "arrsize", arrsize, "coarsize", coarsize,                            &
!  "artsize", artsize, "arstart", arstart, "stype", stype

call MPI_Type_create_subarray( arrdim, artsize, arrsize, arstart,      &
 MPI_ORDER_FORTRAN, MPI_INTEGER, filetype, ierr )
if ( ierr .ne. MPI_SUCCESS ) then
 write (*,'(a,i0)') 'ERROR: cgca_m2mpiio/cgca_pswci2:' //              &
   " MPI_type_create_subarray filetype: rank: ", rank
 error stop
end if

call MPI_Type_commit( filetype, ierr )

! Define subarray for this process, ie what portion of the local array
! is to be written (excludes halos); starting positions use C-indexing.

arstart(:) = 1

call MPI_Type_create_subarray( arrdim, asizehal, arrsize, arstart,     &
 MPI_ORDER_FORTRAN, MPI_INTEGER, mpi_subarray, ierr)
if ( ierr .ne. MPI_SUCCESS ) then
 write (*,'(a,i0)') 'ERROR: cgca_m2mpiio/cgca_pswci2:' //              &
   " MPI_type_create_subarray mpi_subarray: rank: ", rank
 error stop
end if

call MPI_Type_commit(mpi_subarray, ierr)

!  "Striping information cannot be changed on an existing
! file, so to set the stripe count (and stripe size) for the amount of
! parallelism you want to achieve, the file must be deleted if it exists."
! From: Cray Getting Started on MPI I/O manual, S-2490-40 - Dec 2009:
! http://docs.cray.com/books/S-2490-40/

!if ( rank .eq. 0 ) then
!  call MPI_File_delete( fname, MPI_INFO_NULL, ierr )
!  if ( ierr .ne. MPI_SUCCESS )                                         &
!    error stop "ERROR: cgca_pswci2: MPI_file_delete: rank 0"
!end if

! All ranks wait till rank 0 deletes the file
!call MPI_Barrier( comm, ierr )
!if ( ierr .ne. MPI_SUCCESS ) then
!  write (*,*) 'ERROR: cgca_pswci2: MPI_file_open: rank ', rank
!  error stop
!end if

! Overwriting MPI/IO files does not involve erasing the file first.
! So if the old file was bigger, the new smaller file will still
! be sized on disk as the old file, with only a part of it overwritten
! with new data. That would be bad. To avoid this
! image 1 removes all previous contents of this file, if it exists.
rm: if ( img .eq. 1 ) then

  ! this should not be necessary, but Cray issues a caution otherwise
  funit = 0

  open( newunit=funit, file=fname, status="replace", iostat=ierr,      &
        iomsg=iomsg )
  if ( ierr .ne. 0 ) then
    write (*,*) "ERROR: cgca_m2mpiio/cgca_pswci2: open( fname )," //    &
      " iostat=", ierr, "iomsg:", iomsg
    error stop
  end if

  write( funit, * , iostat=ierr) ""
  if ( ierr .ne. 0 ) then
    write (*,*) "ERROR: cgca_m2mpiio/cgca_pswci2: write( fname )," //   &
      " iostat=", ierr
    error stop
  end if

!  flush( funit, iostat=ierr )
!  if ( ierr .ne. 0 ) then
!    write (*,*) "ERROR: cgca_m2mpiio: flush( fname ), iostat=",ierr
!    error stop
!  end if

  close( funit, iostat=ierr )
  if ( ierr .ne. 0 ) then
    write (*,*) "ERROR: cgca_m2mpiio/cgca_pswci2: close( fname )," //   &
      " iostat=", ierr
    error stop
  end if

end if rm

! all images wait till image 1 erases the previous file
sync all

!  Open the file for writing only and attach to file handle fh
!  No IO hints are passed since MPI_INFO_NULL is specified
call MPI_File_open( comm, fname, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE),&
                    MPI_INFO_NULL, fh, ierr )
if ( ierr .ne. MPI_SUCCESS ) then
  write (*,'(a,i0)') "ERROR: cgca_m2mpiio/cgca_pswci2:" //             &
    " MPI_file_open: rank: ", rank
  error stop
end if

!  Set view for this process using appropriate datatype
call MPI_File_set_view(                                                &
 fh, disp, MPI_INTEGER, filetype, 'native', MPI_INFO_NULL, ierr)
if ( ierr .ne. MPI_SUCCESS ) then
  write (*,'(a,i0)') 'ERROR: cgca_m2mpiio/cgca_pswci2:' //             &
    " MPI_file_set_view: rank: ", rank
  error stop
end if

! Write all the data for this process.
! Remove halo data by passing an explicit Fortran subarray
call MPI_File_write_all( fh, coarray(:,:,:,stype), 1, mpi_subarray,    &
                         mpistat, ierr )
if ( ierr .ne. MPI_SUCCESS ) then
  write (*,'(a,i0)') 'ERROR: cgca_m2mpiio/cgca_pswci2:' //             &
    " MPI_file_write_all: rank: ", rank
  error stop
end if

!  Close file
call MPI_File_close( fh, ierr )
if ( ierr .ne. MPI_SUCCESS ) then
 write (*,'(a,i0)') 'ERROR: cgca_m2mpiio/cgca_pswci2:' //              &
   " MPI_file_close: rank: ", rank
 error stop
end if

call MPI_Type_free( filetype, ierr )
call MPI_Type_free( mpi_subarray, ierr )

end subroutine cgca_pswci2

!*roboend*

end module cgca_m2mpiio

max_stars_repo_name: lcebaman/casup

max_issues_repo_name: lcebaman/casup

max_forks_repo_name: lcebaman/casup


text:
       SUBROUTINE SPLIFT (X,Y,YP,YPP,N,W,IERR,ISX,A1,B1,AN,BN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      PARAMETER (FOUR=4.D0)
CRAY      PARAMETER (FOUR=4.0)
C
C  NJTJ
C  ###  CRAY CONVERSIONS
C  ###    1)Comment out the implicit double precision.
C  ###    2)Switch double precision parameter
C  ###      to single precision parameter
C  ###  CRAY CONVERSIONS
C  NJTJ
C
C     SANDIA MATHEMATICAL PROGRAM LIBRARY
C     APPLIED MATHEMATICS DIVISION 2613
C     SANDIA LABORATORIES
C     ALBUQUERQUE, NEW MEXICO  87185
C     CONTROL DATA 6600/7600  VERSION 7.2  MAY 1978
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                    ISSUED BY SANDIA LABORATORIES
C  *                   A PRIME CONTRACTOR TO THE
C  *                UNITED STATES DEPARTMENT OF ENERGY
C  * * * * * * * * * * * * * * * NOTICE  * * * * * * * * * * * * * * *
C  * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
C  * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
C  * UNITED STATES DEPARTMENT OF ENERGY NOR ANY OF THEIR EMPLOYEES,
C  * NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR EMPLOYEES
C  * MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL
C  * LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS OR
C  * USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT OR PROCESS
C  * DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE
C  * OWNED RIGHTS.
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS
C  * PART IS SAND77-1441.
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     WRITTEN BY RONDALL E. JONES
C
C     ABSTRACT
C         SPLIFT FITS AN INTERPOLATING CUBIC SPLINE TO THE N DATA POINT
C         GIVEN IN X AND Y AND RETURNS THE FIRST AND SECOND DERIVATIVES
C         IN YP AND YPP.  THE RESULTING SPLINE (DEFINED BY X, Y, AND
C         YPP) AND ITS FIRST AND SECOND DERIVATIVES MAY THEN BE
C         EVALUATED USING SPLINT.  THE SPLINE MAY BE INTEGRATED USING
C         SPLIQ.  FOR A SMOOTHING SPLINE FIT SEE SUBROUTINE SMOO.
C
C     DESCRIPTION OF ARGUMENTS
C         THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST
C         E.G.   X(N), Y(N), YP(N), YPP(N), W(3N)
C
C       --INPUT--
C
C         X    - ARRAY OF ABSCISSAS OF DATA (IN INCREASING ORDER)
C         Y    - ARRAY OF ORDINATES OF DATA
C         N    - THE NUMBER OF DATA POINTS.  THE ARRAYS X, Y, YP, AND
C                YPP MUST BE DIMENSIONED AT LEAST N.  (N .GE. 4)
C         ISX  - MUST BE ZERO ON THE INITIAL CALL TO SPLIFT.
C                IF A SPLINE IS TO BE FITTED TO A SECOND SET OF DATA
C                THAT HAS THE SAME SET OF ABSCISSAS AS A PREVIOUS SET,
C                AND IF THE CONTENTS OF W HAVE NOT BEEN CHANGED SINCE
C                THAT PREVIOUS FIT WAS COMPUTED, THEN ISX MAY BE
C                SET TO ONE FOR FASTER EXECUTION.
C         A1,B1,AN,BN - SPECIFY THE END CONDITIONS FOR THE SPLINE WHICH
C                ARE EXPRESSED AS CONSTRAINTS ON THE SECOND DERIVATIVE
C                OF THE SPLINE AT THE END POINTS (SEE YPP).
C                THE END CONDITION CONSTRAINTS ARE
C                        YPP(1) = A1*YPP(2) + B1
C                AND
C                        YPP(N) = AN*YPP(N-1) + BN
C                WHERE
C                        ABS(A1).LT. 1.0  AND  ABS(AN).LT. 1.0.
C
C                THE SMOOTHEST SPLINE (I.E., LEAST INTEGRAL OF SQUARE
C                OF SECOND DERIVATIVE) IS OBTAINED BY A1=B1=AN=BN=0.
C                IN THIS CASE THERE IS AN INFLECTION AT X(1) AND X(N).
C                IF THE DATA IS TO BE EXTRAPOLATED (SAY, BY USING SPLIN
C                TO EVALUATE THE SPLINE OUTSIDE THE RANGE X(1) TO X(N))
C                THEN TAKING A1=AN=0.5 AND B1=BN=0 MAY YIELD BETTER
C                RESULTS.  IN THIS CASE THERE IS AN INFLECTION
C                AT X(1) - (X(2)-X(1)) AND AT X(N) + (X(N)-X(N-1)).
C                IN THE MORE GENERAL CASE OF A1=AN=A  AND B1=BN=0,
C                THERE IS AN INFLECTION AT X(1) - (X(2)-X(1))*A/(1.0-A)
C                AND AT X(N) + (X(N)-X(N-1))*A/(1.0-A).
C
C                A SPLINE THAT HAS A GIVEN FIRST DERIVATIVE YP1 AT X(1)
C                AND YPN AT Y(N) MAY BE DEFINED BY USING THE
C                FOLLOWING CONDITIONS.
C
C                A1=-0.5
C
C                B1= 3.0*((Y(2)-Y(1))/(X(2)-X(1))-YP1)/(X(2)-X(1))
C
C                AN=-0.5
C
C                BN=-3.0*((Y(N)-Y(N-1))/(X(N)-X(N-1))-YPN)/(X(N)-X(N-1)
C
C       --OUTPUT--
C
C         YP   - ARRAY OF FIRST DERIVATIVES OF SPLINE (AT THE X(I))
C         YPP  - ARRAY OF SECOND DERIVATIVES OF SPLINE (AT THE X(I))
C         IERR - A STATUS CODE
C              --NORMAL CODE
C                 1 MEANS THAT THE REQUESTED SPLINE WAS COMPUTED.
C              --ABNORMAL CODES
C                 2 MEANS THAT N, THE NUMBER OF POINTS, WAS .LT. 4.
C                 3 MEANS THE ABSCISSAS WERE NOT STRICTLY INCREASING.
C
C       --WORK--
C
C         W    - ARRAY OF WORKING STORAGE DIMENSIONED AT LEAST 3N.
       DIMENSION X(N),Y(N),YP(N),YPP(N),W(N,3)
C
       IF (N.LT.4) THEN
         IERR = 2
         RETURN
       ENDIF
       NM1  = N-1
       NM2  = N-2
       IF (ISX.GT.0) GO TO 40
       DO 5 I=2,N
         IF (X(I)-X(I-1) .LE. 0) THEN
           IERR = 3
           RETURN
         ENDIF
 5     CONTINUE
C
C     DEFINE THE TRIDIAGONAL MATRIX
C
       W(1,3) = X(2)-X(1)
       DO 10 I=2,NM1
         W(I,2) = W(I-1,3)
         W(I,3) = X(I+1)-X(I)
 10      W(I,1) = 2*(W(I,2)+W(I,3))
       W(1,1) = FOUR
       W(1,3) =-4*A1
       W(N,1) = FOUR
       W(N,2) =-4*AN
C
C     L U DECOMPOSITION
C
       DO 30 I=2,N
         W(I-1,3) = W(I-1,3)/W(I-1,1)
 30    W(I,1) = W(I,1) - W(I,2)*W(I-1,3)
C
C     DEFINE *CONSTANT* VECTOR
C
 40   YPP(1) = 4*B1
      DOLD = (Y(2)-Y(1))/W(2,2)
      DO 50 I=2,NM2
        DNEW   = (Y(I+1) - Y(I))/W(I+1,2)
        YPP(I) = 6*(DNEW - DOLD)
        YP(I)  = DOLD
 50   DOLD = DNEW
      DNEW = (Y(N)-Y(N-1))/(X(N)-X(N-1))
      YPP(NM1) = 6*(DNEW - DOLD)
      YPP(N) = 4*BN
      YP(NM1)= DOLD
      YP(N) = DNEW
C
C     FORWARD SUBSTITUTION
C
      YPP(1) = YPP(1)/W(1,1)
      DO 60 I=2,N
 60   YPP(I) = (YPP(I) - W(I,2)*YPP(I-1))/W(I,1)
C
C     BACKWARD SUBSTITUTION
C
       DO 70 J=1,NM1
         I = N-J
   70 YPP(I) = YPP(I) - W(I,3)*YPP(I+1)
C
C     COMPUTE FIRST DERIVATIVES
C
      YP(1) = (Y(2)-Y(1))/(X(2)-X(1)) - (X(2)-X(1))*(2*YPP(1)
     1  + YPP(2))/6
      DO 80 I=2,NM1
 80   YP(I) = YP(I) + W(I,2)*(YPP(I-1) + 2*YPP(I))/6
      YP(N) = YP(N) + (X(N)-X(NM1))*(YPP(NM1) + 2*YPP(N))/6
C
      IERR = 1
      RETURN
      END

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

module pnetcdf_test
! Aim to read a netcdf restart file.
! the local directory. This is because you don't want to allocate
! state_ens_handle%vars because you will run out of memory

use types_mod,            only : r8, missing_r8
use mpi_utilities_mod,    only : initialize_mpi_utilities,     &
                                 finalize_mpi_utilities,       &
                                 task_count, my_task_id
use utilities_mod,        only : nc_check
use model_mod,            only : fill_variable_list,           &
                                 info_file_name,               &
                                 construct_file_name_in
use ensemble_manager_mod, only : ensemble_type, print_ens_handle
use io_filenames_mod,     only : restart_files_in
use assim_tools_mod,      only : test_state_copies

use mpi
use pnetcdf

implicit none

public :: parallel_read_transpose

contains

subroutine parallel_read_transpose(state_ens_handle, restart_in_file_name, domain, dart_index)

type(ensemble_type), intent(inout) :: state_ens_handle
character(len=129),  intent(in)    :: restart_in_file_name
integer,             intent(in)    :: domain
integer,             intent(inout) :: dart_index !< This is for mulitple domains

integer(KIND=MPI_OFFSET_KIND), allocatable :: stride(:)

real(r8), allocatable ::  data_block(:)

! ncfile
integer            :: ncfile, ndims, nvars, ngatts, unlimited
integer            :: var_ndims
integer            :: var_size
integer            :: var_dimids(3) !dimensions ids (i,j,k)
integer            :: var_ids(5)
character(len=256) :: var_names(5)

integer :: j, i, k

integer(KIND=MPI_OFFSET_KIND),allocatable :: dim_sizes(:)
character(len=256)           ,allocatable :: dim_names(:)

integer :: total_var_size(5)
integer :: numdims
integer :: numvars

integer(KIND=MPI_OFFSET_KIND),allocatable :: start(:)
integer(KIND=MPI_OFFSET_KIND),allocatable :: count(:)

!character(len=256)  :: dir
character(len=1024) :: filename
character(len=256)  :: variable_list(5)

real(r8)                      :: begin

integer :: num_copies
integer :: my_pe
integer :: icopy
integer :: istart
integer :: iend

!call initialize_mpi_utilities('test io')
if(my_task_id()==0) then
call print_ens_handle(state_ens_handle)
endif

num_copies = state_ens_handle%num_copies
my_pe      = state_ens_handle%my_pe

!dir      = '/glade/p/work/hendric/test_cases/pop_large_test_case/'
filename = trim(info_file_name(1))

! specific for POP
! SATL_CUR, TEMP,CUR, UVEL_VUR, VVEL_CUR, PSURC_CUR
numvars   = 5
var_names = fill_variable_list(numvars)

istart = 1

COPIES: do icopy = 1,num_copies
   begin = MPI_WTIME()

   filename = restart_files_in(icopy, domain)

   if(my_task_id()==0) write(*,*) "opening : " , trim(filename)

   ! Open netcdf file
   call nc_check(nf90mpi_open(MPI_COMM_WORLD,trim(filename),       &
                 NF90_NOWRITE + NF90_64BIT_OFFSET, MPI_INFO_NULL,  &
                 ncfile),'open error')

   ! inquire general information
   call nc_check(nfmpi_inq(ncfile, ndims=ndims, nvars=nvars,       &
                 ngatts=ngatts,unlimdimid=unlimited), 'nf90mpi_inq')

   if(my_task_id()==0) then
      121 format('ndims [',I2,'], nvars [',I3,'], ngatts[',I4,']')
      print 121, ndims, nvars, ngatts
      print *, " "
   endif

   ! inquire dimension information
   allocate(dim_sizes(ndims))
   allocate(dim_names(ndims))

   do i = 1,ndims
     call nc_check(nf90mpi_inquire_dimension(ncfile, i,            &
                   name=dim_names(i), len=dim_sizes(i)),'dim_len')
     if(my_task_id()==0) then
       120 format('dim ',A,' size = ',I4)
       print 120, trim(dim_names(i)), dim_sizes(i)
     endif
   enddo

   if(my_task_id() == 0) print*, " "

   ! get the var_names ids
   do i = 1,numvars
     call nc_check(nf90mpi_inq_varid(ncfile, var_names(i), &
                   varid=var_ids(i)), 'inq_varid')
     if(my_task_id() == 0) then
       123 format('var_id[',I2,'] variable_name = ',A9)
       print 123, var_ids(i), trim(var_names(i))
     endif
   enddo

   if(my_task_id() == 0) print*, " "

   do i = 1,numvars
     call nc_check(nf90mpi_inquire_variable(ncfile, var_ids(i), &
                   ndims=var_ndims, dimids=var_dimids),'inq_var')

     allocate(start(var_ndims))
     allocate(count(var_ndims))
     allocate(stride(var_ndims))

     !! hyperslab read
     !k = 2
     !start(k) = (dim_sizes(var_dimids(k))/task_count())*(my_task_id())+1
     !count(k) = (dim_sizes(var_dimids(k))/task_count())
     !var_size = count(k)

     !do j = 1,var_ndims
     !  if (j /= k) then
     !    stride(j) = task_count()
     !    start(j) = 1
     !    count(j) = dim_sizes(var_dimids(j))
     !    var_size = var_size*count(j)
     !  endif
     !enddo

     ! strided read
     k = 1
     start(k)  = my_task_id()+1
     count(k)  = (dim_sizes(var_dimids(k))/task_count())
     stride(k) = task_count()
     var_size  = count(k)

     do j = 1,var_ndims
       if (j /= k) then
         stride(j) = 1
         start(j)  = 1
         count(j)  = dim_sizes(var_dimids(j))
         var_size  = var_size*count(j)
       endif
     enddo

     iend = istart+var_size

     allocate(data_block(var_size))
     call nc_check(nfmpi_get_vars_double_all(ncfile,var_ids(i), &
                   start=start,count=count,stride=stride,       &
                   dvals=data_block),'get')

     state_ens_handle%copies(icopy,istart:iend) = data_block
     !if (my_task_id()==0) then
     !  open(unit=12,FILE="x2.txt",FORM="FORMATTED",STATUS="OLD", &
     !       ACTION="WRITE")
     !  do j = 1, var_size
     !     write(12,*) data_block(j)
     !  end do
     !endif
     istart = iend

     deallocate(start)
     deallocate(count)
     deallocate(stride)
     deallocate(data_block)

   enddo

   if(my_task_id() == 0) then
     888 format('MPI Wall Time = ',F5.2, 'sec')
     print 888, MPI_WTIME()-begin
     print *, " "
   endif

   ! Get the variable dimensions
   do i = 1,numvars
     call nc_check(nf90mpi_inquire_variable(ncfile, var_ids(i), &
                   ndims=numdims), 'inq_dimid')

     total_var_size(i) = product(dim_sizes(1:numdims))

     if(my_task_id() == 0) then
       124 format('var_id[',I2,'] num_dims = ',I1,' total_var_size = ',I9)
       print 124, var_ids(i), numdims, total_var_size(i)
     endif
   enddo

   deallocate(dim_sizes)
   deallocate(dim_names)

enddo COPIES

  i = my_task_id()+1
  write(filename,'(A,I1.1)') 'p.',i
  open(unit=i,FILE=filename,FORM="FORMATTED",STATUS="OLD", &
       ACTION="WRITE")
  do j = 1, state_ens_handle%my_num_vars
     write(i,*) state_ens_handle%copies(1,j)
  end do
  close(unit=i)

end subroutine parallel_read_transpose

end module pnetcdf_test

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


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

20110319 21:27:27 nbsp Welcome to the Wiki. The wiki is here so that we can all explore, discuss and compile anything and everything about Davis — especially the little, enjoyable things. http://daviswiki.org/Yolo_Fruit_Stand?actiondiff&version288&version187 Your comment about the Yolo Fruit Stand contains an accusation about how they conduct their business. How should the wiki react to that comment? How would you feel if you were the owner of a business and someone made a statement like that about you? Users/JasonAller


text:
Subroutine crkkpi(i1, i2, xsk1, xsk2, xsk3, xsk4, xsk5, xsk6, xsk7, xsk8, xsk9, xsk10, xsk11, sigk, iblock, lbp1, lbp2, emm1, emm2)
  Parameter (maxstr=150001, maxr=1, amn=0.939457, amp=0.93828, ap1=0.13496, amrho=0.769, amomga=0.782, ameta=0.5473, ap2=0.13957, am0=1.232, pi=3.1415926, cutoff=1.8966, avmass=0.9383)
  Parameter (aka=0.498, ala=1.1157, asa=1.1974)
  Parameter (mx=4, my=4, mz=8, mpx=4, mpy=4, mpz=10, mpzp=10)
  Common /aa/r(3, maxstr)
  Common /bb/p(3, maxstr)
  Common /cc/e(maxstr)
  Common /ee/id(maxstr), lb(maxstr)
  Common /input1/masspr, massta, iseed, iavoid, dt
  Common /rndf77/nseed
  Save
  iblock = 1907
  x1 = ranart(nseed)*sigk
  xsk2 = xsk1 + xsk2
  xsk3 = xsk2 + xsk3
  xsk4 = xsk3 + xsk4
  xsk5 = xsk4 + xsk5
  xsk6 = xsk5 + xsk6
  xsk7 = xsk6 + xsk7
  xsk8 = xsk7 + xsk8
  xsk9 = xsk8 + xsk9
  xsk10 = xsk9 + xsk10
  If (x1<=xsk1) Then
    lb(i1) = 3 + int(3*ranart(nseed))
    lb(i2) = 3 + int(3*ranart(nseed))
    e(i1) = ap2
    e(i2) = ap2
    Goto 100
  Else If (x1<=xsk2) Then
    lb(i1) = 3 + int(3*ranart(nseed))
    lb(i2) = 25 + int(3*ranart(nseed))
    e(i1) = ap2
    e(i2) = amrho
    Goto 100
  Else If (x1<=xsk3) Then
    lb(i1) = 3 + int(3*ranart(nseed))
    lb(i2) = 28
    e(i1) = ap2
    e(i2) = amomga
    Goto 100
  Else If (x1<=xsk4) Then
    lb(i1) = 3 + int(3*ranart(nseed))
    lb(i2) = 0
    e(i1) = ap2
    e(i2) = ameta
    Goto 100
  Else If (x1<=xsk5) Then
    lb(i1) = 25 + int(3*ranart(nseed))
    lb(i2) = 25 + int(3*ranart(nseed))
    e(i1) = amrho
    e(i2) = amrho
    Goto 100
  Else If (x1<=xsk6) Then
    lb(i1) = 25 + int(3*ranart(nseed))
    lb(i2) = 28
    e(i1) = amrho
    e(i2) = amomga
    Goto 100
  Else If (x1<=xsk7) Then
    lb(i1) = 25 + int(3*ranart(nseed))
    lb(i2) = 0
    e(i1) = amrho
    e(i2) = ameta
    Goto 100
  Else If (x1<=xsk8) Then
    lb(i1) = 28
    lb(i2) = 28
    e(i1) = amomga
    e(i2) = amomga
    Goto 100
  Else If (x1<=xsk9) Then
    lb(i1) = 28
    lb(i2) = 0
    e(i1) = amomga
    e(i2) = ameta
    Goto 100
  Else If (x1<=xsk10) Then
    lb(i1) = 0
    lb(i2) = 0
    e(i1) = ameta
    e(i2) = ameta
  Else
    iblock = 222
    Call rhores(i1, i2)
    lb(i1) = 29
    e(i2) = 0.
  End If
  100 Continue
  lbp1 = lb(i1)
  lbp2 = lb(i2)
  emm1 = e(i1)
  emm2 = e(i2)
  Return
End Subroutine crkkpi

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
      double precision bbsqmin,bbsqmax,wsqmin,wsqmax,m3456min,m3456max
      common/limits/bbsqmin,bbsqmax,wsqmin,wsqmax,m3456min,m3456max

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


!***********************************/
!	Name:変位量の入力
!	Alias:OGetDisplacementMK2
!	Description:
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.01.10
!	Update:
!	Other:
!***********************************/
subroutine OGetDisplacement(MG,MW)
    use StructVar_Mod
    use ConstantVar_Mod
    implicit none

    type(MoveGrid), intent(inout) :: MG
    type(MotionWait), intent(in) :: MW

    MG%GCD%Translation = MW%Accumulate%Translation

    MG%GCD%Rotation = MW%Accumulate%Rotation
    !ここでの入力変位量は局所座標における変位量になっているため， 逆回転行列を掛けてグローバル座標系での表記に戻してやる必要がある
    MG%GCD%Translation(:) = matmul(MG%TM%Top2Next(1:3,1:3),MG%GCD%Translation)

return
end subroutine OGetDisplacement


text:
MODULE nrutil
  USE nrtype
  IMPLICIT NONE
  INTEGER(I4B), PARAMETER :: NPAR_ARTH=16,NPAR2_ARTH=8
  INTEGER(I4B), PARAMETER :: NPAR_GEOP=4,NPAR2_GEOP=2
  INTEGER(I4B), PARAMETER :: NPAR_CUMSUM=16
  INTEGER(I4B), PARAMETER :: NPAR_CUMPROD=8
  INTEGER(I4B), PARAMETER :: NPAR_POLY=8
  INTEGER(I4B), PARAMETER :: NPAR_POLYTERM=8
  INTERFACE array_copy
     MODULE PROCEDURE array_copy_r, array_copy_d, array_copy_i
  END INTERFACE
  INTERFACE swap
     MODULE PROCEDURE swap_i,swap_r,swap_rv,swap_c, &
          swap_cv,swap_cm,swap_z,swap_zv,swap_zm, &
          masked_swap_rs,masked_swap_rv,masked_swap_rm
  END INTERFACE
  INTERFACE reallocate
     MODULE PROCEDURE reallocate_rv,reallocate_rm,&
          reallocate_iv,reallocate_im,reallocate_hv
  END INTERFACE
  INTERFACE imaxloc
     MODULE PROCEDURE imaxloc_r,imaxloc_i
  END INTERFACE
  INTERFACE assert
     MODULE PROCEDURE assert1,assert2,assert3,assert4,assert_v
  END INTERFACE
  INTERFACE assert_eq
     MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn
  END INTERFACE
  INTERFACE arth
     MODULE PROCEDURE arth_r, arth_d, arth_i
  END INTERFACE
  INTERFACE geop
     MODULE PROCEDURE geop_r, geop_d, geop_i, geop_c, geop_dv
  END INTERFACE
  INTERFACE cumsum
     MODULE PROCEDURE cumsum_r,cumsum_i
  END INTERFACE
  INTERFACE poly
     MODULE PROCEDURE poly_rr,poly_rrv,poly_dd,poly_ddv,&
          poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv
  END INTERFACE
  INTERFACE poly_term
     MODULE PROCEDURE poly_term_rr,poly_term_cc
  END INTERFACE
  INTERFACE outerprod
     MODULE PROCEDURE outerprod_r,outerprod_d
  END INTERFACE
  INTERFACE outerdiff
     MODULE PROCEDURE outerdiff_r,outerdiff_d,outerdiff_i
  END INTERFACE
  INTERFACE scatter_add
     MODULE PROCEDURE scatter_add_r,scatter_add_d
  END INTERFACE
  INTERFACE scatter_max
     MODULE PROCEDURE scatter_max_r,scatter_max_d
  END INTERFACE
  INTERFACE diagadd
     MODULE PROCEDURE diagadd_rv,diagadd_r
  END INTERFACE
  INTERFACE diagmult
     MODULE PROCEDURE diagmult_rv,diagmult_r
  END INTERFACE
  INTERFACE get_diag
     MODULE PROCEDURE get_diag_rv, get_diag_dv
  END INTERFACE
  INTERFACE put_diag
     MODULE PROCEDURE put_diag_rv, put_diag_r
  END INTERFACE
CONTAINS
  !BL
  SUBROUTINE array_copy_r(src,dest,n_copied,n_not_copied)
    REAL(SP), DIMENSION(:), INTENT(IN) :: src
    REAL(SP), DIMENSION(:), INTENT(OUT) :: dest
    INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied
    n_copied=min(size(src),size(dest))
    n_not_copied=size(src)-n_copied
    dest(1:n_copied)=src(1:n_copied)
  END SUBROUTINE array_copy_r
  !BL
  SUBROUTINE array_copy_d(src,dest,n_copied,n_not_copied)
    REAL(DP), DIMENSION(:), INTENT(IN) :: src
    REAL(DP), DIMENSION(:), INTENT(OUT) :: dest
    INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied
    n_copied=min(size(src),size(dest))
    n_not_copied=size(src)-n_copied
    dest(1:n_copied)=src(1:n_copied)
  END SUBROUTINE array_copy_d
  !BL
  SUBROUTINE array_copy_i(src,dest,n_copied,n_not_copied)
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: src
    INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: dest
    INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied
    n_copied=min(size(src),size(dest))
    n_not_copied=size(src)-n_copied
    dest(1:n_copied)=src(1:n_copied)
  END SUBROUTINE array_copy_i
  !BL
  !BL
  SUBROUTINE swap_i(a,b)
    INTEGER(I4B), INTENT(INOUT) :: a,b
    INTEGER(I4B) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_i
  !BL
  SUBROUTINE swap_r(a,b)
    REAL(SP), INTENT(INOUT) :: a,b
    REAL(SP) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_r
  !BL
  SUBROUTINE swap_rv(a,b)
    REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b
    REAL(SP), DIMENSION(SIZE(a)) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_rv
  !BL
  SUBROUTINE swap_c(a,b)
    COMPLEX(SPC), INTENT(INOUT) :: a,b
    COMPLEX(SPC) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_c
  !BL
  SUBROUTINE swap_cv(a,b)
    COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: a,b
    COMPLEX(SPC), DIMENSION(SIZE(a)) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_cv
  !BL
  SUBROUTINE swap_cm(a,b)
    COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: a,b
    COMPLEX(SPC), DIMENSION(size(a,1),size(a,2)) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_cm
  !BL
  SUBROUTINE swap_z(a,b)
    COMPLEX(DPC), INTENT(INOUT) :: a,b
    COMPLEX(DPC) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_z
  !BL
  SUBROUTINE swap_zv(a,b)
    COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: a,b
    COMPLEX(DPC), DIMENSION(SIZE(a)) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_zv
  !BL
  SUBROUTINE swap_zm(a,b)
    COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: a,b
    COMPLEX(DPC), DIMENSION(size(a,1),size(a,2)) :: dum
    dum=a
    a=b
    b=dum
  END SUBROUTINE swap_zm
  !BL
  SUBROUTINE masked_swap_rs(a,b,mask)
    REAL(SP), INTENT(INOUT) :: a,b
    LOGICAL(LGT), INTENT(IN) :: mask
    REAL(SP) :: swp
    if (mask) then
       swp=a
       a=b
       b=swp
    end if
  END SUBROUTINE masked_swap_rs
  !BL
  SUBROUTINE masked_swap_rv(a,b,mask)
    REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b
    LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask
    REAL(SP), DIMENSION(size(a)) :: swp
    where (mask)
       swp=a
       a=b
       b=swp
    end where
  END SUBROUTINE masked_swap_rv
  !BL
  SUBROUTINE masked_swap_rm(a,b,mask)
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b
    LOGICAL(LGT), DIMENSION(:,:), INTENT(IN) :: mask
    REAL(SP), DIMENSION(size(a,1),size(a,2)) :: swp
    where (mask)
       swp=a
       a=b
       b=swp
    end where
  END SUBROUTINE masked_swap_rm
  !BL
  !BL
  FUNCTION reallocate_rv(p,n)
    REAL(SP), DIMENSION(:), POINTER :: p, reallocate_rv
    INTEGER(I4B), INTENT(IN) :: n
    INTEGER(I4B) :: nold,ierr
    allocate(reallocate_rv(n),stat=ierr)
    if (ierr /= 0) call &
         nrerror('reallocate_rv: problem in attempt to allocate memory')
    if (.not. associated(p)) RETURN
    nold=size(p)
    reallocate_rv(1:min(nold,n))=p(1:min(nold,n))
    deallocate(p)
  END FUNCTION reallocate_rv
  !BL
  FUNCTION reallocate_iv(p,n)
    INTEGER(I4B), DIMENSION(:), POINTER :: p, reallocate_iv
    INTEGER(I4B), INTENT(IN) :: n
    INTEGER(I4B) :: nold,ierr
    allocate(reallocate_iv(n),stat=ierr)
    if (ierr /= 0) call &
         nrerror('reallocate_iv: problem in attempt to allocate memory')
    if (.not. associated(p)) RETURN
    nold=size(p)
    reallocate_iv(1:min(nold,n))=p(1:min(nold,n))
    deallocate(p)
  END FUNCTION reallocate_iv
  !BL
  FUNCTION reallocate_hv(p,n)
    CHARACTER(1), DIMENSION(:), POINTER :: p, reallocate_hv
    INTEGER(I4B), INTENT(IN) :: n
    INTEGER(I4B) :: nold,ierr
    allocate(reallocate_hv(n),stat=ierr)
    if (ierr /= 0) call &
         nrerror('reallocate_hv: problem in attempt to allocate memory')
    if (.not. associated(p)) RETURN
    nold=size(p)
    reallocate_hv(1:min(nold,n))=p(1:min(nold,n))
    deallocate(p)
  END FUNCTION reallocate_hv
  !BL
  FUNCTION reallocate_rm(p,n,m)
    REAL(SP), DIMENSION(:,:), POINTER :: p, reallocate_rm
    INTEGER(I4B), INTENT(IN) :: n,m
    INTEGER(I4B) :: nold,mold,ierr
    allocate(reallocate_rm(n,m),stat=ierr)
    if (ierr /= 0) call &
         nrerror('reallocate_rm: problem in attempt to allocate memory')
    if (.not. associated(p)) RETURN
    nold=size(p,1)
    mold=size(p,2)
    reallocate_rm(1:min(nold,n),1:min(mold,m))=&
         p(1:min(nold,n),1:min(mold,m))
    deallocate(p)
  END FUNCTION reallocate_rm
  !BL
  FUNCTION reallocate_im(p,n,m)
    INTEGER(I4B), DIMENSION(:,:), POINTER :: p, reallocate_im
    INTEGER(I4B), INTENT(IN) :: n,m
    INTEGER(I4B) :: nold,mold,ierr
    allocate(reallocate_im(n,m),stat=ierr)
    if (ierr /= 0) call &
         nrerror('reallocate_im: problem in attempt to allocate memory')
    if (.not. associated(p)) RETURN
    nold=size(p,1)
    mold=size(p,2)
    reallocate_im(1:min(nold,n),1:min(mold,m))=&
         p(1:min(nold,n),1:min(mold,m))
    deallocate(p)
  END FUNCTION reallocate_im
  !BL
  FUNCTION ifirstloc(mask)
    LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask
    INTEGER(I4B) :: ifirstloc
    INTEGER(I4B), DIMENSION(1) :: loc
    loc=maxloc(merge(1,0,mask))
    ifirstloc=loc(1)
    if (.not. mask(ifirstloc)) ifirstloc=size(mask)+1
  END FUNCTION ifirstloc
  !BL
  FUNCTION imaxloc_r(arr)
    REAL(SP), DIMENSION(:), INTENT(IN) :: arr
    INTEGER(I4B) :: imaxloc_r
    INTEGER(I4B), DIMENSION(1) :: imax
    imax=maxloc(arr(:))
    imaxloc_r=imax(1)
  END FUNCTION imaxloc_r
  !BL
  FUNCTION imaxloc_i(iarr)
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr
    INTEGER(I4B), DIMENSION(1) :: imax
    INTEGER(I4B) :: imaxloc_i
    imax=maxloc(iarr(:))
    imaxloc_i=imax(1)
  END FUNCTION imaxloc_i
  !BL
  FUNCTION iminloc(arr)
    REAL(SP), DIMENSION(:), INTENT(IN) :: arr
    INTEGER(I4B), DIMENSION(1) :: imin
    INTEGER(I4B) :: iminloc
    imin=minloc(arr(:))
    iminloc=imin(1)
  END FUNCTION iminloc
  !BL
  SUBROUTINE assert1(n1,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    LOGICAL, INTENT(IN) :: n1
    if (.not. n1) then
       write (*,*) 'nrerror: an assertion failed with this tag:', &
            string
       STOP 'program terminated by assert1'
    end if
  END SUBROUTINE assert1
  !BL
  SUBROUTINE assert2(n1,n2,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    LOGICAL, INTENT(IN) :: n1,n2
    if (.not. (n1 .and. n2)) then
       write (*,*) 'nrerror: an assertion failed with this tag:', &
            string
       STOP 'program terminated by assert2'
    end if
  END SUBROUTINE assert2
  !BL
  SUBROUTINE assert3(n1,n2,n3,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    LOGICAL, INTENT(IN) :: n1,n2,n3
    if (.not. (n1 .and. n2 .and. n3)) then
       write (*,*) 'nrerror: an assertion failed with this tag:', &
            string
       STOP 'program terminated by assert3'
    end if
  END SUBROUTINE assert3
  !BL
  SUBROUTINE assert4(n1,n2,n3,n4,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    LOGICAL, INTENT(IN) :: n1,n2,n3,n4
    if (.not. (n1 .and. n2 .and. n3 .and. n4)) then
       write (*,*) 'nrerror: an assertion failed with this tag:', &
            string
       STOP 'program terminated by assert4'
    end if
  END SUBROUTINE assert4
  !BL
  SUBROUTINE assert_v(n,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    LOGICAL, DIMENSION(:), INTENT(IN) :: n
    if (.not. all(n)) then
       write (*,*) 'nrerror: an assertion failed with this tag:', &
            string
       STOP 'program terminated by assert_v'
    end if
  END SUBROUTINE assert_v
  !BL
  FUNCTION assert_eq2(n1,n2,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    INTEGER, INTENT(IN) :: n1,n2
    INTEGER :: assert_eq2
    if (n1 == n2) then
       assert_eq2=n1
    else
       write (*,*) 'nrerror: an assert_eq failed with this tag:', &
            string
       STOP 'program terminated by assert_eq2'
    end if
  END FUNCTION assert_eq2
  !BL
  FUNCTION assert_eq3(n1,n2,n3,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    INTEGER, INTENT(IN) :: n1,n2,n3
    INTEGER :: assert_eq3
    if (n1 == n2 .and. n2 == n3) then
       assert_eq3=n1
    else
       write (*,*) 'nrerror: an assert_eq failed with this tag:', &
            string
       STOP 'program terminated by assert_eq3'
    end if
  END FUNCTION assert_eq3
  !BL
  FUNCTION assert_eq4(n1,n2,n3,n4,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    INTEGER, INTENT(IN) :: n1,n2,n3,n4
    INTEGER :: assert_eq4
    if (n1 == n2 .and. n2 == n3 .and. n3 == n4) then
       assert_eq4=n1
    else
       write (*,*) 'nrerror: an assert_eq failed with this tag:', &
            string
       STOP 'program terminated by assert_eq4'
    end if
  END FUNCTION assert_eq4
  !BL
  FUNCTION assert_eqn(nn,string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    INTEGER, DIMENSION(:), INTENT(IN) :: nn
    INTEGER :: assert_eqn
    if (all(nn(2:) == nn(1))) then
       assert_eqn=nn(1)
    else
       write (*,*) 'nrerror: an assert_eq failed with this tag:', &
            string
       STOP 'program terminated by assert_eqn'
    end if
  END FUNCTION assert_eqn
  !BL
  SUBROUTINE nrerror(string)
    CHARACTER(LEN=*), INTENT(IN) :: string
    write (*,*) 'nrerror: ',string
    STOP 'program terminated by nrerror'
  END SUBROUTINE nrerror
  !BL
  FUNCTION arth_r(first,increment,n)
    REAL(SP), INTENT(IN) :: first,increment
    INTEGER(I4B), INTENT(IN) :: n
    REAL(SP), DIMENSION(n) :: arth_r
    INTEGER(I4B) :: k,k2
    REAL(SP) :: temp
    if (n > 0) arth_r(1)=first
    if (n <= NPAR_ARTH) then
       do k=2,n
          arth_r(k)=arth_r(k-1)+increment
       end do
    else
       do k=2,NPAR2_ARTH
          arth_r(k)=arth_r(k-1)+increment
       end do
       temp=increment*NPAR2_ARTH
       k=NPAR2_ARTH
       do
          if (k >= n) exit
          k2=k+k
          arth_r(k+1:min(k2,n))=temp+arth_r(1:min(k,n-k))
          temp=temp+temp
          k=k2
       end do
    end if
  END FUNCTION arth_r
  !BL
  FUNCTION arth_d(first,increment,n)
    REAL(DP), INTENT(IN) :: first,increment
    INTEGER(I4B), INTENT(IN) :: n
    REAL(DP), DIMENSION(n) :: arth_d
    INTEGER(I4B) :: k,k2
    REAL(DP) :: temp
    if (n > 0) arth_d(1)=first
    if (n <= NPAR_ARTH) then
       do k=2,n
          arth_d(k)=arth_d(k-1)+increment
       end do
    else
       do k=2,NPAR2_ARTH
          arth_d(k)=arth_d(k-1)+increment
       end do
       temp=increment*NPAR2_ARTH
       k=NPAR2_ARTH
       do
          if (k >= n) exit
          k2=k+k
          arth_d(k+1:min(k2,n))=temp+arth_d(1:min(k,n-k))
          temp=temp+temp
          k=k2
       end do
    end if
  END FUNCTION arth_d
  !BL
  FUNCTION arth_i(first,increment,n)
    INTEGER(I4B), INTENT(IN) :: first,increment,n
    INTEGER(I4B), DIMENSION(n) :: arth_i
    INTEGER(I4B) :: k,k2,temp
    if (n > 0) arth_i(1)=first
    if (n <= NPAR_ARTH) then
       do k=2,n
          arth_i(k)=arth_i(k-1)+increment
       end do
    else
       do k=2,NPAR2_ARTH
          arth_i(k)=arth_i(k-1)+increment
       end do
       temp=increment*NPAR2_ARTH
       k=NPAR2_ARTH
       do
          if (k >= n) exit
          k2=k+k
          arth_i(k+1:min(k2,n))=temp+arth_i(1:min(k,n-k))
          temp=temp+temp
          k=k2
       end do
    end if
  END FUNCTION arth_i
  !BL
  !BL
  FUNCTION geop_r(first,factor,n)
    REAL(SP), INTENT(IN) :: first,factor
    INTEGER(I4B), INTENT(IN) :: n
    REAL(SP), DIMENSION(n) :: geop_r
    INTEGER(I4B) :: k,k2
    REAL(SP) :: temp
    if (n > 0) geop_r(1)=first
    if (n <= NPAR_GEOP) then
       do k=2,n
          geop_r(k)=geop_r(k-1)*factor
       end do
    else
       do k=2,NPAR2_GEOP
          geop_r(k)=geop_r(k-1)*factor
       end do
       temp=factor**NPAR2_GEOP
       k=NPAR2_GEOP
       do
          if (k >= n) exit
          k2=k+k
          geop_r(k+1:min(k2,n))=temp*geop_r(1:min(k,n-k))
          temp=temp*temp
          k=k2
       end do
    end if
  END FUNCTION geop_r
  !BL
  FUNCTION geop_d(first,factor,n)
    REAL(DP), INTENT(IN) :: first,factor
    INTEGER(I4B), INTENT(IN) :: n
    REAL(DP), DIMENSION(n) :: geop_d
    INTEGER(I4B) :: k,k2
    REAL(DP) :: temp
    if (n > 0) geop_d(1)=first
    if (n <= NPAR_GEOP) then
       do k=2,n
          geop_d(k)=geop_d(k-1)*factor
       end do
    else
       do k=2,NPAR2_GEOP
          geop_d(k)=geop_d(k-1)*factor
       end do
       temp=factor**NPAR2_GEOP
       k=NPAR2_GEOP
       do
          if (k >= n) exit
          k2=k+k
          geop_d(k+1:min(k2,n))=temp*geop_d(1:min(k,n-k))
          temp=temp*temp
          k=k2
       end do
    end if
  END FUNCTION geop_d
  !BL
  FUNCTION geop_i(first,factor,n)
    INTEGER(I4B), INTENT(IN) :: first,factor,n
    INTEGER(I4B), DIMENSION(n) :: geop_i
    INTEGER(I4B) :: k,k2,temp
    if (n > 0) geop_i(1)=first
    if (n <= NPAR_GEOP) then
       do k=2,n
          geop_i(k)=geop_i(k-1)*factor
       end do
    else
       do k=2,NPAR2_GEOP
          geop_i(k)=geop_i(k-1)*factor
       end do
       temp=factor**NPAR2_GEOP
       k=NPAR2_GEOP
       do
          if (k >= n) exit
          k2=k+k
          geop_i(k+1:min(k2,n))=temp*geop_i(1:min(k,n-k))
          temp=temp*temp
          k=k2
       end do
    end if
  END FUNCTION geop_i
  !BL
  FUNCTION geop_c(first,factor,n)
    COMPLEX(SP), INTENT(IN) :: first,factor
    INTEGER(I4B), INTENT(IN) :: n
    COMPLEX(SP), DIMENSION(n) :: geop_c
    INTEGER(I4B) :: k,k2
    COMPLEX(SP) :: temp
    if (n > 0) geop_c(1)=first
    if (n <= NPAR_GEOP) then
       do k=2,n
          geop_c(k)=geop_c(k-1)*factor
       end do
    else
       do k=2,NPAR2_GEOP
          geop_c(k)=geop_c(k-1)*factor
       end do
       temp=factor**NPAR2_GEOP
       k=NPAR2_GEOP
       do
          if (k >= n) exit
          k2=k+k
          geop_c(k+1:min(k2,n))=temp*geop_c(1:min(k,n-k))
          temp=temp*temp
          k=k2
       end do
    end if
  END FUNCTION geop_c
  !BL
  FUNCTION geop_dv(first,factor,n)
    REAL(DP), DIMENSION(:), INTENT(IN) :: first,factor
    INTEGER(I4B), INTENT(IN) :: n
    REAL(DP), DIMENSION(size(first),n) :: geop_dv
    INTEGER(I4B) :: k,k2
    REAL(DP), DIMENSION(size(first)) :: temp
    if (n > 0) geop_dv(:,1)=first(:)
    if (n <= NPAR_GEOP) then
       do k=2,n
          geop_dv(:,k)=geop_dv(:,k-1)*factor(:)
       end do
    else
       do k=2,NPAR2_GEOP
          geop_dv(:,k)=geop_dv(:,k-1)*factor(:)
       end do
       temp=factor**NPAR2_GEOP
       k=NPAR2_GEOP
       do
          if (k >= n) exit
          k2=k+k
          geop_dv(:,k+1:min(k2,n))=geop_dv(:,1:min(k,n-k))*&
               spread(temp,2,size(geop_dv(:,1:min(k,n-k)),2))
          temp=temp*temp
          k=k2
       end do
    end if
  END FUNCTION geop_dv
  !BL
  !BL
  RECURSIVE FUNCTION cumsum_r(arr,seed) RESULT(ans)
    REAL(SP), DIMENSION(:), INTENT(IN) :: arr
    REAL(SP), OPTIONAL, INTENT(IN) :: seed
    REAL(SP), DIMENSION(size(arr)) :: ans
    INTEGER(I4B) :: n,j
    REAL(SP) :: sd
    n=size(arr)
    if (n == 0_i4b) RETURN
    sd=0.0_sp
    if (present(seed)) sd=seed
    ans(1)=arr(1)+sd
    if (n < NPAR_CUMSUM) then
       do j=2,n
          ans(j)=ans(j-1)+arr(j)
       end do
    else
       ans(2:n:2)=cumsum_r(arr(2:n:2)+arr(1:n-1:2),sd)
       ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2)
    end if
  END FUNCTION cumsum_r
  !BL
  RECURSIVE FUNCTION cumsum_i(arr,seed) RESULT(ans)
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: arr
    INTEGER(I4B), OPTIONAL, INTENT(IN) :: seed
    INTEGER(I4B), DIMENSION(size(arr)) :: ans
    INTEGER(I4B) :: n,j,sd
    n=size(arr)
    if (n == 0_i4b) RETURN
    sd=0_i4b
    if (present(seed)) sd=seed
    ans(1)=arr(1)+sd
    if (n < NPAR_CUMSUM) then
       do j=2,n
          ans(j)=ans(j-1)+arr(j)
       end do
    else
       ans(2:n:2)=cumsum_i(arr(2:n:2)+arr(1:n-1:2),sd)
       ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2)
    end if
  END FUNCTION cumsum_i
  !BL
  !BL
  RECURSIVE FUNCTION cumprod(arr,seed) RESULT(ans)
    REAL(SP), DIMENSION(:), INTENT(IN) :: arr
    REAL(SP), OPTIONAL, INTENT(IN) :: seed
    REAL(SP), DIMENSION(size(arr)) :: ans
    INTEGER(I4B) :: n,j
    REAL(SP) :: sd
    n=size(arr)
    if (n == 0_i4b) RETURN
    sd=1.0_sp
    if (present(seed)) sd=seed
    ans(1)=arr(1)*sd
    if (n < NPAR_CUMPROD) then
       do j=2,n
          ans(j)=ans(j-1)*arr(j)
       end do
    else
       ans(2:n:2)=cumprod(arr(2:n:2)*arr(1:n-1:2),sd)
       ans(3:n:2)=ans(2:n-1:2)*arr(3:n:2)
    end if
  END FUNCTION cumprod
  !BL
  !BL
  FUNCTION poly_rr(x,coeffs)
    REAL(SP), INTENT(IN) :: x
    REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs
    REAL(SP) :: poly_rr
    REAL(SP) :: pow
    REAL(SP), DIMENSION(:), ALLOCATABLE :: vec
    INTEGER(I4B) :: i,n,nn
    n=size(coeffs)
    if (n <= 0) then
       poly_rr=0.0_sp
    else if (n < NPAR_POLY) then
       poly_rr=coeffs(n)
       do i=n-1,1,-1
          poly_rr=x*poly_rr+coeffs(i)
       end do
    else
       allocate(vec(n+1))
       pow=x
       vec(1:n)=coeffs
       do
          vec(n+1)=0.0_sp
          nn=ishft(n+1,-1)
          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
          if (nn == 1) exit
          pow=pow*pow
          n=nn
       end do
       poly_rr=vec(1)
       deallocate(vec)
    end if
  END FUNCTION poly_rr
  !BL
  FUNCTION poly_dd(x,coeffs)
    REAL(DP), INTENT(IN) :: x
    REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs
    REAL(DP) :: poly_dd
    REAL(DP) :: pow
    REAL(DP), DIMENSION(:), ALLOCATABLE :: vec
    INTEGER(I4B) :: i,n,nn
    n=size(coeffs)
    if (n <= 0) then
       poly_dd=0.0_dp
    else if (n < NPAR_POLY) then
       poly_dd=coeffs(n)
       do i=n-1,1,-1
          poly_dd=x*poly_dd+coeffs(i)
       end do
    else
       allocate(vec(n+1))
       pow=x
       vec(1:n)=coeffs
       do
          vec(n+1)=0.0_dp
          nn=ishft(n+1,-1)
          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
          if (nn == 1) exit
          pow=pow*pow
          n=nn
       end do
       poly_dd=vec(1)
       deallocate(vec)
    end if
  END FUNCTION poly_dd
  !BL
  FUNCTION poly_rc(x,coeffs)
    COMPLEX(SPC), INTENT(IN) :: x
    REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs
    COMPLEX(SPC) :: poly_rc
    COMPLEX(SPC) :: pow
    COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec
    INTEGER(I4B) :: i,n,nn
    n=size(coeffs)
    if (n <= 0) then
       poly_rc=0.0_sp
    else if (n < NPAR_POLY) then
       poly_rc=coeffs(n)
       do i=n-1,1,-1
          poly_rc=x*poly_rc+coeffs(i)
       end do
    else
       allocate(vec(n+1))
       pow=x
       vec(1:n)=coeffs
       do
          vec(n+1)=0.0_sp
          nn=ishft(n+1,-1)
          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
          if (nn == 1) exit
          pow=pow*pow
          n=nn
       end do
       poly_rc=vec(1)
       deallocate(vec)
    end if
  END FUNCTION poly_rc
  !BL
  FUNCTION poly_cc(x,coeffs)
    COMPLEX(SPC), INTENT(IN) :: x
    COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: coeffs
    COMPLEX(SPC) :: poly_cc
    COMPLEX(SPC) :: pow
    COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec
    INTEGER(I4B) :: i,n,nn
    n=size(coeffs)
    if (n <= 0) then
       poly_cc=0.0_sp
    else if (n < NPAR_POLY) then
       poly_cc=coeffs(n)
       do i=n-1,1,-1
          poly_cc=x*poly_cc+coeffs(i)
       end do
    else
       allocate(vec(n+1))
       pow=x
       vec(1:n)=coeffs
       do
          vec(n+1)=0.0_sp
          nn=ishft(n+1,-1)
          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
          if (nn == 1) exit
          pow=pow*pow
          n=nn
       end do
       poly_cc=vec(1)
       deallocate(vec)
    end if
  END FUNCTION poly_cc
  !BL
  FUNCTION poly_rrv(x,coeffs)
    REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x
    REAL(SP), DIMENSION(size(x)) :: poly_rrv
    INTEGER(I4B) :: i,n,m
    m=size(coeffs)
    n=size(x)
    if (m <= 0) then
       poly_rrv=0.0_sp
    else if (m < n .or. m < NPAR_POLY) then
       poly_rrv=coeffs(m)
       do i=m-1,1,-1
          poly_rrv=x*poly_rrv+coeffs(i)
       end do
    else
       do i=1,n
          poly_rrv(i)=poly_rr(x(i),coeffs)
       end do
    end if
  END FUNCTION poly_rrv
  !BL
  FUNCTION poly_ddv(x,coeffs)
    REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x
    REAL(DP), DIMENSION(size(x)) :: poly_ddv
    INTEGER(I4B) :: i,n,m
    m=size(coeffs)
    n=size(x)
    if (m <= 0) then
       poly_ddv=0.0_dp
    else if (m < n .or. m < NPAR_POLY) then
       poly_ddv=coeffs(m)
       do i=m-1,1,-1
          poly_ddv=x*poly_ddv+coeffs(i)
       end do
    else
       do i=1,n
          poly_ddv(i)=poly_dd(x(i),coeffs)
       end do
    end if
  END FUNCTION poly_ddv
  !BL
  FUNCTION poly_msk_rrv(x,coeffs,mask)
    REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x
    LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask
    REAL(SP), DIMENSION(size(x)) :: poly_msk_rrv
    poly_msk_rrv=unpack(poly_rrv(pack(x,mask),coeffs),mask,0.0_sp)
  END FUNCTION poly_msk_rrv
  !BL
  FUNCTION poly_msk_ddv(x,coeffs,mask)
    REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x
    LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask
    REAL(DP), DIMENSION(size(x)) :: poly_msk_ddv
    poly_msk_ddv=unpack(poly_ddv(pack(x,mask),coeffs),mask,0.0_dp)
  END FUNCTION poly_msk_ddv
  !BL
  !BL
  RECURSIVE FUNCTION poly_term_rr(a,b) RESULT(u)
    REAL(SP), DIMENSION(:), INTENT(IN) :: a
    REAL(SP), INTENT(IN) :: b
    REAL(SP), DIMENSION(size(a)) :: u
    INTEGER(I4B) :: n,j
    n=size(a)
    if (n <= 0) RETURN
    u(1)=a(1)
    if (n < NPAR_POLYTERM) then
       do j=2,n
          u(j)=a(j)+b*u(j-1)
       end do
    else
       u(2:n:2)=poly_term_rr(a(2:n:2)+a(1:n-1:2)*b,b*b)
       u(3:n:2)=a(3:n:2)+b*u(2:n-1:2)
    end if
  END FUNCTION poly_term_rr
  !BL
  RECURSIVE FUNCTION poly_term_cc(a,b) RESULT(u)
    COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a
    COMPLEX(SPC), INTENT(IN) :: b
    COMPLEX(SPC), DIMENSION(size(a)) :: u
    INTEGER(I4B) :: n,j
    n=size(a)
    if (n <= 0) RETURN
    u(1)=a(1)
    if (n < NPAR_POLYTERM) then
       do j=2,n
          u(j)=a(j)+b*u(j-1)
       end do
    else
       u(2:n:2)=poly_term_cc(a(2:n:2)+a(1:n-1:2)*b,b*b)
       u(3:n:2)=a(3:n:2)+b*u(2:n-1:2)
    end if
  END FUNCTION poly_term_cc
  !BL
  !BL
  FUNCTION zroots_unity(n,nn)
    INTEGER(I4B), INTENT(IN) :: n,nn
    COMPLEX(SPC), DIMENSION(nn) :: zroots_unity
    INTEGER(I4B) :: k
    REAL(SP) :: theta
    zroots_unity(1)=1.0
    theta=TWOPI/n
    k=1
    do
       if (k >= nn) exit
       zroots_unity(k+1)=cmplx(cos(k*theta),sin(k*theta),SPC)
       zroots_unity(k+2:min(2*k,nn))=zroots_unity(k+1)*&
            zroots_unity(2:min(k,nn-k))
       k=2*k
    end do
  END FUNCTION zroots_unity
  !BL
  FUNCTION outerprod_r(a,b)
    REAL(SP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(SP), DIMENSION(size(a),size(b)) :: outerprod_r
    outerprod_r = spread(a,dim=2,ncopies=size(b)) * &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerprod_r
  !BL
  FUNCTION outerprod_d(a,b)
    REAL(DP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(DP), DIMENSION(size(a),size(b)) :: outerprod_d
    outerprod_d = spread(a,dim=2,ncopies=size(b)) * &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerprod_d
  !BL
  FUNCTION outerdiv(a,b)
    REAL(SP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(SP), DIMENSION(size(a),size(b)) :: outerdiv
    outerdiv = spread(a,dim=2,ncopies=size(b)) / &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerdiv
  !BL
  FUNCTION outersum(a,b)
    REAL(SP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(SP), DIMENSION(size(a),size(b)) :: outersum
    outersum = spread(a,dim=2,ncopies=size(b)) + &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outersum
  !BL
  FUNCTION outerdiff_r(a,b)
    REAL(SP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(SP), DIMENSION(size(a),size(b)) :: outerdiff_r
    outerdiff_r = spread(a,dim=2,ncopies=size(b)) - &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerdiff_r
  !BL
  FUNCTION outerdiff_d(a,b)
    REAL(DP), DIMENSION(:), INTENT(IN) :: a,b
    REAL(DP), DIMENSION(size(a),size(b)) :: outerdiff_d
    outerdiff_d = spread(a,dim=2,ncopies=size(b)) - &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerdiff_d
  !BL
  FUNCTION outerdiff_i(a,b)
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: a,b
    INTEGER(I4B), DIMENSION(size(a),size(b)) :: outerdiff_i
    outerdiff_i = spread(a,dim=2,ncopies=size(b)) - &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerdiff_i
  !BL
  FUNCTION outerand(a,b)
    LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: a,b
    LOGICAL(LGT), DIMENSION(size(a),size(b)) :: outerand
    outerand = spread(a,dim=2,ncopies=size(b)) .and. &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION outerand
  !BL
  SUBROUTINE scatter_add_r(dest,source,dest_index)
    REAL(SP), DIMENSION(:), INTENT(OUT) :: dest
    REAL(SP), DIMENSION(:), INTENT(IN) :: source
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index
    INTEGER(I4B) :: m,n,j,i
    n=assert_eq2(size(source),size(dest_index),'scatter_add_r')
    m=size(dest)
    do j=1,n
       i=dest_index(j)
       if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j)
    end do
  END SUBROUTINE scatter_add_r
  SUBROUTINE scatter_add_d(dest,source,dest_index)
    REAL(DP), DIMENSION(:), INTENT(OUT) :: dest
    REAL(DP), DIMENSION(:), INTENT(IN) :: source
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index
    INTEGER(I4B) :: m,n,j,i
    n=assert_eq2(size(source),size(dest_index),'scatter_add_d')
    m=size(dest)
    do j=1,n
       i=dest_index(j)
       if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j)
    end do
  END SUBROUTINE scatter_add_d
  SUBROUTINE scatter_max_r(dest,source,dest_index)
    REAL(SP), DIMENSION(:), INTENT(OUT) :: dest
    REAL(SP), DIMENSION(:), INTENT(IN) :: source
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index
    INTEGER(I4B) :: m,n,j,i
    n=assert_eq2(size(source),size(dest_index),'scatter_max_r')
    m=size(dest)
    do j=1,n
       i=dest_index(j)
       if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j))
    end do
  END SUBROUTINE scatter_max_r
  SUBROUTINE scatter_max_d(dest,source,dest_index)
    REAL(DP), DIMENSION(:), INTENT(OUT) :: dest
    REAL(DP), DIMENSION(:), INTENT(IN) :: source
    INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index
    INTEGER(I4B) :: m,n,j,i
    n=assert_eq2(size(source),size(dest_index),'scatter_max_d')
    m=size(dest)
    do j=1,n
       i=dest_index(j)
       if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j))
    end do
  END SUBROUTINE scatter_max_d
  !BL
  SUBROUTINE diagadd_rv(mat,diag)
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    REAL(SP), DIMENSION(:), INTENT(IN) :: diag
    INTEGER(I4B) :: j,n
    n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagadd_rv')
    do j=1,n
       mat(j,j)=mat(j,j)+diag(j)
    end do
  END SUBROUTINE diagadd_rv
  !BL
  SUBROUTINE diagadd_r(mat,diag)
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    REAL(SP), INTENT(IN) :: diag
    INTEGER(I4B) :: j,n
    n = min(size(mat,1),size(mat,2))
    do j=1,n
       mat(j,j)=mat(j,j)+diag
    end do
  END SUBROUTINE diagadd_r
  !BL
  SUBROUTINE diagmult_rv(mat,diag)
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    REAL(SP), DIMENSION(:), INTENT(IN) :: diag
    INTEGER(I4B) :: j,n
    n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagmult_rv')
    do j=1,n
       mat(j,j)=mat(j,j)*diag(j)
    end do
  END SUBROUTINE diagmult_rv
  !BL
  SUBROUTINE diagmult_r(mat,diag)
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    REAL(SP), INTENT(IN) :: diag
    INTEGER(I4B) :: j,n
    n = min(size(mat,1),size(mat,2))
    do j=1,n
       mat(j,j)=mat(j,j)*diag
    end do
  END SUBROUTINE diagmult_r
  !BL
  FUNCTION get_diag_rv(mat)
    REAL(SP), DIMENSION(:,:), INTENT(IN) :: mat
    REAL(SP), DIMENSION(size(mat,1)) :: get_diag_rv
    INTEGER(I4B) :: j
    j=assert_eq2(size(mat,1),size(mat,2),'get_diag_rv')
    do j=1,size(mat,1)
       get_diag_rv(j)=mat(j,j)
    end do
  END FUNCTION get_diag_rv
  !BL
  FUNCTION get_diag_dv(mat)
    REAL(DP), DIMENSION(:,:), INTENT(IN) :: mat
    REAL(DP), DIMENSION(size(mat,1)) :: get_diag_dv
    INTEGER(I4B) :: j
    j=assert_eq2(size(mat,1),size(mat,2),'get_diag_dv')
    do j=1,size(mat,1)
       get_diag_dv(j)=mat(j,j)
    end do
  END FUNCTION get_diag_dv
  !BL
  SUBROUTINE put_diag_rv(diagv,mat)
    REAL(SP), DIMENSION(:), INTENT(IN) :: diagv
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    INTEGER(I4B) :: j,n
    n=assert_eq2(size(diagv),min(size(mat,1),size(mat,2)),'put_diag_rv')
    do j=1,n
       mat(j,j)=diagv(j)
    end do
  END SUBROUTINE put_diag_rv
  !BL
  SUBROUTINE put_diag_r(scal,mat)
    REAL(SP), INTENT(IN) :: scal
    REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat
    INTEGER(I4B) :: j,n
    n = min(size(mat,1),size(mat,2))
    do j=1,n
       mat(j,j)=scal
    end do
  END SUBROUTINE put_diag_r
  !BL
  SUBROUTINE unit_matrix(mat)
    REAL(SP), DIMENSION(:,:), INTENT(OUT) :: mat
    INTEGER(I4B) :: i,n
    n=min(size(mat,1),size(mat,2))
    mat(:,:)=0.0_sp
    do i=1,n
       mat(i,i)=1.0_sp
    end do
  END SUBROUTINE unit_matrix
  !BL
  FUNCTION upper_triangle(j,k,extra)
    INTEGER(I4B), INTENT(IN) :: j,k
    INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra
    LOGICAL(LGT), DIMENSION(j,k) :: upper_triangle
    INTEGER(I4B) :: n
    n=0
    if (present(extra)) n=extra
    upper_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) < n)
  END FUNCTION upper_triangle
  !BL
  FUNCTION lower_triangle(j,k,extra)
    INTEGER(I4B), INTENT(IN) :: j,k
    INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra
    LOGICAL(LGT), DIMENSION(j,k) :: lower_triangle
    INTEGER(I4B) :: n
    n=0
    if (present(extra)) n=extra
    lower_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) > -n)
  END FUNCTION lower_triangle
  !BL
  FUNCTION vabs(v)
    REAL(SP), DIMENSION(:), INTENT(IN) :: v
    REAL(SP) :: vabs
    vabs=sqrt(dot_product(v,v))
  END FUNCTION vabs
  !BL
END MODULE nrutil

max_stars_repo_name: aasensio/DeepLearning

max_issues_repo_name: aasensio/DeepLearning

max_forks_repo_name: aasensio/DeepLearning


text:
*DECK CSCAL
      SUBROUTINE CSCAL (N, CA, CX, INCX)
C***BEGIN PROLOGUE  CSCAL
C***PURPOSE  Multiply a vector by a constant.
C***LIBRARY   SLATEC (BLAS)
C***CATEGORY  D1A6
C***TYPE      COMPLEX (SSCAL-S, DSCAL-D, CSCAL-C)
C***KEYWORDS  BLAS, LINEAR ALGEBRA, SCALE, VECTOR
C***AUTHOR  Lawson, C. L., (JPL)
C           Hanson, R. J., (SNLA)
C           Kincaid, D. R., (U. of Texas)
C           Krogh, F. T., (JPL)
C***DESCRIPTION
C
C                B L A S  Subprogram
C    Description of Parameters
C
C     --Input--
C        N  number of elements in input vector(s)
C       CA  complex scale factor
C       CX  complex vector with N elements
C     INCX  storage spacing between elements of CX
C
C     --Output--
C       CX  complex result (unchanged if N .LE. 0)
C
C     Replace complex CX by complex CA*CX.
C     For I = 0 to N-1, replace CX(IX+I*INCX) with CA*CX(IX+I*INCX),
C     where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
C
C***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C                 Krogh, Basic linear algebra subprograms for Fortran
C                 usage, Algorithm No. 539, Transactions on Mathematical
C                 Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   791001  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900821  Modified to correct problem with a negative increment.
C           (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CSCAL
      COMPLEX CA, CX(*)
      INTEGER I, INCX, IX, N
C***FIRST EXECUTABLE STATEMENT  CSCAL
      IF (N .LE. 0) RETURN
C
      IF (INCX .EQ. 1) GOTO 20
C
C     Code for increment not equal to 1.
C
      IX = 1
      IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
      DO 10 I = 1,N
        CX(IX) = CA*CX(IX)
        IX = IX + INCX
   10 CONTINUE
      RETURN
C
C     Code for increment equal to 1.
C
   20 DO 30 I = 1,N
        CX(I) = CA*CX(I)
   30 CONTINUE
      RETURN
      END


!
!     --------------------------------------------------------------
!     R L S P 3 2
!     --------------------------------------------------------------
!                                                                  *
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville             October 1996   *
!
      SUBROUTINE RLSP32(K, JA1, JA2, JA3, K1, K2, K3, K4, KA, IRE, IAT, REC)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONSTS_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  12:24:23  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE dlsa3_I
      USE dlsa4_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: K
      INTEGER  :: JA1
      INTEGER  :: JA2
      INTEGER  :: JA3
      INTEGER  :: K1
      INTEGER  :: K2
      INTEGER  :: K3
      INTEGER  :: K4
      INTEGER  :: KA
      INTEGER  :: IRE
      INTEGER  :: IAT
      REAL(DOUBLE) , INTENT(OUT) :: REC
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      REAL(DOUBLE) :: RE
!-----------------------------------------------
!
      REC = ONE
      IF (JA3 - JA2 > 1) THEN
         IAT = 0
         CALL DLSA3 (K, JA2, JA3, K3, IRE, IAT, RE)
         IF (IAT == 0) RETURN
         REC = RE*REC
      ENDIF
!
      IAT = 0
      CALL DLSA4 (K, JA1, JA2, K1, K2, K3, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
!
      IAT = 0
      CALL DLSA4 (K, JA2, JA3, K3, K4, KA, IRE, IAT, RE)
      IF (IAT == 0) RETURN
      REC = RE*REC
      RETURN
      END SUBROUTINE RLSP32


text:
      subroutine fd_calc_heat(iflg)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Calculate fd coefficients for use with heat conduction terms of
!D1 anisotropic perms
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 3.0
!D2
!D2 Initial implementation: Programmer: George Zyvoloski
!D2 1-May-2010
!D2
!***********************************************************************
!D3
!D3 REQUIREMENTS TRACEABILITY
!D3
!D3 2.3.7 Sources and sinks
!D3 2.6   Provide Input/Output Data Files
!D3 3.0   INPUT AND OUTPUT REQUIREMENTS
!D3
!***********************************************************************
!D4
!D4  SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use comai
      use combi
      use comci
      use comdi
      use comji
      use comdti
      use comki
      use davidi
      implicit none

      integer i,icode,iflg,neqp1,i1,i2,jj,kb,kc,ic,ik,kd
      integer i3,i4,kk,idir,imodel,j,kcmin,nsize,nsize_sx
      integer kb_adv_x, kb_adv_y, kb_adv_z
      real*8 area_i,area_kb,sx2c,dism
      real*8 cord1,cord2,cord3,cord1j,cord2j,cord3j
      real*8 cord1k,cord2k,cord3k
      real*8 dis2,dis2i,disx2, disy2, disz2
      real*8 cosxi0,cosyi0,coszi0,cosx,cosy,cosz
      real*8 cosmin,cosikc,areat,dis2min
      real*8 disikb,dis_totalkb,dis_totali

      neqp1 = neq +1
      if(iflg.eq.1) then
c find areas for diffrence terms to augment  anisotropy
       ic = 0
       nsize = nelm(neqp1)-neqp1
       if(.not.allocated(istrw))then
        allocate(istrw(nsize))
        isoy = 1
        isoz = 1
calculate storage requirements
        nsize_sx = 3*neq
        allocate(sx(nsize_sx,3))
        istrw = 0
       else
c error condition
        if(iout.ne.0)
     &   write(iout,*) 'error condition:istrw previously allocated'
        if(iptty.ne.0)
     &   write(iptty,*) 'error condition:istrw previously allocated'
       endif
       do i = 1, n0
            areat = 0.0
            cord1 = cord(i,1)
            cord2 = cord(i,2)
            cord3 = 0.0
            if(icnl.eq.0) then
             cord3 = cord(i,3)
            endif
            do ik = 1,3
             kd = ncon_adv(i,ik)
             if(kd.gt.0) then
             i1 = nelmdg(i)+1
             i2 = nelm(i+1)
             do jj = i1,i2
              kb = nelm(jj)
              if(kb.eq.kd) go to 1400
             enddo
1400         continue
             if(kb.ne.i) then
              cord1j = cord(kb,1)
              cord2j = cord(kb,2)
              cord3j = 0.0
              if(icnl.eq.0) then
               cord3j = cord(kb,3)
              endif
c calculate distance between i-kb
              disx2 = (cord1-cord1j)**2
              disy2 = (cord2-cord2j)**2
              disz2 = (cord3-cord3j)**2
              dis2i = disx2 + disy2 + disz2
              disikb =sqrt(dis2i)
c find cosine squared
              cosxi0 = disx2/dis2i
              cosyi0 = disy2/dis2i
              coszi0 = disz2/dis2i
c find maximum length in i-kb  (different from i-kb) direction for each gridblock
c gridblock i  (use dot product)
            cosmin = 0.1
            kcmin = 0
            i3 = i1
            i4 = i2
            do kk = i3,i4
             kc = nelm(kk)
             if(kc.ne.i.and.kc.ne.kb) then
              cord1k = cord(kc,1)
              cord2k = cord(kc,2)
              cord3k = 0.0
              if(icnl.eq.0) then
               cord3k = cord(kc,3)
              endif
              disx2 = (cord1-cord1k)**2
	        disy2 = (cord2-cord2k)**2
	        disz2 = (cord3-cord3k)**2
	        dis2 = disx2 + disy2 + disz2
c find cosine squared
	       cosx = disx2/dis2
	       cosy = disy2/dis2
             cosz = disz2/dis2
c compare cosines
             cosikc = abs(cosx-cosxi0)+abs(cosy-cosyi0)+abs(cosz-coszi0)
               if(cosikc.lt.cosmin) then
                cosmin = cosikc
                kcmin = kc
                dis2min = dis2
               endif
             endif
             enddo
c found distance calculate total distance
             dis_totali = 0.5*(sqrt(dis2min) +sqrt(dis2i))
c gridblock kb  (use dot product)
            cosmin = 0.1
            kcmin = 0
            i3 = nelm(kb)+1
            i4 = nelm(kb+1)
            do kk = i3,i4
             kc = nelm(kk)
             if(kc.ne.i.and.kc.ne.kb) then
              cord1k = cord(kc,1)
              cord2k = cord(kc,2)
              cord3k = 0.0
              if(icnl.eq.0) then
               cord3k = cord(kc,3)
              endif
            disx2 = (cord1j-cord1k)**2
	      disy2 = (cord2j-cord2k)**2
	      disz2 = (cord3j-cord3k)**2
	      dis2 = disx2 + disy2 + disz2
c find cosine squared
	       cosx = disx2/dis2
	       cosy = disy2/dis2
             cosz = disz2/dis2
c compare cosines
             cosikc = abs(cosx-cosxi0)+abs(cosy-cosyi0)+abs(cosz-coszi0)
               if(cosikc.lt.cosmin) then
                cosmin = cosikc
                kcmin = kc
                dis2min = dis2
               endif
             endif
c found distance calculate total distance
             dis_totalkb = 0.5*(sqrt(dis2min) +sqrt(dis2i))
            enddo
c with distances and volumes can calculate areas for control volumes
            area_i = sx1(i)/dis_totali/disikb
            area_kb = sx1(kb)/dis_totalkb/disikb
c
            if(kb.gt.i) then
c the "divide by 3" is for isotropic coefficients
             ic = ic +1
             istrw(jj-neqp1) = ic
             sx(ic,1) = -max(area_i,area_kb)/3.
            endif
           endif
           endif
           enddo
       enddo
      endif
      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM