submodule (PDEelliptic) elliptic2d

implicit none (type, external)

contains

module procedure elliptic2D_polarization

!------------------------------------------------------------
!-------SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS
!-------INCLUDES FULL OF POLARIZATION CURRENT, INCLUDING CONVECTIVE
!-------TERMS.  VELOCITIES SHOULD BE TRIMMED (WITHOUT GHOST CELLS).
!-------THIS VERSION OF THE *INTEGRATED* POTENTIAL SOLVER OBVIATES
!-------ALL OTHERS SINCE A PURELY ELECTRSTATIC FORM CAN BE RECOVERED
!-------BY ZEROING OUT THE INERTIAL CAPACITANCE.
!-------
!-------THIS FORM IS INTENDED TO  WORK WITH CURVILINEAR MESHES.
!-------NOTE THAT THE FULL GRID VARIABLES (X%DX3ALL, ETC.) MUST
!-------BE USED HERE!!!
!-------The equation solved by this subroutine is:
!-------
!-------  d/dx2(A dV/dx2) + d/dx3(A' dV/dx3) + B dV/dx2 - C dV/x3 + ...
!-------  d/dx2(D d/dt(dV/dx2)) + d/dx3(D d/dt(dV/dx3)) + ...
!-------  d/dx2( D*v2 d^V/dx2^2 + D*v3*d^2V/dx3/dx2 ) + ...
!-------  d/dx3( D*v2 d^2V/dx2/dx3 + D*v3 d^2V/dx3^2 ) = srcterm
!-------
!-------  for GEMINI:  A=SigP2, A'=SigP3, B=d/dx3(SigH), C=d/dx2(SigH),
!-------               D=Cm
!------------------------------------------------------------

real(wp), dimension(1:size(SigP2,1),1:size(SigP2,2)) :: SigPh2    !I'm too lazy to recode these as SigP2h2, etc.
real(wp), dimension(1:size(SigP2,1),1:size(SigP2,2)) :: SigPh3
real(wp), dimension(1:size(SigP2,1),1:size(SigP2,2)) :: Cmh2
real(wp), dimension(1:size(SigP2,1),1:size(SigP2,2)) :: Cmh3

real(wp) :: coeff    !coefficient for calculating polarization terms
integer :: ix2,ix3,lx2,lx3    !this overwrites the
integer :: lPhi,lent
integer :: iPhi,ient
integer, dimension(:), allocatable :: ir,ic
real(wp), dimension(:), allocatable :: M
real(wp), dimension(:), allocatable :: b
real(wp) :: tstart,tfin

integer :: utrace

type(MUMPS_STRUC) :: mumps_par

!ONLY ROOT NEEDS TO ASSEMBLE THE MATRIX
!if (myid==0) then
lx2=size(SigP2,1)    !note that these are full-grid sizes since grid module globals are not in scope
lx3=size(SigP2,2)
lPhi=lx2*lx3
!      lent=5*(lx2-2)*(lx3-2)+2*lx2+2*(lx3-2)    !static model; left here as a reference
lent=17*(lx2-2)*(lx3-2)+2*lx2+2*(lx3-2)-3*2*(lx2-2)-3*2*(lx3-2)
!! interior+boundary-x3_adj-x2_adj.  Note that are 3 sets of entries for each adjacent point
allocate(ir(lent),ic(lent),M(lent),b(lPhi))
if (debug) print *, 'MUMPS will attempt a solve of size:  ',lx2,lx3
if (debug) print *, 'Total unknowns and nonzero entries in matrix:  ',lPhi,lent

!PREP INPUT DATA FOR SOLUTION OF SYSTEM
SigPh2(1,:)= 0
SigPh2(2:lx2,:)=0.5_wp*(SigP2(1:lx2-1,:)+SigP2(2:lx2,:))
!! note the different conductiances here ot be associated with derivatives in different directions
SigPh3(:,1)= 0
SigPh3(:,2:lx3)=0.5_wp*(SigP3(:,1:lx3-1)+SigP3(:,2:lx3))
Cmh2(1,:)= 0
Cmh2(2:lx2,:)=0.5_wp*(Cm(1:lx2-1,:)+Cm(2:lx2,:))
Cmh3(:,1)= 0
Cmh3(:,2:lx3)=0.5_wp*(Cm(:,1:lx3-1)+Cm(:,2:lx3))

!------------------------------------------------------------
!-------DEFINE A MATRIX USING SPARSE STORAGE (CENTRALIZED
!-------ASSEMBLED MATRIX INPUT, SEE SECTION 4.5 OF MUMPS USER
!-------GUIDE).
!------------------------------------------------------------
!LOAD UP MATRIX ELEMENTS
M(:)= 0
b = pack(srcterm,.true.)           !boundaries overwritten later, polarization terms also added later.
ient=1

loopx3: do ix3=1,lx3
  loopx2: do ix2=1,lx2
    iPhi=lx2*(ix3-1)+ix2     !linear index referencing Phi(ix2,ix3) as a column vector.  Also row of big matrix

    if (ix2==1) then
      !! BOTTOM GRID POINTS + CORNER
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vminx2(ix3)
      ient=ient+1
      cycle
    elseif (ix2==lx2) then
      !! TOP GRID POINTS + CORNER
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vmaxx2(ix3)
      ient=ient+1
      cycle
    elseif (ix3==1) then
      !! LEFT BOUNDARY
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vminx3(ix2)
      ient=ient+1
      cycle
    elseif (ix3==lx3) then
      !! RIGHT BOUNDARY
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vmaxx3(ix2)
      ient=ient+1
      cycle
    endif

    !! INTERIOR LOCATION
    !> ix2-1,ix3-2 grid point
    coeff=-Cm(ix2,ix3-1)*v2(ix2,ix3-1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    if (ix3==2) then    !out of bounds, use nearest BC, and add to known vector
      b(iPhi)=b(iPhi)-coeff*Vminx3(ix2-1)
    else    !in bounds, add to matrix
      ir(ient)=iPhi
      ic(ient)=iPhi-2*lx2-1

      M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

      ient=ient+1
    end if

    !> ix2,ix3-2 grid point
    coeff=-Cm(ix2,ix3-1)*v3(ix2,ix3-1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)*dx3iall(ix3-1)) )
    if (ix3==2) then
    !! bit of intentional code duplication here and in the following sections to keep things organized in a way I can debug...
      b(iPhi)=b(iPhi)-coeff*Vminx3(ix2)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-2*lx2

      M(ient)=coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

      ient=ient+1
    end if

    !> ix2+1,ix3-2 grid point
    coeff=Cm(ix2,ix3-1)*v2(ix2,ix3-1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    if (ix3==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx3(ix2+1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-2*lx2+1

      M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

      ient=ient+1
    end if

    !> ix2-2,ix3-1
    coeff=-Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(ix3-1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-lx2-2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !!!ix2,ix3-1 grid point in ix2,ix3 equation
    ir(ient)=iPhi
    ic(ient)=iPhi-lx2

    M(ient)=SigPh3(ix2,ix3)/(dx3iall(ix3)*dx3all(ix3))+gradSigH2(ix2,ix3)/(dx3all(ix3)+dx3all(ix3+1))    !static terms

    coeff=Cmh3(ix2,ix3)/(dt*dx3iall(ix3)*dx3all(ix3))
    M(ient)=M(ient)+coeff   !polarization time derivative terms
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,ix3-1)
    !! add in polarziation terms that include previous time step potential at this grid point

    coeff=Cm(ix2,ix3-1)*v3(ix2,ix3-1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3)*dx3iall(ix3-1)) )+ &
    Cm(ix2,ix3-1)*v3(ix2,ix3-1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)*dx3iall(ix3-1)) )
    M(ient)=M(ient)+coeff   !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    coeff=Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) )+ &
    Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    M(ient)=M(ient)+coeff   !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

    ient=ient+1

    !> ix2+2,ix3-1 grid point
    coeff=-Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2+1)+dx2all(ix2+2))*(dx2all(ix2)+dx2all(ix2+1))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(ix3-1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-lx2+2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !> ix2-2,ix3 grid point
    coeff=-Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)*dx2iall(ix2-1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(ix3)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-2

      M(ient)=coeff
      !! d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

      ient=ient+1
    end if

    !> ix2-1,ix3 grid point
    ir(ient)=iPhi
    ic(ient)=iPhi-1

    M(ient)=SigPh2(ix2,ix3)/(dx2iall(ix2)*dx2all(ix2))-gradSigH3(ix2,ix3)/(dx2all(ix2)+dx2all(ix2+1))    !static

    coeff=Cmh2(ix2,ix3)/(dt*dx2iall(ix2)*dx2all(ix2))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2-1,ix3)    !BC's and pol. time deriv.

    coeff=Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2)*dx2iall(ix2-1)) )+ &
    Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)*dx2iall(ix2-1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) )+ &
    Cm(ix2,ix3-1)*v2(ix2,ix3-1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    M(ient)=M(ient)+coeff     !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

    ient=ient+1

    !!!ix2,ix3 grid point (main diagonal)
    ir(ient)=iPhi
    ic(ient)=iPhi

    M(ient)=-SigPh2(ix2+1,ix3)/(dx2iall(ix2)*dx2all(ix2+1)) &
    -SigPh2(ix2,ix3)/(dx2iall(ix2)*dx2all(ix2)) &
    -SigPh3(ix2,ix3+1)/(dx3iall(ix3)*dx3all(ix3+1)) &
    -SigPh3(ix2,ix3)/(dx3iall(ix3)*dx3all(ix3))    !static

    coeff=-Cmh2(ix2+1,ix3)/(dt*dx2iall(ix2)*dx2all(ix2+1)) &
    -Cmh2(ix2,ix3)/(dt*dx2iall(ix2)*dx2all(ix2)) &
    -Cmh3(ix2,ix3+1)/(dt*dx3iall(ix3)*dx3all(ix3+1)) &
    -Cmh3(ix2,ix3)/(dt*dx3iall(ix3)*dx3all(ix3))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,ix3)    !BC's and pol. time deriv.

    coeff=Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)*dx2iall(ix2+1)) ) &
    -Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2)*dx2iall(ix2-1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=Cm(ix2,ix3+1)*v3(ix2,ix3+1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)*dx3all(ix3+1)) ) &
    -Cm(ix2,ix3-1)*v3(ix2,ix3-1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3)*dx3iall(ix3-1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    ient=ient+1

    !!!ix2+1,ix3 grid point
    ir(ient)=iPhi
    ic(ient)=iPhi+1

    M(ient)=SigPh2(ix2+1,ix3)/(dx2iall(ix2)*dx2all(ix2+1))+gradSigH3(ix2,ix3)/(dx2all(ix2)+dx2all(ix2+1))    !static

    coeff=Cmh2(ix2+1,ix3)/(dt*dx2iall(ix2)*dx2all(ix2+1))
    M(ient)=M(ient)+coeff    !pol. time deriv. terms
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2+1,ix3)    !BC's and pol. time deriv.

    coeff=-Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+2)*dx2iall(ix2+1)) ) &
    -Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)*dx2iall(ix2+1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=-Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) ) &
    -Cm(ix2,ix3-1)*v2(ix2,ix3-1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

    ient=ient+1

    !ix2+2,ix3 grid point
    coeff=Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+2)*dx2iall(ix2+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(ix3)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+2

      M(ient)=coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

      ient=ient+1
    end if

    !ix2-2,ix3+1 grid point
    coeff=Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(ix3+1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+lx2-2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !!!ix2,ix3+1 grid point
    ir(ient)=iPhi
    ic(ient)=iPhi+lx2

    M(ient)=SigPh3(ix2,ix3+1)/(dx3iall(ix3)*dx3all(ix3+1))-gradSigH2(ix2,ix3)/(dx3all(ix3)+dx3all(ix3+1))    !static

    coeff=Cmh3(ix2,ix3+1)/(dt*dx3iall(ix3)*dx3all(ix3+1))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,ix3+1)    !BC's and pol. time deriv.

    coeff=-Cm(ix2,ix3+1)*v3(ix2,ix3+1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+2)*dx3iall(ix3+1)) ) &
    -Cm(ix2,ix3+1)*v3(ix2,ix3+1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)*dx3iall(ix3+1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    coeff=-Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) ) &
    -Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

    ient=ient+1

    !ix2+2,ix3+1 grid point
    coeff=Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(ix3+1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+lx2+2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2-1,ix3+2 grid point
    coeff=-Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    ( (dx3all(ix3+1)+dx3all(ix3+2))*(dx3all(ix3)+dx3all(ix3+1))*(dx2all(ix2)+dx2all(ix2+1)) )
    if (ix3==lx3-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx3(ix2-1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+2*lx2-1

      M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

      ient=ient+1
    end if

    !ix2,ix3+2 grid point
    coeff=Cm(ix2,ix3+1)*v3(ix2,ix3+1)/( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+2)*dx3iall(ix3+1)) )
    if (ix3==lx3-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx3(ix2)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+2*lx2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2+1,ix3+2 grid point
    coeff=Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) )
    if (ix3==lx3-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx3(ix2+1)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+2*lx2+1

      M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

      ient=ient+1
    end if

  end do loopx2
end do loopx3
!end if

!FIRE UP MUMPS
!if (myid == 0) then
if (debug) print *, 'Filled ',ient-1,' matrix entries.  Initializing MUMPS...'
!end if
mumps_par%COMM = MPI_COMM_WORLD
mumps_par%JOB = -1
mumps_par%SYM = 0
mumps_par%PAR = 1

call MUMPS_exec(mumps_par)

call quiet_mumps(mumps_par)

!LOAD OUR PROBLEM
!if ( myid==0 ) then
mumps_par%N=lPhi
mumps_par%NZ=lent
allocate( mumps_par%IRN ( mumps_par%NZ ) )
allocate( mumps_par%JCN ( mumps_par%NZ ) )
allocate( mumps_par%A( mumps_par%NZ ) )
allocate( mumps_par%RHS ( mumps_par%N  ) )
mumps_par%IRN=ir
mumps_par%JCN=ic
mumps_par%A=M
mumps_par%RHS=b
deallocate(ir,ic,M,b)     !clear memory before solve begins!!!

if (perflag .and. it/=1) then       !used cached permutation
  allocate(mumps_par%PERM_IN(mumps_par%N))
  mumps_par%PERM_IN=mumps_perm
  mumps_par%ICNTL(7)=1
end if

!may solve some memory allocation issues, uncomment if MUMPS throws errors
!about not having enough memory
!mumps_par%ICNTL(14)=50
!end if

!SOLVE (ALL WORKERS NEED TO SEE THIS CALL)
mumps_par%JOB = 6

call MUMPS_exec(mumps_par)

call check_mumps_status(mumps_par, 'elliptic2D_polarization')

!STORE PERMUTATION USED, SAVE RESULTS, CLEAN UP MUMPS ARRAYS
!(can save ~25% execution time and improves scaling with openmpi
! ~25% more going from 1-2 processors)
!if ( myid==0 ) then
if (debug) print *, 'Now organizing results...'

if (perflag .and. it==1) then
  allocate(mumps_perm(mumps_par%N))     !we don't have a corresponding deallocate statement
  mumps_perm=mumps_par%SYM_PERM
end if

elliptic2D_polarization=reshape(mumps_par%RHS,[lx2,lx3])

if (debug) print *, 'Now attempting deallocations...'

deallocate( mumps_par%IRN )
deallocate( mumps_par%JCN )
deallocate( mumps_par%A   )
deallocate( mumps_par%RHS )
!end if

mumps_par%JOB = -2

call MUMPS_exec(mumps_par)

end procedure elliptic2D_polarization

module procedure elliptic2D_polarization_periodic

!------------------------------------------------------------
!-------SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS
!-------INCLUDES FULL OF POLARIZATION CURRENT, INCLUDING CONVECTIVE
!-------TERMS.  VELOCITIES SHOULD BE TRIMMED (WITHOUT GHOST CELLS).
!-------THIS VERSION OF THE *INTEGRATED* POTENTIAL SOLVER OBVIATES
!-------ALL OTHERS SINCE A PURELY ELECTRSTATIC FORM CAN BE RECOVERED
!-------BY ZEROING OUT THE INERTIAL CAPACITANCE.
!-------
!-------THIS FORM IS INTENDED TO  WORK WITH CARTESIAN MESHES ONLY.
!-------NOTE THAT THE FULL GRID VARIABLES (X%DX3ALL, ETC.) MUST
!-------BE USED HERE!!!
!-------
!-------THIS FUNCTION WORKS ON A PERIODIC MESH BY USING A CIRCULANT MATRIX
!-------The equation solved by this subroutine is:
!-------
!-------  d/dx2(A dV/dx2) + d/dx3(A dV/dx3) + B dV/dx2 + C dV/x3 + ...
!-------  d/dx2(D d/dt(dV/dx2)) + d/dx3(D d/dt(dV/dx3)) + ...
!-------  d/dx2( D*v2 d^V/dx2^2 + D*v3*d^2V/dx3/dx2 ) + ...
!-------  d/dx3( D*v2 d^2V/dx2/dx3 + D*v3 d^2V/dx3^2 ) = srcterm
!------------------------------------------------------------

real(wp), dimension(1:size(SigP,1),1:size(SigP,2)) :: SigPh2
real(wp), dimension(1:size(SigP,1),1:size(SigP,2)) :: SigPh3
real(wp), dimension(1:size(SigP,1),1:size(SigP,2)) :: Cmh2
real(wp), dimension(1:size(SigP,1),1:size(SigP,2)) :: Cmh3

real(wp) :: coeff    !coefficient for calculating polarization terms
integer :: ix2,ix3,lx2,lx3    !this overwrites the values stored in the grid module, which is fine, but perhaps redundant
integer :: lPhi,lent
integer :: iPhi,ient
integer, dimension(:), allocatable :: ir,ic
real(wp), dimension(:), allocatable :: M
real(wp), dimension(:), allocatable :: b
real(wp) :: tstart,tfin

type(MUMPS_STRUC) :: mumps_par

integer :: lcount,ix2tmp,ix3tmp

real(wp), dimension(size(SigP,1),size(SigP,2)) :: tmpresults

!ONLY ROOT NEEDS TO ASSEMBLE THE MATRIX
!if (myid==0) then
lx2=size(SigP,1)    !these are full-grid sizes since grid module globals are not in scope
lx3=size(SigP,2)
lPhi=lx2*lx3
!      lent=5*(lx2-2)*(lx3-2)+2*lx2+2*(lx3-2)
!! static model
!      lent=17*(lx2-2)*(lx3A-2)+2*lx2+2*(lx3-2)-3*2*(lx2-2)-3*2*(lx3-2)
!! interior+boundary-x3_adj-x2_adj.  Note that are 3 sets of entries for each adjacent point.
!! The x3 adjacent points do not need to be removed in teh case of periodic boundary conditions.
!! This is technicall correct, but a bit misleading, I think.
!! Shouldn't it be lent=17*(lx2-2)*(lx3-2)+2*(lx2-2)+2*lx3-3*2*(lx2-2)-3*2*(lx3-2)???????
!      lent=17*(lx2-2)*(lx3+1-2)+2*(lx3+1)+3*(lx2-2)-3*2*(lx3+1-2)
!! true interior with x3 boundaries which are not treated as interior in periodic solves
!! + add x2 boundaries (note that these are now size lx3+1) + 3 entries for each x3 ghost cell that we are adding
!! - x2_adj (x2 is not periodici, two sets of three points each, note again the larger x3 size as compared to aperiodic solutions).
lent=17*(lx2-2)*(lx3-2)+2*(lx3)+17*2*(lx2-2)-3*2*lx3
!! true interior with x3 boundaries which are not treated as interior in periodic solves
!! + add x2 boundaries (note that these are now size lx3+1) + x3 edge cells (treated here as interior)
!! - x2_adj (x2 is not periodici, two sets of three points each, note again the larger x3 size as compared to aperiodic solutions).
allocate(ir(lent),ic(lent),M(lent),b(lPhi))

if (debug) print *, 'MUMPS will attempt a solve of size:  ',lx2,lx3
if (debug) print *, 'Total unknowns and nonzero entries in matrix:  ',lPhi,lent

!NOTE THAT THESE NEED TO BE PERIODIC IN X3
SigPh2(1,:)= 0
SigPh2(2:lx2,:)=0.5_wp*(SigP(1:lx2-1,:)+SigP(2:lx2,:))
SigPh3(:,1)=0.5_wp*(SigP(:,lx3)+SigP(:,1))    !needs to be left interface value so average of first and last grid point
SigPh3(:,2:lx3)=0.5_wp*(SigP(:,1:lx3-1)+SigP(:,2:lx3))
Cmh2(1,:)= 0
Cmh2(2:lx2,:)=0.5_wp*(Cm(1:lx2-1,:)+Cm(2:lx2,:))
Cmh3(:,1)=0.5_wp*(Cm(:,lx3)+Cm(:,1))
Cmh3(:,2:lx3)=0.5_wp*(Cm(:,1:lx3-1)+Cm(:,2:lx3))

!------------------------------------------------------------
!-------DEFINE A MATRIX USING SPARSE STORAGE (CENTRALIZED
!-------ASSEMBLED MATRIX INPUT, SEE SECTION 4.5 OF MUMPS USER
!-------GUIDE).
!------------------------------------------------------------
if (debug) print *, 'Loading up matrix entries...'

!LOAD UP MATRIX ELEMENTS
lcount=0
M(:)= 0
b=pack(srcterm,.true.)           !boundaries overwritten later, polarization terms also added later.
ient=1
do ix3=1,lx3
!! note that we have one extra ghost cell now to deal with due to the way we've chosen to implement periodic boundary conditions
  do ix2=1,lx2
    iPhi=lx2*(ix3-1)+ix2
    !! linear index referencing Phi(ix2,ix3) as a column vector.  Also row of big matrix

    if (ix2==1) then
      !! BOTTOM GRID POINTS + CORNER
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vminx2(ix3)
      ient=ient+1
      cycle
    elseif (ix2==lx2) then
      !! TOP GRID POINTS + CORNER
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1
      b(iPhi)=Vmaxx2(ix3)
      ient=ient+1
      cycle
    endif

    !! TREAT AS AN INTERIOR LOCATION, THIS INCLUDE X3 EDGES NOW SINCE PERIODIC,CIRCULANT

    !! ZZZ - NEED TO WRAP INDICES AROUND:  X 1)
    !! matrix row/column entries; 2) references to dx3i*(anything but ix3);  3) references to conductances/bcs/etc.

    !ix2-1,ix3-2 grid point
    coeff=-Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v2(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    ir(ient)=iPhi
    ix3tmp=mod(ix3-2-1+lx3,lx3)+1
    ix2tmp=ix2-1
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    !              ic(ient)=iPhi-2*lx2-1+lPhi-1
    !! add the grid size to wrap the index around (would be negative otherwise),
    !! -1 at end because the last grid opint is actually the same as the first for our implementation
    !            else
    !              ic(ient)=iPhi-2*lx2-1
    !            end if
    M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )
    ient=ient+1

    !ix2,ix3-2 grid point
    coeff=-Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v3(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)* &
    dx3iall(mod(ix3-1-1+lx3,lx3)+1)) )
    ir(ient)=iPhi
    ix3tmp=mod(ix3-2-1+lx3,lx3)+1
    ix2tmp=ix2
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    !              ic(ient)=iPhi-2*lx2+lPhi-1    !again add grid size to wrap to end
    !            else
    !              ic(ient)=iPhi-2*lx2
    !            end if
    M(ient)=coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term
    ient=ient+1

    !ix2+1,ix3-2 grid point
    coeff=Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v2(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    ir(ient)=iPhi
    ix3tmp=mod(ix3-2-1+lx3,lx3)+1
    ix2tmp=ix2+1
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )
    ient=ient+1

    !ix2-2,ix3-1
    coeff=-Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(mod(ix3-1-1+lx3,lx3)+1)
    else
      ir(ient)=iPhi
      ix3tmp=mod(ix3-1-1+lx3,lx3)+1
      ix2tmp=ix2-2
      ic(ient)=lx2*(ix3tmp-1)+ix2tmp
      !              ic(ient)=iPhi-lx2-2

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2,ix3-1 grid point in ix2,ix3 equation
    ir(ient)=iPhi
    !            ic(ient)=iPhi-lx2
    ix3tmp=mod(ix3-1-1+lx3,lx3)+1
    ix2tmp=ix2
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp

    M(ient)=SigPh3(ix2,ix3)/(dx3iall(ix3)*dx3all(ix3))+gradSigH2(ix2,ix3)/(dx3all(ix3)+dx3all(ix3+1))    !static terms

    coeff=Cmh3(ix2,ix3)/(dt*dx3iall(ix3)*dx3all(ix3))
    M(ient)=M(ient)+coeff   !polarization time derivative terms
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,mod(ix3-1-1+lx3,lx3)+1)
    !! add in polarziation terms that include previous time step potential at this grid point

    coeff=Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v3(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3)*dx3iall(mod(ix3-1-1+lx3,lx3)+1)) )+ &
    Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v3(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)*dx3iall(mod(ix3-1-1+lx3,lx3)+1)) )
    M(ient)=M(ient)+coeff   !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    coeff=Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) )+ &
    Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    M(ient)=M(ient)+coeff   !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

    ient=ient+1

    !ix2+2,ix3-1 grid point
    coeff=-Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2+1)+dx2all(ix2+2))*(dx2all(ix2)+dx2all(ix2+1))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(mod(ix3-1-1+lx3,lx3)+1)
    else
      ir(ient)=iPhi
      !              ic(ient)=iPhi-lx2+2
      ix3tmp=mod(ix3-1-1+lx3,lx3)+1
      ix2tmp=ix2+2
      ic(ient)=lx2*(ix3tmp-1)+ix2tmp

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2-2,ix3 grid point
    coeff=-Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)*dx2iall(ix2-1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(ix3)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi-2

      M(ient)=coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

      ient=ient+1
    end if

    !ix2-1,ix3 grid point
    ir(ient)=iPhi
    ic(ient)=iPhi-1

    M(ient)=SigPh2(ix2,ix3)/(dx2iall(ix2)*dx2all(ix2))-gradSigH3(ix2,ix3)/(dx2all(ix2)+dx2all(ix2+1))    !static

    coeff=Cmh2(ix2,ix3)/(dt*dx2iall(ix2)*dx2all(ix2))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2-1,ix3)    !BC's and pol. time deriv.

    coeff=Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2)*dx2iall(ix2-1)) )+ &
    Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)*dx2iall(ix2-1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=Cm(ix2,mod(ix3+1-1,lx3)+1)*v2(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) )+ &
    Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v2(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    M(ient)=M(ient)+coeff     !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

    ient=ient+1

    !ix2,ix3 grid point (main diagonal)
    ir(ient)=iPhi
    ic(ient)=iPhi

    M(ient)=-SigPh2(ix2+1,ix3)/(dx2iall(ix2)*dx2all(ix2+1)) &
    -SigPh2(ix2,ix3)/(dx2iall(ix2)*dx2all(ix2)) &
    -SigPh3(ix2,mod(ix3+1-1,lx3)+1)/(dx3iall(ix3)*dx3all(ix3+1)) &
    -SigPh3(ix2,ix3)/(dx3iall(ix3)*dx3all(ix3))    !static

    coeff=-Cmh2(ix2+1,ix3)/(dt*dx2iall(ix2)*dx2all(ix2+1)) &
    -Cmh2(ix2,ix3)/(dt*dx2iall(ix2)*dx2all(ix2)) &
    -Cmh3(ix2,mod(ix3+1-1,lx3)+1)/(dt*dx3iall(ix3)*dx3all(ix3+1)) &
    -Cmh3(ix2,ix3)/(dt*dx3iall(ix3)*dx3all(ix3))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,ix3)    !BC's and pol. time deriv.

    coeff=Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)*dx2iall(ix2+1)) ) &
    -Cm(ix2-1,ix3)*v2(ix2-1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2)*dx2iall(ix2-1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=Cm(ix2,mod(ix3+1-1,lx3)+1)*v3(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)*dx3all(ix3+1)) ) &
    -Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v3(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3)*dx3iall(mod(ix3-1-1+lx3,lx3)+1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    ient=ient+1

    !ix2+1,ix3 grid point
    ir(ient)=iPhi
    ic(ient)=iPhi+1

    M(ient)=SigPh2(ix2+1,ix3)/(dx2iall(ix2)*dx2all(ix2+1))+gradSigH3(ix2,ix3)/(dx2all(ix2)+dx2all(ix2+1))    !static

    coeff=Cmh2(ix2+1,ix3)/(dt*dx2iall(ix2)*dx2all(ix2+1))
    M(ient)=M(ient)+coeff    !pol. time deriv. terms
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2+1,ix3)    !BC's and pol. time deriv.

    coeff=-Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+2)*dx2iall(ix2+1)) ) &
    -Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)*dx2iall(ix2+1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

    coeff=-Cm(ix2,mod(ix3+1-1,lx3)+1)*v2(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) ) &
    -Cm(ix2,mod(ix3-1-1+lx3,lx3)+1)*v2(ix2,mod(ix3-1-1+lx3,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3-1)+dx3all(ix3))*(dx2all(ix2)+dx2all(ix2+1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )

    ient=ient+1

    !ix2+2,ix3 grid point
    coeff=Cm(ix2+1,ix3)*v2(ix2+1,ix3)/( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+2)*dx2iall(ix2+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(ix3)
    else
      ir(ient)=iPhi
      ic(ient)=iPhi+2

      M(ient)=coeff    !d/dx2( Cm*v2*d^2/dx2^2(Phi) ) term

      ient=ient+1
    end if

    !ix2-2,ix3+1 grid point
    coeff=Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==2) then
      b(iPhi)=b(iPhi)-coeff*Vminx2(mod(ix3+1-1,lx3)+1)
    else
      ir(ient)=iPhi
      !              ic(ient)=iPhi+lx2-2
      ix3tmp=mod(ix3+1-1,lx3)+1
      ix2tmp=ix2-2
      ic(ient)=lx2*(ix3tmp-1)+ix2tmp

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2,ix3+1 grid point
    ir(ient)=iPhi
    !            ic(ient)=iPhi+lx2
    ix3tmp=mod(ix3+1-1,lx3)+1
    ix2tmp=ix2
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp

    M(ient)=SigPh3(ix2,mod(ix3+1-1,lx3)+1)/ &
    (dx3iall(ix3)*dx3all(ix3+1))-gradSigH2(ix2,ix3)/(dx3all(ix3)+dx3all(ix3+1))    !static

    coeff=Cmh3(ix2,mod(ix3+1-1,lx3)+1)/(dt*dx3iall(ix3)*dx3all(ix3+1))
    M(ient)=M(ient)+coeff    !pol. time deriv.
    b(iPhi)=b(iPhi)+coeff*Phi0(ix2,mod(ix3+1-1,lx3)+1)    !BC's and pol. time deriv.

    coeff=-Cm(ix2,mod(ix3+1-1,lx3)+1)*v3(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+2)*dx3iall(mod(ix3+1-1,lx3)+1)) ) &
    -Cm(ix2,mod(ix3+1-1,lx3)+1)*v3(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)*dx3iall(mod(ix3+1-1,lx3)+1)) )
    M(ient)=M(ient)+coeff    !d/dx3( Cm*v3*d^2/dx3^2(Phi) ) term

    coeff=-Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) ) &
    -Cm(ix2-1,ix3)*v3(ix2-1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2-1)+dx2all(ix2))*(dx3all(ix3)+dx3all(ix3+1)) )
    M(ient)=M(ient)+coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

    ient=ient+1

    !ix2+2,ix3+1 grid point
    coeff=Cm(ix2+1,ix3)*v3(ix2+1,ix3)/ &
    ( (dx2all(ix2)+dx2all(ix2+1))*(dx2all(ix2+1)+dx2all(ix2+2))*(dx3all(ix3)+dx3all(ix3+1)) )
    if (ix2==lx2-1) then
      b(iPhi)=b(iPhi)-coeff*Vmaxx2(mod(ix3+1-1,lx3)+1)
    else
      ir(ient)=iPhi
      !              ic(ient)=iPhi+lx2+2
      ix3tmp=mod(ix3+1-1,lx3)+1
      ix2tmp=ix2+2
      ic(ient)=lx2*(ix3tmp-1)+ix2tmp

      M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )

      ient=ient+1
    end if

    !ix2-1,ix3+2 grid point
    !            coeff=-Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    !                  ( (dx3all(ix3+1)+dx3all(ix3+2))*(dx3all(ix3)+dx3all(ix3+1))*(dx2all(ix2)+dx2all(ix2+1)) )
    coeff=-Cm(ix2,mod(ix3+1-1,lx3)+1)*v2(ix2,mod(ix3+1-1,lx3)+1)/ &      !mods to wrap indices around
    ( (dx3all(ix3+1)+dx3all(ix3+2))*(dx3all(ix3)+dx3all(ix3+1))*(dx2all(ix2)+dx2all(ix2+1)) )
    ir(ient)=iPhi
    !            if (ix3>=lx3-1) then    !this needs to also handle the case where ix3=lx3!!!  Likewise for statements that follow...
    ix3tmp=mod(ix3+2-1,lx3)+1
    ix2tmp=ix2-1
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    !              ic(ient)=iPhi+2*lx2-1-lPhi+1    !wrap to beginning
    !            else
    !              ic(ient)=iPhi+2*lx2-1
    !            end if
    M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )
    ient=ient+1

    !ix2,ix3+2 grid point
    !            coeff=Cm(ix2,ix3+1)*v3(ix2,ix3+1)/ &
    !                  ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+2)*dx3iall(ix3+1)) )
    coeff=Cm(ix2,mod(ix3+1-1,lx3)+1)*v3(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+2)*dx3iall(mod(ix3+1-1,lx3)+1)) )
    ir(ient)=iPhi
    !            if (ix3>=lx3-1) then
    ix3tmp=mod(ix3+2-1,lx3)+1
    ix2tmp=ix2
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    !              ic(ient)=iPhi+2*lx2-lPhi+1    !substract grid size to wrap around to the beginning
    !            else
    !              ic(ient)=iPhi+2*lx2
    !            end if
    M(ient)=coeff    !d/dx2( Cm*v3*d^2/dx2dx3(Phi) )
    ient=ient+1

    !ix2+1,ix3+2 grid point
    !            coeff=Cm(ix2,ix3+1)*v2(ix2,ix3+1)/ &
    !                  ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) )
    coeff=Cm(ix2,mod(ix3+1-1,lx3)+1)*v2(ix2,mod(ix3+1-1,lx3)+1)/ &
    ( (dx3all(ix3)+dx3all(ix3+1))*(dx3all(ix3+1)+dx3all(ix3+2))*(dx2all(ix2)+dx2all(ix2+1)) )
    ir(ient)=iPhi
    !            if (ix3>=lx3-1) then     !this should actually work for any case...
    ix3tmp=mod(ix3+2-1,lx3)+1
    ix2tmp=ix2+1
    ic(ient)=lx2*(ix3tmp-1)+ix2tmp
    !              ic(ient)=iPhi+2*lx2+1-lPhi+1
    !            else
    !              ic(ient)=iPhi+2*lx2+1
    !            end if
    M(ient)=coeff    !d/dx3( Cm*v2*d^2/dx3dx2(Phi) )
    ient=ient+1
  end do
