text:
module div_and_grad_module

  use multifab_module
  use ml_layout_module
  use define_bc_module
  use bc_module

  implicit none

  private

  public :: compute_grad, compute_div

contains

  subroutine compute_grad(mla,phi,gradp,dx, &
                          start_incomp,start_bccomp,start_outcomp,num_comp, &
                          the_bc_level,increment_bccomp_in)

    ! compute the face-centered gradient of a cell-centered field

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(in   ) :: phi(:)
    type(multifab) , intent(inout) :: gradp(:,:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    integer        , intent(in   ) :: start_incomp,start_bccomp,start_outcomp,num_comp
    type(bc_level) , intent(in   ) :: the_bc_level(:)
    logical,  intent(in), optional :: increment_bccomp_in

    ! local
    integer :: n,i,dm,nlevs,ng_p,ng_g,comp,bccomp,outcomp
    integer :: lo(mla%dim),hi(mla%dim)
    logical :: increment_bccomp

    real(kind=dp_t), pointer :: pp(:,:,:,:)
    real(kind=dp_t), pointer :: gpx(:,:,:,:)
    real(kind=dp_t), pointer :: gpy(:,:,:,:)
    real(kind=dp_t), pointer :: gpz(:,:,:,:)

    type(mfiter) :: mfi
    type(box) :: xnodalbox, ynodalbox, znodalbox
    integer :: xlo(mla%dim), xhi(mla%dim)
    integer :: ylo(mla%dim), yhi(mla%dim)
    integer :: zlo(mla%dim), zhi(mla%dim)

    type(bl_prof_timer),save :: bpt

    call build(bpt,"compute_grad")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_p = phi(1)%ng
    ng_g = gradp(1,1)%ng

    increment_bccomp = .true.
    if (present(increment_bccomp_in)) then
       increment_bccomp = increment_bccomp_in
    end if

    !$omp parallel private(mfi,n,i,xnodalbox,ynodalbox,znodalbox,xlo,ylo,zlo) &
    !$omp private(xhi,yhi,zhi,pp,gpx,gpy,gpz,lo,hi,comp,bccomp,outcomp)
    do n=1,nlevs
       call mfiter_build(mfi, phi(n), tiling=.true.)

     do while (more_tile(mfi))
        i = get_fab_index(mfi)

        xnodalbox = get_nodaltilebox(mfi,1)
        xlo = lwb(xnodalbox)
        xhi = upb(xnodalbox)
        ynodalbox = get_nodaltilebox(mfi,2)
        ylo = lwb(ynodalbox)
        yhi = upb(ynodalbox)
        znodalbox = get_nodaltilebox(mfi,3)
        zlo = lwb(znodalbox)
        zhi = upb(znodalbox)

 !      do i=1,nfabs(phi(n))
          pp  => dataptr(phi(n), i)
          gpx => dataptr(gradp(n,1), i)
          gpy => dataptr(gradp(n,2), i)
          lo = lwb(get_box(phi(n), i))
          hi = upb(get_box(phi(n), i))

          do comp=start_incomp,start_incomp+num_comp-1
             if (increment_bccomp) then
                bccomp = start_bccomp + (comp-start_incomp)
             else
                bccomp = start_bccomp
             end if
             outcomp = start_outcomp + (comp-start_incomp)
             select case (dm)
             case (2)
                call compute_grad_2d(pp(:,:,1,comp), ng_p, &
                                     gpx(:,:,1,outcomp), gpy(:,:,1,outcomp), ng_g, &
                                     lo, hi, dx(n,:), &
                                     the_bc_level(n)%adv_bc_level_array(i,:,:,bccomp),xlo,xhi,ylo,yhi)
             case (3)
                gpz => dataptr(gradp(n,3), i)
                call compute_grad_3d(pp(:,:,:,comp), ng_p, &
                                     gpx(:,:,:,outcomp), gpy(:,:,:,outcomp), gpz(:,:,:,outcomp), ng_g, &
                                     lo, hi, dx(n,:), &
                                     the_bc_level(n)%adv_bc_level_array(i,:,:,bccomp),xlo,xhi,ylo,yhi,zlo,zhi)
             end select
          end do
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  contains

    subroutine compute_grad_2d(phi,ng_p,gpx,gpy,ng_g,lo,hi,dx,bc,xlo,xhi,ylo,yhi)

      integer        , intent(in   ) :: ng_p,ng_g,lo(:),hi(:)
      integer        , intent(in   ) :: xlo(:),xhi(:),ylo(:),yhi(:)
      real(kind=dp_t), intent(in   ) :: phi(lo(1)-ng_p:,lo(2)-ng_p:)
      real(kind=dp_t), intent(inout) :: gpx(lo(1)-ng_g:,lo(2)-ng_g:)
      real(kind=dp_t), intent(inout) :: gpy(lo(1)-ng_g:,lo(2)-ng_g:)
      real(kind=dp_t), intent(in   ) :: dx(:)
      integer        , intent(in   ) :: bc(:,:)

      ! local
      integer :: i,j
      real(kind=dp_t) :: dxinv,twodxinv

      dxinv = 1.d0/dx(1)
      twodxinv = 2.d0*dxinv

      ! x-faces
      do j=xlo(2),xhi(2)
         do i=xlo(1),xhi(1)
            gpx(i,j) = ( phi(i,j)-phi(i-1,j) ) * dxinv
         end do
      end do

      ! alter stencil at boundary since ghost value represents value at boundary
      if (xlo(1) .eq. lo(1)) then
      if (bc(1,1) .eq. FOEXTRAP .or. bc(1,1) .eq. HOEXTRAP .or. bc(1,1) .eq. EXT_DIR) then
         i=lo(1)
         do j=lo(2),hi(2)
            gpx(i,j) = ( phi(i,j)-phi(i-1,j) ) * twodxinv
         end do
      end if
      end if
      if (xhi(1) .eq. hi(1)+1) then
      if (bc(1,2) .eq. FOEXTRAP .or. bc(1,2) .eq. HOEXTRAP .or. bc(1,2) .eq. EXT_DIR) then
         i=hi(1)+1
         do j=lo(2),hi(2)
            gpx(i,j) = ( phi(i,j)-phi(i-1,j) ) * twodxinv
         end do
      end if
      end if

      ! y-faces
      do j=ylo(2),yhi(2)
         do i=ylo(1),yhi(1)
            gpy(i,j) = ( phi(i,j)-phi(i,j-1) ) * dxinv
         end do
      end do

      ! alter stencil at boundary since ghost value represents value at boundary
      if (ylo(2) .eq. lo(2)) then
      if (bc(2,1) .eq. FOEXTRAP .or. bc(2,1) .eq. HOEXTRAP .or. bc(2,1) .eq. EXT_DIR) then
         j=lo(2)
         do i=lo(1),hi(1)
            gpy(i,j) = ( phi(i,j)-phi(i,j-1) ) * twodxinv
         end do
      end if
      end if

      ! alter stencil at boundary since ghost value represents value at boundary
      if (yhi(2) .eq. hi(2)+1) then
      if (bc(2,2) .eq. FOEXTRAP .or. bc(2,2) .eq. HOEXTRAP .or. bc(2,2) .eq. EXT_DIR) then
         j=hi(2)+1
         do i=lo(1),hi(1)
            gpy(i,j) = ( phi(i,j)-phi(i,j-1) ) * twodxinv
         end do
      end if
      end if

    end subroutine compute_grad_2d

    subroutine compute_grad_3d(phi,ng_p,gpx,gpy,gpz,ng_g,lo,hi,dx,bc,xlo,xhi,ylo,yhi,zlo,zhi)

      integer        , intent(in   ) :: ng_p,ng_g,lo(:),hi(:)
      integer        , intent(in   ) :: xlo(:),xhi(:),ylo(:),yhi(:),zlo(:),zhi(:)
      real(kind=dp_t), intent(in   ) :: phi(lo(1)-ng_p:,lo(2)-ng_p:,lo(3)-ng_p:)
      real(kind=dp_t), intent(inout) :: gpx(lo(1)-ng_g:,lo(2)-ng_g:,lo(3)-ng_g:)
      real(kind=dp_t), intent(inout) :: gpy(lo(1)-ng_g:,lo(2)-ng_g:,lo(3)-ng_g:)
      real(kind=dp_t), intent(inout) :: gpz(lo(1)-ng_g:,lo(2)-ng_g:,lo(3)-ng_g:)
      real(kind=dp_t), intent(in   ) :: dx(:)
      integer        , intent(in   ) :: bc(:,:)

      ! local
      integer :: i,j,k
      real(kind=dp_t) :: dxinv,twodxinv

      dxinv = 1.d0/dx(1)
      twodxinv = 2.d0*dxinv

      ! x-faces
      do k=xlo(3),xhi(3)
         do j=xlo(2),xhi(2)
            do i=xlo(1),xhi(1)
               gpx(i,j,k) = ( phi(i,j,k)-phi(i-1,j,k) ) * dxinv
            end do
         end do
      end do

      ! alter stencil at boundary since ghost value represents value at boundary
      if (xlo(1) .eq. lo(1)) then
      if (bc(1,1) .eq. FOEXTRAP .or. bc(1,1) .eq. HOEXTRAP .or. bc(1,1) .eq. EXT_DIR) then
         i=xlo(1)
         do k=xlo(3),xhi(3)
            do j=xlo(2),xhi(2)
               gpx(i,j,k) = ( phi(i,j,k)-phi(i-1,j,k) ) * twodxinv
            end do
         end do
      end if
      end if
      if (xhi(1) .eq. hi(1)+1) then
      if (bc(1,2) .eq. FOEXTRAP .or. bc(1,2) .eq. HOEXTRAP .or. bc(1,2) .eq. EXT_DIR) then
         i=xhi(1)
         do k=xlo(3),xhi(3)
            do j=xlo(2),xhi(2)
               gpx(i,j,k) = ( phi(i,j,k)-phi(i-1,j,k) ) * twodxinv
            end do
         end do
      end if
      end if

      ! y-faces
      do k=ylo(3),yhi(3)
         do j=ylo(2),yhi(2)
            do i=ylo(1),yhi(1)
               gpy(i,j,k) = ( phi(i,j,k)-phi(i,j-1,k) ) * dxinv
            end do
         end do
      end do

      ! alter stencil at boundary since ghost value represents value at boundary
      if (ylo(2) .eq. lo(2)) then
      if (bc(2,1) .eq. FOEXTRAP .or. bc(2,1) .eq. HOEXTRAP .or. bc(2,1) .eq. EXT_DIR) then
         j=ylo(2)
         do k=ylo(3),yhi(3)
            do i=ylo(1),yhi(1)
               gpy(i,j,k) = ( phi(i,j,k)-phi(i,j-1,k) ) * twodxinv
            end do
         end do
      end if
      end if
      if (yhi(2) .eq. hi(2)+1) then
      if (bc(2,2) .eq. FOEXTRAP .or. bc(2,2) .eq. HOEXTRAP .or. bc(2,2) .eq. EXT_DIR) then
         j=yhi(2)
         do k=ylo(3),yhi(3)
            do i=ylo(1),yhi(1)
               gpy(i,j,k) = ( phi(i,j,k)-phi(i,j-1,k) ) * twodxinv
            end do
         end do
      end if
      end if

      ! z-faces
      do k=zlo(3),zhi(3)
         do j=zlo(2),zhi(2)
            do i=zlo(1),zhi(1)
               gpz(i,j,k) = ( phi(i,j,k)-phi(i,j,k-1) ) * dxinv
            end do
         end do
      end do

      ! alter stencil at boundary since ghost value represents value at boundary
      if (zlo(3) .eq. lo(3)) then
      if (bc(3,1) .eq. FOEXTRAP .or. bc(3,1) .eq. HOEXTRAP .or. bc(3,1) .eq. EXT_DIR) then
         k=zlo(3)
         do j=zlo(2),zhi(2)
            do i=zlo(1),zhi(1)
               gpz(i,j,k) = ( phi(i,j,k)-phi(i,j,k-1) ) * twodxinv
            end do
         end do
      end if
      end if

      ! alter stencil at boundary since ghost value represents value at boundary
      if (zhi(3) .eq. hi(3)+1) then
      if (bc(3,2) .eq. FOEXTRAP .or. bc(3,2) .eq. HOEXTRAP .or. bc(3,2) .eq. EXT_DIR) then
         k=zhi(3)
         do j=zlo(2),zhi(2)
            do i=zlo(1),zhi(1)
               gpz(i,j,k) = ( phi(i,j,k)-phi(i,j,k-1) ) * twodxinv
            end do
         end do
      end if
      end if

    end subroutine compute_grad_3d

  end subroutine compute_grad

  subroutine compute_div(mla,phi_fc,div,dx,start_incomp,start_outcomp,num_comp,increment_in)

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(in   ) :: phi_fc(:,:)
    type(multifab) , intent(inout) :: div(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:)
    integer        , intent(in   ) :: start_incomp, start_outcomp, num_comp
    logical,  intent(in), optional :: increment_in

    real(kind=dp_t), pointer :: pxp(:,:,:,:)
    real(kind=dp_t), pointer :: pyp(:,:,:,:)
    real(kind=dp_t), pointer :: pzp(:,:,:,:)
    real(kind=dp_t), pointer :: dp(:,:,:,:)
    integer :: i,n,nlevs,dm,ng_p,ng_d,comp,outcomp
    integer :: lo(mla%dim),hi(mla%dim)
    logical :: increment

    type(mfiter) :: mfi
    type(box) :: tilebox
    integer :: tlo(mla%dim), thi(mla%dim)

    type(bl_prof_timer),save :: bpt

    call build(bpt,"compute_div")

    ! do we increment or overwrite div?
    increment = .false.
    if (present(increment_in)) increment = increment_in

    dm = mla%dim
    nlevs = mla%nlevel

    ng_p = phi_fc(1,1)%ng
    ng_d = div(1)%ng

    !$omp parallel private(mfi,n,i,tilebox,tlo,thi,pxp,pyp,pzp,dp,lo,hi,comp,outcomp)
    do n = 1,nlevs
       call mfiter_build(mfi, div(n), tiling=.true.)

     do while (more_tile(mfi))
       i = get_fab_index(mfi)

       tilebox = get_tilebox(mfi)
       tlo = lwb(tilebox)
       thi = upb(tilebox)

!       do i = 1, nfabs(div(n))
          pxp => dataptr(phi_fc(n,1), i)
          pyp => dataptr(phi_fc(n,2), i)
          dp => dataptr(div(n), i)
          lo =  lwb(get_box(div(n), i))
          hi =  upb(get_box(div(n), i))
          do comp=start_incomp,start_incomp+num_comp-1
             outcomp = start_outcomp + (comp-start_incomp)
             select case (dm)
             case (2)
                call compute_div_2d(pxp(:,:,1,comp), pyp(:,:,1,comp), ng_p, &
                                    dp(:,:,1,outcomp), ng_d, dx(n,:),lo, hi, increment,tlo,thi)
             case (3)
                pzp => dataptr(phi_fc(n,3), i)
                call compute_div_3d(pxp(:,:,:,comp), pyp(:,:,:,comp), pzp(:,:,:,comp), ng_p, &
                                    dp(:,:,:,outcomp), ng_d, dx(n,:),lo, hi, increment,tlo,thi)
             end select
          end do
       end do
    end do
    !$omp end parallel

    call destroy(bpt)

  contains

    subroutine compute_div_2d(phix,phiy,ng_p,div,ng_d,dx,lo,hi,increment,tlo,thi)

      integer        , intent(in   ) :: lo(:),hi(:),ng_p,ng_d,tlo(:),thi(:)
      real(kind=dp_t), intent(in   ) :: phix(lo(1)-ng_p:,lo(2)-ng_p:)
      real(kind=dp_t), intent(in   ) :: phiy(lo(1)-ng_p:,lo(2)-ng_p:)
      real(kind=dp_t), intent(inout) ::  div(lo(1)-ng_d:,lo(2)-ng_d:)
      real(kind=dp_t), intent(in   ) ::   dx(:)
      logical        , intent(in   ) :: increment

      integer :: i,j
      real(kind=dp_t) :: dxinv

      dxinv = 1.d0/dx(1)

      if (increment) then

         do j = tlo(2),thi(2)
         do i = tlo(1),thi(1)
            div(i,j) = div(i,j) + &
                 (phix(i+1,j) - phix(i,j)) * dxinv  + &
                 (phiy(i,j+1) - phiy(i,j)) * dxinv
         end do
         end do

      else

         do j = tlo(2),thi(2)
         do i = tlo(1),thi(1)
            div(i,j) = &
                 (phix(i+1,j) - phix(i,j)) * dxinv + &
                 (phiy(i,j+1) - phiy(i,j)) * dxinv
         end do
         end do

      end if

    end subroutine compute_div_2d

    subroutine compute_div_3d(phix,phiy,phiz,ng_p,div,ng_d,dx,lo,hi,increment,tlo,thi)

      integer        , intent(in   ) :: lo(:),hi(:),ng_p,ng_d,tlo(:),thi(:)
      real(kind=dp_t), intent(in   ) :: phix(lo(1)-ng_p:,lo(2)-ng_p:,lo(3)-ng_p:)
      real(kind=dp_t), intent(in   ) :: phiy(lo(1)-ng_p:,lo(2)-ng_p:,lo(3)-ng_p:)
      real(kind=dp_t), intent(in   ) :: phiz(lo(1)-ng_p:,lo(2)-ng_p:,lo(3)-ng_p:)
      real(kind=dp_t), intent(inout) ::  div(lo(1)-ng_d:,lo(2)-ng_d:,lo(3)-ng_d:)
      real(kind=dp_t), intent(in   ) :: dx(:)
      logical        , intent(in   ) :: increment

      integer :: i,j,k
      real(kind=dp_t) :: dxinv

      dxinv = 1.d0/dx(1)

      if (increment) then
         do k = tlo(3),thi(3)
         do j = tlo(2),thi(2)
         do i = tlo(1),thi(1)
            div(i,j,k) = div(i,j,k) + &
                 (phix(i+1,j,k) - phix(i,j,k)) * dxinv + &
                 (phiy(i,j+1,k) - phiy(i,j,k)) * dxinv + &
                 (phiz(i,j,k+1) - phiz(i,j,k)) * dxinv
         end do
         end do
         end do
      else
         do k = tlo(3),thi(3)
         do j = tlo(2),thi(2)
         do i = tlo(1),thi(1)
            div(i,j,k) = &
                 (phix(i+1,j,k) - phix(i,j,k)) * dxinv + &
                 (phiy(i,j+1,k) - phiy(i,j,k)) * dxinv + &
                 (phiz(i,j,k+1) - phiz(i,j,k)) * dxinv
         end do
         end do
         end do
      end if

    end subroutine compute_div_3d

  end subroutine compute_div

end module div_and_grad_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
!
!	Find the Gilula-Haberman correction.
!	eacov is the estimated asymptotic covariance matrix under the model.
!	inform is the estimated information obtained by the Louis approach.
!	loglik is the log likelihood.
!	totalitems is the total weighted number of presented items.
!
!
real(kind=8) function gh(eacov,inform,loglik,totalitems)
implicit none
real(kind=8),intent(in)::eacov(:,:),inform(:,:),loglik,totalitems
!	col counts columns.
!	obs is an observation counter.
!	row counts rows.
integer::col,obs,row
gh=-loglik
do row=1,size(eacov,1)
	do col=1,row
		if(col<row)then
			gh=gh+2.0_8*eacov(row,col)*inform(col,row)
		else
			gh=gh+eacov(row,row)*inform(row,row)
		end if
	end do
end do
gh=gh/totalitems
return
end function gh

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
Describe Users/SarahS here.

20100304 11:43:22 nbsp I just noticed the wave, sorry it took so long to get back to you... Users/JoePomidor

20100521 17:18:10 nbsp Hello, I am unsure if you know this or not, but everyone is able to see whenever a user makes any edit. It seems like you wanted http://daviswiki.org/Target?actiondiff&version2192&version1191 this edit to be anonymous (or maybe you let someone else post from your account), but just in case, I thought I should let you know. And welcome to the Wiki! Users/hankim

20110117 18:24:59 nbsp You ought to write a letter to A) the apartment office; B) the management company; and C) the owner, if you can get the contact info.  A basic demand letter laying out what youre owed and the law, and explaining that youll contact BBB and file suit in small claims court if you arent paid in a timely manner.

