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

module evaluate

!----------------------------------------------------
! This module contains routines for evaluation of the solution.
!
! communication tags in this module are of the form 9xx
!----------------------------------------------------

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

use global
use message_passing
use gridtype_mod
use grid_util
use basis_functions
!----------------------------------------------------

implicit none
private
public evaluate_soln, evaluate_soln_slave, evaluate_soln_local, &
       evaluate_oldsoln_local, copy_old, set_grid_for_old_soln, &
       find_containing_leaf

!----------------------------------------------------
! The following variables are defined:

type(grid_type), pointer, save :: grid_for_old_soln
integer, save :: roundoff_fudge = 1000

!----------------------------------------------------

contains

!          -------------
subroutine evaluate_soln(procs,x,y,u,ux,uy,uxx,uyy,comp,eigen)
!          -------------

!----------------------------------------------------
! This routine evaluates the solution at the points in the arrays (x,y)
! and returns them in the array soln.  The return value is 0.0 for
! points that are outside the domain.  This would only be called by
! the master or a slave requesting an evaluation by another slave.
! If comp and/or eigen is present, it tells which component and/or which
! eigenfunction to return.  The default is 1 for each.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(proc_info), intent(in) :: procs
real(my_real), intent(in) :: x(:),y(:)
real(my_real), optional, intent(out) :: u(:),ux(:),uy(:),uxx(:),uyy(:)
integer, optional, intent(in) :: comp,eigen

!----------------------------------------------------
! Local variables:

integer, allocatable :: send_int(:)
real(my_real), allocatable :: send_real(:)
integer, pointer :: recv_int(:)
real(my_real), pointer :: recv_real(:)
integer :: i, j, ni, nr, proc, first_int, first_real, allocstat

!----------------------------------------------------
! Begin executable code

! Invoke the slaves to evaluate the solution at the points in
! elements they own, and merge the results into a single array

! send the message to evaluate the solution and the points at which to evaluate

call pack_procs_size(this_processors_procs,ni,nr)
allocate(send_int(8+ni),send_real(nr+2*size(x)),stat=allocstat)
if (allocstat /= 0) then
   ierr = ALLOC_FAILED
   call fatal("memory allocation failed in evaluate_soln",procs=procs)
   return
endif
send_int(1) = 3 ! code for "evaluate"
if (present(comp)) then
   send_int(2) = comp
else
   send_int(2) = 1
endif
if (present(eigen)) then
   send_int(3) = eigen
else
   send_int(3) = 1
endif
send_int(4:8) = 0
if (present(u  )) send_int(4) = 1
if (present(ux )) send_int(5) = 1
if (present(uy )) send_int(6) = 1
if (present(uxx)) send_int(7) = 1
if (present(uyy)) send_int(8) = 1
ni = ni+8
first_int = 9
first_real = 1
call pack_procs(this_processors_procs,send_int,first_int,send_real,first_real)
send_real(nr+1:nr+size(x)) = x
send_real(nr+size(x)+1:nr+2*size(x)) = y
nr = nr+2*size(x)
do proc=1,num_proc(procs)
   call phaml_send(procs,proc,send_int,ni,send_real,nr,101)
end do
deallocate(send_real,stat=allocstat)

! receive the solution and flags indicating which points were evaluated
! from each slave

if (present(u  )) u   = 0.0_my_real
if (present(ux )) ux  = 0.0_my_real
if (present(uy )) uy  = 0.0_my_real
if (present(uxx)) uxx = 0.0_my_real
if (present(uyy)) uyy = 0.0_my_real

do i=1,num_proc(procs)
   call phaml_recv(procs,proc,recv_int,ni,recv_real,nr,910)
   j = 0
   if (present(u  )) then
      where (recv_int == 1) u   = recv_real(j+1:j+size(x))
      j = j+size(x)
   endif
   if (present(ux )) then
      where (recv_int == 1) ux  = recv_real(j+1:j+size(x))
      j = j+size(x)
   endif
   if (present(uy )) then
      where (recv_int == 1) uy  = recv_real(j+1:j+size(x))
      j = j+size(x)
   endif
   if (present(uxx)) then
      where (recv_int == 1) uxx = recv_real(j+1:j+size(x))
      j = j+size(x)
   endif
   if (present(uyy)) then
      where (recv_int == 1) uyy = recv_real(j+1:j+size(x))
      j = j+size(x)
   endif
   deallocate(recv_int,recv_real,stat=allocstat)
end do

end subroutine evaluate_soln

!          -------------------
subroutine evaluate_soln_slave(grid,invoker_procs,x,y,comp,eigen, &
                               u,ux,uy,uxx,uyy,which)
!          -------------------

!----------------------------------------------------
! This is the slaves version of evaluate_soln.  If the u's are present, the
! solution is returned in them.  Otherwise, it sends the solutions indicated
! by which to the invoker.  Use of u's or which should be exclusive.
! comp and eigen tell which component and eigenfunction to return.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), intent(in) :: grid
type(proc_info), intent(in) :: invoker_procs
real(my_real), intent(in) :: x(:),y(:)
integer, intent(in) :: comp, eigen
real(my_real), intent(out), optional :: u(:),ux(:),uy(:),uxx(:),uyy(:)
integer, optional, intent(in) :: which(:)

!----------------------------------------------------
! Local variables:

real(my_real) :: loc_u(size(x)), loc_ux(size(x)), loc_uy(size(x)), &
                 loc_uxx(size(x)),loc_uyy(size(x)),xc(3), yc(3)
integer :: have_it(size(x))
integer :: i,j,k,elem,nbasis,astat,isub,nr
real(my_real), allocatable :: basis(:),basisx(:),basisy(:),basisxx(:),basisyy(:)
real(my_real), allocatable :: send_real(:)
integer :: loc_which(5)

!----------------------------------------------------
! Begin executable code

! determine which solutions to evaluate

if (present(which)) then
   loc_which = which
else
   loc_which = 0
   if (present(u  )) loc_which(1) = 1
   if (present(ux )) loc_which(2) = 1
   if (present(uy )) loc_which(3) = 1
   if (present(uxx)) loc_which(4) = 1
   if (present(uyy)) loc_which(5) = 1
endif

have_it = 0

! look for points in my partition

! For each point ...

do i=1,size(x)

! find the element it is in.  have_it==1 if I own the element.
! elem is 0 if I do not own it.

   call find_containing_element(x(i),y(i),have_it(i),elem,grid)

   if (have_it(i)==1) then