end do
!end if

!FIRE UP MUMPS
!if (myid == 0) then
if (debug) print *,  'Debug count:  ',lcount
if (debug) print *, 'Filled ',ient-1,' out of ',lent,' matrix entries for solving ',iPhi,' of ',lPhi, &
' unknowns.  Initializing MUMPS...'
if (ient-1 /= lent) error stop 'Incorrect number of matrix entries filled in potential solve!!!'

!end if
mumps_par%COMM = MPI_COMM_WORLD
mumps_par%JOB = -1
mumps_par%SYM = 0
mumps_par%PAR = 1

call MUMPS_exec(mumps_par)

call quiet_mumps(mumps_par)

!LOAD OUR PROBLEM
!if ( myid==0 ) then
mumps_par%N=lPhi
mumps_par%NZ=lent
allocate( mumps_par%IRN ( mumps_par%NZ ) )
allocate( mumps_par%JCN ( mumps_par%NZ ) )
allocate( mumps_par%A( mumps_par%NZ ) )
allocate( mumps_par%RHS ( mumps_par%N  ) )
mumps_par%IRN=ir
mumps_par%JCN=ic
mumps_par%A=M
mumps_par%RHS=b
deallocate(ir,ic,M,b)     !clear memory before solve begins!!!

if (perflag .and. it/=1) then       !used cached permutation
  allocate(mumps_par%PERM_IN(mumps_par%N))
  mumps_par%PERM_IN=mumps_perm
  mumps_par%ICNTL(7)=1