City of Davis also has a conflict resolution resource that can help you get your deposit back without having to resort to small claims court (hopefully).  http://cityofdavis.org/cs/mediation/index.cfm Users/TomGarberson

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module mod_xdmf

  use iso_fortran_env, only: real32
  use mod_functional, only: operator(.reverse.)
  use mod_kinds, only: ik, rk
  use mod_io, only: contour_writer_t
  ! use mod_grid, only : grid_t, grid_2d_t, grid_3d_t
  use hdf5_interface, only: hdf5_file

  implicit none

  private
  public :: xdmf_contour_t

  type, extends(contour_writer_t) :: xdmf_contour_t
    private
    type(hdf5_file) :: h5_file

    character(:), allocatable :: xml_filename
    character(:), allocatable :: h5_filename

    character(:), allocatable :: dimension_str
    character(:), allocatable :: topology_type
    character(:), allocatable :: time_str
    integer(ik) :: xml_file_id
  contains
    procedure, public :: initialize => init_xdmf
    procedure, public :: finalize => final_xdmf
    procedure, public :: add_header => write_xdmf_header
    procedure, private :: add_int_1d => xdmf_add_int_1d
    procedure, private :: add_int_2d => xdmf_add_int_2d
    procedure, private :: add_int_3d => xdmf_add_int_3d
    procedure, private :: add_real64_1d => xdmf_add_real64_1d
    procedure, private :: add_real64_2d => xdmf_add_real64_2d
    procedure, private :: add_real64_3d => xdmf_add_real64_3d
    procedure, private :: init_xml
    procedure, private :: final_xml
    procedure, private :: write_xml_header
    procedure, private :: write_xml_grid
    procedure, private :: add_xml_variable_1d
    procedure, private :: add_xml_variable_2d
    procedure, private :: add_xml_variable_3d
  end type xdmf_contour_t