! evaluate the basis functions at this point

      nbasis = VERTICES_PER_ELEMENT
      do j=1,EDGES_PER_ELEMENT
         if (grid%edge(grid%element(elem)%edge(j))%degree > 1) then
            nbasis = nbasis + grid%edge(grid%element(elem)%edge(j))%degree - 1
         endif
      end do
      if (grid%element(elem)%degree > 2) then
         nbasis = nbasis + &
               ((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
      endif
      allocate(basis(nbasis),basisx(nbasis),basisy(nbasis),basisxx(nbasis), &
               basisyy(nbasis),stat=astat)
      if (astat /= 0) then
         ierr = ALLOC_FAILED
         call fatal("memory allocation failed in evaluate_soln_slave")
         stop
      endif
      xc = grid%vertex(grid%element(elem)%vertex)%coord%x
      yc = grid%vertex(grid%element(elem)%vertex)%coord%y
      if (loc_which(4)==1 .or. loc_which(5)==1) then
         call p_hier_basis_func(x(i),y(i),xc,yc, &
                                (/grid%edge(grid%element(elem)%edge)%degree, &
                                  grid%element(elem)%degree /), "a", &
                                  basis,basisx,basisy,basisxx,basisyy)
      elseif (loc_which(2)==1 .or. loc_which(3)==1) then
         call p_hier_basis_func(x(i),y(i),xc,yc, &
                                (/grid%edge(grid%element(elem)%edge)%degree, &
                                  grid%element(elem)%degree /), "a", &
                                  basis,basisx,basisy)
      else
         call p_hier_basis_func(x(i),y(i),xc,yc, &
                                (/grid%edge(grid%element(elem)%edge)%degree, &
                                  grid%element(elem)%degree /), "a", &
                                  basis)
      endif

! compute the solution at this point

      loc_u  (i) = 0.0_my_real
      loc_ux (i) = 0.0_my_real
      loc_uy (i) = 0.0_my_real
      loc_uxx(i) = 0.0_my_real
      loc_uyy(i) = 0.0_my_real

      do j=1,VERTICES_PER_ELEMENT
         if (loc_which(1)==1) loc_u  (i) = loc_u  (i) + &
              basis  (j)*grid%vertex_solution(grid%element(elem)%vertex(j), &
              comp,eigen)
         if (loc_which(2)==1) loc_ux (i) = loc_ux (i) + &
              basisx (j)*grid%vertex_solution(grid%element(elem)%vertex(j), &
              comp,eigen)
         if (loc_which(3)==1) loc_uy (i) = loc_uy (i) + &
              basisy (j)*grid%vertex_solution(grid%element(elem)%vertex(j), &
              comp,eigen)
         if (loc_which(4)==1) loc_uxx(i) = loc_uxx(i) + &
              basisxx(j)*grid%vertex_solution(grid%element(elem)%vertex(j), &
              comp,eigen)
         if (loc_which(5)==1) loc_uyy(i) = loc_uyy(i) + &
              basisyy(j)*grid%vertex_solution(grid%element(elem)%vertex(j), &
              comp,eigen)
      end do
      isub = VERTICES_PER_ELEMENT
      do j=1,EDGES_PER_ELEMENT
         if (grid%edge(grid%element(elem)%edge(j))%degree < 2) cycle
         do k=1,grid%edge(grid%element(elem)%edge(j))%degree-1
            isub = isub + 1
            if (loc_which(1)==1) loc_u  (i) = loc_u  (i) + &
               basis  (isub)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
               comp,eigen)
            if (loc_which(2)==1) loc_ux (i) = loc_ux (i) + &
               basisx (isub)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
               comp,eigen)
            if (loc_which(3)==1) loc_uy (i) = loc_uy (i) + &
               basisy (isub)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
               comp,eigen)
            if (loc_which(4)==1) loc_uxx(i) = loc_uxx(i) + &
               basisxx(isub)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
               comp,eigen)
            if (loc_which(5)==1) loc_uyy(i) = loc_uyy(i) + &
               basisyy(isub)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
               comp,eigen)
         end do
      end do
      if (grid%element(elem)%degree > 2) then
         do k=1,((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
            isub = isub + 1
            if (loc_which(1)==1) loc_u  (i) = loc_u  (i) + &
               basis  (isub)*grid%element(elem)%solution(k,comp,eigen)
            if (loc_which(2)==1) loc_ux (i) = loc_ux (i) + &
               basisx (isub)*grid%element(elem)%solution(k,comp,eigen)
            if (loc_which(3)==1) loc_uy (i) = loc_uy (i) + &
               basisy (isub)*grid%element(elem)%solution(k,comp,eigen)
            if (loc_which(4)==1) loc_uxx(i) = loc_uxx(i) + &
               basisxx(isub)*grid%element(elem)%solution(k,comp,eigen)
            if (loc_which(5)==1) loc_uyy(i) = loc_uyy(i) + &
               basisyy(isub)*grid%element(elem)%solution(k,comp,eigen)
         end do
      endif

   else

      loc_u  (i) = 0.0_my_real
      loc_ux (i) = 0.0_my_real
      loc_uy (i) = 0.0_my_real
      loc_uxx(i) = 0.0_my_real
      loc_uyy(i) = 0.0_my_real

   endif
   deallocate(basis,basisx,basisy,basisxx,basisyy,stat=astat)
end do

! Copy the solution to the u's or send the solution to the invoker.
! Note that my_proc(invoker_procs) is not me, it is the invoker.

if (present(which)) then
   nr = size(x)*count(loc_which==1)
   allocate(send_real(nr),stat=astat)
   if (astat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in evaluate_soln_slave")
      stop
   endif
   j = 0
   if (loc_which(1) == 1) then
      send_real(j+1:j+size(x)) = loc_u
      j = j + size(x)
   endif
   if (loc_which(2) == 1) then
      send_real(j+1:j+size(x)) = loc_ux
      j = j + size(x)
   endif
   if (loc_which(3) == 1) then
      send_real(j+1:j+size(x)) = loc_uy
      j = j + size(x)
   endif
   if (loc_which(4) == 1) then
      send_real(j+1:j+size(x)) = loc_uxx
      j = j + size(x)
   endif
   if (loc_which(5) == 1) then
      send_real(j+1:j+size(x)) = loc_uyy
      j = j + size(x)
   endif
   call phaml_send(invoker_procs,my_proc(invoker_procs),have_it,size(x), &
                   send_real,nr,910)
   deallocate(send_real)
else
   if (present(u  )) u   = loc_u
   if (present(ux )) ux  = loc_ux
   if (present(uy )) uy  = loc_uy
   if (present(uxx)) uxx = loc_uxx
   if (present(uyy)) uyy = loc_uyy
endif

end subroutine evaluate_soln_slave

!          -----------------------
subroutine find_containing_element(x,y,have_it,elem,grid)
!          -----------------------

!----------------------------------------------------
! This routine looks for an element containing the point (x,y).  If it
! determines that (x,y) is not in an element owned by this processor, it
! quits looking and returns (have_it=0,elem=0).  Otherwise it returns
! have_it=1 and the element index in elem.  If the point falls on an
! element boundary, it is indeterminant as to which of the containing
! elements it returns.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
integer, intent(out) :: have_it,elem
type(grid_type), intent(in) :: grid

!----------------------------------------------------
! Local variables:

integer :: e,root,i,neigh(EDGES_PER_ELEMENT)
real(my_real) :: bc(VERTICES_PER_ELEMENT)
!----------------------------------------------------
! Begin executable code

! find an element in the initial grid that contains (x,y)

! first look for it by moving from a triangle to a neighbor in the direction
! of the point, which is characterized by a negative barycentric coordinate

root = -10
e = grid%head_level_elem(1)
do
! if all barycentric coordinates are positive, it's in there
   bc = barycentric(x,y,e,grid,no_det=.true.)
   if (all(bc >= -roundoff_fudge*epsilon(0.0_my_real))) then
      root = e
      exit
   endif
   neigh = grid%initial_neighbor(:,e)
   do i=1,VERTICES_PER_ELEMENT
      if (bc(i) < 0) then
         e = neigh(i)
         if (e /= BOUNDARY) exit
      endif
   end do
   if (e == BOUNDARY) exit
end do

! if that failed, go through all the initial elements

if (root == -10) then
   e = grid%head_level_elem(1)
   do while (e /= END_OF_LIST)
      if (all(barycentric(x,y,e,grid,no_det=.true.) >= &
              -roundoff_fudge*epsilon(0.0_my_real))) then
         root = e
         exit
      endif
      e = grid%element(e)%next
   end do
endif

! go down the refinement tree to find the leaf element that contains (x,y)

if (root /= -10) then
   call find_containing_leaf(x,y,root,have_it,elem,grid)
else
   have_it = 0
   elem = 0
endif

end subroutine find_containing_element

!                    --------------------
recursive subroutine find_containing_leaf(x,y,root,have_it,elem,grid)
!                    --------------------

!----------------------------------------------------
! This routine recursively goes down the refinement tree, starting at root,
! to find a leaf element containing (x,y) as long as (x,y) may be in an
! element I own.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
integer, intent(in) :: root
integer, intent(out) :: have_it,elem
type(grid_type), intent(in) :: grid

!----------------------------------------------------
! Local variables:

integer :: children(MAX_CHILD), i, allc(MAX_CHILD)
!----------------------------------------------------
! Begin executable code

! if root is a leaf, we've found it; but verify ownership

allc = ALL_CHILDREN
children = get_child_lid(grid%element(root)%gid,allc,grid%elem_hash)
if (children(1) == NO_CHILD) then
   if (grid%element(root)%iown) then
      have_it = 1
      elem = root
   else
      have_it = 0
      elem = 0
   endif
   return
endif

! otherwise, look at the barycentric coordinates of (x,y) in each child,
! until one is found where they are all positive

have_it = 0
elem = 0
do i=1,MAX_CHILD
   if (all(barycentric(x,y,children(i),grid,no_det=.true.) >= &
           -roundoff_fudge*epsilon(0.0_my_real))) then
      call find_containing_leaf(x,y,children(i),have_it,elem,grid)
      exit
   endif
end do

end subroutine find_containing_leaf

!        -----------
function barycentric(x,y,elem,grid,no_det) result(zeta)
!        -----------

!----------------------------------------------------
! This function returns the barycentric coordinates of (x,y) in element elem.
! If no_det is present and .true., the coordinates are not divided by the
! determinant (area of the triangle).  This is useful if you are only looking
! at the signs of the coordinates, where it can be wrong if the point is
! very close to an edge and the area is very small so roundoff can be very bad.
! RESTRICTION triangles
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
integer, intent(in) :: elem
type(grid_type), intent(in) :: grid
real(my_real), dimension(VERTICES_PER_ELEMENT) :: zeta
logical, optional, intent(in) :: no_det

!----------------------------------------------------
! Local variables:

real(quad_real) :: x1,x2,x3,y1,y2,y3,det
real(quad_real) :: xy1,xy2,xy3,yx1,yx2,yx3,x1y2,x2y3,x3y1,x1y3,x2y1,x3y2
real(quad_real) :: sump,summ
logical :: local_no_det

!----------------------------------------------------
! Begin executable code

! check for request for no determinant

if (present(no_det)) then
   local_no_det = no_det
else
   local_no_det = .false.
endif

! local variables for the vertices, to make the code easier to read

x1 = real(grid%vertex(grid%element(elem)%vertex(1))%coord%x,quad_real)
x2 = real(grid%vertex(grid%element(elem)%vertex(2))%coord%x,quad_real)
x3 = real(grid%vertex(grid%element(elem)%vertex(3))%coord%x,quad_real)
y1 = real(grid%vertex(grid%element(elem)%vertex(1))%coord%y,quad_real)
y2 = real(grid%vertex(grid%element(elem)%vertex(2))%coord%y,quad_real)
y3 = real(grid%vertex(grid%element(elem)%vertex(3))%coord%y,quad_real)

! compute the barycentric coordinates of the point

! reduce roundoff by summing all the positive parts and negative parts
! separately, and then adding the two partial sums

! products needed for the sums

xy1 = real(x,quad_real)*y1
xy2 = real(x,quad_real)*y2
xy3 = real(x,quad_real)*y3
yx1 = real(y,quad_real)*x1
yx2 = real(y,quad_real)*x2
yx3 = real(y,quad_real)*x3
x1y2 = x1*y2
x2y3 = x2*y3
x3y1 = x3*y1
x1y3 = x1*y3
x2y1 = x2*y1
x3y2 = x3*y2

! compute the determinant

! det = x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)

sump=0.0_quad_real; summ=0.0_quad_real
if (x1y2 > 0.0_quad_real) then; sump=sump+x1y2; else; summ=summ+x1y2; endif
if (x1y3 < 0.0_quad_real) then; sump=sump-x1y3; else; summ=summ-x1y3; endif
if (x2y3 > 0.0_quad_real) then; sump=sump+x2y3; else; summ=summ+x2y3; endif
if (x2y1 < 0.0_quad_real) then; sump=sump-x2y1; else; summ=summ-x2y1; endif
if (x3y1 > 0.0_quad_real) then; sump=sump+x3y1; else; summ=summ+x3y1; endif
if (x3y2 < 0.0_quad_real) then; sump=sump-x3y2; else; summ=summ-x3y2; endif
det = sump + summ

! if the request is to not divide by the determinant, only use its sign

if (local_no_det) then
   det = sign(1.0_quad_real,det)
endif

! compute the coordinates

! zeta(1) = (x*(y2-y3) + y*(x3-x2) + (x2*y3-x3*y2))/det

sump=0.0_quad_real; summ=0.0_quad_real
if (xy2 > 0.0_quad_real) then; sump=sump+xy2; else; summ=summ+xy2; endif
if (xy3 < 0.0_quad_real) then; sump=sump-xy3; else; summ=summ-xy3; endif
if (yx3 > 0.0_quad_real) then; sump=sump+yx3; else; summ=summ+yx3; endif
if (yx2 < 0.0_quad_real) then; sump=sump-yx2; else; summ=summ-yx2; endif
if (x2y3 > 0.0_quad_real) then; sump=sump+x2y3; else; summ=summ+x2y3; endif
if (x3y2 < 0.0_quad_real) then; sump=sump-x3y2; else; summ=summ-x3y2; endif
zeta(1) = sump/det + summ/det

! zeta(2) = (x*(y3-y1) + y*(x1-x3) + (x3*y1-x1*y3))/det

sump=0.0_quad_real; summ=0.0_quad_real
if (xy3 > 0.0_quad_real) then; sump=sump+xy3; else; summ=summ+xy3; endif
if (xy1 < 0.0_quad_real) then; sump=sump-xy1; else; summ=summ-xy1; endif
if (yx1 > 0.0_quad_real) then; sump=sump+yx1; else; summ=summ+yx1; endif
if (yx3 < 0.0_quad_real) then; sump=sump-yx3; else; summ=summ-yx3; endif
if (x3y1 > 0.0_quad_real) then; sump=sump+x3y1; else; summ=summ+x3y1; endif
if (x1y3 < 0.0_quad_real) then; sump=sump-x1y3; else; summ=summ-x1y3; endif
zeta(2) = sump/det + summ/det

! zeta(3) = (x*(y1-y2) + y*(x2-x1) + (x1*y2-x2*y1))/det

sump=0.0_quad_real; summ=0.0_quad_real
if (xy1 > 0.0_quad_real) then; sump=sump+xy1; else; summ=summ+xy1; endif
if (xy2 < 0.0_quad_real) then; sump=sump-xy2; else; summ=summ-xy2; endif
if (yx2 > 0.0_quad_real) then; sump=sump+yx2; else; summ=summ+yx2; endif
if (yx1 < 0.0_quad_real) then; sump=sump-yx1; else; summ=summ-yx1; endif
if (x1y2 > 0.0_quad_real) then; sump=sump+x1y2; else; summ=summ+x1y2; endif
if (x2y1 < 0.0_quad_real) then; sump=sump-x2y1; else; summ=summ-x2y1; endif
zeta(3) = sump/det + summ/det

if (VERTICES_PER_ELEMENT /= 3) then
   zeta(4:VERTICES_PER_ELEMENT) = 0.0_my_real
   call warning("function barycentric needs to be changed for nontriangles")
endif

end function barycentric

!          -------------------
subroutine evaluate_soln_local(grid,x,y,elem,comp,eigen,u,ux,uy,uxx,uyy,uxy)
!          -------------------

!----------------------------------------------------
! This routine evaluates the solution and/or derivatives of the solution
! on this processor.  The points (x,y) must all be in element elem.
! comp and eigen tell which components and eigenfunctions to return.
! The three dimensions of u correspond to comp, eigen and point.
! If ux or uy is present, then both must be present.  If uxx or uyy is
! present, then all 5 must be present.  This is only to reduce
! the number of forms of the call to basis_function, so this restriction can
! easily be removed.  Likewise, if uxy is present then u, ux and uy must be
! present and uxx and uyy either both present or not present.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), intent(in) :: grid
real(my_real), intent(in) :: x(:),y(:)
integer, intent(in) :: elem
integer, intent(in) :: comp(:),eigen(:)
real(my_real), optional, intent(out) :: u(:,:,:),ux(:,:,:),uy(:,:,:), &
                                        uxx(:,:,:),uyy(:,:,:),uxy(:,:,:)
! dimensions are (comp,eigen,point)

!----------------------------------------------------
! Local variables:

real(my_real) :: xc(3), yc(3)
real(my_real), allocatable :: basis(:,:),basisx(:,:),basisy(:,:),basisxx(:,:), &
                              basisyy(:,:),basisxy(:,:),solnvect(:,:)
integer :: i,j,k,l,p,nbasis,astat,isub
logical :: useblas

!----------------------------------------------------
! Begin executable code

! evaluate the basis functions at these points

nbasis = VERTICES_PER_ELEMENT
do j=1,EDGES_PER_ELEMENT
   if (grid%edge(grid%element(elem)%edge(j))%degree > 1) then
      nbasis = nbasis + grid%edge(grid%element(elem)%edge(j))%degree - 1
   endif
end do
if (grid%element(elem)%degree > 2) then
   nbasis = nbasis + &
         ((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
endif
allocate(basis(nbasis,size(x)),basisx(nbasis,size(x)),basisy(nbasis,size(x)), &
         basisxx(nbasis,size(x)),basisyy(nbasis,size(x)), &
         basisxy(nbasis,size(x)),stat=astat)
if (astat /= 0) then
   ierr = ALLOC_FAILED
   call fatal("memory allocation failed in evaluate_soln_local")
   stop
endif
if (present(ux) .neqv. present(uy)) then
   ierr = PHAML_INTERNAL_ERROR
   call fatal("need both or neither of ux and uy in evaluate_soln_local")
   stop
endif
if (present(uxx) .or. present(uyy)) then
   if (.not. present(u) .or. .not. present(ux) .or. .not. present(uy) .or. &
       .not. present(uxx) .or. .not. present(uyy)) then
      ierr = PHAML_INTERNAL_ERROR
      call fatal("if uxx or uyy is present in evaluate_soln_local then all args except uxy must be present")
      stop
   endif
endif
if (present(uxy)) then
   if (.not. present(u) .or. .not. present(ux) .or. .not. present(uy) .or. &
       (present(uxx) .neqv. present(uyy))) then
      ierr = PHAML_INTERNAL_ERROR
      call fatal("evaluate_soln_local: uxy present requires u, ux and uy present and uxx and uyy both or neither present")
      stop
   endif
endif
xc = grid%vertex(grid%element(elem)%vertex)%coord%x
yc = grid%vertex(grid%element(elem)%vertex)%coord%y
if (present(u) .and. .not. present(ux)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basis)
elseif (.not. present(u) .and. present(ux)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basisx=basisx,basisy=basisy)
elseif (present(uxy) .and. present(uxx)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy,basisxx,basisyy,basisxy)
elseif (present(uxy) .and. .not. present(uxx)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy,basisxy=basisxy)
elseif (present(uxx)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy,basisxx,basisyy)
else
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                           grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy)
endif

! compute the solution at this point

if (present(u)) u = 0.0_my_real
if (present(ux)) ux = 0.0_my_real
if (present(uy)) uy = 0.0_my_real
if (present(uxx)) uxx = 0.0_my_real
if (present(uyy)) uyy = 0.0_my_real
if (present(uxy)) uxy = 0.0_my_real

! special BLAS code for size(comp)==size(eigen)==1.  Copy solution values to
! an array and use GEMM.  TEMP Probably can do something similar with /=1.
! Also need the first 2 dimensions of the u's to be 1 since we're
! (illegally!) passing a rank 3 array to a rank 2 dummy argument.

useblas = size(comp)==1 .and. size(eigen)==1
if (present(u)) useblas = useblas .and. size(u,dim=1) == 1 .and. &
                          size(u,dim=2) == 1
if (present(ux)) useblas = useblas .and. size(ux,dim=1) == 1 .and. &
                          size(ux,dim=2) == 1
if (present(uy)) useblas = useblas .and. size(uy,dim=1) == 1 .and. &
                          size(uy,dim=2) == 1
if (present(uxx)) useblas = useblas .and. size(uxx,dim=1) == 1 .and. &
                          size(uxx,dim=2) == 1
if (present(uyy)) useblas = useblas .and. size(uyy,dim=1) == 1 .and. &
                          size(uyy,dim=2) == 1
if (present(uxy)) useblas = useblas .and. size(uxy,dim=1) == 1 .and. &
                          size(uxy,dim=2) == 1

if (useblas) then
   allocate(solnvect(1,nbasis),stat=astat)
   if (astat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in evaluate_soln_local")
      stop
   endif

   do j=1,VERTICES_PER_ELEMENT
      solnvect(1,j) = grid%vertex_solution(grid%element(elem)%vertex(j), &
                                           comp(1),eigen(1))
   end do
   isub = VERTICES_PER_ELEMENT
   do j=1,EDGES_PER_ELEMENT
      if (grid%edge(grid%element(elem)%edge(j))%degree < 2) cycle
      do k=1,grid%edge(grid%element(elem)%edge(j))%degree-1
         isub = isub + 1
         solnvect(1,isub) = grid%edge(grid%element(elem)%edge(j))%solution(k, &
                            comp(1),eigen(1))
      end do
   end do
   if (grid%element(elem)%degree > 2) then
      do k=1,((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
         isub = isub + 1
         solnvect(1,isub) = grid%element(elem)%solution(k,comp(1),eigen(1))
      end do
   endif

   if (my_real == kind(1.0)) then
      if (present(u)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                 basis,nbasis,1.0,u,1)
      if (present(ux)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                  basisx,nbasis,1.0,ux,1)
      if (present(uy)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                  basisy,nbasis,1.0,uy,1)
      if (present(uxx)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                  basisxx,nbasis,1.0,uxx,1)
      if (present(uyy)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                  basisyy,nbasis,1.0,uyy,1)
      if (present(uxy)) call sgemm("N","N",1,size(x),nbasis,1.0,solnvect,1, &
                                  basisxy,nbasis,1.0,uxy,1)
   elseif (my_real == kind(1.0d0)) then
      if (present(u)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                 basis,nbasis,1.0d0,u,1)
      if (present(ux)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                  basisx,nbasis,1.0d0,ux,1)
      if (present(uy)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                  basisy,nbasis,1.0d0,uy,1)
      if (present(uxx)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                  basisxx,nbasis,1.0d0,uxx,1)
      if (present(uyy)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                  basisyy,nbasis,1.0d0,uyy,1)
      if (present(uxy)) call dgemm("N","N",1,size(x),nbasis,1.0d0,solnvect,1, &
                                  basisxy,nbasis,1.0d0,uxy,1)
   else
      ierr = PHAML_INTERNAL_ERROR
      call fatal("my_real is neither single nor double precision. Can't call GEMM")
      stop
   endif

   deallocate(solnvect,stat=astat)

else ! general case with size(comp)/=1 .or. size(eigen)/=1

   do j=1,VERTICES_PER_ELEMENT
      do p=1,size(x)
       do i=1,size(comp)
        do l=1,size(eigen)
         if (present(u)) u(i,l,p) = u(i,l,p) + &
            basis(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                       comp(i),eigen(l))
         if (present(ux)) ux(i,l,p) = ux(i,l,p) + &
            basisx(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                        comp(i),eigen(l))
         if (present(uy)) uy(i,l,p) = uy(i,l,p) + &
            basisy(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                        comp(i),eigen(l))
         if (present(uxx)) uxx(i,l,p) = uxx(i,l,p) + &
            basisxx(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                        comp(i),eigen(l))
         if (present(uyy)) uyy(i,l,p) = uyy(i,l,p) + &
            basisyy(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                        comp(i),eigen(l))
         if (present(uxy)) uxy(i,l,p) = uxy(i,l,p) + &
            basisxy(j,p)*grid%vertex_solution(grid%element(elem)%vertex(j), &
                        comp(i),eigen(l))
        end do
       end do
      end do
   end do
   isub = VERTICES_PER_ELEMENT
   do j=1,EDGES_PER_ELEMENT
      if (grid%edge(grid%element(elem)%edge(j))%degree < 2) cycle
      do k=1,grid%edge(grid%element(elem)%edge(j))%degree-1
         isub = isub + 1
         do p=1,size(x)
          do i=1,size(comp)
           do l=1,size(eigen)
            if (present(u)) u(i,l,p) = u(i,l,p) + &
               basis(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k, &
                             comp(i),eigen(l))
            if (present(ux)) ux(i,l,p) = ux(i,l,p) + &
               basisx(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k,&
                              comp(i),eigen(l))
            if (present(uy)) uy(i,l,p) = uy(i,l,p) + &
               basisy(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k,&
                              comp(i),eigen(l))
            if (present(uxx)) uxx(i,l,p) = uxx(i,l,p) + &
              basisxx(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k,&
                              comp(i),eigen(l))
            if (present(uyy)) uyy(i,l,p) = uyy(i,l,p) + &
              basisyy(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k,&
                              comp(i),eigen(l))
            if (present(uxy)) uxy(i,l,p) = uxy(i,l,p) + &
              basisxy(isub,p)*grid%edge(grid%element(elem)%edge(j))%solution(k,&
                              comp(i),eigen(l))
           end do
          end do
         end do
      end do
   end do
   if (grid%element(elem)%degree > 2) then
      do k=1,((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
         isub = isub + 1
         do p=1,size(x)
          do i=1,size(comp)
           do l=1,size(eigen)
            if (present(u)) u(i,l,p) = u(i,l,p) + &
               basis(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
            if (present(ux)) ux(i,l,p) = ux(i,l,p) + &
               basisx(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
            if (present(uy)) uy(i,l,p) = uy(i,l,p) + &
               basisy(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
            if (present(uxx)) uxx(i,l,p) = uxx(i,l,p) + &
               basisxx(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
            if (present(uyy)) uyy(i,l,p) = uyy(i,l,p) + &
               basisyy(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
            if (present(uxy)) uxy(i,l,p) = uxy(i,l,p) + &
               basisxy(isub,p)*grid%element(elem)%solution(k,comp(i),eigen(l))
           end do
          end do
         end do
      end do
   endif

endif ! size(comp)==size(eigen)==1

deallocate(basis,basisx,basisy,basisxx,basisyy,basisxy,stat=astat)

end subroutine evaluate_soln_local

!          ----------------------
subroutine evaluate_oldsoln_local(x,y,u,ux,uy,uxx,uyy,comp,eigen)
!          ----------------------

!----------------------------------------------------
! This routine evaluates the old solution and/or derivatives of the old solution
! on this processor.
! If present, comp and/or eigen tell which component and eigenfunction to return
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
real(my_real), optional, intent(out) :: u,ux,uy,uxx,uyy
integer, optional, intent(in) :: comp, eigen

!----------------------------------------------------
! Local variables:

integer :: j,k,nbasis,astat,isub,elem,loc_comp,loc_eigen
real(my_real), allocatable :: basis(:),basisx(:),basisy(:),basisxx(:),basisyy(:)
real(my_real) :: xc(3), yc(3), solnval
type(grid_type), pointer :: grid

!----------------------------------------------------
! Begin executable code

! make sure the old solution exists

grid => grid_for_old_soln
if (.not. grid%oldsoln_exists) then
   call warning("Old solution has not been saved, and hence cannot be evaluated", &
                "Returning 0.0 as old solution.")
   if (present(u)) u=0.0_my_real
   if (present(ux)) ux=0.0_my_real
   if (present(uy)) uy=0.0_my_real
   if (present(uxx)) uxx=0.0_my_real
   if (present(uyy)) uyy=0.0_my_real
   return
endif

! determine the component and eigenfunction

if (present(comp)) then
   if (comp <= 0) then
      call fatal("phaml_evaluate_oldsoln: comp must be a positive integer.  It is ", &
                  intlist=(/comp/))
      stop
   endif
   if (comp > grid%system_size) then
      call fatal("phaml_evaluate_oldsoln: comp must be no larger than system size.  They are ", &
                 intlist=(/comp,grid%system_size/))
      stop
   endif
   loc_comp = comp
else
   loc_comp = 1
endif
if (present(eigen)) then
   if (eigen <= 0) then
      call fatal("phaml_evaluate_oldsoln: eigen must be a positive integer.  It is ", &
                 intlist=(/eigen/))
      stop
   endif
   if (eigen > max(1,grid%num_eval)) then
      call fatal("phaml_evaluate_oldsoln: eigen must be no larger than the number of eigenvalues computed.  They are ", &
                 intlist=(/eigen,grid%num_eval/))
      stop
   endif
   loc_eigen = eigen
else
   loc_eigen = 1
endif

! find the element containing (x,y)

call find_old_containing_element(x,y,elem,grid)
if (elem == 0) then
   if (present(u)) u=0.0_my_real
   if (present(ux)) ux=0.0_my_real
   if (present(uy)) uy=0.0_my_real
   if (present(uxx)) uxx=0.0_my_real
   if (present(uyy)) uyy=0.0_my_real
   return
endif

! evaluate the basis functions at the point

nbasis = VERTICES_PER_ELEMENT
do j=1,EDGES_PER_ELEMENT
   if (grid%edge(grid%element(elem)%edge(j))%degree > 1) then
      nbasis = nbasis + grid%edge(grid%element(elem)%edge(j))%degree - 1
   endif
end do
if (grid%element(elem)%degree > 2) then
   nbasis = nbasis + &
         ((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
endif
allocate(basis(nbasis),basisx(nbasis),basisy(nbasis),basisxx(nbasis), &
         basisyy(nbasis),stat=astat)
if (astat /= 0) then
   ierr = ALLOC_FAILED
   call fatal("memory allocation failed in evaluate_oldsoln_local")
   stop
endif
xc = grid%vertex(grid%element(elem)%vertex)%coord%x
yc = grid%vertex(grid%element(elem)%vertex)%coord%y
if (present(uxx) .or. present(uyy)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                           grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy,basisxx,basisyy)
elseif (present(ux) .or. present(uy)) then
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                            grid%element(elem)%degree /),"a", &
                          basis,basisx,basisy)
else
   call p_hier_basis_func(x,y,xc,yc, &
                          (/grid%edge(grid%element(elem)%edge)%degree, &
                           grid%element(elem)%degree /),"a", &
                          basis)
endif

! compute the solution at this point

if (present(u)) u = 0.0_my_real
if (present(ux)) ux = 0.0_my_real
if (present(uy)) uy = 0.0_my_real
if (present(uxx)) uxx = 0.0_my_real
if (present(uyy)) uyy = 0.0_my_real

do j=1,VERTICES_PER_ELEMENT
   solnval = grid%vertex_oldsoln(grid%element(elem)%vertex(j),loc_comp,loc_eigen)
   if (present(u)) u = u + solnval*basis(j)
   if (present(ux)) ux = ux + solnval*basisx(j)
   if (present(uy)) uy = uy + solnval*basisy(j)
   if (present(uxx)) uxx = uxx + solnval*basisxx(j)
   if (present(uyy)) uyy = uyy + solnval*basisyy(j)
end do
isub = VERTICES_PER_ELEMENT
do j=1,EDGES_PER_ELEMENT
   if (grid%edge(grid%element(elem)%edge(j))%degree < 2) cycle
   do k=1,grid%edge(grid%element(elem)%edge(j))%degree-1
      isub = isub + 1
      if (.not. associated(grid%edge(grid%element(elem)%edge(j))%oldsoln)) cycle
      if (k > size(grid%edge(grid%element(elem)%edge(j))%oldsoln,dim=1)) cycle
      solnval = grid%edge(grid%element(elem)%edge(j))%oldsoln(k,loc_comp,loc_eigen)
      if (present(u)) u = u + solnval*basis(isub)
      if (present(ux)) ux = ux + solnval*basisx(isub)
      if (present(uy)) uy = uy + solnval*basisy(isub)
      if (present(uxx)) uxx = uxx + solnval*basisxx(isub)
      if (present(uyy)) uyy = uyy + solnval*basisyy(isub)
   end do
end do
if (grid%element(elem)%degree > 2) then
   do k=1,((grid%element(elem)%degree-2)*(grid%element(elem)%degree-1))/2
      isub = isub + 1
      if (.not. associated(grid%element(elem)%oldsoln)) cycle
      solnval = grid%element(elem)%oldsoln(k,loc_comp,loc_eigen)
      if (present(u)) u = u + solnval*basis(isub)
      if (present(ux)) ux = ux + solnval*basisx(isub)
      if (present(uy)) uy = uy + solnval*basisy(isub)
      if (present(uxx)) uxx = uxx + solnval*basisxx(isub)
      if (present(uyy)) uyy = uyy + solnval*basisyy(isub)
   end do
endif

deallocate(basis,basisx,basisy,basisxx,basisyy,stat=astat)

end subroutine evaluate_oldsoln_local

!          ---------------------------
subroutine find_old_containing_element(x,y,elem,grid)
!          ---------------------------

!----------------------------------------------------
! This routine looks for an "effective old leaf" element containing the point
! (x,y).  An element is an effective old leaf if oldsoln is allocated.
! If the point falls on an element boundary, it is indeterminant as to which
! of the containing elements it returns.
! If the point does not fall in any element, 0 is returned.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
integer, intent(out) :: elem
type(grid_type), intent(in) :: grid

!----------------------------------------------------
! Local variables:

integer :: e,root,i,neigh(EDGES_PER_ELEMENT)
real(my_real) :: bc(VERTICES_PER_ELEMENT)
!----------------------------------------------------
! Begin executable code

! find an element in the initial grid that contains (x,y)

! first look for it by moving from a triangle to a neighbor in the direction
! of the point, which is characterized by a negative barycentric coordinate

root = -10
e = grid%head_level_elem(1)
do
! if all barycentric coordinates are positive, it's in there
   bc = barycentric(x,y,e,grid,no_det=.true.)
   if (all(bc >= -roundoff_fudge*epsilon(0.0_my_real))) then
      root = e
      exit
   endif
   neigh = grid%initial_neighbor(:,e)
   do i=1,VERTICES_PER_ELEMENT
      if (bc(i) < 0) then
         e = neigh(i)
         if (e /= BOUNDARY) exit
      endif
   end do
   if (e == BOUNDARY) exit
end do

! if that failed, go through all the initial elements

if (root == -10) then
   e = grid%head_level_elem(1)
   do while (e /= END_OF_LIST)
      if (all(barycentric(x,y,e,grid,no_det=.true.) >= &
              -roundoff_fudge*epsilon(0.0_my_real))) then
         root = e
         exit
      endif
      e = grid%element(e)%next
   end do
endif

! go down the refinement tree to find the effective leaf element that
! contains (x,y)

if (root /= -10) then
   call find_old_containing_leaf(x,y,root,elem,grid)
else
   elem = 0
endif

end subroutine find_old_containing_element

!                    ------------------------
recursive subroutine find_old_containing_leaf(x,y,root,elem,grid)
!                    ------------------------

!----------------------------------------------------
! This routine recursively goes down the refinement tree, starting at root,
! to find a element containing (x,y) that is marked as oldleaf
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

real(my_real), intent(in) :: x,y
integer, intent(in) :: root
integer, intent(out) :: elem
type(grid_type), intent(in) :: grid

!----------------------------------------------------
! Local variables:

integer :: children(MAX_CHILD), i, allc(MAX_CHILD)
!----------------------------------------------------
! Begin executable code

! if root is oldleaf, we've found it

if (grid%element(root)%oldleaf) then
   elem = root
   return
endif

! otherwise, check the children

allc = ALL_CHILDREN
children = get_child_lid(grid%element(root)%gid,allc,grid%elem_hash)
if (children(1) == NO_CHILD) then
   ierr = PHAML_INTERNAL_ERROR
   call fatal("evaluate_oldsoln could not find a containing element marked as oldleaf.")
   stop
endif

! look at the barycentric coordinates of (x,y) in each child,
! until one is found where they are all positive

elem = 0
do i=1,MAX_CHILD
   if (all(barycentric(x,y,children(i),grid,no_det=.true.) >= &
           -roundoff_fudge*epsilon(0.0_my_real))) then
      call find_old_containing_leaf(x,y,children(i),elem,grid)
      exit
   endif
end do

end subroutine find_old_containing_leaf

!          --------
subroutine copy_old(grid)
!          --------

!----------------------------------------------------
! This routine copies solution into oldsoln
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), intent(inout) :: grid
!----------------------------------------------------
! Local variables:

integer :: lev, elem, astat, i, edge
!----------------------------------------------------
! Begin executable code

do lev=1,grid%nlev

! for each element

   elem = grid%head_level_elem(lev)
   do while (elem /= END_OF_LIST)
      if (grid%element(elem)%isleaf) then
         if (associated(grid%element(elem)%solution)) then

! if the element is a leaf and has a solution, then allocate oldsoln to the
! right size and copy solution to it

            if (associated(grid%element(elem)%oldsoln)) then
               if (size(grid%element(elem)%oldsoln,dim=1) /= &
                   size(grid%element(elem)%solution,dim=1) .or. &
                   size(grid%element(elem)%oldsoln,dim=2) /= &
                   size(grid%element(elem)%solution,dim=2) .or. &
                   size(grid%element(elem)%oldsoln,dim=3) /= &
                   size(grid%element(elem)%solution,dim=3)) then
                  deallocate(grid%element(elem)%oldsoln,stat=astat)
               endif
            endif
            if (.not. associated(grid%element(elem)%oldsoln)) then
               allocate(grid%element(elem)%oldsoln( &
                          size(grid%element(elem)%solution,dim=1), &
                          size(grid%element(elem)%solution,dim=2), &
                          size(grid%element(elem)%solution,dim=3)),stat=astat)
               if (astat /= 0) then
                  ierr = ALLOC_FAILED
                  call fatal("memory allocation failed in copy_old")
                  stop
               endif
            endif
            grid%element(elem)%oldsoln = grid%element(elem)%solution
            grid%element(elem)%oldleaf = .true.

         else

! if element is a leaf and does not have a solution, make sure it does not
! have oldsoln either

            if (associated(grid%element(elem)%oldsoln)) then
               deallocate(grid%element(elem)%oldsoln,stat=astat)
            endif
            grid%element(elem)%oldleaf = .true.

         endif
      else

! if element is not a leaf, make sure it does not have an oldsoln or say
! that it is an oldleaf

         if (associated(grid%element(elem)%oldsoln)) then
            deallocate(grid%element(elem)%oldsoln,stat=astat)
         endif
         grid%element(elem)%oldleaf = .false.

      endif

! for each edge

      do i=1,EDGES_PER_ELEMENT
         edge = grid%element(elem)%edge(i)

         if (grid%element(elem)%isleaf) then
            if (associated(grid%edge(edge)%solution)) then

! if the element is a leaf and the edge solution is allocated, make sure the
! edge oldsoln is the same size and copy solution to it

               if (associated(grid%edge(edge)%oldsoln)) then
                  if (size(grid%edge(edge)%oldsoln,dim=1) /= &
                      size(grid%edge(edge)%solution,dim=1) .or. &
                      size(grid%edge(edge)%oldsoln,dim=2) /= &
                      size(grid%edge(edge)%solution,dim=2) .or. &
                      size(grid%edge(edge)%oldsoln,dim=3) /= &
                      size(grid%edge(edge)%solution,dim=3)) then
                     deallocate(grid%edge(edge)%oldsoln,stat=astat)
                  endif
               endif
               if (.not. associated(grid%edge(edge)%oldsoln)) then
                  allocate(grid%edge(edge)%oldsoln( &
                             size(grid%edge(edge)%solution,dim=1), &
                             size(grid%edge(edge)%solution,dim=2), &
                             size(grid%edge(edge)%solution,dim=3)),stat=astat)
                  if (astat /= 0) then
                     ierr = ALLOC_FAILED
                     call fatal("memory allocation failed in copy_old")
                     stop
                  endif
               endif
               grid%edge(edge)%oldsoln = grid%edge(edge)%solution

            else

! if element is a leaf and edge does not have a solution, make sure it does not
! have oldsoln either

               if (associated(grid%edge(edge)%oldsoln)) then
                  deallocate(grid%edge(edge)%oldsoln,stat=astat)
               endif
            endif

! if element is not a leaf, don't change the edge because the neighbor might
! be a leaf and want the edge set, and it won't cause any damage to leave it

         endif
      end do ! next edge

      elem = grid%element(elem)%next
   end do ! next element

end do ! next level

! make sure the oldsoln is the same size as solution and copy

if (associated(grid%vertex_oldsoln)) then
   if (size(grid%vertex_oldsoln,dim=1) /= &
       size(grid%vertex_solution,dim=1) .or. &
       size(grid%vertex_oldsoln,dim=2) /= &
       size(grid%vertex_solution,dim=2) .or. &
       size(grid%vertex_oldsoln,dim=3) /= &
       size(grid%vertex_solution,dim=3)) then
      deallocate(grid%vertex_oldsoln,stat=astat)
   endif
endif
if (.not. associated(grid%vertex_oldsoln)) then
   allocate(grid%vertex_oldsoln(size(grid%vertex_solution,1), &
            size(grid%vertex_solution,2),size(grid%vertex_solution,3)), &
            stat=astat)
   if (astat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in copy_old")
      stop
   endif
endif
grid%vertex_oldsoln = grid%vertex_solution

grid%oldsoln_exists = .true.

end subroutine copy_old

!          ---------------------
subroutine set_grid_for_old_soln(grid)
!          ---------------------

!----------------------------------------------------
! This routine sets grid_for_old_soln so that the grid is accessible even
! though it cannot be passed through the user routines pdecoef, etc.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), target :: grid

!----------------------------------------------------
! Begin executable code

grid_for_old_soln => grid

end subroutine set_grid_for_old_soln

end module evaluate

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
 MODULE curv_db

   ! stores information of curves used for Loads & Prescribed Velocities
   ! created/modified in RDLOAD & RDVEL0
   ! used by APLFOL EXPLIT
   ! by DUMPIN RESTAR

   USE param_db,ONLY: mnam
   IMPLICIT NONE

   SAVE

   TYPE curpar    !Definition: List of nodes (edges, level structure, etc)
     CHARACTER (len=mnam) :: name ! curve name (label)
     INTEGER(kind=4):: dim
     REAL(kind=8),POINTER:: vec(:)
     TYPE (curpar),POINTER:: next
   END TYPE curpar

   INTEGER (kind=4) ::  nlcur   !number of curves
   INTEGER (kind=4) ::  numpac  !Total number of parameters defining curves
   CHARACTER (len=mnam), POINTER :: curnam(:) !(numpac) curve names (labels)
   REAL (kind=8),POINTER :: curpa(:)          !(numpac) parameters defining curves
   TYPE(curpar),POINTER:: headcp, tailcp     !pointers to curve data

 CONTAINS

   FUNCTION getcun (cnam)

     !   returns curve number

     USE lispa0, ONLY :  lures
     IMPLICIT NONE
     INTEGER (kind=4)  :: getcun               ! curve number
     CHARACTER (len=*), INTENT(IN) :: cnam     ! curve name

     LOGICAL          :: found
     INTEGER (kind=4) :: i

     found = .FALSE.
     DO i = 1,nlcur
       found = (TRIM(cnam) == TRIM(curnam(i))) ! compare searched name with the name in db
       IF (found) THEN
         getcun = i                ! returned number of the curve found
         RETURN
       END IF
     END DO

     WRITE(lures,"('Curve label',A)",ERR=9999) TRIM(cnam)
     WRITE(*    ,"('Curve label',A)") TRIM(cnam)
     CALL runend('GETCUN: Curve Label not found')

   RETURN
   9999 CALL runen2('')
   END FUNCTION getcun

   SUBROUTINE dump_curv

     !   dumps curve database
     IMPLICIT NONE

     WRITE (50,ERR=9999) nlcur, numpac   !number of curves  & size of curpa
     WRITE (50,ERR=9999) curnam, curpa   !curve names & curve parameters

     RETURN
     9999 CALL runen2('')
   END SUBROUTINE dump_curv

   SUBROUTINE rest_curv

     !   restores curve database

     IMPLICIT NONE
     TYPE(curpar),POINTER:: aux
     INTEGER :: i,ib,ie

     READ (51) nlcur, numpac     !number of curves  & size of curpa
     ALLOCATE ( curnam(nlcur), curpa(numpac) ) !get memory
     READ (51) curnam, curpa     !curve names & curve parameters
     ! Construct data base from data
     DO i=1,nlcur                !for each curve
       ALLOCATE (aux)            !get memory
       aux%name= curnam(i)       !pass curve name
       ib = curpa(i)             !first position in curpa
       IF( i == nlcur )THEN      !for the last curve
         ie = numpac             !last position in curpa
       ELSE
         ie = curpa(i+1) - 1     !last position for this curve
       END IF
       aux%dim = ie - ib + 1     !number of parameters for this curve
       ALLOCATE( aux%vec(aux%dim) )   !get memory
       aux%vec = curpa(ib:ie)    !pass curve parameters
       CALL add_cur(aux,headcp,tailcp)   !add in list
     END DO

     RETURN

   END SUBROUTINE rest_curv

   SUBROUTINE gen_curpa

     !   generates array curpa

     IMPLICIT NONE
     TYPE(curpar),POINTER:: aux
     INTEGER :: i,ib,ie

     IF( ASSOCIATED(curnam) )DEALLOCATE(curnam) !remove if array exists
     IF( ASSOCIATED(curpa ) )DEALLOCATE(curpa ) !remove if array exists
     nlcur  = 0          !initializes number of curves
     numpac = 0          !initializes number of parameters
     aux => headcp       !point to first curve
     DO              !first loop for each curve
       IF( .NOT. ASSOCIATED(aux) )EXIT  !exit if last curve
       nlcur  = nlcur  + 1              !increase number of curves
       numpac = numpac + aux%dim + 1    !increase number of parameters
       aux => aux%next                  !point to next curve
     END DO

     ALLOCATE (curnam(nlcur), curpa(numpac) ) !get memory
     aux => headcp       !point to first curve
     ib = nlcur+1        !pointer to first curva parameters
     DO i=1,nlcur    !second loop for each curve
       curnam(i) = aux%name     !curname
       curpa(i)  = ib           !pointer to curve parameters in curpa
       ie = ib + aux%dim - 1    !last position
       curpa(ib:ie) = aux%vec   !pass arguments
       IF( INT(curpa(ib)) == 3 )curpa(ib+4) = ib+5 !initializes actual position
       ib = ie + 1              !new first position
       aux => aux%next          !point to next curve
     END DO

     RETURN

   END SUBROUTINE gen_curpa

!-----------------------------------------------------------------------
  SUBROUTINE add_cur(new,head,tail)
!-----------------------------------------------------------------------
  !This subroutine adds data to the end of the list
  TYPE(curpar),POINTER:: new, head, tail

  IF (.NOT.ASSOCIATED(head)) THEN
    !list is empty, start it
    head => new
    tail => new
  ELSE
    !Check if a list is empty
    tail%next => new
    tail => new
  END IF
  NULLIFY (tail%next)

  RETURN
  END SUBROUTINE add_cur

!-----------------------------------------------------------------------
  SUBROUTINE del_cur(curname)
!-----------------------------------------------------------------------
  !Delete a curve from the list
  IMPLICIT NONE

  CHARACTER (len=mnam), INTENT(IN) :: curname     ! curve name
  TYPE(curpar),POINTER :: aux,prev

  aux => headcp
  NULLIFY(prev) !first curve in the list
  DO
    IF( .NOT.ASSOCIATED (aux) )THEN
      WRITE(55,"(' curve ',A,' does not exists ')",ERR=9999) TRIM(curname)
      RETURN
    END IF
    IF( TRIM(curname) == TRIM(aux%name) )THEN
      IF( ASSOCIATED(prev)) THEN  !not the first curve in the list
        prev%next => aux%next    !
        IF( .NOT.ASSOCIATED (prev%next) )tailcp => prev
      ELSE
        headcp => aux%next
      END IF
      DEALLOCATE(aux%vec)
      DEALLOCATE(aux)
      EXIT
    END IF
    prev => aux
    aux => aux%next
  END DO

  RETURN
  9999 CALL runen2('')
  END SUBROUTINE del_cur

!!-----------------------------------------------------------------------
!  SUBROUTINE dalloc_cur(head,tail)
!!-----------------------------------------------------------------------
!  !Delete a list
!  TYPE(curpar),POINTER:: head,tail
!  !Local variables
!  TYPE(curpar),POINTER:: aux
!
!  IF (.NOT.ASSOCIATED(head)) RETURN
!
!  !Search the starting point
!  NULLIFY(tail)
!  DO
!    IF (.NOT.ASSOCIATED(head)) EXIT
!    aux => head%next
!    CALL del_cur(head)
!    head => aux
!  END DO
!
!  RETURN
!  END SUBROUTINE dalloc_cur

 END MODULE curv_db

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
c
c
c=========================================================
      subroutine out1(meqn,mbc,mx,xlower,dx,q,t,iframe,aux,maux,outaux)
c=========================================================
c
c     # Output the results for a general system of conservation laws
c     # in 1 dimension
c
c     # Write the results to the file fort.q<iframe>
c     # Use format required by matlab script  plotclaw1.m
c     # The same format is used by the amrclaw package.
c     # Here it's adapted to output just the single grid.
c
c     # set outaux = .true. to also output the aux arrays to fort.a<iframe>
c
c
      implicit double precision (a-h,o-z)
      dimension q(meqn,1-mbc:mx+mbc)
      dimension aux(maux,1-mbc:mx+mbc)
      character*10 fname1, fname2, fname3
      logical outaux

c
c     # Write the results to the file fort.q<iframe>
c     # Use format required by matlab script  plotclaw1.m
c     # The same format is used by the amrclaw package.
c     # Here it's adapted to output just the single grid.
c
c     # first create the file name and open file
c
      fname1 = 'fort.qxxxx'
      fname2 = 'fort.txxxx'
      fname3 = 'fort.axxxx'
      nstp = iframe
      do 55 ipos = 10, 7, -1
         idigit = mod(nstp,10)
         fname1(ipos:ipos) = char(ichar('0') + idigit)
         fname2(ipos:ipos) = char(ichar('0') + idigit)
         fname3(ipos:ipos) = char(ichar('0') + idigit)
         nstp = nstp / 10
 55   continue

      open(unit=50,file=fname1,status='unknown',form='formatted')
      open(unit=60,file=fname2,status='unknown',form='formatted')

c
c     # the following parameters are used in amrclaw where there are
c     # multiple grids.  Here they are all set to 1:
      ngrids = 1
      mptr = 1
      level = 1

      write(50,1001) mptr,level,mx
 1001 format(i5,'                 grid_number',/,
     &     i5,'                 AMR_level',/,
     &     i5,'                 mx')

      write(50,1002) xlower,dx
 1002 format(e18.8,'    xlow', /,
     &     e18.8,'    dx', /)
c
      do 10 i=1,mx
         do m=1,meqn
c     # exponents with more than 2 digits cause problems reading
c     # into matlab... reset tiny values to zero:
            if (dabs(q(m,i)) .lt. 1d-99) q(m,i) = 0.d0
         enddo
c
         write(50,1005) (q(m,i), m=1,meqn)
 1005    format(50e26.16)
c
 10   continue
      write(50,*) ' '
      write(50,*) ' '

      if (outaux) then
c     # also output the aux arrays:
         open(unit=70,file=fname3,status='unknown',form='formatted')
         write(70,1001) mptr,level,mx
         write(70,1002) xlower,dx
         do 110 i=1,mx
            do m=1,maux
c     # exponents with more than 2 digits cause problems reading
c     # into matlab... reset tiny values to zero:
               if (dabs(aux(m,i)) .lt. 1d-99) aux(m,i) = 0.d0
            enddo
c
            write(70,1005) (aux(m,i), m=1,maux)
c
 110     continue
         write(70,*) ' '
         close(unit=70)
      endif

      write(60,1000) t,meqn,ngrids,maux,1

 1000 format(e26.16,'    time', /,
     &     i5,'                 meqn'/,
     &     i5,'                 ngrids'/,
     &     i5,'                 maux'/,
     &     i5,'                 ndim'/,/)
c

      close(unit=50)
      close(unit=60)

      return
      end

ext: f

max_stars_repo_name: ClimateImpactLab/classic

max_issues_repo_name: ClimateImpactLab/classic

max_forks_repo_name: cjvogl/classic


text:
module m2
  use m3
  use m4
  implicit none
contains
  subroutine proc2
    print *, 'proc2'
    call proc3
    call proc4
    include 'i1.f90'
  end subroutine
end module

max_stars_repo_name: gronki/fortdep

max_issues_repo_name: gronki/fortdep

max_forks_repo_name: gronki/fortdep


! -----------------------------------------------------------------------------
!
! Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young,
! Peter Jimack, Mike Pilling
!
! Copyright (c) 2017 Sam Cox, Roberto Sommariva
!
! This file is part of the AtChem2 software package.
!
! This file is covered by the MIT license which can be found in the file
! LICENSE.md at the top level of the AtChem2 distribution.
!
! -----------------------------------------------------------------------------

! ******************************************************************** !
! ATCHEM2 -- MODULE constraintFunctions
!
! This module handles the addition and removal of constrained species
! from vectors of concentrations, and calculation of jFac. It also
! handles calculation of environment variabe values at a given time t.
! Also implements getEnvVarNum() as a helper function to return the
! number of the environment variable using the name as a key.
! ******************************************************************** !
module constraint_functions_mod
contains

  function calcPhotolysisRaw( l, m, n, tf ) result ( photolysis )
    use types_mod
    use zenith_data_mod, only : cosx, secx
    implicit none

    real(kind=DP), intent(in) :: l, m, n, tf
    real(kind=DP) :: photolysis

    photolysis = l * cosx ** m * exp( -n * secx ) * tf

    return
  end function calcPhotolysisRaw

  ! ----------------------------------------------------------------- !
  ! Calculate the photolysis rate of photolysis number i from the
  ! photolysis equations, using the current zenith data values
  function calcPhotolysis( i ) result ( photolysis )
    use types_mod
    use photolysis_rates_mod, only : cl, cmm, cnn, transmissionFactor
    implicit none

    integer(kind=NPI), intent(in) :: i
    real(kind=DP) :: photolysis

    photolysis = calcPhotolysisRaw( cl(i), cmm(i), cnn(i), transmissionFactor(i) )

    return
  end function calcPhotolysis

  ! ----------------------------------------------------------------- !
  ! Calculate the value of jFac by comparing the constrained value of
  ! jFacSpecies to the value calculated from the photolysis equations.
  ! If the constrained rate is zero, or the value of cosx is below the
  ! threshold (represented by cosx_below_threshold==.true.) then set
  ! jFac to zero.
  subroutine calcJFac( t, jFac )
    use types_mod
    use zenith_data_mod
    use photolysis_rates_mod, only : jFacSpecies, numConstrainedPhotoRates, &
                                     usePhotolysisConstants, constrainedPhotoNames, &
                                     jFacL, jFacM, jFacN, jFacTransmissionFactor
    use interpolation_functions_mod, only : getConstrainedPhotoRatesAtT
    use interpolation_method_mod, only : getConditionsInterpMethod
    implicit none

    real(kind=DP), intent(in) :: t
    real(kind=DP), intent(out) :: jFac
    real(kind=DP) :: JFacSpeciesAtT
    integer(kind=NPI) :: basePhotoRateNum, i
    logical :: firstTime = .true.

    if ( firstTime .eqv. .true. ) then
      write (*, '(2A)') ' basePhotoRate: ', trim( jFacSpecies )
      firstTime = .false.
    end if

    ! GET INDEX OF basePhotoRate SPECIES IN PHOTO CONSTRAINT ARRAY
    basePhotoRateNum = 0
    do i = 1, numConstrainedPhotoRates
      if ( trim( constrainedPhotoNames(i) ) == trim( jFacSpecies ) ) then
        basePhotoRateNum = i
      end if
    end do

    if ( basePhotoRateNum == 0 ) then
      write (*, '(2A)') ' Error! Missing constrained photo rate data for the JFAC species provided: ', trim( jFacSpecies )!
      stop 2
    end if

    ! GET CURRENT VALUE OF basePhotoRate

    call getConstrainedPhotoRatesAtT( t, basePhotoRateNum, JFacSpeciesAtT )

    if ( JFacSpeciesAtT == 0 ) then
      jFac = 0
    else
      if ( usePhotolysisConstants .eqv. .false. ) then
        if ( cosx_below_threshold .eqv. .false. ) then
          jFac = JFacSpeciesAtT / calcPhotolysisRaw( jFacL, jFacM, jFacN, jFacTransmissionFactor )
        else
          jFac = 0
        end if
      else
        write (*, '(A)') ' Error! JFAC should not be used, as constant photolysis rates have been provided.'
        stop 2
      end if
    end if
    return
  end subroutine calcJFac

  ! TODO: As the numbers of the constrained species will be constant, it may be possible
  ! TODO: to speed this up by computing a vector of bools once then apply this as a mask
  ! TODO: each time. Same for the removal.
  ! ----------------------------------------------------------------- !
  ! Take in z, the vector of concentrations of unconstrained species,
  ! and add the concentrations of the constrained species. Return this
  ! in vector x.
  subroutine addConstrainedSpeciesToProbSpec( z, constrainedConcentrations, constrainedSpecs, x )
    use types_mod
    implicit none

    real(kind=DP), intent(in) :: z(*), constrainedConcentrations(:)
    integer(kind=NPI), intent(in) :: constrainedSpecs(:)
    real(kind=DP), intent(out) :: x(:)
    integer(kind=NPI) :: zCounter, i, j, speciesConstrained

    ! This fills x with the contents of z, plus the contents of
    ! constrainedConcentrations, using constrainedSpecs as the key to
    ! which species are constrained.
    if ( size( constrainedConcentrations ) /= size( constrainedSpecs ) ) then
      stop 'size( constrainedConcentrations ) /= size( constrainedSpecs ) in addConstrainedSpeciesToProbSpec().'
    end if
    zCounter = 1
    do i = 1, size( x )
      speciesConstrained = 0
      do j = 1, size( constrainedSpecs )
        if ( i == constrainedSpecs(j) ) then
          speciesConstrained = j
          exit
        end if
      end do
      if ( speciesConstrained > 0 ) then
        x(i) = constrainedConcentrations(speciesConstrained)
      else if ( speciesConstrained == 0 ) then
        x(i) = z(zCounter)
        zCounter = zCounter + 1
      else
        stop 'Error adding constrained values to measured values'
      end if
    end do
    return
  end subroutine addConstrainedSpeciesToProbSpec

  ! ----------------------------------------------------------------- !
  ! Take in x, the vector of concentrations of all species, and remove
  ! the concentrations of the constrained species. Return the remainder
  ! in vector z.
  subroutine removeConstrainedSpeciesFromProbSpec( x, constrainedSpecs, z )
    use types_mod
    implicit none

    real(kind=DP), intent(in) :: x(:)
    integer(kind=NPI), intent(in) :: constrainedSpecs(:)
    real(kind=DP), intent(inout) :: z(*)
    integer(kind=NPI) :: zCounter, i, k
    logical :: speciesConstrained

    zCounter = 1
    ! loop through x(), check if its items are in constrainedSpecs
    do i = 1, size( x )
      speciesConstrained = .false.
      do k = 1, size( constrainedSpecs )
        if ( i == constrainedSpecs(k) ) then
          ! exit loop if in constrainedSpecs
          speciesConstrained = .true.
          exit
        end if
      end do
      ! if item is not in constrainedSpecs, then add species to z
      if ( speciesConstrained .eqv. .false. ) then
        z(zCounter) = x(i)
        zCounter = zCounter + 1
      end if
    end do
    return
  end subroutine removeConstrainedSpeciesFromProbSpec

  ! ----------------------------------------------------------------- !
  ! Return the values of all environment variables at time t, taking
  ! into account whether each variable is constrained, fixed value,
  ! calculated from a formula, or given a default value. Also take
  ! into account that some quantities are dependent on others, so
  ! enforce an order to the calculations.
  subroutine getEnvVarsAtT( t )
    use, intrinsic :: iso_fortran_env, only : stderr => error_unit
    use types_mod
    use storage_mod, only : maxEnvVarNameLength
    use env_vars_mod
    use interpolation_functions_mod, only : getConstrainedEnvVarAtT
    use interpolation_method_mod, only : getConditionsInterpMethod
    use atmosphere_functions_mod, only : calcAirDensity, convertRHtoH2O
    use solar_functions_mod, only : calcDec, calcZenith
    implicit none

    real(kind=DP), intent(in) :: t
    real(kind=DP) :: this_env_val
    integer(kind=NPI) :: envVarNum, orderedEnvVarNum
    character(len=maxEnvVarNameLength) :: this_env_var_name, orderedEnvVarNames(size( envVarNames ))
    logical :: pressure_set, rh_set, temp_set

    ! loop over each environment variable, in a defined order, rather
    ! than just the order of envVarNames, so we can ensure the right
    ! ones are calculated before each other.
    !
    ! Specifically, handle that some variables need to be set before
    ! others (e.g., pressure, temperature and RH before H2O, and DEC
    ! before JFAC). Currently, this relies on
    ! environmentVariables.config having exactly the lines relating to
    ! these 10 variables.
    !
    ! To add another environment variable, the user would need to add
    ! that line to environmentVariables.config, and then add this as
    ! element 11 in the orderedEnvVarNames initialisation below.  Its
    ! treatment needs defining in each of cases 1-3 and default below.

    if ( size( envVarNames ) /= 10 ) then
      write(stderr,*) 'size( envVarNames ) /= 10 in getEnvVarsAtT().'
    end if
    orderedEnvVarNames(1) = 'PRESS'
    orderedEnvVarNames(2) = 'TEMP'
    orderedEnvVarNames(3) = 'M'
    orderedEnvVarNames(4) = 'RH'
    orderedEnvVarNames(5) = 'H2O'
    orderedEnvVarNames(6) = 'BLHEIGHT'
    orderedEnvVarNames(7) = 'DEC'
    orderedEnvVarNames(8) = 'JFAC'
    orderedEnvVarNames(9) = 'DILUTE'
    orderedEnvVarNames(10) = 'ROOF'

    pressure_set = .false.
    rh_set = .false.
    temp_set = .false.

    do orderedEnvVarNum = 1, size( envVarNames )
      ! loop over in a defined order, then find which number that is
      ! in the unordered list that comes from the input file
      this_env_var_name = orderedEnvVarNames(orderedEnvVarNum)
      envVarNum = getEnvVarNum( trim( this_env_var_name ) )

      ! Find which type it is (calc, constrained, fixed, other)
      select case ( envVarTypesNum(envVarNum) )
        case ( 1 ) ! CALC
          select case ( this_env_var_name )
            case ( 'PRESS', 'TEMP', 'RH', 'BLHEIGHT', 'DILUTE', 'ROOF' )
              write (stderr,*) 'getEnvVarsAtT(): No calculation available for ' // trim( this_env_var_name )
              stop
            case ( 'M' )
              if ( ( temp_set .eqv. .true. ) .and. ( pressure_set .eqv. .true. ) ) then
                this_env_val = calcAirDensity( currentEnvVarValues( getEnvVarNum( 'PRESS' ) ), &
                                               currentEnvVarValues( getEnvVarNum( 'TEMP' ) ) )
              else
                write (stderr,*) 'getEnvVarsAtT(): calcAirDensity() called, but no value is yet given to either TEMP, or PRESS.'
                stop
              end if
            case ( 'H2O' )
              if ( ( rh_set .eqv. .true. ) .and. ( temp_set .eqv. .true. ) .and. ( pressure_set .eqv. .true. ) ) then
                this_env_val = convertRHtoH2O( currentEnvVarValues( getEnvVarNum( 'RH' ) ), &
                                               currentEnvVarValues( getEnvVarNum( 'TEMP' ) ), &
                                               currentEnvVarValues( getEnvVarNum( 'PRESS' ) ) )
              else
                write (stderr,*) 'getEnvVarsAtT(): convertRHtoH2O() called, but no value is yet given to either RH, TEMP, or PRESS.'
                stop
              end if
            case ( 'DEC' )
              this_env_val = calcDec()
            case ( 'JFAC' )
              call calcJFac( t, this_env_val )
            case default
              write (stderr,*) 'getEnvVarsAtT(): invalid environment name ' // trim( this_env_var_name )
              stop
          end select

        case ( 2 ) ! CONSTRAINED
          if ( this_env_var_name == 'ROOF') then
            write (stderr,*) 'ROOF should not be constrained.'
            stop
          end if
          call getConstrainedEnvVarAtT( t, envVarNum, this_env_val )

          if (this_env_var_name == 'PRESS') pressure_set = .true.
          if (this_env_var_name == 'TEMP') temp_set = .true.
          if (this_env_var_name == 'RH') rh_set = .true.

        case ( 3 ) ! FIXED VALUE
          this_env_val = envVarFixedValues(envVarNum)

          if (this_env_var_name == 'PRESS') pressure_set = .true.
          if (this_env_var_name == 'TEMP') temp_set = .true.
          if (this_env_var_name == 'RH') rh_set = .true.

        case default ! DEFAULT VALUES
          select case ( this_env_var_name )
            case ( 'TEMP' )
              this_env_val = 298.15_DP
              temp_set = .true.
            case ( 'H2O' )
              this_env_val = 3.91e+17_DP
            case ( 'PRESS' )
              this_env_val = 1013.25_DP
            case ( 'BLHEIGHT', 'RH', 'DILUTE' )
              this_env_val = -1.0_DP
            case ( 'DEC' )
              this_env_val = 0.41_DP
            case ( 'JFAC')
              this_env_val = 1.0_DP
            case ( 'ROOF' )
              if ( this_env_var_name == 'ROOF') then
                write (stderr,*) 'ROOF should not be a default value.'
                stop
              end if
            case default
              write (stderr,*) 'getEnvVarsAtT(): invalid environment name ' // trim( this_env_var_name )
              stop
          end select
      end select
      currentEnvVarValues(envVarNum) = this_env_val

      ! Copy this_env_var_name to the correct output variable
      select case ( this_env_var_name )
        case ( 'TEMP', 'RH', 'H2O', 'PRESS', 'M', 'BLHEIGHT', 'DILUTE', 'JFAC', 'ROOF' )
        case ( 'DEC' )
          call calcZenith( t, this_env_val )
        case default
          write(stderr,*) 'getEnvVarsAtT(): invalid environment name ' // trim( this_env_var_name )
          stop
      end select
    end do

    return
  end subroutine getEnvVarsAtT

  ! ----------------------------------------------------------------- !
  ! Set envVarNum to the index of name within enVarNames
  function getEnvVarNum( name ) result ( envVarNum )
    use, intrinsic :: iso_fortran_env, only : stderr => error_unit
    use types_mod
    use env_vars_mod, only : envVarNames, numEnvVars
    implicit none

    character(len=*), intent(in) :: name
    integer(kind=SI) :: envVarNum
    integer(kind=SI) :: i

    envVarNum = 0
    do i = 1, numEnvVars
      if ( name == trim( envVarNames(i) ) ) then
        envVarNum = i
        exit
      end if
    end do
    if ( envVarNum == 0 ) then
      write (stderr,*) 'The name ' // trim( name ) // ' is not found in getEnvVarNum().'
      stop
    end if

    return
  end function getEnvVarNum

end module constraint_functions_mod


program main
    use PreprocessingClass

    type(Preprocessing_) :: leaf

    call leaf%import(Name="Tutorial/InputData/grass_leaf")
    call leaf%removeBC(Initial=.true.)
    call leaf%SetSizeOfBC(Initial=.true. , NumOfValue=1)
    call leaf%SetBC(Initial=.true.,   val=0.0d0,val_id=1)
    call leaf%Export(Name="Tutorial/InputData/grass_leaf_BC_changed")

end program


text:
Describe ColeJohnson here.

Im in the BandUh Banduh!, its badass!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module points_basic

implicit none

interface
module real function point_dist(ax, ay, bx, by)
real, intent(in) :: ax, ay, bx, by
end function point_dist
end interface

end module points_basic

program test_points_basic
use, intrinsic :: iso_fortran_env, only: stderr=>error_unit
use points_basic, only: point_dist
implicit none

real :: dist

dist = point_dist(1.,1., 3.,5.)

print *,'distance', dist

if (abs(dist - 4.47213602) >= 1e-5) error stop 'excessive error in computation'

end program

max_stars_repo_name: scivision/fortran-submodule

max_issues_repo_name: scivision/fortran-submodule

max_forks_repo_name: scivision/fortran-submodule


text:
      program testhelmbhrouts
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     This program tests the difference function
c     evaluation routines in helmbhrouts.f
c     against high-precision computations
c     produced in Mathematica.
c
c     The test data files are formatted as
c     follows
c
c     test1pars.txt contains the parameters
c
c     ntheta, nrscale, nterms
c
c     The data loops over the zk, r, j values
c     zk is outermost, then r, then j is innermost
c
c     zk = e^(i*k*2*pi/ntheta), k = 1,ntheta
c     r = 10^(offset+k), k = 1,nrscale (offset = -15)
c     j = 0,nterms
c
c     each row of test1vals.txt is
c
c     j Re[zk] Im[zk] r rscale Re[diffout_j[zk,r]] Im[diffout_j[zk,r]
c     Re[diffin_j[zk,r]] Im[diffin_j[zk,r]
c
c     Where
c
c     diffout_j[zk,r] = H_j[zk*r]+i*2^(j-1)*(j-1)!*2/(pi*(zk*r)^j)
c
c     diffin_j[zk,r] = J_j[zk*r]-(zk*r)^j/(2^j*j!)
c

      implicit real *8 (a-h,o-z)
      real *8 dtoobig, dtoosmall
      complex *16 eye
      data eye / (0.0d0,1.0d0) /
      complex *16, allocatable :: zks(:,:)
      complex *16, allocatable :: diffout(:,:,:), diffin(:,:,:)
      real *8, allocatable :: rs(:,:), rscales(:,:)
      integer, allocatable :: iftest(:,:,:)
      complex *16, allocatable :: hvec(:)
      complex *16, allocatable :: ders(:)
      complex *16, allocatable :: diffs(:)
      real *8, allocatable :: errmax(:), rsmax(:)
      complex *16, allocatable :: zksmax(:)
      real *8, allocatable :: errmaxin(:), rsmaxin(:)
      complex *16, allocatable :: zksmaxin(:)

c     parameters
      dtoobig = 1.0d300
      dtoosmall = 1.0d-300

      open(unit=20,file="../test/test-data/test1pars.txt")
      open(unit=21,file="../test/test-data/test1vals.txt")

      read(20,*) ntheta
      read(20,*) nrscale
      read(20,*) nterms

      write(*,*) "PARAMETERS (ntheta, nrscale, nterms)"
      write(*,*) ntheta, nrscale, nterms

      allocate(zks(nrscale,ntheta),rs(nrscale,ntheta),
     1     rscales(nrscale,ntheta))
      allocate(diffout(0:nterms,nrscale,ntheta),
     1     diffin(0:nterms,nrscale,ntheta))
      allocate(hvec(0:nterms),ders(0:nterms),diffs(0:nterms))
      allocate(zksmax(0:nterms),errmax(0:nterms),rsmax(0:nterms))
      allocate(zksmaxin(0:nterms),errmaxin(0:nterms),
     1     rsmaxin(0:nterms))

      write(*,*) "READING IN DATA ..."
      do i = 1,ntheta
         do j = 1,nrscale
            do k = 0,nterms
               read(21,*,iostat=ios) ktemp, dkr, dki, dr, drsc,
     1              dor, doi, dir, dii
               if (ios .ne. 0) then
                  write(*,*) "bad read"
                  stop
               endif
               if (ktemp .ne. k) write(*,*) "something went wrong"
               zks(j,i) = dkr+eye*dki
               rs(j,i) = dr
               rscales(j,i) = drsc
               diffout(k,j,i) = dor+eye*doi
               diffin(k,j,i) = dir+eye*dii
            enddo
         enddo
      enddo

      do i = 0,nterms
         errmax(i) = 0.0d0
         zksmax(i) = 0.0d0
         rsmax(i) = 0.0d0
         errmaxin(i) = 0.0d0
         zksmaxin(i) = 0.0d0
         rsmaxin(i) = 0.0d0
      enddo

      write(*,*) "RUNNING TEST (ON FIRST AND SECOND QUADRANT)..."

      do i = 1,ntheta/2-1
         do j = 1,nrscale-3

            call diffsloghank(rs(j,i),zks(j,i),rscales(j,i),
     1           diffs,ifders,ders,hvec,nterms)
            do k = 0,nterms
               err = abs(diffs(k)-diffout(k,j,i))/abs(diffout(k,j,i))
               if (err .gt. errmax(k)) then
                  errmax(k) = err
                  zksmax(k) = zks(j,i)
                  rsmax(k) = rs(j,i)
               endif
            enddo

            call diffszkjfun(rs(j,i),zks(j,i),rscales(j,i),
     1           diffs,ifders,ders,hvec,nterms)
            do k = 0,nterms
               err = abs(diffs(k)-diffin(k,j,i))/abs(diffin(k,j,i))
               if (err .gt. errmaxin(k)) then
                  errmaxin(k) = err
                  zksmaxin(k) = zks(j,i)
                  rsmaxin(k) = rs(j,i)
               endif
            enddo

         enddo
      enddo

      write(*,*) "WORST CASE ERROR BY MODE"

      do i = 0,nterms
         write(*,*) i, errmax(i), errmaxin(i)
      enddo

      stop
      end

max_stars_repo_name: askhamwhat/biharm-evals

max_issues_repo_name: askhamwhat/biharm-evals

max_forks_repo_name: askhamwhat/biharm-evals


text:
module calc_temp
    implicit none

    contains

    subroutine compute_profile(a_param, damping_param, pressure, nlayers, result)
        real(8), intent(in) :: a_param
        real(8), intent(in) :: damping_param
        real(8), intent(in) :: pressure(nlayers)
        integer, intent(in) :: nlayers
        real(8), intent(out) :: result(nlayers)

        result = a_param*exp(damping_param*log10(pressure))
    end subroutine
end module

max_stars_repo_name: ucl-exoplanets/taurex-plugin-examples

max_issues_repo_name: ucl-exoplanets/taurex-plugin-examples

max_forks_repo_name: ucl-exoplanets/taurex-plugin-examples


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! check for constant expression : 1.0q0.

program test
  real*16 :: r
  !{error "PGF90-S-0084-Illegal use of symbol  - KIND parameter"}
  r = 1.23456789123456789123456789123456789q0_16

end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


! { dg-do run }
! { dg-options " -std=f2003" }
! Checks internal file read/write of namelists
! (Fortran 2003 feature)
! PR fortran/28224
program nml_internal
  integer   :: i, j
  real      :: r
  namelist /nam/ i, j, r
  character(len=250) :: str

  i = 42
  j = -718
  r = exp(1.0)
  write(str,nml=nam)
  i = -33
  j = 10
  r = sin(1.0)
  read(str,nml=nam)
  if(i /= 42 .or. j /= -718 .or. abs(r-exp(1.0)) > 1e-5) STOP 1
end program nml_internal


!
! accerelation
!

module accerelation
  implicit none
  public aitken

contains

! Aitken's delta-squared accerelation
!
!   y = x0 - (Delta x0)^2 / Delta^2 x0
!     = x2 - (x2 - x1)^2 / (x2 - 2 x1 + x0)
!     = x0 - (x1 - x0)^2 / (x2 - 2 x1 + x0)
!
!         x0: a vector at the n-th step
!         x1: a vector at the n+1-st step
!         x2: a vector at the n+2-nd step
!         y: accelerated vector
!

  ! subroutine aitken(n, x0, incx0, x1, incx1, x2, incx2, y, incy)
  !   integer, intent(in) :: n, incx0, incx1, incx2, incy
  !   double precision, intent(in) :: x0(1:incx0,1:n), x1(1:incx1,1:n), x2(1:incx2,1:n)
  !   double precision, intent(out) :: y(1:incy,1:n)
  !
  !   ! double precision :: g(1:n), h(1:n)
  !   !
  !   ! g(1:n) = x0(1,1:n)
  !   ! daxpy(n, -1.0d0, x1, incx1, g, 1)
  !   !
  !   ! h(1:n) = x2(1,1:n)
  !   ! daxpy(n, -2.0d0, x1, incx1, h, 1)
  !   ! daxpy(n, 1.0d0, x0, incx0, h, 1)
  !   !
  !   ! g(1:n) = g(1:n) * g(1:n) / h(1:n)
  !   ! y(1,1:n) = x(1,1:n)
  !   ! daxpy(n, -1.0d0, g, 1, y, incy)
  !
  !   y(1,1:n) = x0(1,1:n) - (x1(1,1:n) - x0(1,1:n))^2 / (x2(1,1:n) - 2.0d0 * x1(1,1:n) + x0(1,1:n))
  ! end subroutine

  subroutine aitken(n, x0, incx0, x1, incx1, x2, incx2, y, incy)
    integer, intent(in) :: n, incx0, incx1, incx2, incy
    double precision, intent(in) :: x0(1:incx0,1:n), x1(1:incx1,1:n), x2(1:incx2,1:n)
    double precision, intent(out) :: y(1:incy,1:n)

    ! double precision :: g(1:n), h(1:n)
    !
    ! g(1:n) = x0(1,1:n)
    ! daxpy(n, -1.0d0, x1, incx1, g, 1)
    !
    ! h(1:n) = x2(1,1:n)
    ! daxpy(n, -2.0d0, x1, incx1, h, 1)
    ! daxpy(n, 1.0d0, x0, incx0, h, 1)
    !
    ! g(1:n) = g(1:n) * g(1:n) / h(1:n)
    ! y(1,1:n) = x(1,1:n)
    ! daxpy(n, -1.0d0, g, 1, y, incy)

    y(1,1:n) = x0(1,1:n) - (x1(1,1:n) - x0(1,1:n))^2 / (x2(1,1:n) - 2.0d0 * x1(1,1:n) + x0(1,1:n))
  end subroutine

end module accerelation


text:
      subroutine corfree(newcor,nstr,nmon,ipstr,ipestr,ipmon,ipemon)
      use tfstk
      use ffs
      use tffitcode
c      Release memories
      parameter (kfiles=3)
      parameter (kstack=3)
      integer*8 isa,isb,isc,ipyv,ippv,ipobs,ipiv,ipbckup1,ipemibuf,
     $     ipmon,ipstr,ipestr,ipemon,k,ip,ipdefv,itcoupste
      integer*8 iifnam,iidata
      common /mcfiles/icomf,nof(kfiles),iifnam(kfiles),iidata(kfiles),
     1                memax(kfiles)
      common /mcstack/icoma,ipnt(kstack),iistck(kstack),mstkmx(kstack)
      common /pstatis/na,nb,nc,isa,isb,isc
      common /cordefv/ ipdefv,idefvar
      include 'inc/common.inc'
      include 'inc/coroper.inc'
      common /corbtune/ipyv,ippv,ipobs,ipiv,ipbckup1,ipemibuf,nemitd,
     $     nememax,iter
      common /corcoup/ itcoupste
      newcor=0
      if(ibckup.ne.0) call tfree(ibckup)
      if(itmon.ne.0) then
        call tfree(ipmon)
        call tfree(ipemon)
        nmon=0
        nmona=0
        nmonact=0
        ipmon=0
        ipemon=0
        itmon=0
        itemon=0
      endif
      if(itstr.ne.0) then
        call tfree(ipstr)
        call tfree(ipestr)
        nstr=0
        nstra=0
        ipstr=0
        ipestr=0
        itstr=0
        itestr=0
      endif
      icomf=0
      do 11 i=1,kfiles
        do 10 j=1,nof(i)
          k=ilist(mod(j-1,2)+1,iifnam(i)+(j-1)/2)
          if(k.ne.0) call tfree(k)
          k=ilist(mod(j-1,2)+1,iidata(i)+(j-1)/2)
          if(k.ne.0) call tfree(k)
   10   continue
        if(iifnam(i).ne.0) call tfree(iifnam(i))
        if(iidata(i).ne.0) call tfree(iidata(i))
   11 continue
      icoma=0
      do 13 i=1,kstack
        do 12 j=1,ipnt(i)
          k=ilist(mod(j-1,2)+1,iistck(i)+(j-1)/2)
          if(k.ne.0) call tfree(k)
   12   continue
        if(iistck(i).ne.0) call tfree(int8(iistck(i)))
   13 continue
      na=0
      nb=0
      nc=0
      if(istope.ne.0) call tfree(istope)
      istope=0
      nstope=0
      if(ipdefv.ne.0) then
        do 14 i=1,idefvar
          ip=ilist(mod(i-1,2)+1,ipdefv+(i-1)/2)
          if(ip.ne.0) call tfree(ip)
 14     continue
        call tfree(ipdefv)
      endif
      if(ipemibuf.ne.0) then
        call tfree(ipemibuf)
        nemitd=0
      endif
      if(ipyv.ne.0) call tfree(ipyv)
      if(ippv.ne.0) call tfree(ippv)
      if(ipobs.ne.0) call tfree(ipobs)
      if(ipiv.ne.0) call tfree(ipiv)
      if(ipbckup1.ne.0) call tfree(ipbckup1)
      if(itcoupste.ne.0) call tfree(itcoupste)
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
C     ===========================================
      SUBROUTINE MY_GPXYZ(CPX,CPY,CPZ,CET,RAD,XV)
C     ===========================================
C     Auxiliary routine for angle calculation.
C
C     A.Quadt, May 1997.
C     -----------------------------------------

      IMPLICIT NONE

#include "caltru.inc"
#include "partap.inc"

C
C -------> Caltru is Fettab'd
C
      INTEGER    IERR
      REAL       ZP,YP,CPX,CPY,CPZ,CET,RAD,XV(3),THETA,PHI,
     1           SINPHI,SINTH,COSTH,COSPHI,XP

      CALL CCCXYZ (CALTRU_CELLNR,XP,YP,ZP,IERR)
      CALL CCCAPO (XP-XV(1),YP-XV(2),ZP-XV(3),RAD,THETA,PHI)
      RAD    = SQRT((XP-XV(1))**2+(YP-XV(2))**2+(ZP-XV(3))**2)
      SINPHI = SIN(PHI)
      COSPHI = COS(PHI)
      SINTH  = SIN(THETA)
      COSTH  = COS(THETA)

      CPX    = CALTRU_E*SINTH*COSPHI
      CPY    = CALTRU_E*SINTH*SINPHI
      CPZ    = CALTRU_E*COSTH
      CET    = CALTRU_E*SINTH

      RETURN
      END

ext: fpp


text:
	FUNCTION PR_HCDM  ( hh )
C************************************************************************
C* PR_HCDM								*
C*									*
C* This function returns the height in meters which corresponds to a	*
C* given code figure from WMO Code Table 1677.  If the given code	*
C* figure is found to be invalid, then RMISSD is returned.		*
C*									*
C* REAL PR_HCDM  ( HH )							*
C*									*
C* Input parameters:							*
C*	HH		CHAR* 		Code figure from WMO table 1677 *
C*									*
C* Output parameters:							*
C*	PR_HCDM		REAL		Height in meters  		*
C*									*
C**									*
C* Log:									*
C* J. Ator/NCEP		01/02						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	hh
C-----------------------------------------------------------------------
C
C*	Note that code figures 90 through 99 will also be mapped
C*	to RMISSD, since there is no clear way to set a single output
C*	value for the ranges represented by such code figures.
C
	PR_HCDM = RMISSD
C
	CALL ST_INTG  ( hh (1:2), ihh, ier )
	IF  ( ier .eq. 0 )  THEN
	    IF  ( ( ihh .ge. 0 ) .and. ( ihh .le. 50 ) )  THEN
		PR_HCDM = FLOAT ( ihh ) * 30.
	    ELSE IF  ( ( ihh .ge. 56 ) .and. ( ihh .le. 80 ) )  THEN
		PR_HCDM = FLOAT ( ihh - 50 ) * 300.
	    ELSE IF  ( ( ihh .ge. 81 ) .and. ( ihh .le. 89 ) )  THEN
		PR_HCDM = ( FLOAT ( ihh - 80 ) * 1500. ) + 9000.
	    END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!/*****************************************************************************
!*
!* Rokko: Integrated Interface for libraries of eigenvalue decomposition
!*
!* Copyright (C) 2012-2013 by Tatsuya Sakashita <t-sakashita@issp.u-tokyo.ac.jp>,
!*                            Synge Todo <wistaria@comp-phys.org>,
!*                            Tsuyoshi Okubo <t-okubo@issp.u-tokyo.ac.jp>
!*
!* Distributed under the Boost Software License, Version 1.0. (See accompanying
!* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
!*
!*****************************************************************************/

program frank_matrix
  use MPI
  use rokko
  implicit none
  integer :: dim
  type(rokko_distributed_matrix) :: mat, Z
  type(rokko_grid) :: grid
  type(rokko_parallel_dense_solver) :: solver
  type(rokko_localized_vector) :: w
  character(len=100) :: solver_name, tmp_str
  integer arg_len, status
  double precision, allocatable, dimension(:,:) :: array, array_tmp

  integer :: provided,ierr, myrank, nprocs
  integer :: i, j, proc

  call MPI_init_thread(MPI_THREAD_MULTIPLE, provided, ierr)
  call MPI_comm_rank(MPI_COMM_WORLD, myrank, ierr)
  call MPI_comm_size(MPI_COMM_WORLD, nprocs, ierr)

  if (command_argument_count().eq.2) then
     call get_command_argument(1, tmp_str, arg_len, status)
     solver_name = trim(tmp_str)
     call get_command_argument(2, tmp_str, arg_len, status)
     read(tmp_str, *) dim
  else
     write(*,'(A)') "Error: frank_array solver_name dimension"
     stop
  endif

  write(*,*) "solver name = ", trim(solver_name)
  write(*,*) "matrix dimension = ", dim

  call rokko_parallel_dense_solver_construct(solver, solver_name)
  call rokko_grid_construct(grid, MPI_COMM_WORLD, rokko_grid_row_major)

  call rokko_distributed_matrix_construct(mat, dim, dim, grid, solver, rokko_matrix_col_major)
  call rokko_distributed_matrix_construct(Z, dim, dim, grid, solver, rokko_matrix_col_major)
  call rokko_localized_vector_construct(w, dim)

  ! generate frank matrix as a localized matrix
  allocate(array(dim, dim))
  do i=1, dim
     do j=1, dim
        array(i,j) = dble(dim + 1 - max(i, j))
     end do
  end do
  allocate(array_tmp(dim, dim))

  call rokko_distributed_matrix_generate_array(mat,array)
  call rokko_distributed_matrix_print(mat)
  call rokko_parallel_dense_solver_diagonalize_distributed_matrix(solver, mat, w, Z)

  array = 0.0
  call rokko_all_gather(Z, array_tmp)
  array = matmul(transpose(array_tmp), array_tmp)
  call mpi_barrier(mpi_comm_world, ierr)

  !  print*, "array=", array
  do proc = 0, nprocs-1
     if (proc == myrank) then
        print*, "fmyrank=", myrank
     do i = 1, dim
        write(*,'(10f8.4)') (array(i, j), j=1, dim)
     end do
     endif
   call mpi_barrier(mpi_comm_world, ierr)
!   call sleep(0.1)
  end do

  if (myrank.eq.0) then
     write(*,*) "Computed Eigenvalues = "
     do i = 1, dim
        write(*,"(f30.20)") rokko_localized_vector_get(w ,i)
     enddo
  endif

  call mpi_barrier(mpi_comm_world, ierr)

  call rokko_distributed_matrix_destruct(mat)
  call rokko_distributed_matrix_destruct(Z)
  call rokko_localized_vector_destruct(w)
  call rokko_parallel_dense_solver_destruct(solver)
  call rokko_grid_destruct(grid)
  deallocate(array,array_tmp)

  call MPI_finalize(ierr)
end program frank_matrix

max_stars_repo_name: wistaria/rokko

max_issues_repo_name: wistaria/rokko

max_forks_repo_name: wistaria/rokko


text:
c Compute distance and azimuth on a sphere
c changes by Felix Waldhauser (fw)

	subroutine delaz2(alat, alon, blat, blon, del, dist, az)

	implicit none

	doubleprecision	alat, alon	! Coordinates of first point
	real		blat, blon	! Coordinates of second point
	real		del		! Sentral angle (degrees)
	real		dist		! Distance (km)
	real		az		! Azimuth from a to b (degrees)

C	Local variables:
	doubleprecision	acol, bcol
	doubleprecision	azr
	doubleprecision	blatr, blonr
	doubleprecision	colat
	doubleprecision	cosdel
	doubleprecision	delr
	doubleprecision	flat
	doubleprecision	geoa
	doubleprecision	geob
	doubleprecision	rad
	doubleprecision	radius
	doubleprecision alatr, alonr
	doubleprecision diflon
	doubleprecision pi2
	doubleprecision tana, tanb
cfw	real		xtop, xden	! see comments below
	doubleprecision top, den

c	Built-in functions: Declarations not needed
	doubleprecision dtan
	doubleprecision	datan
	doubleprecision	dsin
	doubleprecision	dcos
	doubleprecision	dacos

c	doubleprecision top,den

	data pi2/1.570796d0/
	data rad/1.745329d-02/
	data flat/.993231d0/

c-----convert to radians
	alatr=alat*rad
	alonr=alon*rad
	blatr=blat*rad
	blonr=blon*rad
c-----convert latitudes to geocentric colatitudes
	tana=flat*dtan(alatr)
	geoa=datan(tana)
	acol=pi2-geoa
	tanb=flat*dtan(blatr)
	geob=datan(tanb)
	bcol=pi2-geob
c-----calculate delta
	diflon=blonr-alonr
	cosdel=dsin(acol)*dsin(bcol)*dcos(diflon)+dcos(acol)*
     &	dcos(bcol)
	delr=dacos(cosdel)
c-----calculate azimuth from a to b

c*****	Note the use of single precision xtop and xden instead
c	of the double precision top and den in the original
c	program.
c*****	Note also the call to atan2 instead of datan2.
c	Both of these changes were made so that dyn.load
c	would work in Splus.  For some reason, the ld command
c	ld -r -d didn't find _d_atan2
c						WLE 10/16/91
cfw	xtop = dsin(diflon)
cfw	xden=(dsin(acol)/dtan(bcol))-dcos(diflon)*dcos(acol)
cfw	azr=atan2(xtop,xden)
c changes reversed: fw 02/05/17
	top = dsin(diflon)
	den=(dsin(acol)/dtan(bcol))-dcos(diflon)*dcos(acol)
	azr=datan2(top,den)

c----- convert to degrees
	del=delr/rad
	az=azr/rad
	if(az.lt.0.0) az=360.+az
c-----compute distance in kilometers
	colat=pi2-(alatr+blatr)/2.d0
cfw the equatorial radius of the Earth is 6378.137 km (IUGG value)
cfw the mean equatorial radius from Bott, 1982, is 6378.140 km
cfw	radius=6378.163d0*
cfw     & 	(1.d0+3.35278d-3*((1.d0/3.d0)-(dcos(colat)**2)))
	radius=6378.140*(1.0+3.37853d-3*((1/3)-((dcos(colat))**2)))
	dist=delr*radius
	return
c  ***** end of subroutine delaz *****
	end


text:
C
C
      SUBROUTINE ANEOSD(TQQ,RHOQQ,PQQ,EQQ,SQQ,CVQQ,DPDTQQ,DPDRQQ,
     1                  FKROQQ,CSQQ,KPA,MATQQ)
C
C***********************************************************************
C
C     ANEOS PACKAGE - RUNNING ENTRY POINT - DOUBLE PRECISION ENTRY
C     See ANEOSV for argument list
C
C************************************** 9/89 version slt ***************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CALL ANEOSV (1,TQQ,RHOQQ,MATQQ,
     &  PQQ,EQQ,SQQ,CVQQ,DPDTQQ,DPDRQQ,FKROQQ,CSQQ,KPA,R2PLQQ,R2PHQQ,
     &  ZBARQQ)
      RETURN
      END

max_stars_repo_name: chreinhardt/M-ANEOS

max_issues_repo_name: chreinhardt/M-ANEOS

max_forks_repo_name: chreinhardt/M-ANEOS