end if

!mumps_par%ICNTL(14)=50
!end if

!SOLVE (ALL WORKERS NEED TO SEE THIS CALL)
mumps_par%JOB = 6

call MUMPS_exec(mumps_par)

call check_mumps_status(mumps_par, 'elliptic2D_polarization_periodic')

!STORE PERMUTATION USED, SAVE RESULTS, CLEAN UP MUMPS ARRAYS
!(can save ~25% execution time and improves scaling with openmpi
! ~25% more going from 1-2 processors)
!if ( myid==0 ) then
if (debug) print *, 'Now organizing results...'

if (perflag .and. it==1) then
  allocate(mumps_perm(mumps_par%N))     !we don't have a corresponding deallocate statement
  mumps_perm=mumps_par%SYM_PERM
end if

!IF WE HAVE DONE A PERIODIC SOLVE, THE LAST GRID POINT NEEDS TO BE IGNORED WHEN WE RESHAPE THE POTENTIAL ARRAY.

tmpresults=reshape(mumps_par%RHS,[lx2,lx3])
elliptic2D_polarization_periodic=tmpresults(1:lx2,1:lx3)    !sort of superfluous now that hte solve size is the same as the grid

if (debug) print *, 'Now attempting deallocations...'

deallocate( mumps_par%IRN )
deallocate( mumps_par%JCN )
deallocate( mumps_par%A   )
deallocate( mumps_par%RHS )
!end if

mumps_par%JOB = -2

call MUMPS_exec(mumps_par)

end procedure elliptic2D_polarization_periodic

module procedure elliptic2D_cart
!! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS.
!!
!! ASSUME THAT:
!! * WE ARE RESOLVING THE POTENTIAL ALONG THE FIELD LINE
!! * POTENTIAL VARIES IN X1 AND X3. X2 IS NOMINALLY JUST ONE ELEMENT.
!! * LEFT AND RIGHT BOUNDARIES (IN X3) USE DIRICHLET BOUNARY CONDITIONS
!! * TOP (ALTITUDE) CAN BE NEUMANN OR DIRICHLET.
!! * BOTTOM (ALTITUDE) IS ALWAYS DIRICHLET.
!!
!! This subroutine solves equations of the form:
!!
!!    d/dx1(sig0 dV/dx1) + d/dx3(sigP dV/dx3) = srcterm