contains

  subroutine init_xml(self)
    class(xdmf_contour_t), intent(inout) :: self
    character(50) :: char_buffer
    integer(ik) :: ios

    open(newunit=self%xml_file_id, file=self%xml_filename, &
         iostat=ios, status="replace", action="write")

    if(ios /= 0) then
      print *, "Error opening xdmf file: '"//self%xml_filename//"'"
      error stop "File I/O Problem"
    end if

    select case(self%ndim)
    case(2)
      self%topology_type = '2DRectMesh'

      char_buffer = ''
      write(char_buffer, '(2(i0, 1x))') self%shape
      self%dimension_str = trim(char_buffer)
    case(3)
      self%topology_type = '3DRectMesh'

      char_buffer = ''
      write(char_buffer, '(3(i0, 1x))') self%shape
      self%dimension_str = trim(char_buffer)
    case default
      error stop 'Error: mod_xdmf only supports 2D and 3D grids'
    end select

    char_buffer = ''
    write(char_buffer, '(f0.4)') self%time
    self%time_str = trim(char_buffer)
  end subroutine init_xml

  subroutine final_xml(self)
    class(xdmf_contour_t), intent(inout) :: self

    write(self%xml_file_id, '(2(2x),a)') '</Grid>'
    write(self%xml_file_id, '(1(2x),a)') '</Domain>'
    write(self%xml_file_id, '(a)') '</Xdmf>'
    close(self%xml_file_id)
  end subroutine final_xml

  subroutine write_xml_header(self)
    class(xdmf_contour_t), intent(inout) :: self

    write(self%xml_file_id, '(a)') '<?xml version="1.0"?>'
    write(self%xml_file_id, '(a)') '<Xdmf version="2.2">'
    write(self%xml_file_id, '(1(2x),a)') '<Domain>'
    write(self%xml_file_id, '(2(2x),a)') '<Grid GridType="Uniform" Name="grid">'
    write(self%xml_file_id, '(3(2x),a)') '<Time Value="'//self%time_str//'"/>'
    write(self%xml_file_id, '(3(2x),a,3(i0,1x),a)') '<Topology NumberOfElements="', &
      .reverse.self%shape, '" TopologyType="'//self%topology_type//'"/>'
  end subroutine write_xml_header

  subroutine write_xml_grid(self)
    class(xdmf_contour_t), intent(inout) :: self

    if(self%ndim == 2) then
      write(self%xml_file_id, '(3(2x),a)') '<Geometry GeometryType="VXVYVZ">'

      write(self%xml_file_id, '(4(2x),a,i0,a)') '<DataItem Dimensions="', self%shape(1), '"'// &
        ' Format="HDF" NumberType="Float"'// &
        ' Precision="4">'//self%h5_filename//':/x</DataItem>'

      write(self%xml_file_id, '(4(2x),a,i0,a)') '<DataItem Dimensions="', self%shape(2), '"'// &
        ' Format="HDF" NumberType="Float"'// &
        ' Precision="4">'//self%h5_filename//':/y</DataItem>'

      write(self%xml_file_id, '(3(2x),a)') '</Geometry>'

    else if(self%ndim == 3) then
      write(self%xml_file_id, '(3(2x),a)') '<Geometry GeometryType="VXVYVZ">'

      write(self%xml_file_id, '(4(2x),a,i0,a)') '<DataItem Dimensions="', self%shape(1), '"'// &
        ' Format="HDF" NumberType="Float"'// &
        ' Precision="4">'//self%h5_filename//':/x</DataItem>'

      write(self%xml_file_id, '(4(2x),a,i0,a)') '<DataItem Dimensions="', self%shape(2), '"'// &
        ' Format="HDF" NumberType="Float"'// &
        ' Precision="4">'//self%h5_filename//':/y</DataItem>'

      write(self%xml_file_id, '(4(2x),a,i0,a)') '<DataItem Dimensions="', self%shape(3), '"'// &
        ' Format="HDF" NumberType="Float"'// &
        ' Precision="4">'//self%h5_filename//':/z</DataItem>'

      write(self%xml_file_id, '(3(2x),a)') '</Geometry>'
    end if

  end subroutine write_xml_grid

  subroutine add_xml_variable_1d(self, name)

    class(xdmf_contour_t), intent(inout) :: self
    character(*), intent(in) :: name

    write(self%xml_file_id, '(3(2x),a)') '<Attribute AttributeType="Scalar" Center="Node" Name="'//name//'">'
    write(self%xml_file_id, '(4(2x),a,1(i0,1x),a)') '<DataItem Dimensions="', 1, '" Format="HDF" '// &
      'NumberType="Float" Precision="4">'//self%h5_filename//':/'//name//'</DataItem>'
    write(self%xml_file_id, '(3(2x),a)') '</Attribute>'

  end subroutine add_xml_variable_1d

  subroutine add_xml_variable_2d(self, name, data_shape)

    class(xdmf_contour_t), intent(inout) :: self
    character(*), intent(in) :: name
    integer(ik), dimension(2), intent(in) :: data_shape

    write(self%xml_file_id, '(3(2x),a)') '<Attribute AttributeType="Scalar" Center="Node" Name="'//name//'">'
    write(self%xml_file_id, '(4(2x),a,2(i0,1x),a)') '<DataItem Dimensions="', .reverse.data_shape, '" Format="HDF" '// &
      'NumberType="Float" Precision="4">'//self%h5_filename//':/'//name//'</DataItem>'
    write(self%xml_file_id, '(3(2x),a)') '</Attribute>'

  end subroutine add_xml_variable_2d

  subroutine add_xml_variable_3d(self, name, data_shape)

    class(xdmf_contour_t), intent(inout) :: self
    character(*), intent(in) :: name
    integer(ik), dimension(3), intent(in) :: data_shape

    write(self%xml_file_id, '(3(2x),a)') '<Attribute AttributeType="Scalar" Center="Node" Name="'//name//'">'
    write(self%xml_file_id, '(4(2x),a,3(i0,1x),a)') '<DataItem Dimensions="', .reverse.data_shape, '" Format="HDF" '// &
      'NumberType="Float" Precision="4">'//self%h5_filename//':/'//name//'</DataItem>'
    write(self%xml_file_id, '(3(2x),a)') '</Attribute>'

  end subroutine add_xml_variable_3d

  subroutine init_xdmf(self, timestep, time, nx_ny_nz, ndim)
    !! Initialize the xdmf contour writer

    class(xdmf_contour_t), intent(inout) :: self
    real(rk), intent(in) :: time
    integer(ik), intent(in) :: timestep
    integer(ik), intent(in) :: nx_ny_nz(3)
    integer(ik), intent(in) :: ndim

    ! Locals
    character(50) :: char_buff

    self%format = 'xdmf'
    self%time = time
    self%ndim = ndim
    self%timestep = timestep

    write(char_buff, '(a,i0.7)') self%filename_prefix, self%timestep
    self%filename = trim(char_buff)
    self%h5_filename = self%filename//'.h5'
    self%xml_filename = self%filename//'.xdmf'

    allocate(self%shape(ndim))

    if(this_image() /= 1) then
      error stop "Error, init_xdmf in "//__FILE__//" should only be called by image 1"
    endif

    call self%init_xml()

    call self%h5_file%initialize(filename=self%h5_filename, &
                                 status='new', action='w', comp_lvl=6)

  end subroutine init_xdmf

  subroutine final_xdmf(self)
    class(xdmf_contour_t), intent(inout) :: self

    call self%final_xml()
    call self%h5_file%finalize()
  end subroutine final_xdmf

  subroutine write_xdmf_header(self)
    class(xdmf_contour_t), intent(inout) :: self
    call self%write_xml_header()
    call self%write_xml_grid()
  end subroutine

  subroutine xdmf_add_int_1d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    integer(ik), dimension(:), intent(in):: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_1d(name)
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_int_1d

  subroutine xdmf_add_int_2d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    integer(ik), dimension(:, :), intent(in):: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_2d(name, shape(data))
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_int_2d

  subroutine xdmf_add_int_3d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    integer(ik), dimension(:, :, :), intent(in):: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_3d(name, shape(data))
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_int_3d

  subroutine xdmf_add_real64_1d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    real(rk), dimension(:), intent(in) :: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_1d(name)
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_real64_1d

  subroutine xdmf_add_real64_2d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    real(rk), dimension(:, :), intent(in) :: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_2d(name, shape(data))
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_real64_2d

  subroutine xdmf_add_real64_3d(self, name, data, units, description)
    class(xdmf_contour_t), intent(inout) :: self
    real(rk), dimension(:, :, :), intent(in) :: data
    character(*), intent(in) :: name
    character(*), intent(in) :: units
    character(*), intent(in) :: description

    character(:), allocatable :: dset_name
    dset_name = '/'//trim(name)
    call self%add_xml_variable_3d(name, shape(data))
    call self%h5_file%add(dset_name, data)
    call self%h5_file%writeattr(dset_name, 'description', trim(description))
    call self%h5_file%writeattr(dset_name, 'units', trim(units))
  end subroutine xdmf_add_real64_3d

end module mod_xdmf

! <?xml version="1.0" ?>
! <Xdmf version="2.2">
!   <Domain>
!     <Grid GridType="Uniform" Name="grid">
!       <Time Value="2.3 nanosecond"/>
!       <Topology NumberOfElements="266 601 31" TopologyType="3DSMesh"/>
!       <Geometry GeometryType="X_Y_Z">
!         <DataItem Dimensions="266 601 31" Format="HDF" NumberType="Float" Precision="4">profile_126_2.300ns.h5:/x</DataItem>
!         <DataItem Dimensions="266 601 31" Format="HDF" NumberType="Float" Precision="4">profile_126_2.300ns.h5:/y</DataItem>
!         <DataItem Dimensions="266 601 31" Format="HDF" NumberType="Float" Precision="4">profile_126_2.300ns.h5:/z</DataItem>
!       </Geometry>
!       <Attribute AttributeType="Scalar" Center="Cell" Name="mass_density">
! <DataItem Dimensions="265 600 30" Format="HDF" NumberType="Float" Precision="4">profile_126_2.300ns.h5:/mass_density</DataItem>
!       </Attribute>
!       <Attribute AttributeType="Scalar" Center="Cell" Name="ion_temperature">
!         <DataItem Dimensions="265 600 30" Format="HDF" NumberType="Float" Precision="4">profile_126_2.300ns.h5:/ion_temperature</DataItem>
!       </Attribute>
!     </Grid>
!   </Domain>
! </Xdmf>

max_stars_repo_name: smillerc/lbm

max_issues_repo_name: smillerc/lbm

max_forks_repo_name: smillerc/lbm


text:
The Varsity Facebook Team is for any of those select, obsessed, perverted people who have an obscene number of friends on Facebook. Membership in the team requires a minimum of 500 UC Davis friends (note the change from total friends). If you or anyone you know is in this select group, add their names to the list. This list was initially compiled from the http://ucdavis.facebook.com/group.php?gid2200601328 500 Club facebook group on 4/19/06.

 The Team
    Chris Dietrich
    Users/EddieKim
    Samantha Wang
    Scott Sablan
    Lula AhmedFalol
    Daniel Thorne
    Robert Kenzo Onaga
    Nafeh Malik
    Simone Davis
    Pat Weathers
    Rangashree Varadarajan
    Katie Webber
    Tre Clifton
    Danny Bishop
    Caroline Jordan
    Tolu Wusu
    Daryl Suyat
    Annie Chung
    Nina Beckwith
    X Steve OstrowskiX (banned from facebook)
    Users/MaxMikalonis Max Mikalonis
    Joe Chatham
    Users/MilesThomas Miles Thomas
    Users/RomanRivilis Roman Rivilis

Alumni

    Eduardo De La Torre
    Yoshi Akaha
    Andrew Bianchi
    Scott Judson
    Madisyn Coniglio
    Alexandra Frick
    George Ajlouny
    Ami Vora
    Emily Tung
    Rebecca Lovell Becca Lovell
    Rebecca Schwartz
    Users/StevenLee Steven D. Lee
    Ivan Carrillo
    Molly Fluet
    Paul Harms
    Kalen Gallagher
    Shaneil Sharma
    Spencer Higgins
    Chris Perry
    James Schwab
    Albert Parnell
    Eric Friedman
    Dwaimy RosasRomero
    Darnell Holloway
    Andrew Peake
    Don Ho
    Rob Roy
    Chris Herold
    Users/WencheMolenaar Wenche Molenaar
    Parisa Manteghi

The UC Davis Students student who has the most Facebook friends was: James Schwab (1,702 as of 8/15/07).
The non UC Davis Students student who has (tied for) the least Facebook friends is: Users/Jabberwokky (0 as of 10/12/07).

20060419 21:00:16 nbsp So I just looked up Katie Webber on facebook because I was like who is this uber popular person Ive never heard of, and unless I am totally tripping, it says she has 582 friends.  Granted thats a lot, but thats half as many as shes listed here as having.  It doesnt include from other schools and I certainly didnt bother to tally up those numbers but I doubt that adds more than a hundred if that .... so whats the dealio? Users/KrisFricke

20060419 21:12:22 nbsp If you add up the not UCD friends, then that number alone is ~432 (I just did that really fast, but it shouldnt be too far off), and that places her with 1014 Facebook Friends. Users/MaxMikalonis

20060419 21:15:12 nbsp I started this page looking at peoples total number of friends. i figured that while its a little harder to figure out, its probably a more accurate way to measure how much of a facebook slut someone is (myself included). if you look at katies profile you can see that she has a ton of friends from other campuses. if you click on the link part of Katie has 582 UC Davis friends and then refresh the list to show all friends, it gives a handy little count of her friends from everywhere. when i first saw her name, i thought the same thing. how does this girl know everyone else but not me? i guess old guys like us are of no interest to upstart student focus freshmen girls like katie. sigh... Users/AdamGerber

20060419 21:15:17 nbsp Are you for reals?  theres only seven schools with more then ten.. shes really nickle and diming it up Users/KrisFricke

20060419 21:22:05 nbsp I doubt she specifically was trying to be on this page or anything, but yeah, Some people I just dont understand. Some of them, they have to have met the person to add them. Others are looser, and have their whole High School as friends. Users/MaxMikalonis

20060419 21:22:24 nbsp The question I supose is... Do she and Schwab know each other? :) Users/JosephBleckman

20060419 21:24:57 nbsp The answer is yes! She ran for Senate on Student Focus...... Users/KyleFlick

20060419 21:25:18 nbsp http://ucdavis.facebook.com/friends.php?search_querywebber&category_id8&id3208416 they do know each other! Users/AdamGerber

20060419 21:29:28 nbsp But around 1/2 way through the summer, James started adding every new Freshman on Fb, so in that respect, it doesnt count, but in terms of Kyle Flicks comment, they might. Users/MaxMikalonis

20060419 21:32:30 nbsp I think I remember that. I was kinda weirded out by it. Users/JosephBleckman

20060419 21:48:23 nbsp I usually dont add or approve random peeps... hence my normal friend count. Users/TusharRawat

20060419 21:51:47 nbsp yeah me, I get a friends request, and Im like yeah I know them, but do I know them that WELL? ...  pretty much the opposite of the friendwhoring types Users/KrisFricke
   20060419 22:03:21 nbsp kris, does that make you a friend prude? Users/AdamGerber

20060419 22:10:29 nbsp I only have 475 friends. I feel so very unpopular. Users/RobRoy
  The party king of davis party king, and 2006 City Council Election next mayor of davis has only 475 facebook friends? something about this lonley monarch seems strangely sad and a little Dickensian. i guess your just gonna have to settle with having more real friends while i have to promote my large number of imaginary friends.  Adam Gerber Adam
  http://ucdavis.facebook.com/friends.php?id3201803&category_id0 Your actually at 490, btw. so very close.  Adam Gerber Adam

20060420 00:51:47 nbsp I actually know two of these folks and have had multiple classes with them and yet are not friends with them on facebook.  I am awesome, or incredibly lame because of this. Users/BenWebber

20060420 11:03:13 nbsp I think these should be reverted to just UCD friends.  Hell, I could go back and add everyone from my highschool I ever talked to and double my number too. It doesnt really matter though.  I love James but when he runs into half of that list he doesnt recognize them I know its the same for the rest of the people listed here as well.  I tend to delete people when their profiles show up in the recently updated and I cant remember who they were.  However, Ive noticed that other schools, particularly Harvard and the Ivy Leagues, the numbers average much higher and dont seem to subscribe to my same facebook philosophy. Users/ThomasLloyd

20060422 21:44:52 nbsp I agree about that. Being overly popular in high school doesnt equal Davis relevancy. Users/JosephBleckman

20060423 01:21:47 nbsp Additionally, it is much easier to find people w/500+ Facebook UCD friends than 500+Facebook friends total. However, 500 is quite exclusive, but I guess that is the point. Users/MaxMikalonis

20060426 13:47:31 nbsp So, is there anyone who still wants to keep it at the total friend count? Users/JosephBleckman

20060426 16:42:20 nbsp i do. i think that having friends at schools other than ucdavis is still a good sign of how much of a facebook slut someone is. i guess we could switch it over, but i still feel like total friends is a better measurement. Users/AdamGerber

20060426 17:51:15 nbsp Friends at Davis is a really obvious number (its listed directly), while friends at other school needs to be added up by someone with a morbid interest in the subject; as long as we have data for the latter, we may as well list the former as well, and ohmygosh list it by BOTH systems... Users/KrisFricke

20060430 11:42:06 nbsp okay, so i changed it over for all of you. now its just uc davis friends. its alot of the same people as before but just a little less. sadly, im one of the few people who didnt make the cut :( Users/AdamGerber

20060822 08:13:47 nbsp Congrats to Rob Roy for making the list! Users/AndrewBianchi

20070224 14:31:31 nbsp I have 66 friends on facebook. About two in real life. Users/AliceChoe

20070312 18:37:45 nbsp This team isnt so 1337 anymore, the bar should be raised. Users/GregWebb

20070322 01:25:41 nbsp I think that people should only be counted if they actually know and/or have met their facebook friends.  There are a few people who added me as a facebook friend about 2 weeks after I got here Freshman year.  I do realize there is absolutely no way to enforce this.. it was just a passing thought. Users/KirstenOnell

20071011 22:21:12 nbsp As a member of this stupid team I really think it is ridiculous to tally up these popularity statistics. Users/AliceChoe has the best comment on this page. I have 857 facebook friends  how obscene. I dont know 90% of the people. There are no perks to being in this team. Why are good people wasting their time documenting a pointless popularity contest. Users/RobRoy

20071011 22:23:35 nbsp facebookers of the world unite! you have nothing to lose but your friends! delete your account! Users/JessicaRockwell

20071202 23:48:54 nbsp with all those names, I wouldnt call us a team at all P Users/EricFriedman

20080214 10:58:25 nbsp Ive been at UCD for over five years and Ive yet to get 200 friends.  I think the bar is high enough. Users/BrentLaabs


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine vortex(node1,node2,nodem,nelem,lakon,kon,ipkon,
     &     nactdog,identity,ielprop,prop,iflag,v,xflow,f,
     &     nodef,idirf,df,cp,R,numf,set,mi,ttime,time,iaxial,iplausi)
!
!     vortex element
!
!     author: Yannick Muller
!
      implicit none
!
      logical identity
      character*8 lakon(*)
      character*81 set(*)
!
      integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
     &     ielprop(*),nodef(*),idirf(*),index,iflag,iaxial,
     &     inv,ipkon(*),kon(*),t_chang,nelemswirl,mi(*),iplausi
!
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,cp,
     &     p1,p2,T1,T2,km1,pi,ttime,time,r2d,r1d,eta,U1,
     &     c1u,c2u,cinput,r1,r2,omega,K1,ciu,expon,
     &     Ui,Kr,cte1,cte2,qred_crit,A,xflow_oil
!
!
!
      if(iflag.eq.0) then
         identity=.true.
!
         if(nactdog(2,node1).ne.0)then
            identity=.false.
         elseif(nactdog(2,node2).ne.0)then
            identity=.false.
         elseif(nactdog(1,nodem).ne.0)then
            identity=.false.
         endif
!
      elseif(iflag.eq.1)then
         if(v(1,nodem).ne.0.d0) then
            xflow=v(1,nodem)
            return
         endif
!
         xflow=0.d0
!
      elseif(iflag.eq.2)then
!
         numf=4
         index=ielprop(nelem)
         kappa=(cp/(cp-R))
         km1=kappa-1
         pi=4.d0*datan(1.d0)
!
!     radius downstream
         r2d=prop(index+1)
!
!     radius upstream
         r1d=prop(index+2)
!
!     pressure correction factor
         eta=prop(index+3)
!
         p1=v(2,node1)
         p2=v(2,node2)
!
         xflow=v(1,nodem)*iaxial
!
         if(xflow.gt.0.d0) then
            inv=1.d0
            p1=v(2,node1)
            p2=v(2,node2)
            T1=v(0,node1)
            T2=v(0,node2)
            R1=r1d
            R2=r2d
!
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!
         elseif(xflow.lt.0.d0) then
            inv=-1.d0
            R1=r2d
            R2=r1d
            p1=v(2,node2)
            p2=v(2,node1)
            T1=v(0,node2)
            T2=v(0,node1)
            xflow=-v(1,nodem)*iaxial
!
            nodef(1)=node2
            nodef(2)=node2
            nodef(3)=nodem
            nodef(4)=node1
!
         endif
!
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=2
!
         kappa=(cp/(cp-R))
!
!     FREE VORTEX
!
         if(lakon(nelem)(4:5).eq.'FR')then
!
!     rotation induced loss (correction factor)
           K1= prop(index+4)
!
!     tangential velocity of the disk at vortex entry
             U1=prop(index+5)
!
!     number of the element generating the upstream swirl
             nelemswirl=nint(prop(index+6))
!
!     rotation speed (revolution per minutes)
             omega=prop(index+7)
!
!     Temperature change
            t_chang=prop(index+8)
!
            if(omega.gt.0.d0) then
!
!     rotation speed is given if the swirl comes from a rotating part
!     typically the blade of a coverplate
!

!     C_u is given by radius r1d (see definition of the flow direction)
!     C_u related to radius r2d is a function of r1d
!
               if(inv.gt.0) then
                  c1u=omega*r1
!
!     flow rotation at outlet
                  c2u=c1u*r1/r2
!
               elseif(inv.lt.0) then
                  c2u=omega*r2
!
                  c1u=c2u*r2/r1
               endif
!
            elseif(nelemswirl.gt.0) then
!     preswirl nozzle
               if(lakon(nelemswirl)(2:5).eq.'ORPN') then
                  cinput=prop(ielprop(nelemswirl)+5)
!     rotating orifices
               else if((lakon(nelemswirl)(2:5).eq.'ORMM').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORMA').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORPM').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORPA')) then
                  cinput=prop(ielprop(nelemswirl)+7)
!     forced vortex
               elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
                  cinput=prop(ielprop(nelemswirl)+7)
!     free vortex
               elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
                  cinput=prop(ielprop(nelemswirl)+9)
!     Moehring
               elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
                  cinput=prop(ielprop(nelemswirl)+10)
!     RCAVO
               elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
     &                 (lakon(nelemswirl)(2:4).eq.'ROA'))then
                  cinput=prop(ielprop(nelemswirl)+6)
!     RCAVI
               elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
                  cinput=prop(ielprop(nelemswirl)+5)
              else
                write(*,*) '*ERROR in vortex:'
                write(*,*) ' element',nelemswirl
                write(*,*) ' referred by element',nelem
                write(*,*) ' is not a swirl generating element'
                cinput=0.d0
               endif
!
               cinput=U1+K1*(cinput-U1)
!
               if(inv.gt.0) then
                  c1u=cinput
                  c2u=c1u*R1/R2
               elseif(inv.lt.0) then
                  c2u=cinput
                  c1u=c2u*R2/R1
               endif
            endif
!
!     storing the tengential velocity for later use (wirbel cascade)
            if(inv.gt.0) then
               prop(index+9)=c2u
            elseif(inv.lt.0) then
               prop(index+9)=c1u
            endif
!
!    inner rotation
!
            if(R1.lt.R2) then
               ciu=c1u
            elseif(R1.ge.R2) then
               ciu=c2u
            endif
!
            expon=kappa/km1
!
            if(R2.ge.R1) then
!
               cte1=c1u**2/(2*Cp*T1)
               cte2=1-(R1/R2)**2

               f=p2/p1-1d0-eta*((1+cte1*cte2)**expon-1d0)
!
               df(1)=-p2/p1**2
!
               df(2)=eta*expon*cte1/T1*cte2*
     &              (1+cte1*cte2)**(expon-1)
!
               df(3)=0
!
               df(4)=1/p1
!
            elseif(R2.lt.R1) then
!
               cte1=c2u**2/(2*Cp*T2)
               cte2=1-(R2/R1)**2
!
               f=p1/p2-1d0-eta*((1+cte1*cte2)**expon-1d0)
!
               df(1)=1/p2
!
               df(2)=eta*expon*cte1/T1*cte2*
     &           (1+cte1*cte2)**(expon-1)
!
               df(3)=0
!
               df(4)=-p1/p2**2
!
            endif
!
!     FORCED VORTEX
!
         elseif(lakon(nelem)(4:5).eq.'FO') then
!
!    core swirl ratio
            Kr=prop(index+4)
!
!     rotation speed (revolution per minutes) of the rotating part
!     responsible for the swirl
            omega=prop(index+5)