real(wp), dimension(1:size(sig0,1),1:size(sig0,3)) :: sig0h1
real(wp), dimension(1:size(sigP,1),1:size(sigP,3)) :: sigPh3

integer :: ix1,ix3,lx1,lx3
integer :: lPhi,lent
integer :: iPhi,ient
integer, dimension(:), allocatable :: ir,ic
real(wp), dimension(:), allocatable :: M
real(wp), dimension(:), allocatable :: b
real(wp) :: tstart,tfin

type(MUMPS_STRUC) :: mumps_par

!ONLY ROOT NEEDS TO ASSEMBLE THE MATRIX
!if (myid==0) then
lx1=size(sig0,1)
lx3=size(sig0,3)
lPhi=lx1*lx3
if (flagdirich==0) then
  lent=5*(lx1-2)*(lx3-2)+2*lx1+2*(lx3-2)+2*lx3    !first +1 for Neumann bottom, second for Neumann top
else
  !        lent=5*(lx1-2)*(lx3-2)+2*lx1+2*(lx3-2)+1    !first +1 for Neumann bottom
  lent=5*(lx1-2)*(lx3-2)+2*(lx1-2)+2*lx3+lx3
end if
allocate(ir(lent),ic(lent),M(lent),b(lPhi))

if (debug) print *, 'MUMPS will attempt a solve of size:  ',lx1,lx3
if (debug) print *, 'Total unknowns and nonzero entries in matrix:  ',lPhi,lent

!AVERAGE CONDUCTANCES TO THE GRID INTERFACE POINTS
sig0h1(1,:)=0
sig0h1(2:lx1,:)=0.5_wp*(sig0(1:lx1-1,1,:)+sig0(2:lx1,1,:))
sigPh3(:,1)=0
sigPh3(:,2:lx3)=0.5_wp*(sigP(:,1,1:lx3-1)+sigP(:,1,2:lx3))

!------------------------------------------------------------
!-------DEFINE A MATRIX USING SPARSE STORAGE (CENTRALIZED
!-------ASSEMBLED MATRIX INPUT, SEE SECTION 4.5 OF MUMPS USER
!-------GUIDE).
!------------------------------------------------------------

!LOAD UP MATRIX ELEMENTS
M(:)=0
b=pack(srcterm,.true.)           !boundaries overwritten later
ient=1
do ix3=1,lx3
  do ix1=1,lx1
    iPhi=lx1*(ix3-1)+ix1     !linear index referencing Phi(ix1,ix3) as a column vector.  Also row of big matrix

    if (ix1==1) then
    !! (LOGICAL) BOTTOM GRID POINTS + CORNER, USE NEUMANN HERE, PRESUMABLY ZERO.
    !! ZZZ - potential problem here if we have inverted grid...
      if (gridflag/=1) then
        ir(ient)=iPhi
        ic(ient)=iPhi
        M(ient)=-1
        ient=ient+1
        ir(ient)=iPhi
        ic(ient)=iPhi+1
        M(ient)=1
        !            b(iPhi)=Vminx1(ix3)
        b(iPhi)= 0    !force bottom current to zero
        ient=ient+1
      else
        if (flagdirich/=0) then    !ZZZ - need to check non-inverted???
          ir(ient)=iPhi
          ic(ient)=iPhi
          M(ient)=1
          b(iPhi)=Vmaxx1(1,ix3)
          ient=ient+1
        else
          ir(ient)=iPhi
          ic(ient)=iPhi-1
          M(ient)=-1/dx1(lx1)
          ient=ient+1
          ir(ient)=iPhi
          ic(ient)=iPhi
          M(ient)=1/dx1(lx1)
          b(iPhi)=Vmaxx1(1,ix3)
          ient=ient+1
        end if
      end if
    elseif (ix1==lx1) then    !(LOGICAL) TOP GRID POINTS + CORNER
      if (gridflag/=1) then    !non-inverted?
        if (flagdirich/=0) then    !ZZZ - need to check non-inverted???
          ir(ient)=iPhi
          ic(ient)=iPhi
          M(ient)=1
          b(iPhi)=Vmaxx1(1,ix3)
          ient=ient+1
        else
          ir(ient)=iPhi
          ic(ient)=iPhi-1
          M(ient)=-1/dx1(lx1)
          ient=ient+1
          ir(ient)=iPhi
          ic(ient)=iPhi
          M(ient)=1/dx1(lx1)
          b(iPhi)=Vmaxx1(1,ix3)
          ient=ient+1
        end if
      else
        ir(ient)=iPhi
        ic(ient)=iPhi
        M(ient)=-1
        ient=ient+1
        ir(ient)=iPhi
        ic(ient)=iPhi+1
        M(ient)=1
        !            b(iPhi)=Vminx1(ix3)
        b(iPhi)=0
        !! force bottom current to zero
        ient=ient+1
      end if
    elseif (ix3==1) then      !LEFT BOUNDARY
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1.0
      b(iPhi)=Vminx3(ix1,1)
      ient=ient+1
    elseif (ix3==lx3) then    !RIGHT BOUNDARY
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=1.0
      b(iPhi)=Vmaxx3(ix1,1)
      ient=ient+1
    else                      !INTERIOR
      !ix1,ix3-1 grid point in ix1,ix3 equation
      ir(ient)=iPhi
      ic(ient)=iPhi-lx1
      M(ient)=sigPh3(ix1,ix3)/(dx3iall(ix3)*dx3all(ix3))
      ient=ient+1

      !ix1-1,ix3 grid point
      ir(ient)=iPhi
      ic(ient)=iPhi-1
      M(ient)=sig0h1(ix1,ix3)/(dx1i(ix1)*dx1(ix1))
      ient=ient+1

      !ix1,ix3 grid point
      ir(ient)=iPhi
      ic(ient)=iPhi
      M(ient)=-sig0h1(ix1+1,ix3)/(dx1i(ix1)*dx1(ix1+1)) &
      -sig0h1(ix1,ix3)/(dx1i(ix1)*dx1(ix1)) &
      -sigPh3(ix1,ix3+1)/(dx3iall(ix3)*dx3all(ix3+1)) &
      -sigPh3(ix1,ix3)/(dx3iall(ix3)*dx3all(ix3))
      ient=ient+1

      !ix1+1,ix3 grid point
      ir(ient)=iPhi
      ic(ient)=iPhi+1
      M(ient)=sig0h1(ix1+1,ix3)/(dx1i(ix1)*dx1(ix1+1))
      ient=ient+1

      !ix1,ix3+1 grid point
      ir(ient)=iPhi
      ic(ient)=iPhi+lx1
      M(ient)=sigPh3(ix1,ix3+1)/(dx3iall(ix3)*dx3all(ix3+1))
      ient=ient+1
    end if
  end do
end do
!end if
if (debug) print *, 'Number of entries used:  ',ient-1

!FIRE UP MUMPS
mumps_par%COMM = MPI_COMM_WORLD
mumps_par%JOB = -1
mumps_par%SYM = 0
mumps_par%PAR = 1

call MUMPS_exec(mumps_par)

call quiet_mumps(mumps_par)

!LOAD OUR PROBLEM
!if ( myid==0 ) then
mumps_par%N=lPhi
mumps_par%NZ=lent
allocate( mumps_par%IRN ( mumps_par%NZ ) )
allocate( mumps_par%JCN ( mumps_par%NZ ) )
allocate( mumps_par%A( mumps_par%NZ ) )
allocate( mumps_par%RHS ( mumps_par%N  ) )
mumps_par%IRN=ir
mumps_par%JCN=ic
mumps_par%A=M
mumps_par%RHS=b
deallocate(ir,ic,M,b)     !clear memory before solve begins!!!

if (perflag .and. it/=1) then       !used cached permutation
  if (debug) print *, 'Using a previously stored permutation'
  allocate(mumps_par%PERM_IN(mumps_par%N))
  mumps_par%PERM_IN = mumps_perm
  mumps_par%ICNTL(7) = 1
end if

!may solve some memory allocation issues, uncomment if MUMPS throws errors
!about not having enough memory
! e.g. INFO(1) = -9
! however this error may also mean there is a deeper problem with the code.
!,mumps_par%ICNTL(14) = 50
!end if

!> SOLVE (ALL WORKERS NEED TO SEE THIS CALL)
mumps_par%JOB = 6

call MUMPS_exec(mumps_par)

call check_mumps_status(mumps_par, 'elliptic2D_cart')

!STORE PERMUTATION USED, SAVE RESULTS, CLEAN UP MUMPS ARRAYS
!(can save ~25% execution time and improves scaling with openmpi
! ~25% more going from 1-2 processors).  WOW - this halves execution
! time on some big 2048*2048 solves!!!
!if ( myid==0 ) then
if (debug) print *, 'Now organizing results...'

if (perflag .and. it==1) then
  if (debug) print *, 'Storing ordering for future time step use...'
  allocate(mumps_perm(mumps_par%N))     !we don't have a corresponding deallocate statement
  mumps_perm=mumps_par%SYM_PERM
end if

elliptic2D_cart=reshape(mumps_par%RHS,[lx1,1,lx3])

if (debug) print *, 'Now attempting deallocations...'

deallocate( mumps_par%IRN )
deallocate( mumps_par%JCN )
deallocate( mumps_par%A   )
deallocate( mumps_par%RHS )
!end if

mumps_par%JOB = -2

call MUMPS_exec(mumps_par)

end procedure elliptic2D_cart

end submodule elliptic2d


program square_root
implicit none

real(kind=8) :: A, root, temp_root, int1, int2, dec1, dec2, n_decimal
integer :: i, j, N, integer_part

write(*,*) 'number?'
read(*,*) A

!-------------------------------------------
! For the integer part of the square root:
!-------------------------------------------

N = nint(A) ! just because we need an integer for the do loop below
do i = 0, N
   int1 = dble(i)
   int2 = int1 + 1.0d0
   if (int2*int2 .gt. A) then
      if (int1*int1 .le. A) then
         root = int1 !at this point we have the integer part of the root
         exit
      end if
   end if
end do

!------------------------------------------
! For the decimal part of the square root:
!------------------------------------------

n_decimal = 1.0d0
do j = 1, 16 !level of precision
   n_decimal = n_decimal/10.0 !moving to the right by one decimal place
   write(*,'(A5I4)') 'j', j
   do i = 0, 9 !the digit at the j-th decimal place can have values
                  ! from 0 to 9
      dec1 = root + dble(i)*n_decimal
      dec2 = root + dble(i + 1)*n_decimal
      write(*,'(I3)', advance='no') i
      write(*,*) dec1, dec2
      if (dec2*dec2 .gt. A) then
         if(dec1*dec1 .le. A) then
            root = dec1 ! and now we have the square root with precision j
            exit
         end if
      end if
   end do
   write(*,*)
end do

write(*,*) root

end program square_root