!
!     Temperature change
            t_chang=prop(index+6)
!
            if(R2.ge.R1) then
               Ui=omega*R1
               c1u=Ui*kr
               c2u=c1u*R2/R1
            elseif(R2.lt.R1) then
               Ui=omega*R2
               c2u=Ui*kr
               c1u=c2u*R1/R2
            endif
!
!     storing the tengential velocity for later use (wirbel cascade)
            if(inv.gt.0) then
               prop(index+7)=c2u
            elseif(inv.lt.0) then
               prop(index+7)=c1u
            endif
!
            expon=kappa/km1
!
            if(((R2.ge.R1).and.(xflow.gt.0d0))
     &           .or.((R2.lt.R1).and.(xflow.lt.0d0)))then
!
               cte1=(c1u)**2/(2*Cp*T1)
               cte2=(R2/R1)**2-1
!
                f=p2/p1-1-eta*((1+cte1*cte2)**expon-1)
!
!     pressure node1
               df(1)=-p2/p1**2
!
!     temperature node1
               df(2)=eta*expon*cte1/T1*cte2*(1+cte1*cte2)**(expon-1)
!
!     massflow nodem
               df(3)=0
!
!     pressure node2
               df(4)=1/p1
!
            elseif(((R2.lt.R1).and.(xflow.gt.0d0))
     &              .or.((R2.gt.R1).and.(xflow.lt.0d0)))then
               cte1=(c2u)**2/(2*Cp*T2)
               cte2=(R1/R2)**2-1
!
               f=p1/p2-1-eta*((1+cte1*cte2)**expon-1)
!
!     pressure node1
               df(1)=1/p2
!
!     temperature node1
               df(2)=eta*expon*cte1/T2*cte2*(1+cte1*cte2)**(expon-1)
!
!     massflow nodem
               df(3)=0
!
!     pressure node2
               df(4)=-p1/p2**2
!
            endif
         endif
!
!     outpout
!
         elseif(iflag.eq.3) then
!
         index=ielprop(nelem)
         kappa=(cp/(cp-R))
         km1=kappa-1
         pi=4.d0*datan(1.d0)
!
!     radius downstream
         r2d=prop(index+1)
!
!     radius upstream
         r1d=prop(index+2)
!
!     pressure correction factor
         eta=prop(index+3)
!
         p1=v(2,node1)
         p2=v(2,node2)
!
         xflow=v(1,nodem)*iaxial
!
         if(xflow.gt.0.d0) then
            inv=1.d0
            p1=v(2,node1)
            p2=v(2,node2)
            T1=v(0,node1)
            T2=v(0,node2)
            R1=r1d
            R2=r2d
!
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!
         elseif(xflow.lt.0.d0) then
            inv=-1.d0
            R1=r2d
            R2=r1d
            p1=v(2,node2)
            p2=v(2,node1)
            T1=v(0,node2)
            T2=v(0,node1)
            xflow=v(1,nodem)*iaxial
!
            nodef(1)=node2
            nodef(2)=node2
            nodef(3)=nodem
            nodef(4)=node1
!
         endif
!
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=2
!
         kappa=(cp/(cp-R))
!
!     FREE VORTEX
!
         if(lakon(nelem)(4:5).eq.'FR')then
!
!     rotation induced loss (correction factor)
           K1= prop(index+4)
!
!     tengential velocity of the disk at vortex entry
             U1=prop(index+5)
!
!     number of the element generating the upstream swirl
             nelemswirl=nint(prop(index+6))
!
!     rotation speed (revolution per minutes)
             omega=prop(index+7)
!
!     Temperature change
            t_chang=prop(index+8)
!
            if(omega.gt.0.d0) then
!
!     rotation speed is given if the swirl comes from a rotating part
!     typically the blade of a coverplate
!
!     C_u is given by radius r1d (see definition of the flow direction)
!     C_u related to radius r2d is a function of r1d
!
               if(inv.gt.0) then
                  c1u=omega*r1
!
!     flow rotation at outlet
                  c2u=c1u*r1/r2
!
               elseif(inv.lt.0) then
                  c2u=omega*r2
!
                  c1u=c2u*r2/r1
               endif
!
            elseif(nelemswirl.gt.0) then
!     swirl generating element
!
!     preswirl nozzle
               if(lakon(nelemswirl)(2:5).eq.'ORPN') then
                  cinput=prop(ielprop(nelemswirl)+5)
!     rotating orifices
               elseif((lakon(nelemswirl)(2:5).eq.'ORMM').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORMA').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORPM').or.
     &                 (lakon(nelemswirl)(2:5).eq.'ORPA')) then
                  cinput=prop(ielprop(nelemswirl)+7)
!     forced vortex
               elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
                  cinput=prop(ielprop(nelemswirl)+7)
!     free vortex
               elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
                  cinput=prop(ielprop(nelemswirl)+9)
!     Moehring
               elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
                  cinput=prop(ielprop(nelemswirl)+10)
!     RCAVO
               elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
     &                 (lakon(nelemswirl)(2:4).eq.'ROA'))then
                  cinput=prop(ielprop(nelemswirl)+6)
!     RCAVI
               elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
                  cinput=prop(ielprop(nelemswirl)+5)
              else
                write(*,*) '*ERROR in vortex:'
                write(*,*) ' element',nelemswirl
                write(*,*) ' referred by element',nelem
                write(*,*) ' is not a swirl generating element'
                cinput=0.d0
               endif
!
               cinput=U1+K1*(cinput-U1)
!
               if(inv.gt.0) then
                  c1u=cinput
                  c2u=c1u*R1/R2
               elseif(inv.lt.0) then
                  c2u=cinput
                  c1u=c2u*R2/R1
               endif
            endif
!
!     storing the tengential velocity for later use (wirbel cascade)
            if(inv.gt.0) then
               prop(index+9)=c2u
            elseif(inv.lt.0) then
               prop(index+9)=c1u
            endif
!
!    inner rotation
!
            if(R1.lt.R2) then
               ciu=c1u
            elseif(R1.ge.R2) then
               ciu=c2u
            endif
!
            expon=kappa/km1
!
            if(R2.ge.R1) then
!
               cte1=c1u**2/(2*Cp*T1)
               cte2=1-(R1/R2)**2

               f=p2/p1-1d0-eta*((1+cte1*cte2)**expon-1d0)
!
               df(1)=-p2/p1**2
!
               df(2)=eta*expon*cte1/T1*cte2*
     &              (1+cte1*cte2)**(expon-1)
!
               df(3)=0
!
               df(4)=1/p1
!
            elseif(R2.lt.R1) then
!
               cte1=c2u**2/(2*Cp*T2)
               cte2=1-(R2/R1)**2
!
               f=p1/p2-1d0-eta*((1+cte1*cte2)**expon-1d0)
!
               df(1)=1/p2
!
               df(2)=eta*expon*cte1/T1*cte2*
     &           (1+cte1*cte2)**(expon-1)
!
               df(3)=0
!
               df(4)=-p1/p2**2
!
            endif
!
!     FORCED VORTEX
!
         elseif(lakon(nelem)(4:5).eq.'FO') then
!
!    core swirl ratio
            Kr=prop(index+4)
!
!     rotation speed (revolution per minutes) of the rotating part
!     responsible for the swirl
            omega=prop(index+5)
!
!     Temperature change
            t_chang=prop(index+6)
!
!    no element generating the upstream swirl
             nelemswirl=0
!
            if(R2.ge.R1) then
               Ui=omega*R1
               c1u=Ui*kr
               c2u=c1u*R2/R1
            elseif(R2.lt.R1) then
               Ui=omega*R2
               c2u=Ui*kr
               c1u=c2u*R1/R2
            endif
!
!     storing the tengential velocity for later use (wirbel cascade)
            if(inv.gt.0) then
               prop(index+7)=c2u
            elseif(inv.lt.0) then
               prop(index+7)=c1u
            endif
!
            expon=kappa/km1
         endif
!
         xflow_oil=0.d0
!
         write(1,*) ''
         write(1,55) ' from node ',node1,
     &' to node ', node2,' :   air massflow rate = ',xflow,' ',
     &' , oil massflow rate = ',xflow_oil,' '
!
         if(inv.eq.1) then
            write(1,56)'       Inlet node ',node1,' :     Tt1 = ',T1,
     &           '  , Ts1 = ',T1,'  , Pt1 = ',p1
            write(1,*)'             Element ',nelem,lakon(nelem)
            write(1,57)'             C1u = ',C1u,
     &'  , C2u = ',C2u
            write(1,56)'      Outlet node ',node2,' :    Tt2 = ',T2,
     &           '  , Ts2 = ',T2,'  , Pt2 = ',p2
!
         else if(inv.eq.-1) then
            write(1,56)'       Inlet node ',node2,':     Tt1 = ',T1,
     &           '  , Ts1 = ',T1,'  , Pt1 = ',p1
            write(1,*)'             Element ',nelem,lakon(nelem)
            write(1,57)'             C1u = ',C1u,
     &'  , C2u = ',C2u
            write(1,56)'      Outlet node ',node1,'     Tt2 = ',
     &           T2,'  , Ts2 = ',T2,'  , Pt2 = ',p2
         endif
      endif
!
 55   format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
 56   format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
 57   format(1x,a,e11.4,a,e11.4,a)
!
      xflow=xflow/iaxial
      df(3)=df(3)*iaxial
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c { dg-do compile }
c { dg-options "-Wall" }
c PR42852 -Wall warns about truncated lines when only a continuation character is truncated
      print *, "Hello!"                                                  &   !xxxxx
     & //   " World!"
      print *, "Hello!"                                                  &   xxxxx
     & //   " World!"
      print *, "Hello!"                                                  //
     & //   " World!"
      end
c { dg-warning "Line truncated" " " { target *-*-* } 6 }
c { dg-warning "Line truncated" " " { target *-*-* } 8 }

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
	SUBROUTINE GDBHLO (  hilo, hlsym, grid, kx, ky, iret )
C************************************************************************
C* GDBHLO								*
C*									*
C* This subroutine finds and marks relative extrema for the grid.	*
C*									*
C* GDBHLO  ( HILO, HLSYM, GRID, KX, KY, IRET )				*
C*									*
C* Input parameters:							*
C*	HILO		CHAR*		HILO input string 		*
C*	HLSYM		CHAR*		Text specifier			*
C*	GRID (KX,KY)	REAL		Scaled grid data		*
C*	KX		INTEGER		Grid x dimension		*
C*	KY		INTEGER		Grid y dimension		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* K. Brill/NMC		 6/93						*
C* D. Keiser/GSC	 7/95	Changed DC to DG in call to find extrema*
C* D. Keiser/GSC	10/95	Changed ierr to acct for pos numbers	*
C* S. Jacobs/NCEP	 5/96	Fixed plotting location of the values	*
C* S. Jacobs/NCEP	 9/97	Changed call to GQTEXT and GSTEXT	*
C* D.W.Plummer/NCEP	 6/01	Added special symbols processing	*
C* T. Lee/SAIC		10/01	Grouped symbols and labels		*
C* T. Lee/SAIC		10/02	Reset symbol size			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	PARAMETER	( JHGRP = 5, JLGRP = 6 )
C*
	CHARACTER*(*)	hilo, hlsym
	REAL		grid (*)
C*
	LOGICAL		intflg, valflg (2)
	CHARACTER*32	csymbl (2), text, symtxt, valtxt
	INTEGER		isymbl (2), isymtp (2), iprecn (2), knt (2),
     +			icolr (2)
	REAL		range (4)
	REAL		xmax (128), ymax (128), vmax (128)
	REAL		xmin (128), ymin (128), vmin (128)
	INTEGER		isymwd
C------------------------------------------------------------------------
	iret  = 0
	CALL IN_HILO ( hilo, icolr, csymbl, isymbl, isymtp,
     +		       valflg, iprecn, range, krad, knt, intflg,
     +		       iret )
	IF ( ( icolr (1) + icolr (2) ) .eq. 0 ) RETURN
	IF ( knt (1) .gt. 128) knt (1) = 128
	IF ( knt (2) .gt. 128) knt (2) = 128
C
C*	Find the extrema.
C
	CALL DG_HILO ( grid, kx, ky, krad, intflg, knt (1),
     +		       knt (2), range (1), range (2), range (3),
     +		       range (4), nummax, xmax, ymax, vmax,
     +		       nummin, xmin, ymin, vmin, ierr )
C
	IF ( ierr .ge. 0 ) THEN
C
C*          Write out error message regarding internal buffers
C*	    in DG_HILO.
C
            IF ( ierr .ne. 0 ) CALL ER_WMSG ( 'DG', ierr, ' ', ier2 )
C
C*	    Save the TEXT characteristics.
C
	    CALL GQTEXT ( ifntsv, ihwsv, sizsv, iwidsv,
     +			  ibdrsv, irrtsv, ijstsv, ier )
C
C*	    Set the TEXT for the HILO strings.
C
            IF ( hlsym .eq. ' ' ) THEN
                text = '1.5;1'
            ELSE
                text = hlsym
            END IF
 	    CALL IN_HLSY ( text, symtxt, valtxt, labloc, ier )
C
C*	    Set the marker size to same size as text..
C
	    sizmrk = 0
	    IF ( isymtp (1) .ge. 2 .or. isymtp (2) .ge. 2 ) THEN
		CALL ST_RLST ( symtxt, '/', 0., 1, sizmrk, n, ier )
	    END IF
	    CALL IN_TEXT ( symtxt, ier )
	    CALL GQTEXT ( itxfn, itxhw, siztxt, itxwid, ibrdr, irrotn,
     +			  ijust, ier )
	    isymwd = itxwid
C
C*	    Label the highs.
C
	    siztxt = 0
	    IF ( icolr (1) .gt. 0 .and. nummax .gt. 0 ) THEN
		CALL GSCOLR ( icolr (1), ier )
		IF ( isymtp (1) .ge. 2 ) THEN
		    IF ( isymtp (1) .eq. 3 )  THEN
			CALL GSSPCL ( sizmrk, isymwd, ier )
			code    = isymbl (1)
			iixoff  = 0
			iiyoff  = 0
		    ELSE
		      CALL GSMRKR ( isymbl (1), 0, sizmrk, isymwd, ier )
		    END IF
		END IF
		DO i = 1, nummax
		    CALL GSGRP ( JHGRP, ier )
C
		    IF  ( isymtp (1) .ge. 2 )  THEN
			IF  ( isymtp (1) .eq. 3 )  THEN
			    CALL GSPCL ( 'G', 1, code,
     +					 xmax (i), ymax (i),
     +					 iixoff, iiyoff, ier )
			  ELSE
			    CALL GMARK ( 'G', 1, xmax(i), ymax(i), ier )
			END IF
		    END IF
C
		    IF ( isymtp (1) .eq. 1 ) THEN
			CALL ST_LSTR ( csymbl (1), lens, ier )
			ixoff = -  ( lens - 1 )
			CALL IN_TEXT ( symtxt, ier )
			CALL GQTEXT ( itxfn, itxhw, siztxt,
     +				      itxwid, ibrdr, irrotn,
     +				      ijust, ier )
			CALL GTEXT ( 'G', xmax (i), ymax (i),
     +				         csymbl (1), 0., ixoff, 0,
     +					 ier )
		    END IF
		    size1 = MAX ( siztxt, sizmrk )
		    IF ( valflg (1) ) THEN
			IF ( iprecn (1) .gt. 0 ) THEN
			    CALL ST_RLCH ( vmax (i), iprecn (1),
     +					   text, ier )
			ELSE
			    ivm = NINT ( vmax (i) )
			    CALL ST_INCH ( ivm, text, ier )
			END IF
			CALL ST_LSTR ( text, lens, ier )
			ixoff = - ( lens - 1 )
			IF ( labloc .eq. 1 ) ixoff = (-2) * lens
			IF ( labloc .eq. 3 ) ixoff = 1
			CALL IN_TEXT ( valtxt, ier )
			CALL GQTEXT ( itxfn, itxhw, size2, itxwid,
     +				      ibrdr, irrotn, ijust, ier )
			iyoff = NINT ( size1 / size2 + 1 )
			CALL GTEXT ( 'G', xmax (i), ymax (i), text,
     +				 	0., ixoff, -iyoff, ier )
		    END IF
		    CALL GEGRP ( ier )
		END DO
	    END IF
C
C*	    Label the lows.
C
	    siztxt = 0
	    IF ( icolr (2) .gt. 0 .and. nummin .gt. 0 ) THEN
		CALL GSCOLR ( icolr (2), ier  )
		IF ( isymtp (2) .ge. 2 ) THEN
		    IF ( isymtp (2) .eq. 3 )  THEN
			CALL GSSPCL ( sizmrk, isymwd, ier )
			code   = isymbl (2)
			iixoff = 0
			iiyoff = 0
		    ELSE
		      CALL GSMRKR ( isymbl (2), 0, sizmrk, isymwd, ier )
		    END IF
		END IF
		DO i = 1, nummin
		    CALL GSGRP ( JLGRP, ier )
		    IF  ( isymtp (2) .ge. 2 )  THEN
			IF  ( isymtp (2) .eq. 3 )  THEN
			    CALL GSPCL ( 'G', 1, code,
     +					 xmin (i), ymin (i),
     +					 iixoff, iiyoff, ier )
			  ELSE
			    CALL GMARK ( 'G', 1, xmin(i), ymin(i), ier )
			END IF
		    END IF
C
		    IF ( isymtp (2) .eq. 1 ) THEN
			CALL ST_LSTR ( csymbl (2), lens, ier )
			ixoff = -  ( lens - 1 )
			CALL IN_TEXT ( symtxt, ier )
			CALL GQTEXT ( itxfn, itxhw, siztxt,
     +				      itxwid, ibrdr, irrotn,
     +				      ijust, ier )
			CALL GTEXT ( 'G', xmin (i), ymin (i),
     +				         csymbl (2), 0., ixoff, 0,
     +					 ier )
		    END IF
		    size1 = MAX ( siztxt, sizmrk )
		    IF ( valflg (2) ) THEN
			IF ( iprecn (2) .gt. 0 ) THEN
     			    CALL ST_RLCH ( vmin (i), iprecn (2),
     +					   text, ier )
			ELSE
			    ivm = NINT ( vmin (i) )
			    CALL ST_INCH ( ivm, text, ier )
			END IF
			CALL ST_LSTR ( text, lens, ier )
			ixoff = - ( lens - 1 )
			IF ( labloc .eq. 1 ) ixoff = (-2) * lens
			IF ( labloc .eq. 3 ) ixoff = 1
			CALL IN_TEXT ( valtxt, ier )
			CALL GQTEXT ( itxfn, itxhw, size2, itxwid,
     +				      ibrdr, irrotn, ijust, ier )
			iyoff = NINT ( size1 / size2 + 1 )
			CALL GTEXT ( 'G', xmin (i), ymin (i), text,
     +				 	0., ixoff, -iyoff, ier )
		    END IF
		    CALL GEGRP ( ier )
	    	END DO
	    END IF
C
C*	    Reset text characteristics.
C
	    CALL GSTEXT ( ifntsv, ihwsv, sizsv, iwidsv,
     +			  ibdrsv, irrtsv, ijstsv, ier )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c
c diagsic.f modified from AY's sea ice diag code nre 27/2/4
c
c AY (23/09/04) : upgraded to output sea-ice albedo

      subroutine diagsic(istep)

#include "seaice.cmn"

      real hamax(2),hamin(2),ticemax,ticemin,haaav(2),ticeaav,
     :     albicemax,albicemin,albiceaav,area

      integer i,j,l,icepts,istep

      do l=1,2
         hamax(l) = varice(l,1,1)
         hamin(l) = hamax(l)
         haaav(l) = 0.0
      enddo
      ticemax = tice(1,1)
      ticemin = ticemax
      ticeaav = 0.0
      albicemax = albice(1,1)
      albicemin = albicemax
      albiceaav = 0.0
      icepts = 0
      area = 0.0

      do j=1,jmax
         do i=1,imax
            do l=1,2
               hamax(l) = max(hamax(l),varice(l,i,j))
               hamin(l) = min(hamin(l),varice(l,i,j))
            enddo
            ticemax = max(ticemax,tice(i,j))
            ticemin = min(ticemin,tice(i,j))
            albicemax = max(albicemax,albice(i,j))
            albicemin = min(albicemin,albice(i,j))
            if(varice(2,i,j).gt.0.0)then
               haaav(1) = haaav(1) + varice(1,i,j)*ds(j)
               haaav(2) = haaav(2) + varice(2,i,j)*ds(j)
               ticeaav = ticeaav + tice(i,j)
               albiceaav = albiceaav + albice(i,j)
               icepts = icepts + 1
               area = area + ds(j)
            endif
         enddo
      enddo
      if(icepts.gt.0)then
crma         haaav(1) = haaav(1)/icepts
crma         haaav(2) = haaav(2)/icepts
crma         ticeaav  = ticeaav /icepts
crma         albiceaav = albiceaav / icepts
         haaav(1) = haaav(1)/area
         haaav(2) = haaav(2)/area
         ticeaav  = ticeaav /area
         albiceaav = albiceaav / area

      endif
      print*,'Sea-ice properties at time ',istep
      print*,'Ice height : max = ',hamax(1),'; min = ',hamin(1)
     +     ,'; avg = ',haaav(1)
      print*,'Ice cover  : max = ',hamax(2),'; min = ',hamin(2)
     +     ,'; avg = ',haaav(2)
      print*,'Ice temp.  : max = ',ticemax,'; min = ',ticemin
     +     ,'; avg = ',ticeaav
      print*,'Ice albd.  : max = ',albicemax,'; min = ',albicemin
     +     ,'; avg = ',albiceaav

      end


text:
program i
  logical :: a
end program i

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
!*****************************************************************************
!
! Rokko: Integrated Interface for libraries of eigenvalue decomposition
!
! Copyright (C) 2012-2020 by Rokko Developers https://github.com/t-sakashita/rokko
!
! 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 rokko
  use various_mpi_comm_mod
  implicit none
  integer :: dim
  type(rokko_distributed_matrix) :: mat, Z
  type(rokko_grid) :: grid
  type(rokko_mapping_bc) :: map
  type(rokko_parallel_dense_ev) :: solver
  type(rokko_eigen_vector) :: w
  character(len=:), allocatable :: library, routine
  character(len=:), allocatable :: library_routine, tmp_str
  integer :: comm
  integer :: provided, ierr, myrank, nprocs
  integer :: i

  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)
  comm = create_cart_comm(mpi_comm_world)

  if (comm /= MPI_COMM_NULL) then
     if (command_argument_count() >= 1) then
        call get_command_argument_deferred(1, library_routine)
     else
        library_routine = rokko_parallel_dense_ev_default_solver()
     endif
     call rokko_split_solver_name(library_routine, library, routine)

     if (command_argument_count() == 2) then
        call get_command_argument_deferred(2, tmp_str)
        read(tmp_str, *) dim
     else
        dim = 10
     endif

     if (myrank == 0) then
        print *,"library = ", library
        print *,"routine = ", routine
        print *,"dimension = ", dim
     endif

     call rokko_construct(solver, library)
     call rokko_construct(grid, comm, rokko_grid_row_major)
     call rokko_default_mapping(solver, dim, grid, map)
     call rokko_construct(mat, map)
     call rokko_construct(Z, map)
     call rokko_construct(w, dim)

     if (myrank == 0) then
        write(*,'("Created ", i0, "x", i0, " size communicator with new cartesian topology")') &
             & rokko_get_nprow(grid), rokko_get_npcol(grid)
     endif

     ! generate frank matrix
     call rokko_frank_matrix_generate(mat)
     call rokko_print(mat)

     call rokko_diagonalize(solver, mat, w, Z)

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

     call rokko_destruct(mat)
     call rokko_destruct(Z)
     call rokko_destruct(w)
     call rokko_destruct(solver)
     call rokko_destruct(map)
     call rokko_destruct(grid)

     call mpi_comm_free(comm, ierr)
  endif
  call MPI_finalize(ierr)