text:
! wrapftemplates_std.f
! This file is generated by Shroud nowrite-version. Do not edit.
! Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
! other Shroud Project Developers.
! See the top-level COPYRIGHT file for details.
!
! SPDX-License-Identifier: (BSD-3-Clause)
!
!>
!! \file wrapftemplates_std.f
!! \brief Shroud generated wrapper for std namespace
!<
! splicer begin namespace.std.file_top
! splicer end namespace.std.file_top
module templates_std_mod
    use iso_c_binding, only : C_INT, C_NULL_PTR, C_PTR
    ! splicer begin namespace.std.module_use
    ! splicer end namespace.std.module_use
    implicit none

    ! splicer begin namespace.std.module_top
    ! splicer end namespace.std.module_top

    ! helper capsule_data_helper
    type, bind(C) :: TEM_SHROUD_capsule_data
        type(C_PTR) :: addr = C_NULL_PTR  ! address of C++ memory
        integer(C_INT) :: idtor = 0       ! index of destructor
    end type TEM_SHROUD_capsule_data

    type vector_int
        type(TEM_SHROUD_capsule_data) :: cxxmem
        ! splicer begin namespace.std.class.vector_int.component_part
        ! splicer end namespace.std.class.vector_int.component_part
    contains
        procedure :: dtor => vector_int_dtor
        procedure :: push_back => vector_int_push_back
        procedure :: at => vector_int_at
        procedure :: get_instance => vector_int_get_instance
        procedure :: set_instance => vector_int_set_instance
        procedure :: associated => vector_int_associated
        ! splicer begin namespace.std.class.vector_int.type_bound_procedure_part
        ! splicer end namespace.std.class.vector_int.type_bound_procedure_part
    end type vector_int

    type vector_double
        type(TEM_SHROUD_capsule_data) :: cxxmem
        ! splicer begin namespace.std.class.vector_double.component_part
        ! splicer end namespace.std.class.vector_double.component_part
    contains
        procedure :: dtor => vector_double_dtor
        procedure :: push_back => vector_double_push_back
        procedure :: at => vector_double_at
        procedure :: get_instance => vector_double_get_instance
        procedure :: set_instance => vector_double_set_instance
        procedure :: associated => vector_double_associated
        ! splicer begin namespace.std.class.vector_double.type_bound_procedure_part
        ! splicer end namespace.std.class.vector_double.type_bound_procedure_part
    end type vector_double

    interface operator (.eq.)
        module procedure vector_int_eq
        module procedure vector_double_eq
    end interface

    interface operator (.ne.)
        module procedure vector_int_ne
        module procedure vector_double_ne
    end interface

    interface

        ! ----------------------------------------
        ! Function:  vector
        ! Attrs:     +intent(ctor)
        ! Exact:     c_function_shadow_scalar
        subroutine c_vector_int_ctor(SHT_rv) &
                bind(C, name="TEM_vector_int_ctor")
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(OUT) :: SHT_rv
        end subroutine c_vector_int_ctor

        ! ----------------------------------------
        ! Function:  ~vector
        ! Attrs:     +intent(dtor)
        ! Requested: c_subroutine_void_scalar
        ! Match:     c_subroutine
        subroutine c_vector_int_dtor(self) &
                bind(C, name="TEM_vector_int_dtor")
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(INOUT) :: self
        end subroutine c_vector_int_dtor

        ! ----------------------------------------
        ! Function:  void push_back
        ! Attrs:     +intent(subroutine)
        ! Requested: c_subroutine_void_scalar
        ! Match:     c_subroutine
        ! ----------------------------------------
        ! Argument:  const int & value +intent(in)
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_&
        ! Match:     c_default
        subroutine c_vector_int_push_back(self, value) &
                bind(C, name="TEM_vector_int_push_back")
            use iso_c_binding, only : C_INT
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            integer(C_INT), intent(IN) :: value
        end subroutine c_vector_int_push_back

        ! ----------------------------------------
        ! Function:  int & at
        ! Attrs:     +deref(pointer)+intent(function)
        ! Requested: c_function_native_&_pointer
        ! Match:     c_function_native_&
        ! ----------------------------------------
        ! Argument:  size_type n +value
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_scalar
        ! Match:     c_default
        function c_vector_int_at(self, n) &
                result(SHT_rv) &
                bind(C, name="TEM_vector_int_at")
            use iso_c_binding, only : C_PTR, C_SIZE_T
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            integer(C_SIZE_T), value, intent(IN) :: n
            type(C_PTR) SHT_rv
        end function c_vector_int_at

        ! ----------------------------------------
        ! Function:  int & at
        ! Attrs:     +api(buf)+deref(pointer)+intent(function)
        ! Requested: c_function_native_&_buf_pointer
        ! Match:     c_function_native_&
        ! ----------------------------------------
        ! Argument:  size_type n +value
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_scalar
        ! Match:     c_default
        function c_vector_int_at_bufferify(self, n) &
                result(SHT_rv) &
                bind(C, name="TEM_vector_int_at_bufferify")
            use iso_c_binding, only : C_PTR, C_SIZE_T
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            integer(C_SIZE_T), value, intent(IN) :: n
            type(C_PTR) SHT_rv
        end function c_vector_int_at_bufferify

        ! splicer begin namespace.std.class.vector_int.additional_interfaces
        ! splicer end namespace.std.class.vector_int.additional_interfaces

        ! ----------------------------------------
        ! Function:  vector
        ! Attrs:     +intent(ctor)
        ! Exact:     c_function_shadow_scalar
        subroutine c_vector_double_ctor(SHT_rv) &
                bind(C, name="TEM_vector_double_ctor")
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(OUT) :: SHT_rv
        end subroutine c_vector_double_ctor

        ! ----------------------------------------
        ! Function:  ~vector
        ! Attrs:     +intent(dtor)
        ! Requested: c_subroutine_void_scalar
        ! Match:     c_subroutine
        subroutine c_vector_double_dtor(self) &
                bind(C, name="TEM_vector_double_dtor")
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(INOUT) :: self
        end subroutine c_vector_double_dtor

        ! ----------------------------------------
        ! Function:  void push_back
        ! Attrs:     +intent(subroutine)
        ! Requested: c_subroutine_void_scalar
        ! Match:     c_subroutine
        ! ----------------------------------------
        ! Argument:  const double & value +intent(in)
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_&
        ! Match:     c_default
        subroutine c_vector_double_push_back(self, value) &
                bind(C, name="TEM_vector_double_push_back")
            use iso_c_binding, only : C_DOUBLE
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            real(C_DOUBLE), intent(IN) :: value
        end subroutine c_vector_double_push_back

        ! ----------------------------------------
        ! Function:  double & at
        ! Attrs:     +deref(pointer)+intent(function)
        ! Requested: c_function_native_&_pointer
        ! Match:     c_function_native_&
        ! ----------------------------------------
        ! Argument:  size_type n +value
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_scalar
        ! Match:     c_default
        function c_vector_double_at(self, n) &
                result(SHT_rv) &
                bind(C, name="TEM_vector_double_at")
            use iso_c_binding, only : C_PTR, C_SIZE_T
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            integer(C_SIZE_T), value, intent(IN) :: n
            type(C_PTR) SHT_rv
        end function c_vector_double_at

        ! ----------------------------------------
        ! Function:  double & at
        ! Attrs:     +api(buf)+deref(pointer)+intent(function)
        ! Requested: c_function_native_&_buf_pointer
        ! Match:     c_function_native_&
        ! ----------------------------------------
        ! Argument:  size_type n +value
        ! Attrs:     +intent(in)
        ! Requested: c_in_native_scalar
        ! Match:     c_default
        function c_vector_double_at_bufferify(self, n) &
                result(SHT_rv) &
                bind(C, name="TEM_vector_double_at_bufferify")
            use iso_c_binding, only : C_PTR, C_SIZE_T
            import :: TEM_SHROUD_capsule_data
            implicit none
            type(TEM_SHROUD_capsule_data), intent(IN) :: self
            integer(C_SIZE_T), value, intent(IN) :: n
            type(C_PTR) SHT_rv
        end function c_vector_double_at_bufferify

        ! splicer begin namespace.std.class.vector_double.additional_interfaces
        ! splicer end namespace.std.class.vector_double.additional_interfaces

        ! splicer begin namespace.std.additional_interfaces
        ! splicer end namespace.std.additional_interfaces
    end interface

    interface vector_double
        module procedure vector_double_ctor
    end interface vector_double

    interface vector_int
        module procedure vector_int_ctor
    end interface vector_int

contains

    ! ----------------------------------------
    ! Function:  vector
    ! Attrs:     +intent(ctor)
    ! Exact:     f_ctor
    ! Attrs:     +intent(ctor)
    ! Exact:     c_ctor
    function vector_int_ctor() &
            result(SHT_rv)
        type(vector_int) :: SHT_rv
        ! splicer begin namespace.std.class.vector_int.method.ctor
        call c_vector_int_ctor(SHT_rv%cxxmem)
        ! splicer end namespace.std.class.vector_int.method.ctor
    end function vector_int_ctor

    ! ----------------------------------------
    ! Function:  ~vector
    ! Attrs:     +intent(dtor)
    ! Exact:     f_dtor
    ! Attrs:     +intent(dtor)
    ! Exact:     c_dtor
    subroutine vector_int_dtor(obj)
        class(vector_int) :: obj
        ! splicer begin namespace.std.class.vector_int.method.dtor
        call c_vector_int_dtor(obj%cxxmem)
        ! splicer end namespace.std.class.vector_int.method.dtor
    end subroutine vector_int_dtor

    ! Generated by cxx_template
    ! ----------------------------------------
    ! Function:  void push_back
    ! Attrs:     +intent(subroutine)
    ! Exact:     f_subroutine
    ! Attrs:     +intent(subroutine)
    ! Exact:     c_subroutine
    ! ----------------------------------------
    ! Argument:  const int & value +intent(in)
    ! Attrs:     +intent(in)
    ! Requested: f_in_native_&
    ! Match:     f_default
    ! Attrs:     +intent(in)
    ! Requested: c_in_native_&
    ! Match:     c_default
    subroutine vector_int_push_back(obj, value)
        use iso_c_binding, only : C_INT
        class(vector_int) :: obj
        integer(C_INT), intent(IN) :: value
        ! splicer begin namespace.std.class.vector_int.method.push_back
        call c_vector_int_push_back(obj%cxxmem, value)
        ! splicer end namespace.std.class.vector_int.method.push_back
    end subroutine vector_int_push_back

    ! Generated by cxx_template - arg_to_buffer
    ! ----------------------------------------
    ! Function:  int & at
    ! Attrs:     +deref(pointer)+intent(function)
    ! Exact:     f_function_native_&_buf_pointer
    ! Attrs:     +api(buf)+deref(pointer)+intent(function)
    ! Requested: c_function_native_&_buf_pointer
    ! Match:     c_function_native_&
    ! ----------------------------------------
    ! Argument:  size_type n +value
    ! Attrs:     +intent(in)
    ! Requested: f_in_native_scalar
    ! Match:     f_default
    ! Attrs:     +intent(in)
    ! Requested: c_in_native_scalar
    ! Match:     c_default
    function vector_int_at(obj, n) &
            result(SHT_rv)
        use iso_c_binding, only : C_INT, C_PTR, C_SIZE_T, c_f_pointer
        class(vector_int) :: obj
        integer(C_SIZE_T), value, intent(IN) :: n
        integer(C_INT), pointer :: SHT_rv
        ! splicer begin namespace.std.class.vector_int.method.at
        type(C_PTR) :: SHC_rv_ptr
        SHC_rv_ptr = c_vector_int_at_bufferify(obj%cxxmem, n)
        call c_f_pointer(SHC_rv_ptr, SHT_rv)
        ! splicer end namespace.std.class.vector_int.method.at
    end function vector_int_at

    ! Return pointer to C++ memory.
    function vector_int_get_instance(obj) result (cxxptr)
        use iso_c_binding, only: C_PTR
        class(vector_int), intent(IN) :: obj
        type(C_PTR) :: cxxptr
        cxxptr = obj%cxxmem%addr
    end function vector_int_get_instance

    subroutine vector_int_set_instance(obj, cxxmem)
        use iso_c_binding, only: C_PTR
        class(vector_int), intent(INOUT) :: obj
        type(C_PTR), intent(IN) :: cxxmem
        obj%cxxmem%addr = cxxmem
        obj%cxxmem%idtor = 0
    end subroutine vector_int_set_instance

    function vector_int_associated(obj) result (rv)
        use iso_c_binding, only: c_associated
        class(vector_int), intent(IN) :: obj
        logical rv
        rv = c_associated(obj%cxxmem%addr)
    end function vector_int_associated

    ! splicer begin namespace.std.class.vector_int.additional_functions
    ! splicer end namespace.std.class.vector_int.additional_functions

    ! ----------------------------------------
    ! Function:  vector
    ! Attrs:     +intent(ctor)
    ! Exact:     f_ctor
    ! Attrs:     +intent(ctor)
    ! Exact:     c_ctor
    function vector_double_ctor() &
            result(SHT_rv)
        type(vector_double) :: SHT_rv
        ! splicer begin namespace.std.class.vector_double.method.ctor
        call c_vector_double_ctor(SHT_rv%cxxmem)
        ! splicer end namespace.std.class.vector_double.method.ctor
    end function vector_double_ctor

    ! ----------------------------------------
    ! Function:  ~vector
    ! Attrs:     +intent(dtor)
    ! Exact:     f_dtor
    ! Attrs:     +intent(dtor)
    ! Exact:     c_dtor
    subroutine vector_double_dtor(obj)
        class(vector_double) :: obj
        ! splicer begin namespace.std.class.vector_double.method.dtor
        call c_vector_double_dtor(obj%cxxmem)
        ! splicer end namespace.std.class.vector_double.method.dtor
    end subroutine vector_double_dtor

    ! Generated by cxx_template
    ! ----------------------------------------
    ! Function:  void push_back
    ! Attrs:     +intent(subroutine)
    ! Exact:     f_subroutine
    ! Attrs:     +intent(subroutine)
    ! Exact:     c_subroutine
    ! ----------------------------------------
    ! Argument:  const double & value +intent(in)
    ! Attrs:     +intent(in)
    ! Requested: f_in_native_&
    ! Match:     f_default
    ! Attrs:     +intent(in)
    ! Requested: c_in_native_&
    ! Match:     c_default
    subroutine vector_double_push_back(obj, value)
        use iso_c_binding, only : C_DOUBLE
        class(vector_double) :: obj
        real(C_DOUBLE), intent(IN) :: value
        ! splicer begin namespace.std.class.vector_double.method.push_back
        call c_vector_double_push_back(obj%cxxmem, value)
        ! splicer end namespace.std.class.vector_double.method.push_back
    end subroutine vector_double_push_back

    ! Generated by cxx_template - arg_to_buffer
    ! ----------------------------------------
    ! Function:  double & at
    ! Attrs:     +deref(pointer)+intent(function)
    ! Exact:     f_function_native_&_buf_pointer
    ! Attrs:     +api(buf)+deref(pointer)+intent(function)
    ! Requested: c_function_native_&_buf_pointer
    ! Match:     c_function_native_&
    ! ----------------------------------------
    ! Argument:  size_type n +value
    ! Attrs:     +intent(in)
    ! Requested: f_in_native_scalar
    ! Match:     f_default
    ! Attrs:     +intent(in)
    ! Requested: c_in_native_scalar
    ! Match:     c_default
    function vector_double_at(obj, n) &
            result(SHT_rv)
        use iso_c_binding, only : C_DOUBLE, C_PTR, C_SIZE_T, c_f_pointer
        class(vector_double) :: obj
        integer(C_SIZE_T), value, intent(IN) :: n
        real(C_DOUBLE), pointer :: SHT_rv
        ! splicer begin namespace.std.class.vector_double.method.at
        type(C_PTR) :: SHC_rv_ptr
        SHC_rv_ptr = c_vector_double_at_bufferify(obj%cxxmem, n)
        call c_f_pointer(SHC_rv_ptr, SHT_rv)
        ! splicer end namespace.std.class.vector_double.method.at
    end function vector_double_at

    ! Return pointer to C++ memory.
    function vector_double_get_instance(obj) result (cxxptr)
        use iso_c_binding, only: C_PTR
        class(vector_double), intent(IN) :: obj
        type(C_PTR) :: cxxptr
        cxxptr = obj%cxxmem%addr
    end function vector_double_get_instance

    subroutine vector_double_set_instance(obj, cxxmem)
        use iso_c_binding, only: C_PTR
        class(vector_double), intent(INOUT) :: obj
        type(C_PTR), intent(IN) :: cxxmem
        obj%cxxmem%addr = cxxmem
        obj%cxxmem%idtor = 0
    end subroutine vector_double_set_instance

    function vector_double_associated(obj) result (rv)
        use iso_c_binding, only: c_associated
        class(vector_double), intent(IN) :: obj
        logical rv
        rv = c_associated(obj%cxxmem%addr)
    end function vector_double_associated

    ! splicer begin namespace.std.class.vector_double.additional_functions
    ! splicer end namespace.std.class.vector_double.additional_functions

    ! splicer begin namespace.std.additional_functions
    ! splicer end namespace.std.additional_functions

    function vector_int_eq(a,b) result (rv)
        use iso_c_binding, only: c_associated
        type(vector_int), intent(IN) ::a,b
        logical :: rv
        if (c_associated(a%cxxmem%addr, b%cxxmem%addr)) then
            rv = .true.
        else
            rv = .false.
        endif
    end function vector_int_eq

    function vector_int_ne(a,b) result (rv)
        use iso_c_binding, only: c_associated
        type(vector_int), intent(IN) ::a,b
        logical :: rv
        if (.not. c_associated(a%cxxmem%addr, b%cxxmem%addr)) then
            rv = .true.
        else
            rv = .false.
        endif
    end function vector_int_ne

    function vector_double_eq(a,b) result (rv)
        use iso_c_binding, only: c_associated
        type(vector_double), intent(IN) ::a,b
        logical :: rv
        if (c_associated(a%cxxmem%addr, b%cxxmem%addr)) then
            rv = .true.
        else
            rv = .false.
        endif
    end function vector_double_eq

    function vector_double_ne(a,b) result (rv)
        use iso_c_binding, only: c_associated
        type(vector_double), intent(IN) ::a,b
        logical :: rv
        if (.not. c_associated(a%cxxmem%addr, b%cxxmem%addr)) then
            rv = .true.
        else
            rv = .false.
        endif
    end function vector_double_ne