end program frank_matrix

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE REDUCE_KFFD_TO_KAAD ( PART_VEC_F_AO )

! Call routines to reduce the KFFD differential stiffness matrix from the F-set to the A, O-sets

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  ERR, F04, F06, L2E, LINK2E, L2E_MSG, SC1, WRT_ERR, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FACTORED_MATRIX, FATAL_ERR, NDOFF, NDOFA, NDOFO, NTERM_KFFD, NTERM_KAAD,    &
                                         NTERM_KAOD, NTERM_KOOD, NTERM_KOODs, NTERM_GOA
      USE PARAMS, ONLY                :  EPSIL, KOORAT, SPARSTOR, RCONDK
      USE TIMDAT, ONLY                :  HOUR, MINUTE, SEC, SFRAC, TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  REDUCE_KFFD_TO_KAAD_BEGEND
      USE CONSTANTS_1, ONLY           :  ONE
      USE SPARSE_MATRICES, ONLY       :  I_KFFD, J_KFFD, KFFD, I_KAAD, J_KAAD, KAAD, I_KAOD, J_KAOD, KAOD, I_GOA, J_GOA, GOA,      &
                                         I_KOOD, I2_KOOD, J_KOOD, KOOD, I_KOODs, I2_KOODs, J_KOODs, KOODs

      USE SPARSE_MATRICES, ONLY       :  SYM_GOA, SYM_KFFD, SYM_KAAD, SYM_KAOD, SYM_KOOD
      USE SCRATCH_MATRICES

      USE REDUCE_KFFD_TO_KAAD_USE_IFs

      IMPLICIT NONE

      CHARACTER, PARAMETER            :: CR13 = CHAR(13)   ! This causes a carriage return simulating the "+" action in a FORMAT
      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'REDUCE_KFFD_TO_KAAD'
      CHARACTER(  1*BYTE)             :: SYM_CRS2            ! Storage format for matrix CRS2 (either 'Y' for sym storage or
!                                                              'N' for nonsymmetric storage)

      INTEGER(LONG), INTENT(IN)       :: PART_VEC_F_AO(NDOFF)! Partitioning vector (F set into A and O sets)
      INTEGER(LONG)                   :: AROW_MAX_TERMS      ! Output from MATMULT_SFS_NTERM and input to MATMULT_SFS

      INTEGER(LONG)                   :: I,J                 ! DO loop indices
      INTEGER(LONG)                   :: KAAD_ROW_MAX_TERMS  ! Output from subr PARTITION_SIZE (max terms in any row of matrix)
      INTEGER(LONG)                   :: KAOD_ROW_MAX_TERMS  ! Output from subr PARTITION_SIZE (max terms in any row of matrix)
      INTEGER(LONG)                   :: KOOD_ROW_MAX_TERMS  ! Output from subr PARTITION_SIZE (max terms in any row of matrix)
      INTEGER(LONG)                   :: NTERM_CRS1          ! Number of terms in matrix CRS1
      INTEGER(LONG)                   :: NTERM_CRS2          ! Number of terms in matrix CRS2
      INTEGER(LONG), PARAMETER        :: NUM1        = 1     ! Used in subr's that partition matrices
      INTEGER(LONG), PARAMETER        :: NUM2        = 2     ! Used in subr's that partition matrices
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = REDUCE_KFFD_TO_KAAD_BEGEND

      INTRINSIC                       :: DABS

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Partition KAAD from KFFD (This is KAAD before reduction, or KAAD(bar) )

      IF (NDOFA > 0) THEN

         CALL PARTITION_SS_NTERM ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD,      PART_VEC_F_AO, PART_VEC_F_AO,  &
                                    NUM1, NUM1, KAAD_ROW_MAX_TERMS, 'KAAD', NTERM_KAAD, SYM_KAAD )

         CALL ALLOCATE_SPARSE_MAT ( 'KAAD', NDOFA, NTERM_KAAD, SUBR_NAME )

         IF (NTERM_KAAD > 0) THEN
            CALL PARTITION_SS ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD, KFFD, PART_VEC_F_AO, PART_VEC_F_AO,    &
                                 NUM1, NUM1, KAAD_ROW_MAX_TERMS, 'KAAD', NTERM_KAAD, NDOFA, SYM_KAAD, I_KAAD, J_KAAD, KAAD )
         ENDIF

      ENDIF

! Partition KAOD from KFFD

      IF ((NDOFA > 0) .AND. (NDOFO > 0)) THEN

         CALL PARTITION_SS_NTERM ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD,      PART_VEC_F_AO, PART_VEC_F_AO,  &
                                    NUM1, NUM2, KAOD_ROW_MAX_TERMS, 'KAOD', NTERM_KAOD, SYM_KAOD )

         CALL ALLOCATE_SPARSE_MAT ( 'KAOD', NDOFA, NTERM_KAOD, SUBR_NAME )

         IF (NTERM_KAOD > 0) THEN
            CALL PARTITION_SS ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD, KFFD, PART_VEC_F_AO, PART_VEC_F_AO,    &
                                 NUM1, NUM2, KAOD_ROW_MAX_TERMS, 'KAOD', NTERM_KAOD, NDOFA, SYM_KAOD, I_KAOD, J_KAOD, KAOD )
         ENDIF

      ENDIF

! Partition KOOD from KFFD

      IF (NDOFO > 0) THEN

         CALL PARTITION_SS_NTERM ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD,      PART_VEC_F_AO, PART_VEC_F_AO,  &
                                    NUM2, NUM2, KOOD_ROW_MAX_TERMS, 'KOOD', NTERM_KOOD, SYM_KOOD )

         CALL ALLOCATE_SPARSE_MAT ( 'KOOD', NDOFO, NTERM_KOOD, SUBR_NAME )

         IF (NTERM_KOOD > 0) THEN
            CALL PARTITION_SS ( 'KFFD', NTERM_KFFD, NDOFF, NDOFF, SYM_KFFD, I_KFFD, J_KFFD, KFFD, PART_VEC_F_AO, PART_VEC_F_AO,    &
                                 NUM2, NUM2, KOOD_ROW_MAX_TERMS, 'KOOD', NTERM_KOOD, NDOFO, SYM_KOOD, I_KOOD, J_KOOD, KOOD )
         ENDIF

      ENDIF