end module templates_std_mod

ext: f

max_stars_repo_name: ExternalRepositories/shroud

max_issues_repo_name: ExternalRepositories/shroud

max_forks_repo_name: ExternalRepositories/shroud


text:
      subroutine jacobi(a,n,np,d,v,nrot)
C-----------------------------------------------------------------------
C     Subroutine used to find the eigenvalues and eigenvectors of a
C     SYMMETRIC matrix, a, with physical dimension, np, and logical
C     dimension, n.
C     Eigenvalues stored in d.
C     Normalised eigenvectors formed and stored in v.
C     nrot gives the total number of steps.
C     Taken from Numerical Recipes Online book: www.nr.com.
C     Slightly modified to deal with double precision reals.
C-----------------------------------------------------------------------
      integer n,np,nrot,nmax
      real*8 a(np,np),d(np),v(np,np)
      parameter (nmax=600)
      integer i,ip,iq,j
      real*8 c,g,h,s,sm,t,tau,theta,tresh,b(nmax),z(nmax)
      do ip=1,n
         do iq=1,n
            v(ip,iq)=0.d0
         enddo
         v(ip,ip)=1.d0
      enddo
      do ip=1,n
         b(ip)=a(ip,ip)
         d(ip)=b(ip)
         z(ip)=0.d0
      enddo
      nrot=0
      do i=1,50
         sm=0.d0
         do ip=1,n-1
            do iq=ip+1,n
               sm=sm+abs(a(ip,iq))
            enddo
         enddo
         if(sm.eq.0.d0)return
         if(i.lt.4)then
            tresh=0.2d0*sm/(dble(n))**2
         else
            tresh=0.d0
         endif
         do ip=1,n-1
            do iq=ip+1,n
               g=100.d0*abs(a(ip,iq))
               if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip)))
     .            .and.(abs(d(iq))+g.eq.abs(d(iq))))then
                  a(ip,iq)=0.d0
               elseif(abs(a(ip,iq)).gt.tresh)then
                  h=d(iq)-d(ip)
                  if(abs(h)+g.eq.abs(h))then
                     t=a(ip,iq)/h
                  else
                     theta=0.5d0*h/a(ip,iq)
                     t=1.d0/(abs(theta)+sqrt(1.d0+theta**2))
                     if(theta.lt.0.d0)t=-t
                  endif
                  c=1.d0/sqrt(1.d0+t**2)
                  s=t*c
                  tau=s/(1.d0+c)
                  h=t*a(ip,iq)
                  z(ip)=z(ip)-h
                  z(iq)=z(iq)+h
                  d(ip)=d(ip)-h
                  d(iq)=d(iq)+h
                  a(ip,iq)=0.d0
                  do j=1,ip-1
                     g=a(j,ip)
                     h=a(j,iq)
                     a(j,ip)=g-s*(h+g*tau)
                     a(j,iq)=h+s*(g-h*tau)
                  enddo
                  do j=ip+1,iq-1
                     g=a(ip,j)
                     h=a(j,iq)
                     a(ip,j)=g-s*(h+g*tau)
                     a(j,iq)=h+s*(g-h*tau)
                  enddo
                  do j=iq+1,n
                     g=a(ip,j)
                     h=a(iq,j)
                     a(ip,j)=g-s*(h+g*tau)
                     a(iq,j)=h+s*(g-h*tau)
                  enddo
                  do j=1,n
                     g=v(j,ip)
                     h=v(j,iq)
                     v(j,ip)=g-s*(h+g*tau)
                     v(j,iq)=h+s*(g-h*tau)
                  enddo
                  nrot=nrot+1
               endif
            enddo
         enddo
         do ip=1,n
            b(ip)=b(ip)+z(ip)
            d(ip)=b(ip)
            z(ip)=0.d0
         enddo
      enddo
      write(6,'(''Too many iterations in jacobi'')')
      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
Loshandra Ostrava is a piano and guitar (and other instruments?) area music teachers music instructor/coach. She is also a yoga therapy/instructor/coach, and does shiatsu therapy.

She can be reached at Mailto(loshan9@yahoo.com, mail), (530)4005860, (530)7538356. She accepts Davis Dollars.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE residue_par (gcr, gcz, gcl)
      USE vmec_main, p5 => cp5
      USE vmec_params, ONLY: rss, zcs, rsc, zcc,                               &
                             meven, modd, ntmax, signgs
      USE realspace, ONLY: phip
      USE xstuff
      USE precon2d
      USE parallel_include_module
      USE parallel_vmec_module, ONLY: tlglob_arr, trglob_arr,                  &
                                      lactive, SAXLASTNTYPE
      USE blocktridiagonalsolver, ONLY: L_COLSCALE

      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(0:ntor,0:mpol1,ns,ntmax), INTENT(INOUT) ::           &
         gcr, gcz, gcl
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: n0=0, m0=0, m1=1
      INTEGER, PARAMETER :: n3d=0, nasym=1
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: nsfix, jedge, delIter
      REAL(dp) :: r1, fac, tmp, tmp2(ns), ftotal

      INTEGER  :: i, j, k, l, m, blksize, left, right
      INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
      INTEGER :: MPI_STAT(MPI_STATUS_SIZE)
      REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: send_buf
      REAL(dp), ALLOCATABLE, DIMENSION(:) :: recv_buf
      REAL(dp) :: tredon, tredoff
!-----------------------------------------------
      CALL second0 (treson)

#ifdef _HBANGLE
!FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
      IF (lfreeb .AND. lrfp) THEN
         fac = 0
         IF (ictrl_prec2d .EQ. 0) THEN
            fac = 1.E-1_dp
         END IF
         gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
      END IF
#else
!
!     SYMMETRIC PERTURBATIONS (BASED ON POLAR RELATIONS):
!        Rss(n) = Zcs(n), n != 0
!     ASYMMETRIC PERTURBATIONS:
!        Rsc(n) = Zcc(n), ALL n
!
!     INTERNALLY:
!        XC(rss) = .5*(Rss + Zcs), XC(zcs) = .5*(Rss - Zcs) -> 0
!        XC(rsc) = .5*(Rsc + Zcc), XC(zcc) = .5*(Rsc - Zcc) -> 0
!     THIS IMPLIES THE CONSTRAINT
!        3D ONLY : GC(zcs) = 0;
!        ASYM:     GC(zcc) = 0
!

      IF (lthreed) THEN
         CALL constrain_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
      END IF
      IF (lasym) THEN
         CALL constrain_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
      END IF

      IF (lfreeb .AND. lrfp) THEN
         fac = 0
         IF (ictrl_prec2d .EQ. 0) THEN
            fac = 1.E-1_dp
         END IF
         gcr(0,m0,ns,:) = fac*gcr(0,m0,ns,:)
         gcz(0,m0,ns,:) = fac*gcz(0,m0,ns,:)
      END IF
#endif

!     PRECONDITIONER MUST BE CALCULATED USING RAW (UNPRECONDITIONED) FORCES
      IF (ictrl_prec2d .GE. 2 .OR. ictrl_prec2d .EQ. -1) RETURN

!
!     COMPUTE INVARIANT RESIDUALS
!
      r1 = one/(2*r0scale)**2
      jedge = 0
      delIter = iter2-iter1

      IF (delIter       .lt. 50 .and.                                        &
          (fsqr + fsqz) .LT. 1.E-6_dp) THEN
         jedge = 1
      ENDIF

      CALL getfsq_par (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)

      CALL second0(tredon)
      tmp = SUM(gcl(:,:,tlglob:trglob,:)*gcl(:,:,tlglob:trglob,:))
      CALL MPI_Allreduce(tmp,ftotal,1,MPI_REAL8,MPI_SUM,NS_COMM,MPI_ERR)
      CALL second0(tredoff)
      allreduce_time = allreduce_time + (tredoff - tredon)
      fsql = fnormL*ftotal
      IF(rank .EQ. nranks-1) THEN
         fedge = r1*fnorm*SUM(gcr(:,:,ns,:)**2 + gcz(:,:,ns,:)**2)
      END IF
!
!     PERFORM PRECONDITIONING AND COMPUTE RESIDUES
!
      IF (ictrl_prec2d .EQ. 1) THEN

         IF (l_colscale .AND. lactive) THEN
            CALL SAXLASTNTYPE(pgc, pcol_scale, pgc)
         END IF

         LRESIDUECALL = .TRUE.
         CALL block_precond_par(pgc)
         LRESIDUECALL = .FALSE.

         IF (.NOT.lfreeb .AND. ANY(gcr(:,:,ns,:) .NE. zero)) THEN
            STOP 'gcr(ns) != 0 for fixed boundary in residue'
         END IF
         IF (.NOT.lfreeb .AND. ANY(gcz(:,:,ns,:) .NE. zero)) THEN
            STOP 'gcz(ns) != 0 for fixed boundary in residue'
         END IF
         IF (ANY(gcl(1:,m0,:,zsc) .NE. zero)) THEN
            STOP 'gcl(m=0,n>0,sc) != 0 in residue'
         END IF
         IF (lthreed .AND. ANY(gcl(n0,:,:,zcs) .NE. zero)) THEN
            STOP 'gcl(n=0,m,cs) != 0 in residue'
         END IF

         fsqr1 = SUM(gcr*gcr)
         fsqz1 = SUM(gcz*gcz)
         fsql1 = SUM(gcl*gcl)

      ELSE
!        m = 1 constraint scaling

         IF (lthreed) THEN
            CALL scale_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
         END IF
         IF (lasym) THEN
            CALL scale_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
         END IF

         jedge = 0
         CALL scalfor_par (gcr, arm, brm, ard, brd, crd, jedge)
         jedge = 1
         CALL scalfor_par (gcz, azm, bzm, azd, bzd, crd, jedge)

         CALL getfsq_par (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)

         DO l = tlglob, trglob
            gcl(:,:,l,:) = pfaclam(:,:,l,:)*gcl(:,:,l,:)
            tmp2(l) = SUM(gcl(:,:,l,:)**2)
         END DO
         CALL Gather1XArray(tmp2)
         ftotal = SUM(tmp2(2:ns))
         fsql1 = hs*ftotal

         CALL PadSides(pgc)

      ENDIF

      CALL second0 (tresoff)
      residue_time = residue_time + (tresoff-treson)

      END SUBROUTINE residue_par

      SUBROUTINE constrain_m1_par(gcr, gcz)
      USE vmec_main
      USE parallel_include_module
      USE precon2d, ONLY: ictrl_prec2d
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(0:ntor,ns), INTENT(INOUT) :: gcr, gcz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      REAL(dp), PARAMETER :: FThreshold = 1.E-6_dp
      REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: temp
!-----------------------------------------------
!
!     COMPUTE INTERNAL gr, gz
!     NOTE: internal gz => 0 for both values of lconm1 (although gz is different)
!     FOR lconm1=T, gcr(internal) = gcr+gcz, gcz(internal) = gcr-gcz->0
!
      ALLOCATE(temp(0:ntor,ns))
      IF (lconm1) THEN
         temp(:,tlglob:trglob) = gcr(:,tlglob:trglob)
         gcr(:,tlglob:trglob) = osqrt2*(gcr(:,tlglob:trglob) +                 &
                                        gcz(:,tlglob:trglob))
         gcz(:,tlglob:trglob) = osqrt2*(temp(:,tlglob:trglob) -                &
                                        gcz(:,tlglob:trglob))
      END IF

!v8.50: ADD iter2<2 so reset=<WOUT_FILE> works
      IF (fsqz         .LT. FThreshold .OR.                                    &
     &    iter2        .LT. 2          .OR.                                    &
     &    ictrl_prec2d .NE. 0) THEN
         gcz(:,tlglob:trglob) = 0
      END IF

      DEALLOCATE(temp)
      END SUBROUTINE constrain_m1_par

      SUBROUTINE scale_m1_par(gcr, gcz)
      USE vmec_main
      USE parallel_include_module
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(0:ntor,ns), INTENT(inout) :: gcr, gcz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: nodd=2
      INTEGER :: n
      REAL(dp) :: fac(ns)
!-----------------------------------------------
      IF (.not.lconm1) RETURN

      fac(tlglob:trglob) = (ard(tlglob:trglob,nodd) +                          &
                            brd(tlglob:trglob,nodd))                           &
                         / (ard(tlglob:trglob,nodd) +                          &
                            brd(tlglob:trglob,nodd) +                          &
                            azd(tlglob:trglob,nodd) +                          &
                            bzd(tlglob:trglob,nodd))
      DO n = 0, ntor
         gcr(n,tlglob:trglob) = fac(tlglob:trglob)*gcr(n,tlglob:trglob)
      END DO

      fac(tlglob:trglob) = (azd(tlglob:trglob,nodd) +                          &
                            bzd(tlglob:trglob,nodd))                           &
                         / (ard(tlglob:trglob,nodd) +                          &
                            brd(tlglob:trglob,nodd) +                          &
                            azd(tlglob:trglob,nodd) +                          &
                            bzd(tlglob:trglob,nodd))
      DO n = 0, ntor
         gcz(n,tlglob:trglob) = fac(tlglob:trglob)*gcz(n,tlglob:trglob)
      END DO

      END SUBROUTINE scale_m1_par

      SUBROUTINE residue (gcr, gcz, gcl)
      USE vmec_main, p5 => cp5
      USE vmec_params, ONLY: rss, zcs, rsc, zcc,                               &
                             meven, modd, ntmax, signgs
      USE realspace, ONLY: phip
      USE xstuff
      USE precon2d
#ifdef _HBANGLE
      USE angle_constraints, ONLY: scalfor_rho
#endif
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(ns,0:ntor,0:mpol1,ntmax), INTENT(inout) :: &
        gcr, gcz, gcl
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: n0=0, m0=0, m1=1
      INTEGER, PARAMETER :: n3d=0, nasym=1
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: nsfix, jedge, delIter
      REAL(dp) :: r1, fac
      INTEGER  :: i, j, k, l
!-----------------------------------------------
!
!     IMPOSE M=1 MODE CONSTRAINT TO MAKE THETA ANGLE
!     INVARIANT TO PHI-SHIFTS (AND THETA SHIFTS FOR ASYMMETRIC CASE)
!     (ZCS = RSS, ZSS = RCS ARE THE CORRECT POLAR RELATIONS)
!

#ifdef _HBANGLE
!FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
      IF (lfreeb .AND. lrfp) THEN
         fac = 0
         IF (ictrl_prec2d .EQ. 0) THEN
            fac = 1.E-1_dp
         END IF
         gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
      END IF
#else
!
!     SYMMETRIC PERTURBATIONS (BASED ON POLAR RELATIONS):
!        RSS(n) = ZCS(n), n != 0
!     ASYMMETRIC PERTURBATIONS:
!        RSC(n) = ZCC(n), ALL n
!
!     INTERNALLY:
!        XC(rss) = .5*(Rss + Zcs), XC(zcs) = .5*(Rss - Zcs) -> 0
!        XC(rsc) = .5*(Rsc + Zcc), XC(zcc) = .5*(Rsc - Zcc) -> 0
!     THIS IMPLIES THE CONSTRAINT
!        3D ONLY : GC(zcs) = 0;
!        ASYM:     GC(zcc) = 0
!

      IF (lthreed) THEN
         CALL constrain_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
      END IF
      IF (lasym) THEN
         CALL constrain_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
      END IF

!FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
      IF (lfreeb .AND. lrfp) THEN
         fac = 0
         IF (ictrl_prec2d .EQ. 0) THEN
            fac = 1.E-1_dp
         END IF
         gcr(ns,0,m0,:) = fac*gcr(ns,0,m0,:)
         gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
      END IF
#endif

!     PRECONDITIONER MUST BE CALCULATED USING RAW (UNPRECONDITIONED) FORCES
      IF (ictrl_prec2d .GE. 2 .OR.                                             &
          ictrl_prec2d .EQ. -1) THEN
         RETURN
      END IF

!
!     COMPUTE INVARIANT RESIDUALS
!
      r1 = one/(2*r0scale)**2
      jedge = 0
!SPH-JAH013108: MUST INCLUDE EDGE FORCE (INITIALLY) FOR V3FITA TO WORK
!ADD A V3FIT RELATED FLAG? ADD fsq criterion first
      delIter = iter2 - iter1

      IF (delIter       .lt. 50 .and.                                          &
          (fsqr + fsqz) .lt. 1.E-6_dp) THEN
         jedge = 1
      END IF

      CALL getfsq (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)

      fsql = fnormL*SUM(gcl*gcl)
      fedge = r1*fnorm*SUM(gcr(ns,:,:,:)**2 + gcz(ns,:,:,:)**2)

!
!     PERFORM PRECONDITIONING AND COMPUTE RESIDUES
!

      IF (ictrl_prec2d .EQ. 1) THEN

         CALL block_precond(gc)

         IF (.not.lfreeb .and. ANY(gcr(ns,:,:,:) .ne. zero)) THEN
            STOP 'gcr(ns) != 0 for fixed boundary in residue'
         END IF
         IF (.not.lfreeb .and. ANY(gcz(ns,:,:,:) .ne. zero)) THEN
            STOP 'gcz(ns) != 0 for fixed boundary in residue'
         END IF
         IF (ANY(gcl(:,1:,0,zsc) .ne. zero)) THEN
            STOP 'gcl(m=0,n>0,sc) != 0 in residue'
         END IF
         IF (lthreed .and. ANY(gcl(:,n0,:,zcs) .ne. zero)) THEN
            STOP 'gcl(n=0,m,cs) != 0 in residue'
         END IF

         fsqr1 = SUM(gcr*gcr)
         fsqz1 = SUM(gcz*gcz)
         fsql1 = SUM(gcl*gcl)

      ELSE
#ifdef _HBANGLE
         CALL scalfor_rho(gcr, gcz)
#else
!        m = 1 constraint scaling
         IF (lthreed) THEN
            CALL scale_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
         END IF
         IF (lasym) THEN
            CALL scale_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
         END IF

         jedge = 0
         CALL scalfor (gcr, arm, brm, ard, brd, crd, jedge)
         jedge = 1
         CALL scalfor (gcz, azm, bzm, azd, bzd, crd, jedge)
#endif

!SPH: add fnorm1 ~ 1/R**2, since preconditioned forces gcr,gcz ~ Rmn or Zmn
         CALL getfsq (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)
#ifdef _HBANGLE
!
!     TO IMPROVE CONVERGENCE, REDUCE FORCES INITIALLY IF THEY ARE TOO LARGE
!
         fac = .5_dp
         IF ((iter2 - iter1) .LT. 25 .AND.                                     &
             (fsqr + fsqz)   .GT. 1.E-2_dp) THEN
            fac = fac / SQRT(1.E2_dp*(fsqr+fsqz))
         END IF
         gcr = fac*gcr
         gcz = fac*gcz
#endif

!SPH: THIS IS NOT INVARIANT UNDER PHIP->A*PHIP, AM->A**2*AM IN PROFIL1D
!     (EXTCUR -> A*EXTCUR for FREE BOUNDARY)

         gcl = faclam*gcl
         fsql1 = hs*SUM(gcl*gcl)
!030514      fsql1 = hs*lamscale**2*SUM(gcl*gcl)

      END IF

      END SUBROUTINE residue

      SUBROUTINE constrain_m1(gcr, gcz)
      USE vmec_main
      USE precon2d, ONLY: ictrl_prec2d
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(ns,0:ntor), INTENT(INOUT) :: gcr, gcz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      REAL(dp), PARAMETER :: FThreshold = 1.E-6_dp
      REAL(dp) :: temp(ns,0:ntor)
!-----------------------------------------------
!
!     COMPUTE INTERNAL gr, gz
!     NOTE: internal gz => 0 for both values of lconm1 (although gz is different)
!     FOR lconm1=T, gcr(internal) = gcr+gcz, gcz(internal) = gcr-gcz->0
!
      IF (lconm1) THEN
         temp = gcr
         gcr = osqrt2*(gcr + gcz)
         gcz = osqrt2*(temp - gcz)
      END IF

!v8.50: ADD iter2<2 so reset=<WOUT_FILE> works
      IF (fsqz  .LT. FThreshold .OR.                                           &
     &    iter2 .LT. 2          .OR.                                           &
     &    ictrl_prec2d .NE. 0) THEN
         gcz = 0
      END IF

      END SUBROUTINE constrain_m1

      SUBROUTINE scale_m1(gcr, gcz)
      USE vmec_main
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), DIMENSION(ns,0:ntor), INTENT(inout) :: gcr, gcz
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: nodd=2
      INTEGER :: n
      REAL(dp) :: fac(ns)
!-----------------------------------------------
      IF (.not.lconm1) RETURN

      fac = (ard(:,nodd) + brd(:,nodd))                                        &
          / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
      DO n = 0, ntor
         gcr(:,n) = fac*gcr(:,n)
      END DO

      fac = (azd(:,nodd) + bzd(:,nodd))                                        &
          / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
      DO n = 0, ntor
         gcz(:,n) = fac*gcz(:,n)
      END DO

      END SUBROUTINE scale_m1

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT


text:
! ------------------------
! Thermodynamic Properties
! ------------------------
Module thermo_props  !
   Implicit None