! ----------------------------------------------------------------------------------------------------------------------------------
! Reduce KFFD to KAAD = KAAD(bar) + KAOD*GOA. Note: (KAOD*GOA)' + GOA'*KOOD*GOA = 0 by definition of GOA = -KOOD(-1)*KAOD'
! If PARAM MATSPARS = 'Y', then we use sparse matrix operations (multiply/add/transpose). If not, use full matrix operations

         IF (NTERM_KAOD > 0) THEN                          ! Calc KAOD*GOA & add it orig KAAD

                                                           ! CCS1 will be sparse CCS format version of sparse CRS matrix GOA
            CALL ALLOCATE_SCR_CCS_MAT ( 'CCS1', NDOFA, NTERM_GOA, SUBR_NAME )
            CALL SPARSE_CRS_SPARSE_CCS ( NDOFO, NDOFA, NTERM_GOA, 'GOA', I_GOA, J_GOA, GOA, 'CCS1', J_CCS1, I_CCS1, CCS1, 'Y' )

                                                           ! Sparse multiply to get CRS1 = KAOD*GOA. Use CCS1 for GOA CCS
            CALL MATMULT_SSS_NTERM ( 'KAOD' , NDOFA, NTERM_KAOD, SYM_KAOD, I_KAOD , J_KAOD ,                                       &
                                     'GOA' , NDOFA, NTERM_GOA, SYM_GOA, J_CCS1, I_CCS1, AROW_MAX_TERMS,                            &
                                     'CRS1' ,       NTERM_CRS1 )

            CALL ALLOCATE_SCR_CRS_MAT ( 'CRS1', NDOFA, NTERM_CRS1, SUBR_NAME )

            CALL MATMULT_SSS ( 'KAOD' , NDOFA, NTERM_KAOD , SYM_KAOD, I_KAOD , J_KAOD , KAOD ,                                     &
                               'GOA' , NDOFA, NTERM_GOA , SYM_GOA, J_CCS1, I_CCS1, CCS1, AROW_MAX_TERMS,                           &
                               'CRS1', ONE  , NTERM_CRS1,          I_CRS1, J_CRS1, CRS1 )

            CALL DEALLOCATE_SCR_MAT ( 'CCS1' )

                                                           ! CRS1 = KAOD*GOA has all nonzero terms in it.
            IF      (SPARSTOR == 'SYM   ') THEN            !      If SPARSTOR == 'SYM   ', rewrite CRS1 as sym in CRS2

               CALL SPARSE_CRS_TERM_COUNT ( NDOFA, NTERM_CRS1, 'CRS1 = KAOD*GOA all nonzeros', I_CRS1, J_CRS1, NTERM_CRS2 )
               CALL ALLOCATE_SCR_CRS_MAT ( 'CRS2', NDOFA, NTERM_CRS2, SUBR_NAME )
               CALL CRS_NONSYM_TO_CRS_SYM ( 'CRS1 = KAOD*GOA all nonzeros', NDOFA, NTERM_CRS1, I_CRS1, J_CRS1, CRS1,               &
                                            'CRS2 = KAOD*GOA stored sym'  ,        NTERM_CRS2, I_CRS2, J_CRS2, CRS2 )
               SYM_CRS2 = 'Y'

            ELSE IF (SPARSTOR == 'NONSYM') THEN            !      If SPARSTOR == 'NONSYM', rewrite CRS1 in CRS2 with NTERM_CRS2

               NTERM_CRS2 = NTERM_CRS1
               CALL ALLOCATE_SCR_CRS_MAT ( 'CRS2', NDOFA, NTERM_CRS2, SUBR_NAME )
               DO I=1,NDOFA+1
                  I_CRS2(I) = I_CRS1(I)
               ENDDO
               DO I=1,NTERM_CRS2
                  J_CRS2(I) = J_CRS1(I)
                    CRS2(I) =   CRS1(I)
               ENDDO
               SYM_CRS2 = 'N'

            ELSE                                           ! Error - incorrect SPARSTOR

               WRITE(ERR,932) SUBR_NAME, SPARSTOR
               WRITE(F06,932) SUBR_NAME, SPARSTOR
               FATAL_ERR = FATAL_ERR + 1
               CALL OUTA_HERE ( 'Y' )

            ENDIF

            CALL DEALLOCATE_SCR_MAT ( 'CRS1' )             ! Deallocate CRS1
                                                           ! Sparse add to get CRS1 = KAAD-bar + CRS2
            CALL MATADD_SSS_NTERM ( NDOFA, 'KAAD-bar', NTERM_KAAD, I_KAAD , J_KAAD , SYM_KAAD , 'KOOD*GOA', NTERM_CRS2,            &
                                                                 I_CRS2, J_CRS2, SYM_CRS2, 'CRS1', NTERM_CRS1 )
            CALL ALLOCATE_SCR_CRS_MAT ( 'CRS1', NDOFA, NTERM_CRS1, SUBR_NAME )
            CALL MATADD_SSS (NDOFA, 'KAAD-bar' , NTERM_KAAD , I_KAAD , J_KAAD , KAAD , ONE, 'KOOD*GOA', NTERM_CRS2,                &
                             I_CRS2, J_CRS2, CRS2, ONE, 'CRS1', NTERM_CRS1, I_CRS1, J_CRS1, CRS1 )

            CALL DEALLOCATE_SCR_MAT ( 'CRS2' )             ! Deallocate CRS2

            NTERM_KAAD = NTERM_CRS1                        ! Reallocate KAAD to be size of CRS1
            WRITE(SC1, * ) '    Reallocate KAAD'
            WRITE(SC1,12345,ADVANCE='NO') '       Deallocate KAAD', CR13
            CALL DEALLOCATE_SPARSE_MAT ( 'KAAD' )
            WRITE(SC1,12345,ADVANCE='NO') '       Allocate   KAAD', CR13
            CALL ALLOCATE_SPARSE_MAT ( 'KAAD', NDOFA, NTERM_KAAD, SUBR_NAME )
                                                           ! Set KAAD = CRS1
            DO I=1,NDOFA+1
               I_KAAD(I) = I_CRS1(I)
            ENDDO
            DO J=1,NTERM_KAAD
               J_KAAD(J) = J_CRS1(J)
                 KAAD(J) =   CRS1(J)
            ENDDO

            CALL DEALLOCATE_SCR_MAT ( 'CRS1' )             ! Deallocate CRS1

         ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  932 FORMAT(' *ERROR   932: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' PARAMETER SPARSTOR MUST BE EITHER "SYM" OR "NONSYM" BUT VALUE IS ',A)

 2400 FORMAT(' *ERROR  2400: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' THERE IS AN O-SET BUT GUYAN REDUCTION MATRIX GOA HAS ',I12,' TERMS IN IT. MUST BE > 0')

12345 FORMAT(A,10X,A)

! **********************************************************************************************************************************

      END SUBROUTINE REDUCE_KFFD_TO_KAAD

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


text:
Subroutine PrintTitle(iUnit,nBuild)

Write(iUnit,101)nbuild

101 FORMAT(38x,'********************************************************'/ &
		   38x,'*              ReadExcitedStates v.1                   *'/ &
		   38x,'*               (build number',i6,')                   *'/ &
		   38x,'*                                                      *'/ &
		   38x,'*   Process excited states information from G0X files  *'/ &
		   38x,'*                                                      *'/ &
		   38x,'*     Originally written by Stanislav K. Ignatov       *'/ &
		   38x,'*          Theoretical Chemistry Group                 *'/ &
		   38x,'* N.I.Lobachevsky State University of Nizhny Novgorod  *'/ &
		   38x,'*                 ignatov@ichem.unn.ru                 *'/ &
		   38x,'*               http://www.qchem.unn.ru                *'/ &
		   38x,'*             Nizhny Novgorod, Russia, 2016            *'/ &
		   38x,'********************************************************'/)

End

max_stars_repo_name: skignatov/ReadExcitedStates

max_issues_repo_name: skignatov/ReadExcitedStates

max_forks_repo_name: skignatov/ReadExcitedStates


text:
       PROGRAM POTEN1_3D
C++::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C
C  This program calculates a Schwarzschild's model
C  of a triaxial galaxy using De Zeeuw's approximation (1983)
C   and projects the density on one of the principal planes
C
C JLP Version of 10-10-86
C--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C
	PARAMETER (IDIM=1000)
	PARAMETER (IDIM1=50)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*4 IMAGE(IDIM,IDIM)
	REAL*8 PHI(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*8 RHO(0:IDIM1,0:IDIM1,0:IDIM1)
	CHARACTER NAME*40,COMMENTS*80,ANS*1,AGAIN*1
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ
C
	PI=3.14159265358979323846
C
C GRAV=6.670 E-8 in dyn.cm2.g-2  (C.G.S.), but like QUINN
C we take GRAV=1.
	GRAV=1.D0

	PRINT *,' PROGRAM POTEN VERSION OF JUNE 8th 1986'

C Inquire the format of the files :
	CALL JLP_INQUIFMT

	OPEN(1,FILE='poten.log',STATUS='unknown',
     1	ACCESS='SEQUENTIAL')
C
C  Input of the parameters
C
10	FORMAT(A)

	PRINT *,' NUMBER OF POINTS FOR THE MODEL (LESS THAN 50)?'
	READ(5,*) NPTS
1000	WRITE (6,335)
335	FORMAT(' Enter P, Q : ',$)
	READ (5,*) PP,QQ
	PRINT *,' ENTER PASX, PASY, PASZ :'
	READ(5,*) PASX,PASY,PASZ

	WRITE(1,340),PP,QQ,PASX,PASY,PASZ
340	FORMAT(4X,'P =',F8.3,'  Q =',F8.3,5X,'PASX, PASY, PASZ :',
     1	3(F9.3,2X),/)

C Calculates C1 and C3

	PP2=PP*PP
	QQ2=QQ*QQ
	PP3=PP2*PP
	QQ3=QQ2*QQ
	ASYM2=2.D0*PP2*QQ2-PP2-QQ2
	SYM2=PP2*QQ2+PP2+QQ2
	SYM3=1.D0+PP3+QQ3

	D1=ASYM2/(14.D0*SYM2)
	D2=(1.D0-(3.D0/SYM3))/6.D0
	D3=(PP2-QQ2)/(28.D0*SYM2)
	D4=(PP3-QQ3)/(12.D0*SYM3)

C	C1=(4.D0/35.D0)*ASYM2/RRR
C	C3=(2.D0/35.D0)*(PP2-QQ2)/RRR

	C1=(D1**0.6)*(D2**0.4)
	C2=(D1/D2)**0.4
	C3=(D3**0.6)*(D4**0.4)
	C4=(D3/D4)**0.4

	CTE1=(4.D0/PI)/(4*PI*GRAV)

	WRITE(1,341) C1,C2,C3,C4
341	FORMAT(4X,'C1 =',E12.5,'  C2 =',E12.5,
     1	'C3 =',E12.5,'  C4 =',E12.5,/,/,11X,'X',
     1	13X,'Y',13X,'Z',10X,'PHI(X,Y,Z)',7X,'RHO(X,Y,Z)',/)

	DO 100 IX=0,NPTS
	X=FLOAT(IX)*PASX
	X2=X*X

	DO 100 IY=0,NPTS
	Y=FLOAT(IY)*PASY
	Y2=Y*Y

	DO 100 IZ=0,NPTS
	Z=FLOAT(IZ)*PASZ
	Z2=Z*Z

	RAD2=X2+Y2+Z2
	IF(RAD2.EQ.0.) RAD2=1.E-12
	RAD=DSQRT(RAD2)

C Calculating U, V and W

	URAD=(-1.D0/RAD)*DLOG(RAD+DSQRT(1+RAD2))+1.D0
	VRAD=(-1.D0*C1*RAD2)/((1.D0+C2*RAD2)**1.5)
	WRAD=(-1.D0*C3*RAD2)/((1.D0+C4*RAD2)**1.5)

C Remember :CTE1=(4.D0/PI)/(4*PI*GRAV)
C
	FRAD=1.D0/((1.D0+RAD2)**1.5)
	GRAD=3.D0*C1*C2*(7.D0+2.D0*C2*RAD2)*RAD2/
     1	((1.D0+C2*RAD2)**3.5)
	HRAD=3.D0*C3*C4*(7.D0+2.D0*C4*RAD2)*RAD2/
     1	((1.D0+C4*RAD2)**3.5)

	AA=(2.D0*Z2-X2-Y2)/(2.D0*RAD2)
	BB=3.D0*(X2-Y2)/RAD2

C
	PHI(IX,IY,IZ)=URAD-AA*VRAD+BB*WRAD
	RHO(IX,IY,IZ)=FRAD-AA*GRAD+BB*HRAD

	IF(IY.EQ.1.AND.IZ.EQ.1)THEN
	WRITE(1,342)X,Y,Z,PHI(IX,IY,IZ),RHO(IX,IY,IZ)
	ENDIF
342	FORMAT(2X,3(F12.3,2X),4X,E12.5,4X,E12.5)

100	CONTINUE

	NPL=2*NPTS+1
	NL=2*NPTS+1

	PRINT *,' DO YOU WANT TO STORE THE IMAGE OF THE DENSITY'

	PRINT *,' PROJECTED ON THE XY PLANE ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,1)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

	PRINT *,' PROJECTED ON THE YZ PLANE ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,2)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

	PRINT *,' PROJECTED ON THE ZX PLANE ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,3)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

	PRINT *,' DO YOU WANT TO STORE THE IMAGE OF THE POTENTIAL?'

	PRINT *,' 1: IN THE PLANE X,Y ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,1)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

	PRINT *,' 1: IN THE PLANE Y,Z ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,2)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

	PRINT *,' 1: IN THE PLANE Z,X ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,3)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF

C
999	CLOSE(1)
	PRINT *,' Output in "poten.log".'
	STOP
	END
C
C

	SUBROUTINE PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,IOPT)
C****************************************************
C	Subroutine PROJECT to project the density on one
C	of the planes XY (IOPT=1), YZ (IOPT=2) or ZX (IOPT=3)
C********************************************************
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 RHO(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*4 IMAGE(IDIM,IDIM)
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ
C
C  Projects the density on the plane XY
C
	IF(IOPT.EQ.1) THEN

	DO 101 IX=0,NPTS
	DO 101 IY=0,NPTS
	RRR=0.D0

	DO IZ=0,NPTS
	RRR=2*RHO(IX,IY,IZ)*PASZ+RRR
	END DO
	RRR=RRR*FACT

	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)

101	CONTINUE
	ENDIF

C
C  Projects the density on the plane YZ
C
	IF(IOPT.EQ.2) THEN

	DO 102 IX=0,NPTS
	DO 102 IY=0,NPTS
	RRR=0.D0

	DO IZ=1,NPTS
	RRR=2*RHO(IZ,IX,IY)*PASX+RRR
	END DO
	RRR=RRR*FACT

	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)

102	CONTINUE
	ENDIF

C
C  Projects the density on the plane ZX
C
	IF(IOPT.EQ.3) THEN

	DO 103 IX=0,NPTS
	DO 103 IY=0,NPTS
	RRR=0.D0

	DO IZ=0,NPTS
	RRR=2*RHO(IY,IZ,IX)*PASY+RRR
	END DO
	RRR=RRR*FACT

	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)

103	CONTINUE
	ENDIF

	RETURN
	END

	SUBROUTINE PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,IOPT)
C****************************************************
C	Subroutine PLANE to create an image of the potential
C	in the plane chosen by the user.
C	 If IOPT=1 plane X,Y
C	 If IOPT=2 plane Y,Z
C	 If IOPT=3 plane Z,X
C********************************************************
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 PHI(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*4 IMAGE(IDIM,IDIM)
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ

C Projection on the XY plane
C
	IF(IOPT.EQ.1)THEN
	IZ=0
	DO 101 IX=0,NPTS
	DO 101 IY=0,NPTS
	RRR=PHI(IX,IY,IZ)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
101	CONTINUE
	ENDIF

C Projection on the YZ plane
C
	IF(IOPT.EQ.2)THEN
	IZ=0
	DO 102 IX=0,NPTS
	DO 102 IY=0,NPTS
	RRR=PHI(IZ,IX,IY)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
102	CONTINUE
	ENDIF

C Projection on the XY plane
C
	IF(IOPT.EQ.3)THEN
	IZ=0
	DO 103 IX=0,NPTS
	DO 103 IY=0,NPTS
	RRR=PHI(IY,IZ,IX)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
103	CONTINUE
	ENDIF

	RETURN
	END
C********************************************************************
C         IFAIL=0
C         CALL D01GAF(RAD1,Y,NPTS,S(J),ERR(J),IFAIL)
C D01GAF integrates Y(RAD1). Method of GILL and MILLER
C Comp. Journal 15, pp. 80-83, 1972. An algorithm for the integration
C of unequally spaced data.
C 	RAD1: all distinct in either ascending or descending order
C	S(J) : (output) estimate of the integral
C	ERR(J) : (output) estimate of the uncertainty in the integral
C

ext: for


! { dg-do run }
!
! Test the fix for PR84115.
!
! Contributed by G Steinmetz  <gscfq@t-online.de>
!
  character(:), allocatable :: chr
  allocate (chr, source = "abc")
  call s(chr, "abc")
  chr = "mary had a little lamb"
  call s(chr, "mary had a little lamb")
  deallocate (chr)
contains
  subroutine s(x, carg)
    character(:), allocatable :: x
    character(*) :: carg
    associate (y => x)
      if (y .ne. carg) STOP 1
    end associate
  end
end