Contains
        Subroutine mean_molecular_weight(X, M, MM)
                real(8), intent(in) :: X(21), M(21)
                real(8), intent(out) :: MM
                integer :: i

                MM = 0
                do i = 1, size(X)
                MM = MM + X(i) * M(i)
                end do
                MM = MM / 1000 ! Translate to Kilograms
        End Subroutine mean_molecular_weight

   Subroutine zeta(delta, Ar, z)
      real(8), intent(in) :: delta, Ar(3, 3)
      real(8), intent(out) :: z
      z = 1.d0 + delta * Ar(2, 1)
   End Subroutine zeta

   Subroutine isochoric_heat(tau, R, Ao, Ar, cv)
      real(8), intent(in) :: tau, R, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: cv

      cv = - tau ** 2.d0 * (Ao(3, 2) + Ar(3, 2))
      cv = cv * R
   End Subroutine isochoric_heat

   Subroutine isobaric_heat(delta, tau, R, Ao, Ar, cp)
      real(8), intent(in) :: delta, tau, R, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: cp
      real(8) :: up, down

      up = (1.d0 + delta * Ar(2, 1) - delta * tau * Ar(3, 3)) ** 2.d0
      down = 1.d0 + 2.d0 * delta * Ar(2, 1) + delta ** 2.d0 * Ar(3, 1)

      cp = - tau ** 2.d0 * (Ao(3, 2) + Ar(3, 2)) + up / down
      cp = cp * R
   End Subroutine isobaric_heat

   Subroutine sound_speed(delta, tau, R, T, M, Ao, Ar, w)
      real(8), intent(in) :: delta, tau, R, T, M, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: w
      real(8) :: up, down

      up = (1.d0 + delta * Ar(2, 1) - delta * tau * Ar(3, 3)) ** 2
      down = (tau ** 2 * (Ao(3, 2) + Ar(3, 2)))

      w = 1.d0 + 2.d0 * delta * Ar(2, 1) + delta ** 2 * Ar(3, 1) - up / down
      w = sqrt(w * R * T / M)
   End Subroutine sound_speed

   Subroutine isothermal_thermal_coefficent(delta, tau, rho, Ar, delta_t)
      real(8), intent(in) :: delta, Ar(3, 3), tau, rho
      real(8), intent(out) :: delta_t
      real(8) :: up, down

      up = 1.d0 + delta * Ar(2, 1) - delta * tau * Ar(3, 3)
      down = 1.d0 + 2.d0 * delta * Ar(2, 1) + delta ** 2 * Ar(3, 1)

      delta_t = (1.d0 - up / down) / rho
   End Subroutine isothermal_thermal_coefficent

   Subroutine dp_dt(rho, delta, tau, R, Ar, dpdt)
      real(8), intent(in) :: delta, R, rho, tau
      real(8), intent(out) :: dpdt
      real(8), dimension(3, 3), intent(in) :: Ar

      dpdt = rho * R * (1.d0 + delta * Ar(2, 1) - delta * tau * Ar(3, 3))

   End Subroutine dp_dt

   Subroutine dp_drho(T, delta, R, Ar, dpdrho)
      real(8), intent(in) :: delta, R, T
      real(8), intent(out) :: dpdrho
      real(8), dimension(3, 3), intent(in) :: Ar

      dpdrho = R * T * (1.d0 + 2.d0 * delta * Ar(2, 1) + (delta ** 2) * Ar(3, 1))

   End Subroutine dp_drho

   Subroutine dp_dv(rho, delta, T, R, Ar, dpdv) ! !Recordar que acá hablamos de volumen molar!!
      real(8), intent(in) :: delta, R, rho, T
      real(8), intent(out) :: dpdv
      real(8), dimension(3, 3), intent(in) :: Ar

      dpdv = - (rho ** 2) * R * T * (1.d0 + 2.d0 * delta * Ar(2, 1) + (delta ** 2) * Ar(3, 1))

   End Subroutine dp_dv

   Subroutine pressure(delta, rho, R, T, Ar, p)
      real(8), intent(in) :: delta, rho, R, T, Ar(3, 3)
      real(8), intent(out) :: p
      real(8) :: z

      call zeta(delta, Ar, z)
      p = z * (rho * 1000.d0) * R * T
   End Subroutine pressure

   Subroutine entropy(tau, R, Ao, Ar, s)
      real(8), intent(in) :: tau, R, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: s

      s = (tau * (Ao(2, 2) + Ar(2, 2)) - Ao(1, 1) - Ar(1, 1)) * R
   End Subroutine entropy

   Subroutine internal_energy(tau, R, T, Ao, Ar, u)
      real(8), intent(in) :: tau, R, T, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: u

      u = tau * (Ao(2, 2) + Ar(2, 2)) * R * T
   End Subroutine internal_energy

   Subroutine enthalpy(delta, tau, R, T, Ao, Ar, h)
      real(8), intent(in) :: delta, tau, R, T, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: h

      h = ((1.d0 + tau * (Ao(2, 2) + Ar(2, 2)) + delta * Ar(2, 1))) * R * T
   End Subroutine enthalpy

   Subroutine gibbs_free_energy(delta, R, T, Ao, Ar, g)
      real(8), intent(in) :: delta, R, T, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: g

      g = (1.d0 + Ao(1, 1) + Ar(1, 1) + delta * Ar(2, 1)) * R * T
   End Subroutine gibbs_free_energy

   Subroutine joule_thomson_coeff(delta, tau, rho, R, Ao, Ar, JT)
      real(8), intent(in) :: delta, tau, rho, R, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: JT
      real(8) :: up, down_1, down_2

      up = - (delta * Ar(2, 1) + delta ** 2 * Ar(3, 1) + delta * tau * Ar(3, 3))
      down_1 = (1.d0 + delta * Ar(2, 1) - delta * tau * Ar (3, 3)) ** 2
      down_2 = - tau ** 2 * (Ao(3, 2) + Ar(3, 2)) * (1.d0 + 2.d0 * delta * Ar(2, 1) + delta ** 2 * Ar(3, 1))

      JT = up / (down_1 + down_2) / (R * rho * 1000)

   End Subroutine joule_thomson_coeff

   Subroutine isentropic_exponent(delta, tau, Ao, Ar, k)
      real(8), intent(in) :: delta, tau, Ao(3, 3), Ar(3, 3)
      real(8), intent(out) :: k
      real(8) :: up1, down1, up2, down2

      up1 = 1.d0 + 2 * delta * Ar(2, 1) + delta ** 2 * Ar(3, 1)
      down1 = 1.d0 + delta * Ar(2, 1)
      up2 = (1.d0 + delta * Ar(2, 1) - delta * tau * Ar(2, 3)) ** 2
      down2 = tau ** 2 * (Ao(3, 2) + Ar(2, 2)) * (1.d0 + 2.d0 * delta * Ar(2, 1) + delta ** 2 * Ar(3, 2))

      k = (up1 / down1) * (1.d0 - up2 / down2)
   End Subroutine isentropic_exponent

   Subroutine second_thermal_virial_coeff(rho_r, Ar, B)
      real(8), intent(in) :: rho_r, Ar(3, 3)
      real(8), intent(out) :: B
      real(8) :: delta

      delta = 1d-15
      B = Ar(2, 1) / rho_r
   End Subroutine second_thermal_virial_coeff

   Subroutine third_thermal_virial_coeff(rho_r, Ar, C)
      real(8), intent(in) :: rho_r, Ar(3, 3)
      real(8), intent(out) :: C
      real(8) :: delta

      delta = 1d-15
      C = Ar(3, 1) / (rho_r ** 2)
   End Subroutine third_thermal_virial_coeff

   Subroutine helmholtz_per_mol(x, delta, tau, rho_r, T_r, ar, ar_x, ar_dx, &
                                v_r_x, T_r_x, dar_dn, dadr_dn)
           real(8), intent(in) :: x(21), delta, tau, rho_r, T_r, &
                   ar(3, 3), ar_x(21), ar_dx(21), v_r_x(21), T_r_x(21)
           real(8), dimension(21), intent(out) :: dar_dn, dadr_dn
           real(8), dimension(21) :: drhor_dn, dtr_dn
           integer :: i

           dar_dn = 0.d0
           dadr_dn = 0.d0

           drhor_dn = - rho_r ** 2 * v_r_x
           dtr_dn = v_r_x

           do i = 1, size(dtr_dn)
                drhor_dn = drhor_dn - (- x(i) * (- rho_r ** 2 * v_r_x(i)))
                dtr_dn = dtr_dn - x(i) * t_r_x(i)
           end do

           dar_dn = delta * ar(2, 1) * (1.d0 - 1.d0 / rho_r * dtr_dn)
           dar_dn = dar_dn + tau * ar(2, 2) / T_r * dtr_dn + ar_x

           dadr_dn = delta * ar(3, 1) * (1.d0 - 1.d0 / rho_r * dtr_dn)
           dadr_dn = dadr_dn + tau * ar(3, 3) / T_r * dtr_dn + ar_dx

           do i = 1, size(dtr_dn)
                   dar_dn = dar_dn - X(i) * ar_x(i)
                   dadr_dn = dadr_dn - X(i) * ar_dx(i)
           end do
   End Subroutine helmholtz_per_mol

End Module thermo_props

max_issues_repo_name: fedebenelli/PyForFluids


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine exchange_4(g,h,ibeg,ifin1,jbeg,jfin1)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   compute the right hand side based on exact solution
c---------------------------------------------------------------------

      implicit none

      include 'mpinpb.h'
      include 'applu.incl'

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      double precision  g(0:isiz2+1,0:isiz3+1),
     >        h(0:isiz2+1,0:isiz3+1)
      integer ibeg, ifin1
      integer jbeg, jfin1

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j
      integer ny2
      double precision  dum(1024)

      integer msgid1, msgid3
      integer STATUS(MPI_STATUS_SIZE)
      integer IERROR

      ny2 = ny + 2

c---------------------------------------------------------------------
c   communicate in the east and west directions
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   receive from east
c---------------------------------------------------------------------
      if (jfin1.eq.ny) then
        call MPI_IRECV( dum,
     >                  2*nx,
     >                  dp_type,
     >                  east,
     >                  from_e,
     >                  MPI_COMM_WORLD,
     >                  msgid3,
     >                  IERROR )

        call MPI_WAIT( msgid3, STATUS, IERROR )

        do i = 1,nx
          g(i,ny+1) = dum(i)
          h(i,ny+1) = dum(i+nx)
        end do

      end if

c---------------------------------------------------------------------
c   send west
c---------------------------------------------------------------------
      if (jbeg.eq.1) then
        do i = 1,nx
          dum(i) = g(i,1)
          dum(i+nx) = h(i,1)
        end do

        call MPI_SEND( dum,
     >                 2*nx,
     >                 dp_type,
     >                 west,
     >                 from_e,
     >                 MPI_COMM_WORLD,
     >                 IERROR )

      end if

c---------------------------------------------------------------------
c   communicate in the south and north directions
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   receive from south
c---------------------------------------------------------------------
      if (ifin1.eq.nx) then
        call MPI_IRECV( dum,
     >                  2*ny2,
     >                  dp_type,
     >                  south,
     >                  from_s,
     >                  MPI_COMM_WORLD,
     >                  msgid1,
     >                  IERROR )

        call MPI_WAIT( msgid1, STATUS, IERROR )

        do j = 0,ny+1
          g(nx+1,j) = dum(j+1)
          h(nx+1,j) = dum(j+ny2+1)
        end do

      end if

c---------------------------------------------------------------------
c   send north
c---------------------------------------------------------------------
      if (ibeg.eq.1) then
        do j = 0,ny+1
          dum(j+1) = g(1,j)
          dum(j+ny2+1) = h(1,j)
        end do

        call MPI_SEND( dum,
     >                 2*ny2,
     >                 dp_type,
     >                 north,
     >                 from_s,
     >                 MPI_COMM_WORLD,
     >                 IERROR )

      end if

      return
      end

max_stars_repo_name: jonarbo/KUBE

max_issues_repo_name: jonarbo/KUBE

max_forks_repo_name: jonarbo/KUBE


text:
Also see Museums and Exhibits Museums, Sacramento Galleries and ArtAbout

 Campus Galleries

  The Basement Gallery in the basement of the art building
  Buehler Alumni & Visitors Center several areas
  Carl Gorman Museum a gallery devoted to Native American art
  Craft Center Gallery at the east end of the Craft Center
  Design Gallery and hallway space
  XMemorial Union Art Gallery upstairsX Closed  Now the MU Art Lounge
  Nelson Gallery and back vestibule space
  http://www.lib.ucdavis.edu/ul/about/art/ Shields Library Collection spread throughout the library

 Community Galleries

  The Artery a cooperative gallery
  Boy Scout Cabin open during ArtAbout and by appointment, contact the John Natsoulas Art Gallery to view
  Davis Art Center the Tsao Gallery is at the core
  Gallery 1855 located in the office of the Davis Cemetery, well respected gallery, shows change monthly
  International House the Lounge and the Community Room regularly display works by local artists.
  John Natsoulas Art Gallery a major art gallery, and dont miss the staircase art (really!)
  MichelNovelozo Gallery see the art in the basement too
  Pence Gallery three large gallery spaces
  T.S. Post Studio/Gallery ceramic art by Tom and Sara Post.
  Spirit and Nature Gallery  paintings by Maria Jaoudi
 Establishments with Art
Criteria: public, frequently changing, participate in art events
  Cobalt Salon & Gallery
  CommuniCare Health Center home of the Davis Community Gallery
  Davisville Travel mostly on the west wall
  Design House Furniture Galleries Design House Gallery a large dedicated space, worthwhile
  Nestware
  Primary Concepts a tattoo and piercing establishment with monthly art shows

 Woodland
  wiki:woodland:Centerline Gallery
  wiki:woodland:Studio Artists Gallery  local artists are featured, openings on the first Friday of every month. Classes in various media.

Note: The Davis Art Association printed a very nice Gallery Guide & Public Art Walking Tour brochure in 2004. It lists most of the galleries listed above, along with photographs, addresses, hours, times, etc. Its almost as good as Wiki, plus its portable :). I picked mine up at the Artery. Users/SteveDavison

The Yolo County Visitors Bureau also has the Gallery Guide & Public Art Walking Tour available at their office or they will mail it to you free of charge.