!*****************************************************************************************
!> author: Jacob Williams
!
!  Runge-Kutta integration.

    module rk_module

    use kind_module,       only: wp
    use numbers_module,    only: zero

    implicit none

    private

    type,abstract,public :: rk_class

        !! main integration class:

        integer :: n = 0  !! user specified number of variables
        procedure(deriv_func),pointer  :: f      => null()  !! user-specified derivative function
        procedure(report_func),pointer :: report => null()  !! user-specified report function
        procedure(event_func),pointer  :: g      => null()  !! event function (stop when this is zero)

        contains

        procedure                        :: initialize         !! initialize the class (set n,f, and report)
        procedure,non_overridable,public :: integrate          !! main integration routine
        procedure,non_overridable,public :: integrate_to_event !! integration with event finding
        procedure(step_func),deferred    :: step               !! the step routine for the rk method
        procedure,public                 :: destroy            !! destructor

    end type rk_class

    !extend the abstract class to create an RK4 method:
    ! [all we need to do is set the step function]
    type,extends(rk_class),public :: rk4_class
        !! 4th order Runge-Kutta method.
        contains
        procedure :: step => rk4
    end type rk4_class
    type,extends(rk_class),public :: rk8_10_class
        !! 8th order Runge-Kutta method.
        contains
        procedure :: step => rk8_10
    end type rk8_10_class

    interface

        subroutine deriv_func(me,t,x,xdot)  !! derivative function
        import :: rk_class,wp
        implicit none
            class(rk_class),intent(inout)        :: me
            real(wp),intent(in)                  :: t    !! time
            real(wp),dimension(me%n),intent(in)  :: x    !! state vector
            real(wp),dimension(me%n),intent(out) :: xdot !! derivative of state vector
        end subroutine deriv_func

        subroutine event_func(me,t,x,g)  !! event function
        import :: rk_class,wp
        implicit none
            class(rk_class),intent(inout)        :: me
            real(wp),intent(in)                  :: t !! time
            real(wp),dimension(me%n),intent(in)  :: x !! state vector
            real(wp),intent(out)                 :: g !! g(t,x). The goal is to stop the integration when g=0.
        end subroutine event_func

        subroutine report_func(me,t,x)  !! report function
        import :: rk_class,wp
        implicit none
            class(rk_class),intent(inout)        :: me
            real(wp),intent(in)                  :: t !! time
            real(wp),dimension(me%n),intent(in)  :: x !! state vector
        end subroutine report_func

        subroutine step_func(me,t,x,h,xf)   !! rk step function
        import :: rk_class,wp
        implicit none
            class(rk_class),intent(inout)        :: me
            real(wp),intent(in)                  :: t  !! initial time
            real(wp),dimension(me%n),intent(in)  :: x  !! initial state vector
            real(wp),intent(in)                  :: h  !! time step \( |\Delta t| \)
            real(wp),dimension(me%n),intent(out) :: xf !! final state vector
        end subroutine step_func

    end interface

    public :: rk_test !for testing

    contains
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize the [[rk_class]].

    subroutine initialize(me,n,f,report,g)

    implicit none

    class(rk_class),intent(inout)   :: me
    integer,intent(in)              :: n       !! number of variables
    procedure(deriv_func)           :: f       !! derivative function
    procedure(report_func),optional :: report  !! for reporting the steps
    procedure(event_func),optional  :: g       !! for stopping at an event

    call me%destroy()

    me%n = n
    me%f => f
    if (present(report)) me%report => report
    if (present(g))      me%g      => g

    end subroutine initialize
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[rk_class]].

    subroutine destroy(me)

    implicit none

    class(rk_class),intent(out)   :: me

    end subroutine destroy
!*****************************************************************************************

!*****************************************************************************************
!>
!  Main integration routine for the [[rk_class]].

    subroutine integrate(me,t0,x0,h,tf,xf)

    implicit none

    class(rk_class),intent(inout)        :: me
    real(wp),intent(in)                  :: t0    !! initial time
    real(wp),dimension(me%n),intent(in)  :: x0    !! initial state
    real(wp),intent(in)                  :: h     !! abs(time step)
    real(wp),intent(in)                  :: tf    !! final time
    real(wp),dimension(me%n),intent(out) :: xf    !! final state

    real(wp) :: t,dt,t2
    real(wp),dimension(me%n) :: x
    logical :: last,export

    if (.not. associated(me%f)) error stop 'Error in integrate: f is not associated.'

    export = associated(me%report)

    if (export) call me%report(t0,x0)  !first point

    if (h==zero) then
        xf = x0
    else

        t = t0
        x = x0
        dt = sign(h,tf-t0)  !time step (correct sign)
        do
            t2 = t + dt
            last = ((dt>=zero .and. t2>=tf) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tf))         !
            if (last) dt = tf-t                     !
            call me%step(t,x,dt,xf)
            if (last) exit
            if (export) call me%report(t2,xf)   !intermediate point
            x = xf
            t = t2
        end do

    end if

    if (export) call me%report(tf,xf)   !last point

    end subroutine integrate
!*****************************************************************************************

!*****************************************************************************************
!>
!  Event-finding integration routine for the [[rk_class]].
!  Integrates until g(t,x)=0, or until t=tf (whichever happens first).
!
!@note There are some efficiency improvements that could be made here.
!      This is a work in progress.

    subroutine integrate_to_event(me,t0,x0,h,tmax,tol,tf,xf,gf)

    use brent_module

    implicit none

    class(rk_class),intent(inout)        :: me
    real(wp),intent(in)                  :: t0      !! initial time
    real(wp),dimension(me%n),intent(in)  :: x0      !! initial state
    real(wp),intent(in)                  :: h       !! abs(time step)
    real(wp),intent(in)                  :: tmax    !! max final time if event not located
    real(wp),intent(in)                  :: tol     !! function tolerance for root finding
    real(wp),intent(out)                 :: tf      !! actual final time reached
    real(wp),dimension(me%n),intent(out) :: xf      !! final state (at tf)
    real(wp),intent(out)                 :: gf      !! g value at tf

    !local variables:
    real(wp) :: t,dt,t2,ga,gb,dt_root,dum
    real(wp),dimension(me%n) :: x,g_xf
    logical :: first,last,export
    procedure(report_func),pointer :: report
    type(brent_class) :: solver
    integer :: iflag

    if (.not. associated(me%f)) error stop 'Error in integrate_to_event: f is not associated.'
    if (.not. associated(me%g)) error stop 'Error in integrate_to_event: g is not associated.'
    if (h==zero) error stop 'Error in integrate_to_event: h must not be zero.'

    !If the points are being exported:
    export = associated(me%report)

    !first point:
    if (export) call me%report(t0,x0)

    if (t0==tmax) then
        xf = x0
        tf = t0
        call me%g(t0,x0,gf)
    else

        first = .true.
        t = t0
        x = x0
        call me%g(t0,x0,ga)     !evaluate event function
        dt = sign(h,tmax-t0)    !time step (correct sign)

        do

            t2 = t + dt
            last = ((dt>=zero .and. t2>=tmax) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tmax))         !
            if (last) then
                dt = tmax-t
                t2 = tmax
            end if
            call me%step(t,x,dt,xf)
            call me%g(t2,xf,gb)     !evaluate event function

            if (first .and. abs(ga)<=tol) then

                !we ignore a root at t0 after the first step
                if (abs(gb)<=tol) then !check this one since it could have landed on a root
                    gf = gb
                    tf = t2
                    exit
                else
                    if (last) then  !exiting without having found a root
                        tf = t2
                        gf = gb
                        exit
                    end if
                    if (export) call me%report(t2,xf)   !intermediate point
                    x = xf
                    t = t2
                    ga = gb
                end if

            elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt]

                !find the root:
                call solver%set_function(solver_func)
                call solver%find_zero(zero,dt,tol,dt_root,dum,iflag,ga,gb)
                t2 = t + dt_root
                gf = solver_func(solver,dt_root)
                tf = t2
                xf = g_xf !computed in the solver function
                exit

            else  !no root yet, continue

                if (last) then  !exiting without having found a root
                    tf = t2
                    gf = gb
                    exit
                end if
                if (export) call me%report(t2,xf)   !intermediate point
                x = xf
                t = t2
                ga = gb

            end if

            if (first) first = .false.

        end do

    end if

    if (export) call me%report(t2,xf)   !last point

    contains

        function solver_func(this,delt) result(g)

        !! root solver function. The input is the dt offset from time t.

        implicit none

        class(brent_class),intent(inout) :: this
        real(wp),intent(in) :: delt  !! from [0 to dt]
        real(wp) :: g

        !take a step from t to t+delt and evaluate g function:
        call me%step(t,x,delt,g_xf)
        call me%g(t+delt,g_xf,g)

        end function solver_func

    end subroutine integrate_to_event
!*****************************************************************************************

!*****************************************************************************************
!>
!  Take one Runge Kutta 4 integration step: `t -> t+h (x -> xf)`

    subroutine rk4(me,t,x,h,xf)

    implicit none

    class(rk4_class),intent(inout)       :: me
    real(wp),intent(in)                  :: t   !! initial time
    real(wp),dimension(me%n),intent(in)  :: x   !! initial state
    real(wp),intent(in)                  :: h   !! time step
    real(wp),dimension(me%n),intent(out) :: xf  !! state at time `t+h`

    !local variables:
    real(wp),dimension(me%n) :: f1,f2,f3,f4
    real(wp) :: h2

    !parameters:
    real(wp),parameter :: half = 0.5_wp
    real(wp),parameter :: six  = 6.0_wp

    h2 = half*h

    call me%f(t,x,f1)
    call me%f(t+h2,x+h2*f1,f2)
    call me%f(t+h2,x+h2*f2,f3)
    call me%f(t+h,x+h*f3,f4)

    xf = x + h*(f1+f2+f2+f3+f3+f4)/six

    end subroutine rk4
!*****************************************************************************************

!*****************************************************************************************
!>
!  Take one Runge Kutta 8 integration step: `t -> t+h (x -> xf)`
!  This is Formula (8-10) from Reference [1].
!
!# Reference
!  1. E. B. Shanks, "[Higher Order Approximations of Runge-Kutta Type](http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19650022581.pdf)",
!     NASA Technical Note, NASA TN D-2920, Sept. 1965.

    subroutine rk8_10(me,t,x,h,xf)

    implicit none

    class(rk8_10_class),intent(inout)    :: me
    real(wp),intent(in)                  :: t   !! initial time
    real(wp),dimension(me%n),intent(in)  :: x   !! initial state
    real(wp),intent(in)                  :: h   !! time step
    real(wp),dimension(me%n),intent(out) :: xf  !! state at time `t+h`

    !local variables:
    real(wp),dimension(me%n) :: f0,f1,f2,f3,f4,f5,f6,f7,f8,f9

    !parameters:
    real(wp),parameter :: a1  = 4.0_wp/27.0_wp
    real(wp),parameter :: a2  = 2.0_wp/9.0_wp
    real(wp),parameter :: a3  = 1.0_wp/3.0_wp
    real(wp),parameter :: a4  = 1.0_wp/2.0_wp
    real(wp),parameter :: a5  = 2.0_wp/3.0_wp
    real(wp),parameter :: a6  = 1.0_wp/6.0_wp
    real(wp),parameter :: a8  = 5.0_wp/6.0_wp
    real(wp),parameter :: c   = 1.0_wp/840.0_wp
    real(wp),parameter :: c0  = 41.0_wp
    real(wp),parameter :: c3  = 27.0_wp
    real(wp),parameter :: c4  = 272.0_wp
    real(wp),parameter :: c5  = 27.0_wp
    real(wp),parameter :: c6  = 216.0_wp
    real(wp),parameter :: c8  = 216.0_wp
    real(wp),parameter :: c9  = 41.0_wp
    real(wp),parameter :: aa1 = 4.0_wp/27.0_wp
    real(wp),parameter :: aa2 = 1.0_wp/18.0_wp
    real(wp),parameter :: aa3 = 1.0_wp/12.0_wp
    real(wp),parameter :: aa4 = 1.0_wp/8.0_wp
    real(wp),parameter :: aa5 = 1.0_wp/54.0_wp
    real(wp),parameter :: aa6 = 1.0_wp/4320.0_wp
    real(wp),parameter :: aa7 = 1.0_wp/20.0_wp
    real(wp),parameter :: aa8 = 1.0_wp/288.0_wp
    real(wp),parameter :: aa9 = 1.0_wp/820.0_wp
    real(wp),parameter :: b21 = 3.0_wp
    real(wp),parameter :: b32 = 3.0_wp
    real(wp),parameter :: b43 = 3.0_wp
    real(wp),parameter :: b50 = 13.0_wp
    real(wp),parameter :: b52 = -27.0_wp
    real(wp),parameter :: b53 = 42.0_wp
    real(wp),parameter :: b54 = 8.0_wp
    real(wp),parameter :: b60 = 389.0_wp
    real(wp),parameter :: b62 = -54.0_wp
    real(wp),parameter :: b63 = 966.0_wp
    real(wp),parameter :: b64 = -824.0_wp
    real(wp),parameter :: b65 = 243.0_wp
    real(wp),parameter :: b70 = -231.0_wp
    real(wp),parameter :: b72 = 81.0_wp
    real(wp),parameter :: b73 = -1164.0_wp
    real(wp),parameter :: b74 = 656.0_wp
    real(wp),parameter :: b75 = -122.0_wp
    real(wp),parameter :: b76 = 800.0_wp
    real(wp),parameter :: b80 = -127.0_wp
    real(wp),parameter :: b82 = 18.0_wp
    real(wp),parameter :: b83 = -678.0_wp
    real(wp),parameter :: b84 = 456.0_wp
    real(wp),parameter :: b85 = -9.0_wp
    real(wp),parameter :: b86 = 576.0_wp
    real(wp),parameter :: b87 = 4.0_wp
    real(wp),parameter :: b90 = 1481.0_wp
    real(wp),parameter :: b92 = -81.0_wp
    real(wp),parameter :: b93 = 7104.0_wp
    real(wp),parameter :: b94 = -3376.0_wp
    real(wp),parameter :: b95 = 72.0_wp
    real(wp),parameter :: b96 = -5040.0_wp
    real(wp),parameter :: b97 = -60.0_wp
    real(wp),parameter :: b98 = 720.0_wp

    call me%f(t,x,f0)
    call me%f(t+a1*h,x+aa1*h*f0,f1)
    call me%f(t+a2*h,x+aa2*h*(f0+b21*f1),f2)
    call me%f(t+a3*h,x+aa3*h*(f0+b32*f2),f3)
    call me%f(t+a4*h,x+aa4*h*(f0+b43*f3),f4)
    call me%f(t+a5*h,x+aa5*h*(b50*f0+b52*f2+b53*f3+b54*f4),f5)
    call me%f(t+a6*h,x+aa6*h*(b60*f0+b62*f2+b63*f3+b64*f4+b65*f5),f6)
    call me%f(t+h,x+aa7*h*(b70*f0+b72*f2+b73*f3+b74*f4+b75*f5+b76*f6),f7)
    call me%f(t+a8*h,x+aa8*h*(b80*f0+b82*f2+b83*f3+b84*f4+b85*f5+b86*f6+b87*f7),f8)
    call me%f(t+h,x+aa9*h*(b90*f0+b92*f2+b93*f3+b94*f4+b95*f5+b96*f6+b97*f7+b98*f8),f9)

    xf = x + h*c*(c0*f0+c3*f3+c4*f4+c5*f5+c6*f6+c8*f8+c9*f9)

    end subroutine rk8_10
!*****************************************************************************************

!*****************************************************************************************
!>
!  Unit test of the [[rk_module]].
!  Integrate a two-body orbit around the Earth.

    subroutine rk_test()

    type,extends(rk4_class) :: spacecraft
        !! spacecraft propagation type.
        !! extends the [[rk4_class]] to include data used in the deriv routine
        real(wp) :: mu = zero      !! central body gravitational parameter (km3/s2)
        integer :: fevals = 0      !! number of function evaluations
        logical :: first = .true.  !! first point is being exported
    end type spacecraft

    integer,parameter :: n=6    !! number of state variables
    real(wp),parameter :: tol = 1.0e-12_wp !! event location tolerance

    type(spacecraft) :: s, s2
    real(wp) :: t0,tf,x0(n),dt,xf(n),x02(n),gf,tf_actual

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' rk_test'
    write(*,*) '---------------'
    write(*,*) ''

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

    !constructor (main body is Earth):
    s = spacecraft(n=n,f=twobody,mu=398600.436233_wp,report=twobody_report)

    !initial conditions:
    x0 = [10000.0_wp,10000.0_wp,10000.0_wp,&   !initial state [r,v] (km,km/s)
            1.0_wp,2.0_wp,3.0_wp]
    t0 = zero       !initial time (sec)
    dt = 10.0_wp    !time step (sec)
    tf = 1000.0_wp  !final time (sec)

    s%fevals = 0
    s%first = .true.
    call s%integrate(t0,x0,dt,tf,xf)    !forward
    write(*,*) ''
    write(*,'(A/,*(F15.6/))') 'Final state:',xf

    s%fevals = 0
    s%report => null()    !disable reporting
    call s%integrate(tf,xf,-dt,t0,x02)  !backwards

    write(*,'(A/,*(E20.12/))') 'Error:',x02-x0
    write(*,'(A,I5)') 'Function evaluations:', s%fevals
    write(*,*) ''

    !***************************************************************************
    !event finding test:

    write(*,*) ' Event test - integrate until z = 12,000'
    s2 = spacecraft(n=n,f=twobody,g=twobody_event,mu=398600.436233_wp,report=twobody_report)
    x0 = [10000.0_wp,10000.0_wp,10000.0_wp,&   !initial state [r,v] (km,km/s)
            1.0_wp,2.0_wp,3.0_wp]
    t0 = zero       !initial time (sec)
    dt = 10.0_wp    !time step (sec)
    tf = 1000.0_wp  !final time (sec)
    call s2%integrate_to_event(t0,x0,dt,tf,tol,tf_actual,xf,gf)
    write(*,*) ''
    write(*,'(A/,*(F15.6/))') 'Final time: ',tf_actual
    write(*,'(A/,*(F15.6/))') 'Final state:',xf
    write(*,'(A/,*(F15.6/))') 'Event func :',gf

    contains
!*****************************************************************************************

    !*********************************************************
        subroutine twobody(me,t,x,xdot)

        !! derivative routine for two-body orbit propagation

        implicit none

        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x
        real(wp),dimension(me%n),intent(out) :: xdot

        real(wp),dimension(3) :: r,v,a_grav
        real(wp) :: rmag

        select type (me)
        class is (spacecraft)

            r = x(1:3)
            v = x(4:6)
            rmag = norm2(r)
            a_grav = -me%mu/rmag**3 * r !acceleration due to gravity

            xdot(1:3) = v
            xdot(4:6) = a_grav

            me%fevals = me%fevals + 1

        end select

        end subroutine twobody
    !*********************************************************

    !*********************************************************
        subroutine twobody_report(me,t,x)

        !! report function - write time,state to console

        implicit none

        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x

        select type (me)
        class is (spacecraft)
            if (me%first) then  !print header
                write(*,*) ''
                write(*,'(*(A15,1X))')  'time (sec)','x (km)','y (km)','z (km)',&
                                        'vx (km/s)','vy (km/s)','vz (km/s)'
                me%first = .false.
            end if
        end select

        write(*,'(*(F15.6,1X))') t,x

        end subroutine twobody_report
    !*********************************************************

    !*********************************************************
        subroutine twobody_event(me,t,x,g)

        !! event function (z = 12,000)

        implicit none

        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x
        real(wp),intent(out)                 :: g

        g = x(3) - 12000.0_wp

        end subroutine twobody_event
    !*********************************************************

    end subroutine rk_test
!*****************************************************************************************

!*****************************************************************************************
    end module rk_module
!*****************************************************************************************


text:
!!
!! @HEADER
!!
!!!!**********************************************************************
!!
!!  Zoltan Toolkit for Load-balancing, Partitioning, Ordering and Coloring
!!                  Copyright 2012 Sandia Corporation
!!
!! Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
!! the U.S. Government retains certain rights in this software.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are
!! met:
!!
!! 1. Redistributions of source code must retain the above copyright
!! notice, this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright
!! notice, this list of conditions and the following disclaimer in the
!! documentation and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the Corporation nor the names of the
!! contributors may be used to endorse or promote products derived from
!! this software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
!! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
!! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
!! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!!
!! Questions? Contact Karen Devine	kddevin@sandia.gov
!!                    Erik Boman	egboman@sandia.gov
!!
!!!!**********************************************************************
!!
!! @HEADER
 !!
! ****************************************************************************/
!/*
!   Code imported to Zoltan zdrive from
!
!   zoltanParams_read_file.c
!
!   Read Zoltan parameters from a file, call Zoltan to set the parameters
!
!   zoltanParams library
!
!   Jim Teresco
!
!   Department of Computer Science
!   Williams College
!
!   and
!
!   Computer Science Research Institute
!   Sandia National Laboratories
!
!*/

!  Translated to Fortran by Bill Mitchell, NIST, March 2007

module dr_param_file

!#include <stdio.h>
!#include <mpi.h>
!#include <stdlib.h>
!#include "zoltan.h"

use zoltan
use mpi_h

implicit none
private

! Are other routines public?  If so, add them separated by commas. Use
! an amperstand at the end of the line to continue on another line.

public ztnPrm_read_file

! standard error and standard out are not standardized in Fortran 90, but most
! compilers use units 0 (most) and 6 (universal).  Change here if necessary.

integer, parameter :: stderr = 0, stdout = 6

!/*#define DEBUG 1*/

logical, parameter :: DEBUG = .false.

!struct zoltanParams_list_entry {
!  char *param;
!  char *value;
!  struct zoltanParams_list_entry *next;
!};

! unknown length character strings are tricky in Fortran, just pick a
! length that is hopefully big enough.  Change it here if necessary.

integer, parameter :: MAX_CHAR_LEN = 128

type ztnPrm_list_entry
   character(len=MAX_CHAR_LEN) :: param
   character(len=MAX_CHAR_LEN) :: value
   type (ztnPrm_list_entry), pointer :: next
end type ztnPrm_list_entry

!struct zoltanParams_hier_struct {
!  int partition;
!  struct zoltanParams_list_entry *first;
!};

type ztnPrm_hier_struct
   integer :: partition
   type (ztnPrm_list_entry), pointer :: first
end type ztnPrm_hier_struct

!static struct zoltanParams_hier_struct **zph = NULL;
!static int num_levels = 0;
!static MPI_Comm comm;

! I think the use of **zph is as an allocatable array

type(ztnPrm_hier_struct), save, allocatable :: zph(:)
integer :: num_levels = 0
integer :: comm

contains

!static void check_level(int level) {
!
!  if (!zph) {
!    fprintf(stderr,"check_level: must set number of levels first\n");
!    return;
!  }
!
!  if (level >= num_levels) {
!    fprintf(stderr,"check_level: invalid level\n");
!  }
!}

subroutine check_level(level)
integer :: level

  if (.not. allocated(zph)) then
    write(stderr,*) "check_level: must set number of levels first"
    return
  endif

  if (level >= num_levels) then
    write(stderr,*) "check_level: invalid level"
  endif
end subroutine check_level

!void zoltanParams_hier_free() {
!  int i;
!
!  if (!zph) {
!    fprintf(stderr, "zoltanParams_hier_free warning: not allocated\n");
!    return;
!  }
!
!  for (i=0; i<num_levels; i++) {
!    free(zph[i]);
!  }
!
!  free(zph);
!}

subroutine ztnPrm_hier_free()
  integer :: i

  if (.not. allocated(zph)) then
    write(stderr,*) "ztnPrm_hier_free warning: not allocated"
    return
  endif

  deallocate(zph)
end subroutine ztnPrm_hier_free

!void zoltanParams_hier_set_num_levels(int levels) {
!  int i;

subroutine ztnPrm_hier_set_num_levels(levels)
  integer :: levels
  integer :: i, astat

!#ifdef DEBUG
!  printf("(zoltanParams_hier_set_num_levels) setting to %d\n", levels);
!#endif

  if (DEBUG) then
    write(stdout,*) "(ztnPrm_hier_set_num_levels) setting to ",levels
  endif

!  if (zph) {
!    fprintf(stderr,"zoltanParams_hier_set_num_levels warning: already initialized, reinitializing\n");
!    zoltanParams_hier_free();
!  }

  if (allocated(zph)) then
    write(stderr,*) "ztnPrm_hier_set_num_levels warning: already initialized, reinitializing"
    call ztnPrm_hier_free()
  endif

!  if (levels <= 0) {
!    fprintf(stderr, "(zoltanParams_hier_set_num_levels) num levels must be positive\n");
!    return;
!  }

  if (levels <= 0) then
    write(stderr,*) "(ztnPrm_hier_set_num_levels) num levels must be positive"
    return
  endif

!  num_levels = levels;
!
!  SAFE_MALLOC(zph, struct zoltanParams_hier_struct **,
!	      sizeof(struct zoltanParams_hier_struct *) * levels);
!
!  for (i=0; i<levels; i++) {
!    SAFE_MALLOC(zph[i],  struct zoltanParams_hier_struct *,
!		sizeof (struct zoltanParams_hier_struct));
!    zph[i]->partition = 0;
!    zph[i]->first = NULL;
!  }

  num_levels = levels

  allocate(zph(0:levels-1),stat=astat)
  if (astat /= 0) then
    write(stderr,*) "allocation failed in ztnPrm_hier_set_num_level"
    stop
  endif

  do i=0,levels-1
    zph(i)%partition = 0
    nullify(zph(i)%first)
  end do

!}

end subroutine ztnPrm_hier_set_num_levels

!void zoltanParams_hier_set_partition(int level, int partition) {
!
!#ifdef DEBUG
!  int mypid;
!  MPI_Comm_rank(comm, &mypid);
!
!  printf("[%d] will compute partition %d at level %d\n",
!	 mypid, partition, level);
!#endif
!
!  check_level(level);
!
!  zph[level]->partition = partition;
!}

subroutine ztnPrm_hier_set_partition(level,partition)
integer :: level, partition
integer :: mypid, ierr

  if (DEBUG) then
    call MPI_Comm_rank(comm,mypid,ierr)

    write(stdout,*) "[",mypid,"] will compute partition ",partition," at level ",level
  endif

  call check_level(level)

  zph(level)%partition = partition
end subroutine ztnPrm_hier_set_partition

!void zoltanParams_hier_set_param(int level, char *param, char *value) {
!  struct zoltanParams_list_entry *newparam, *nextparam;

subroutine ztnPrm_hier_set_param(level,param,value)
integer :: level
character(len=*) :: param, value
type(ztnPrm_list_entry), pointer :: newparam, nextparam
integer :: mypid, ierr, astat

!#ifdef DEBUG
!  int mypid;
!  MPI_Comm_rank(comm, &mypid);
!  printf("[%d] will set param <%s> to <%s> at level %d\n",
!	 mypid, param, value, level);
!#endif

  if (DEBUG) then
    call MPI_Comm_rank(comm,mypid,ierr)
    write(stdout,*) "[",mypid,"] will set param ",trim(param)," to ",trim(value)," at level ",level
  endif

!  check_level(level);
!
!  SAFE_MALLOC(newparam, struct zoltanParams_list_entry *,
!	      sizeof(struct zoltanParams_list_entry));

  call check_level(level)

  allocate(newparam,stat=astat)
  if (astat /= 0) then
    write(stderr,*) "allocation failed in ztnPrm_hier_set_param"
    stop
  endif

!  newparam->param = strdup(param);
!  newparam->value = strdup(value);
!  newparam->next = NULL;

  newparam%param = param
  newparam%value = value
  nullify(newparam%next)

!  if (!zph[level]->first) {
!    zph[level]->first = newparam;
!    return;
!  }

  if (.not. associated(zph(level)%first)) then
    zph(level)%first => newparam
    return
  endif

!  nextparam = zph[level]->first;
!  while (nextparam->next) nextparam=nextparam->next;
!  nextparam->next = newparam;

  nextparam => zph(level)%first
  do while (associated(nextparam%next))
    nextparam => nextparam%next
  end do
  nextparam%next => newparam

!}

end subroutine ztnPrm_hier_set_param

!int zoltanParams_hier_get_num_levels() {
!
!  return num_levels;
!}

function ztnPrm_hier_get_num_levels()
integer :: ztnPrm_hier_get_num_levels

  ztnPrm_hier_get_num_levels = num_levels
end function ztnPrm_hier_get_num_levels

!int zoltanParams_hier_get_part(int level) {
!
!  check_level(level);
!
!  return zph[level]->partition;
!}

function ztnPrm_hier_get_part(level)
integer :: level
integer :: ztnPrm_hier_get_part

  call check_level(level)

  ztnPrm_hier_get_part = zph(level)%partition
end function ztnPrm_hier_get_part

!void zoltanParams_hier_use_params(int level, struct Zoltan_Struct *zz, int *ierr) {
!  struct zoltanParams_list_entry *nextparam;
!
!  *ierr = ZOLTAN_OK;
!  check_level(level);
!
!  nextparam = zph[level]->first;
!
!  while (nextparam) {
!    *ierr = Zoltan_Set_Param(zz, nextparam->param, nextparam->value);
!    if (*ierr != ZOLTAN_OK) return;
!    nextparam = nextparam->next;
!  }
!
!}

subroutine ztnPrm_hier_use_params(level,zz,ierr)
integer :: level
type(Zoltan_Struct), pointer :: zz
integer :: ierr
type(ztnPrm_list_entry), pointer :: nextparam

  ierr = ZOLTAN_OK
  call check_level(level)

  nextparam => zph(level)%first

  do while (associated(nextparam))
    ierr = Zoltan_Set_Param(zz, nextparam%param, nextparam%value)
    if (ierr /= ZOLTAN_OK) return
    nextparam => nextparam%next
  end do
end subroutine ztnPrm_hier_use_params

!static int get_num_levels(void *data, int *ierr) {
!
!  *ierr = ZOLTAN_OK;
!  return zoltanParams_hier_get_num_levels();
!}

function get_num_levels(data, ierr)
integer(Zoltan_INT), intent(in) :: data(*)
integer(Zoltan_INT), intent(out) :: ierr
integer(Zoltan_INT) :: get_num_levels

  ierr = ZOLTAN_OK
  get_num_levels = ztnPrm_hier_get_num_levels()
end function get_num_levels

!static int get_part(void *data, int level, int *ierr) {
!
!  *ierr = ZOLTAN_OK;
!
!  return ztnPrm_hier_get_part(level);
!}

function get_part(data, level, ierr)
integer(Zoltan_INT), intent(in) :: data(*)
integer(Zoltan_INT), intent(in) :: level
integer(Zoltan_INT), intent(out) :: ierr
integer(Zoltan_INT) :: get_part

  ierr = ZOLTAN_OK

  get_part = ztnPrm_hier_get_part(level)
end function get_part

!static void get_method(void *data, int level, struct Zoltan_Struct *zz,
!		       int *ierr) {
!
!  zoltanParams_hier_use_params(level, zz, ierr);
!}

subroutine get_method(data,level,azz,ierr)
integer(Zoltan_INT), intent(in) :: data(*)
integer(Zoltan_INT), intent(in) :: level
type(Zoltan_Struct), intent(in), target :: azz
integer(Zoltan_INT), intent(out) :: ierr
type(Zoltan_Struct), pointer :: zz

  zz => azz

  call ztnPrm_hier_use_params(level, zz, ierr)
end subroutine get_method

!void zoltanParams_set_comm(MPI_Comm thecomm) {
!
!  /* remember the comm passed in */
!  MPI_Comm_dup(thecomm, &comm);
!}

subroutine ztnPrm_set_comm(thecomm)
integer :: thecomm
integer :: ierr

! remember the comm passed in

  call MPI_Comm_dup(thecomm, comm, ierr)
end subroutine ztnPrm_set_comm

!void zoltanParams_hier_setup(struct Zoltan_Struct *zz) {
!
!  /* make sure the hierarchical balancing callbacks are in place */
!  if (Zoltan_Set_Fn(zz, ZOLTAN_HIER_NUM_LEVELS_FN_TYPE,
!		    (void (*)()) get_num_levels, NULL) == ZOLTAN_FATAL) {
!    fprintf(stderr,"zoltanParams_hier_setup: set NUM_LEVELS callback failed\n");
!  }
!
!  if (Zoltan_Set_Fn(zz, ZOLTAN_HIER_PARTITION_FN_TYPE,
!		    (void (*)()) get_part, NULL) == ZOLTAN_FATAL) {
!    fprintf(stderr,"zoltanParams_hier_setup: set PARTITION callback failed\n");
!  }
!
!  if (Zoltan_Set_Fn(zz, ZOLTAN_HIER_METHOD_FN_TYPE,
!		    (void (*)()) get_method, NULL) == ZOLTAN_FATAL) {
!    fprintf(stderr,"zoltanParams_hier_setup: set METHOD callback failed\n");
!  }
!}

subroutine ztnPrm_hier_setup(zz)
type(Zoltan_Struct), pointer :: zz
integer(Zoltan_INT) :: dummy(1) = (/0/)

! make sure the hierarchical balancing callbacks are in place

  if (Zoltan_Set_Hier_Num_Levels_Fn(zz, get_num_levels, dummy) == &
      ZOLTAN_FATAL) then
    write(stderr,*) "ztnPrm_hier_setup: set NUM_LEVELS callback failed"
  endif

  if (Zoltan_Set_Hier_Part_Fn(zz, get_part, dummy) == &
      ZOLTAN_FATAL) then
    write(stderr,*) "ztnPrm_hier_setup: set PARTITION callback failed"
  endif

  if (Zoltan_Set_Hier_Method_Fn(zz, get_method, dummy) == &
      ZOLTAN_FATAL) then
    write(stderr,*) "ztnPrm_hier_setup: set METHOD callback failed"
  endif

end subroutine ztnPrm_hier_setup

!/*
!
!  zoltanParams_read_file
!
!  Set up the given Zoltan_Struct with parameters as specified
!  in the given file.
!
!  File format:
!
!  Lines of the format:
!  ZOLTAN_PARAM PARAM_VALUE
!
!  If the parameter is LB_METHOD set to HIER, the next part of the file
!  is interpreted as hierarchical balancing parameters:
!
!  num_levels
!  level 0 partitions for each proc
!  level 0 parameters
!  end with LEVEL END
!  level 1 partitions for each proc
!  level 1 parameters
!  end with LEVEL END
!  ...
!
!  End file with EOF
!
!*/

!void zoltanParams_read_file(struct Zoltan_Struct *lb, char *file,
!			    MPI_Comm thecomm) {
!  FILE *fp;
!  char str1[500], str2[500];
!  int numlevels, level, partition, proc;
!  int ierr;
!  int mypid, numprocs;

subroutine ztnPrm_read_file(lb, file, thecomm)
type(Zoltan_Struct), pointer :: lb
character(len=*) :: file
integer :: thecomm
integer :: fp
character(len=500) :: str1, str2
integer :: numlevels, level, proc
integer :: ierr
integer :: mypid, numprocs
logical :: not2
integer, allocatable :: partition(:)

!  /* remember the comm passed in */
!  MPI_Comm_dup(thecomm, &comm);
!
!  MPI_Comm_rank(comm, &mypid);
!  MPI_Comm_size(comm, &numprocs);

! remember the comm passed in
  call MPI_Comm_dup(thecomm, comm, ierr)

  call MPI_Comm_rank(comm, mypid, ierr)
  call MPI_Comm_size(comm, numprocs, ierr)

!  fp = fopen(file, "r");
!  if (!fp) {
!    fprintf(stderr,"Cannot open file %s for reading", file);
!    return;
!  }

! Assume unit 9 is available.  If it isn't, an error will be reported and
! you can change it to some other positive integer, not too big.

  fp = 9
  open(unit=fp,file=trim(file),action="read",iostat=ierr)
  if (ierr /= 0) then
    write(stderr,*) "cannot open file ",trim(file)," for reading"
    return
  endif

!#ifdef DEBUG
!  if (mypid == 0) {
!    printf("Reading Zoltan parameters from file %s\n", file);
!  }
!#endif

  if (DEBUG) then
    if (mypid == 0) then
      write(stdout,*) "Reading Zoltan parameters from file ",trim(file)
    endif
  endif

!  while (fscanf(fp, "%s %s\n", str1, str2) == 2) {

  do
    call myread(fp, str1, str2, not2)
    if (not2) exit

!    ierr = Zoltan_Set_Param(lb, str1, str2);
!    if (ierr != ZOLTAN_OK) {
!      fprintf(stderr,"Zoltan_Set_Param failed to set param <%s> to <%s>",str1,str2);
!    }
!#ifdef DEBUG
!    else {
!      if (mypid == 0) {
!	printf("Set Zoltan parameter <%s> to <%s>\n", str1, str2);
!      }
!    }
!#endif

    ! get rid of the leading space left on str2
    str2 = adjustl(str2)
    ierr = Zoltan_Set_Param(lb, trim(str1), trim(str2))
    if (ierr /= ZOLTAN_OK) then
      write(stderr,*) "Zoltan_Set_Param failed to set param ",trim(str1)," to ",trim(str2)
    endif
    if (DEBUG) then
      if (ierr == ZOLTAN_OK) then
        if (mypid == 0) then
          write(stdout,*) "Set Zoltan parameter ",trim(str1)," to ",trim(str2)
        endif
      endif
    endif

!    if (strcmp(str1,"LB_METHOD") == 0 && strcmp(str2,"HIER") == 0) {

    if (trim(str1) == "LB_METHOD" .and. trim(str2) == "HIER") then

!      zoltanParams_hier_setup(lb);

      call ztnPrm_hier_setup(lb)

!      /* the rest of the file contains hierarchical balancing parameters */
!      fscanf(fp, "%d", &numlevels);

! the rest of the file contains hierarchical balancing parameters
! The line containing numlevels is already in str1 (NO - it's next in the file)

      read(fp,*) numlevels

!#ifdef DEBUG
!      printf("[%d] read in numlevels=%d\n", mypid, numlevels);
!#endif

      if (DEBUG) then
        write(stdout,*) "[",mypid,"] read in numlevels=",numlevels
      endif

!      zoltanParams_hier_set_num_levels(numlevels);

      call ztnPrm_hier_set_num_levels(numlevels)

!      for (level=0; level<numlevels; level++) {
!	/* first, a list of partitions for each proc should be in the file */
!	for (proc=0; proc<numprocs; proc++) {
!	  fscanf(fp, "%d", &partition);
!	  if (proc == mypid) zoltanParams_hier_set_partition(level, partition);
!	}

      allocate(partition(0:numprocs-1))
      ! probably should check that allocate succeeded
      do level=0,numlevels-1
         read(fp,*) partition ! assumes the line has exactly numprocs numbers
         call ztnPrm_hier_set_partition(level,partition(mypid))

!	/* then parameters until we get LEVEL END */
!	while ((fscanf(fp, "%s %s\n", str1, str2) == 2) &&
!	       (strcmp(str1, "LEVEL") != 0) &&
!	       (strcmp(str2, "END") != 0)) {
!
!	  zoltanParams_hier_set_param(level, str1, str2);
!	}
!      }

! then parameters until we get LEVEL END
        do
          read(fp,*) str1, str2
          str2 = adjustl(str2)
          if (trim(str1) == "LEVEL" .and. trim(str2) == "END") exit
          call ztnPrm_hier_set_param(level, str1, str2)
        end do
      end do
      deallocate(partition)

!    }

    endif

!  }

  end do

!  fclose(fp);

  close(fp)

!}

end subroutine ztnPrm_read_file

! Fortran will generate an error if we try to read 2 strings and there is
! only 1 there.  So we have to read the whole line into a string and
! see if there are 1 or 2 strings in there.  Then read the individual
! strings and return them with a flag indicating if there are 1 or 2.

subroutine myread(runit,str1,str2,not2)
integer :: runit
character(len=*) :: str1, str2
logical :: not2
integer :: iostat

! assume 1000 is plenty long for an input line.

character(len=1000) :: line

! read the whole input line

  read(runit,"(A)",iostat=iostat) line

! end of file?
  if (iostat /= 0) then
     not2 = .true.
  else

! remove leading blanks

  line = adjustl(line)

! read the first string

  read(line,*) str1

! if the length of the whole line with leading and trailing blanks removed
! is the same as the length of the first string, then there is only 1 string

  if (len_trim(line) == len_trim(str1)) then
    not2 = .true.

! otherwise, read the second line

  else
    not2 = .false.
    read(line(len_trim(str1)+1:),"(A)") str2

  endif
  endif

end subroutine myread

end module dr_param_file

max_stars_repo_name: ParticulateFlow/zoltan

max_issues_repo_name: ParticulateFlow/zoltan

max_forks_repo_name: ParticulateFlow/zoltan


text:
program formatF_write

  implicit none
  integer          :: ios
  double precision :: x,y,z

  x =   3.14159d0
  y = -15.137d0
  z = 799.7432d0

  open( unit=1,                   &
        file="formatF_write.out",   &
        form="formatted",         &
        access="sequential",      &
        status="new",             &
        action="write",           &
        position="rewind",        &
        iostat=ios )

  if ( ios .ne. 0 ) then
    print*,"Cannot open file"
  else
   write( unit=1, fmt='(f7.5,f8.3,f9.4)' )  x,y,z
   write( unit=1, fmt='(f6.2,f9.4,f10.5)' ) x,y,z
  end if

  close( unit=1 )

end program formatF_write

max_stars_repo_name: annefou/Fortran


text:
!****************************************************************************
!
!  PURPOSE:  Module for global variables
!
!  VERSION: 0.1, 12-June-2012
!
!  LAST EDITED BY: Kurt
!
!  CHANGES SINCE PREV VERSION:
!****************************************************************************

module Globals
use params

implicit none

!GLOBALS FOR DIRECTORIES
character(len=100)  OutputDir
character(len=100)  OutputDirSims
character(len=100)  InputParamFile
integer :: AdjustCostModel
integer :: NoBankBelief
integer :: OnlyBankBelief
integer :: NoRentalBelief
integer :: nhr
integer :: DoIRF
integer :: DoRf
integer :: DoMqm
integer :: Modify
integer :: beliefshock
integer :: RationInc(ngpAgg)
double precision :: RationIncLevel,omegas
double precision ::Htexu,liquid,InterestCap
!Threadprivate counters
integer :: GeRent,HELOCVal,payextra
integer :: gij,giPh,giPhp,giW,giWp,giR,gif,giz,gizp,ghcc,gih,gil,gim,giAgg,giHD,giAy,giAggp,giCD,giExo,giExop,giC,giZh,girf
double precision :: gliq,gmort,ghouse,gloc,gPh,grf,grm,rf_base,grl,gPr
double precision    :: gvals(2)
integer             :: ginds(2)
integer             :: MtoMsim(nmsim)
integer             :: bcinds(2,nbsim)
double precision    :: adjust,kappa
double precision    :: defwedge
double precision    :: Goldmine(7)
double precision    :: Goldminevals(7)
double precision    :: bcvals(2,nbsim)
double precision    :: bet_h
double precision    :: rent_elas
double precision    :: adjval,hservices
double precision    :: MortDeduct,mortmarkup,helocmarkup
double precision    :: ltax,tautax
double precision    :: beqlev,beqmin
double precision    ::bmin,Ar
integer             :: exurent
double precision    :: scoef1(ngpPh),scoef2(ngpPh),scoef3(ngpPh)
double precision    :: bppcoef
double precision    :: mjprob,mjscale,Gini
double precision    :: HDX,HDY,HDZ
!Steady state simulations
double PRECISION, allocatable, dimension(:,:)::hsimss,msimss,lsimss,bsimss,hconssimss,csimss,beqsimss,psimss
INTEGER, allocatable, dimension(:,:)::hsimIss,buysimss,ownsimss,lsimIss,refisimss,foresimss,sellsimss
integer             :: iAggsteady
integer             :: DoSims
double precision    :: disutil(Jret)
double precision    :: GridBelief(ngpHD)
double precision    :: GridDemand(ngpNu)
double precision    :: GridAY(ngpAY)
double precision    :: Hsteady
double precision    :: GridRf(ngpRf)
double precision    :: GridDTI(ngpC)
double precision    :: Trans(ngpPh)
double precision    :: RfGrid(ngpAgg)
double precision    :: AdjGrid(ngpAgg)
double precision    :: AyGrid(ngpAgg)
integer :: DoInterestCap(ngpAgg)
!Rental cost conversion
double precision    :: phi_h

!Housing Depreciation
double precision    :: deltah,zbar

INTEGER             :: tstart,tend

double precision    :: initashare(ngpzy)
double precision    :: initadist(ngpzy,2)

!Mortgage Fixed Cost
double precision    :: MortFC

!Housing Maintenance Cost
double precision    :: HMRate

!Utility Cost of Foreclosure
double precision    :: epsprph,epsphy
double precision    :: xi

!Utility of homeownership

double precision    :: omega

double precision    :: GSElim
integer             :: GSE
!Economies of scale

double precision    :: phi(Jtot)
double precision    :: nchild(Jtot*2)
double precision    :: adults(Jtot*2)

! Risk aversion

double precision    :: sig

! CES Utility function params

double precision    :: alph
double precision    :: nu

! Discount Factor

double precision    :: bet

!Mortality Shock

double precision    :: delt(Jret)

!Moving Shock

double precision    :: thet(Jwork)

double precision    :: HELOCmax
double precision    :: LTVmax

double precision    :: LTVGrid(ngpAgg)
double precision    :: RmGrid(ngpAgg)
double precision    :: WedgeGrid(ngpAgg)
double precision    :: PhiHGrid(ngpAgg)
double precision    :: Markupgrid(ngpAgg)
double precision    :: RlGrid(ngpAgg)
double precision    :: LTIGrid(ngpAgg)
double precision    :: HELOCGrid(ngpAgg)
double precision    :: FCGrid(ngpAgg)
double precision    :: NuGrid(ngpAgg)
double precision    :: GridPhiH(ngpC)
double precision    :: ZhGrid(ngpAgg)

! DRS For Housing

double precision    :: alpha_h

!Grids for assets etc

double precision    :: GridB(nb)
double precision    :: GridBsim(nbsim)
double precision    :: GridMsim(nmsim)
double precision    :: GridBS(nbs)
double precision    :: GridM(nm)
double precision    :: GridL(nl)
double precision    :: GridH(nh)
double precision    :: PhGrid(ngpPh)
double precision    :: PrGrid(ngpPh,ngpAgg)
double precision    :: GridC(ngpC)
double precision    :: GridZh(ngpZh)

integer             :: WtoRI(ngpW)
integer             :: AtoY(ngpAgg)
integer             :: AtoR(ngpAgg)
integer             :: AtoD(ngpAgg)
integer             :: AtoC(ngpAgg)
integer             :: AtoZ(ngpAgg)
integer             :: CHZYtoA(ngpC,ngpHD,ngpZh,ngpAy)

!Coefficients for the forecast

double precision    :: aa0(ngpAgg,ngpAgg),aa1(ngpAgg,ngpAgg)
double precision    :: aa0p(ngpAgg,ngpAgg),aa1p(ngpAgg,ngpAgg)

double precision    :: PtoV1(ngpPh,ngpAgg,ngpAgg),PtoV2(ngpPh,ngpAgg,ngpAgg)
integer             :: PtoP1(ngpPh,ngpAgg,ngpAgg),PtoP2(ngpPh,ngpAgg,ngpAgg)
double precision :: EPhp(ngpPh,ngpAgg),Phpmat(ngpPh,ngpAgg,ngpAgg)
double precision,dimension(ngpPh)   ::XSHdem

! Mortgage interest rate
double precision    :: rm
double precision    :: rl
double precision    :: ru

double precision    :: rf

! Transition for Aggregate Price
double precision    :: transMatPh(ngpPh,ngpPh)
double precision    :: transMatAgg(ngpAgg,ngpAgg)

! Transition for Individual State
double precision    :: transMatC(ngpC,ngpC)

! Transition for Individual State
double precision    :: transMatRf(ngpRf,ngpRf)

! Transition for Individual State
double precision    :: transMatZh(ngpZh,ngpZh)

! Transition for Individual State
double precision    :: transMatHD(ngpHD,ngpHD)

! Transition for Individual State
double precision    :: transMatAy(ngpAy,ngpAy)

! Transition for Individual State
double precision    :: transMatW(ngpW,ngpW)

TYPE SimulationType
        integer, dimension(nsim)   :: jsim,psimI,eysimI,zysimI,buysim,sellsim,foresim,lsimI,ownsim,refisim,hsimI,diesimI,thetasim
        double precision, dimension(nsim)   :: ysim,ypssim,csim,hsim,hconssim,psim,lsim,taxsim,asim,msim,bsim,qlsim,qmsim,beqsim,hnwsim,delhnwsim,nwsim,delnwsim,mpcsim,csim2
END TYPE SimulationType
type(SimulationType)                    :: simtrans(1:500)    !-1 is previous steady state, Ttransition+1 is to store ass and dur decisions in simulation

double precision   :: Phpath(-Jtot:Tsim),AYpath(-Jtot:Tsim),HDpath(-Jtot:Tsim),Prpath(-Jtot:Tsim),Htpath(-Jtot:Tsim)
integer    :: iPhpath(-Jtot:Tsim),iAYpath(-Jtot:Tsim),iHDpath(-Jtot:Tsim),iAggpath(-Jtot:Tsim)
integer    :: DoExu

!Working Life Value and Policy Functions

double precision, dimension(nbsim,ngpW,Jtot) :: distssr,distssrexu,consdistr,dpconsdistr,hconsdistr,hdemdistr,bdistr,buydistr,dhconsdistr
double precision, allocatable, dimension(:,:,:,:,:,:) :: distsso,distssoexu,consdisto,dpconsdisto,hconsdisto,bdisto,buydisto,foredisto,refidisto,selldisto,dhconsdisto,hdemdisto
double precision, dimension(nbsim,nl,nmsim,nh,ngpW) :: initdistsso

double precision, dimension(nbsim,ngpW,Jtot) :: distr,distrp
double precision, dimension(nbsim,nl,nmsim,nh,ngpW,Jtot) :: disto,distop

double precision, allocatable   :: GridHR(:)
double precision, allocatable   ::    VownStay(:,:,:,:,:,:,:)

double precision, allocatable   ::    VownMove(:,:,:,:,:,:,:)

double precision, allocatable   ::    Vown(:,:,:,:,:,:,:)
double precision, allocatable   ::    EEVown(:,:,:,:,:,:,:)

integer, dimension(nb,ngpExo,ngpPh,Jwork):: hIBpol,lIBpol,mIBpol
integer, dimension(nb,ngpExo,ngpPh,Jret):: hIBpolRet,lIBpolRet,mIBpolRet
integer, dimension(nbs,ngpExo,ngpPh,Jwork):: hIBpolS,lIBpolS,mIBpolS
integer, dimension(nbs,ngpExo,ngpPh,Jret):: hIBpolRetS,lIBpolRetS,mIBpolRetS
double precision, allocatable, dimension(:,:,:,:)   ::  Wrent,consRpol,hconsRpol,bRpol,hRpol,mRpol
double precision, dimension(nbs,ngpExo,ngpPh,Jwork)   ::  WrentS,consRpolS,hconsRpolS,bRpolS,hRpolS,mRpolS
double precision, dimension(nb,ngpExo,ngpPh,Jwork):: Wbuy,consBpol,hconsBpol,bBpol,hBpol,mBpol,lBpol
double precision, dimension(nbs,ngpExo,ngpPh,Jwork):: WbuyS,consBpolS,hconsBpolS,bBpolS,hBpolS,mBpolS,lBpolS
double precision, dimension(nb,ngpExo,ngpPh,Jret)   ::  WrentRet,consRPolRet,hconsRPolRet,bRPolRet,hRPolRet,mRPolRet,WbuyRet,consBPolRet,hconsBPolRet,bBPolRet,hBPolRet,mBPolRet,lBPolRet
double precision, dimension(nbs,ngpExo,ngpPh,Jret)   ::  WrentRetS,consRPolRetS,hconsRPolRetS,bRPolRetS,hRPolRetS,mRPolRetS,WbuyRetS,consBPolRetS,hconsBPolRetS,bBPolRetS,hBPolRetS,mBPolRetS,lBPolRetS

integer, allocatable, dimension(:,:,:,:,:,:,:)::mINpol,lINpol,lIPpol
integer, allocatable, dimension(:,:,:,:,:,:,:)::mINpolRet,lINpolRet,lIPpolRet

double precision, allocatable, dimension(:,:,:,:,:,:,:):: Wsell,buySpol,consSpol,hconsSpol,bSpol,hSpol,mSpol
double precision, allocatable, dimension(:,:,:,:,:,:,:):: Wfore,consFpol,hconsFpol,bFpol,hFpol,mFpol
double precision, allocatable, dimension(:,:,:,:,:,:,:):: Wrefin,consNpol,hconsNpol,bNpol,hNpol,mNpol,lNpol
double precision, allocatable, dimension(:,:,:,:,:,:,:):: Wpay,consPpol,hconsPpol,bPpol,hPpol,mPpol,lPpol

double precision, allocatable, dimension(:,:,:,:,:,:,:):: WsellRet,buySPolRet,consSPolRet,hconsSPolRet,bSPolRet,hSPolRet,mSPolRet
double precision, allocatable, dimension(:,:,:,:,:,:,:):: WforeRet,consFPolRet,hconsFPolRet,bFPolRet,hFPolRet,mFPolRet
double precision, allocatable, dimension(:,:,:,:,:,:,:):: WrefinRet,consNPolRet,hconsNPolRet,bNPolRet,hNPolRet,mNPolRet,lNPolRet
double precision, allocatable, dimension(:,:,:,:,:,:,:):: WpayRet,consPPolRet,hconsPPolRet,bPPolRet,hPPolRet,mPPolRet,lPPolRet

double precision, dimension(nbsim,nl,nmsim,nh,ngpW)::bmatj,mmatj,hmatj,lmatj,ymatj
double precision, dimension(nbsim,nl,nmsim,nh,ngpW,Jtot)::bmat,mmat,hmat,lmat,ymat
double precision, dimension(nbsim,ngpW,Jtot)::bmatr,ymatr
double precision, dimension(nbsim,ngpW)::bmatjr,ymatjr

double precision, dimension(nb,ngpExo,ngpPh,Jwork)    :: Vrent,buyRentpol,consRentpol,hconsRentpol,bRentpol,hRentpol,mRentpol
double precision, dimension(nb,Jwork)    :: EVrent
double precision, dimension(nb,Jwork,ngpExo,ngpPh)    :: EEVrent
double precision, dimension(nb,nl,nm,nh,Jwork)        :: EVown
double precision, dimension(nb,nl,nm,nh,Jret)         :: EVownRet

!Retired Life Value and Policy Functions

double precision, dimension(nb,ngpExo,ngpPh,Jret)    :: VrentRet
double precision, dimension(nb,Jret)    :: EVrentRet
double precision, dimension(nb,Jret,ngpExo,ngpPh)    :: EEVrentRet

double precision, allocatable   ::    VownRet(:,:,:,:,:,:,:)
double precision, allocatable   ::    EEVownRet(:,:,:,:,:,:,:)

!Mortgage Functions

double precision, allocatable   ::  qmret(:,:,:,:,:,:,:)
double precision, allocatable   ::  qm(:,:,:,:,:,:,:)
double precision, allocatable   ::  qmret_e(:,:,:,:,:,:,:)
double precision, allocatable   ::  qmret_e_exo(:,:,:,:,:,:,:)
double precision, allocatable   ::  qm_e(:,:,:,:,:,:,:)
double precision, allocatable   ::  qm_e_exo(:,:,:,:,:,:,:)
double precision, allocatable   ::  ql(:,:,:,:,:,:,:)
double precision, allocatable   ::  qlret(:,:,:,:,:,:,:)
double precision, allocatable   ::  bankgain_mort(:,:,:,:,:,:)
double precision, allocatable   ::  bankgain_loc(:,:,:,:,:,:)
double precision, allocatable   ::  bankgain_mort_move(:,:,:,:,:,:)
double precision, allocatable   ::  bankgain_loc_move(:,:,:,:,:,:)

integer,dimension(nsim)         :: psisimI
integer,dimension(nsim,Jtot)    :: zysimI,thetasim
integer,dimension(psinum,ngpzy)  :: Wind

!GLOBALS FOR GRIDS

double precision, allocatable              	:: fylevgrid(:)        !fixed effect component levels, earnings
double precision, allocatable              	:: fyslogrid(:)        !fixed effect component slopes, earnings

double precision, dimension(Jtot,ngpp)       	:: pgrid,holdpgrid        !pensions grid

double precision, dimension(ngpey)              		:: eygrid        !transitory component, earnings
double precision, dimension(Jwork,ngpzy)             :: zygrid        !permanent component, earnings

double precision, dimension(Jwork,ngpfy,ngpzy,ngpey,ngpAy) :: ygrid 		!earnings
double precision, dimension(Jwork,ngpfy,ngpzy,ngpAy) 		:: ypsgrid 		!persistent component of earnings

!GLOBALS FOR INDICATORS
integer, dimension(ngpW)			:: zWind,fWind
integer, dimension(ngpR)			:: pRind
integer, dimension(ngpp)			:: Rind
integer, dimension(ngpfy)					:: fysloind,fylevind
integer, dimension(ngpfylev,ngpfyslo)		:: fyind

!GLOBALS FOR DISTRIBUTIONS AND TRANSITION MAPRICES
double precision, allocatable		:: initliq(:)
double precision, allocatable		:: initill(:)
double precision, allocatable		:: initnw(:)
double precision, dimension(:,:), allocatable		:: initliqdist
double precision, dimension(:,:), allocatable		:: initnwdist
double precision, dimension(:,:), allocatable		:: initilldist
double precision, allocatable		    :: fylevdist(:)
double precision, allocatable		    :: fyslodist(:)

double precision, dimension(ngpey)		    :: eydist
double precision, dimension(Jwork,ngpzy,ngpzy)	:: zytrans
double precision, dimension(Jwork,ngpzy)			:: zydist

double precision, dimension(Jwork)		    :: kappay
double precision, dimension(Jwork)			:: Vzy

double precision, dimension(psinum)    :: psidist
integer                                 :: inithouse

!PARAMETER GLOBALS
double precision     :: Vfey,Vepsy,Vetay,Vzy0,rhoy,lambday,Vfem,Vepsm,Vetam,rhom,agem,age2m,targetUsualKY,Vphy,phvarygrowth,targetMeanIlliquid,targetMedianIlliquid,targetP75Illiquid,rholevslo
double precision		:: inflation,nomretliq,nomretill,nomborrow,qfree,qill,qborr,taxurcg,mturcg,nomretgov,qgov,taxbase,phaseout,phaseoutamt,servflow,ndweight
double precision     :: gamC,sigC,borrowlim,bondparam,ratematch,debtmatch,maxdeduction,maxdedret
double precision		:: ctax,liqtax,illtax
double precision		:: surprob(Jret),popsize(Jtot),popdist(Jtot)
integer		:: posinitwealth,sameinitwealth,sameborrowlim

double precision, dimension(:,:), allocatable	:: uerand,zyrand,eyrand,dierand,initliqrand,initillrand,thetarand
double precision     :: poprand(nsim)

!GLOBALS TO STORE SIMULATION RESULTS
double precision, dimension(nsim, Jwork)	    :: ysim,yavsim,yavsimforpen,ypssim,eysim,zysim,usim

integer, dimension(nsim)	    	:: fysimI,fyslosimI,fylevsimI,fysimtransI
integer, dimension(nsim, Jtot)	    :: psimI
integer, dimension(nsim, Jtot+1)	:: dedsimI,mpdedsimI
integer, dimension(nsim, Jwork)	    :: eysimI

double precision :: yavall

double precision :: sshdemand
double precision :: ssph(ngpAgg)

double precision, dimension(Jwork,ngpfy)	        :: nbl		! natural borrowing limits (liquid assets)
double precision, dimension(Jwork,ngpfy,ngpzy)	    :: abl		! actual borrowing limits

double precision, dimension(ngpd)					:: dgrid		!deductions grid
double precision, dimension(ngpdret)					:: dgridret		!deductions grid

double precision :: ph_pers,ph_std,ay_std
integer          :: ph_nstd
integer          :: demmodel
double precision :: nwyperc(5),ltvperc(5),nwperc(5),hnwperc(5),liqperc(5),hnwshareperc(5)
!Lifecycle Stats Stuff

real(8), dimension(Jwork)           :: Eyav,Vloglabinc,Eunemp,Elabinc,Euiben
real(8), dimension(Jtot)            :: Econ,Etotcon,Eflow,Vlogtotcon,Vlogflow,Etax,Vlogcon,Vlogcono,Vlogconr,Eadj,Empc,Empadj,Enetlabinc,Eliqpl,CORRil,Ehcon,Ehouse,Emort,Emortpos,Elev,Eass,Eown,Efore,Esell,Epen,Einc,Vloginc,Vloghcon

real(8), dimension(Tsim)     :: pathEage,pathEinc,pathEcon,pathEass,pathEliq,pathEtax,pathEasspos,pathEloc,pathEhouse,pathEmort,pathEmortpos,pathEbeq,pathEbeqpos
real(8), dimension(Tsim)     :: pathEtranpos,pathEpen,pathElocpos,pathEhpos,pathEbuy,pathEmlpos,pathEhnw,pathEhcon
real(8), dimension(Tsim,Jtot)    :: agepathEage,agepathEinc,agepathEcon,agepathEsav,agepathEass,agepathEliq,agepathEtax,agepathEasspos,agepathEliqpos,agepathEdh,agepathEdp
real(8), dimension(Tsim,Jtot)    :: agepathEtranpos,agepathEpen,agepathEincown,agepathEincrent,agepathEnw,AgepathEnwown,agepathEnwrent,agepathEhnw,agepathErent,agepathEhnwshare
real(8), dimension(Tsim,Jtot)    :: agepathEconown,agepathEconrent,agepathEconinclow,agepathEconinchigh,agepathEconlowltv,agepathEconhighltv,agepathEhcon

real(8), dimension(Tsim)     :: pathEown,pathErefi,pathEsell,pathEfore,pathEforem,pathElev
real(8), dimension(Tsim,Jtot)     :: agepathEown,agepathErefi,agepathEsell,agepathEfore,agepathEforem,agepathElev,agepathEloc,agepathEhouse,agepathEmort,agepathEbuy,agepathEmortpos,agepathEbeq
real(8), dimension(Tsim,Jtot)     :: agepathElocpos,agepathEmlpos,agepathEhnwsim,agepathEdelhnwsim
double precision, dimension(Tsim)            :: phrand,aggrand

!THREADPRIVATE GLOBALS
!$OMP THREADPRIVATE(gij,giPh,giPhp,giW,giWp,giR,gliq,gmort,ghouse,gloc,gvals,ginds,gif,giz,ghcc,gih,gil,gim,giAgg,giHD,giAy,giAggp,gPh,EVown,EVownRet,EVRent,EVRentRet,giExo,giExop,giC,giZh,giRf,grf,grm,grl,gPr)

end module globals

max_stars_repo_name: lnsongxf/housing-boom-bust

max_issues_repo_name: floswald/housing-boom-bust

max_forks_repo_name: kurtmitman/housing-boom-bust


text:
c$Header: /data/petsun4/data1/src_solaris/imglin/RCS/stretchout.f,v 1.1 2010/03/03 01:28:53 avi Exp $
c$Log: stretchout.f,v $
c Revision 1.1  2010/03/03  01:28:53  avi
c Initial revision
c
      subroutine stretchout(t)
      real*4 t(4,4)
      real*4 g(3,3),w(3,3),q(3,3),wt(3,3),a(3,3),b(3,3)

      do 1 i=1,3
      do 1 j=1,3
    1 g(i,j)=t(i,j)
      call transpos(g,b,3)
      call matmul(b,g,q,3)
      call eigen(q,w,3)
      call transpos(w,wt,3)
      do 3 i=1,3
    3 q(i,i)=1./sqrt(q(i,i))
      call matmul(w,q,a,3)
      call matmul(a,wt,b,3)
      call matmul(g,b,q,3)
      do 4 i=1,3
      do 4 j=1,3
    4 t(i,j)=q(i,j)

      return
      end

ext: f

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
Address(Bodega Bay Place) is a residential Culdesacs culdesac in West Davis.

 Intersecting Streets

  Westshore Street

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! { dg-do run }
! Test passing character strings by-value.
! PR fortran/32732
program test
  implicit none
  character(len=13) :: chr
  chr =  'Fortran       '
  call sub1(chr)
  if(chr /= 'Fortran       ') call abort()
contains
  subroutine sub1(a)
    character(len=13), VALUE :: a
    a = trim(a)//" rules"
    call sub2(a)
  end subroutine sub1
  subroutine sub2(a)
    character(len=13), VALUE :: a
    print *, a
    if(a /= 'Fortran rules') call abort()
  end subroutine sub2
end program test

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


!>  RK stepper
module pf_my_stepper
  use pf_mod_dtype
  use pf_mod_zndarray
  use pf_mod_rkstepper
  use pf_mod_fftpackage
  use pf_mod_solutions
  use pf_mod_zutils
  use pf_mod_fftops
  implicit none

  ! Define the derived stepper type
  type, extends(pf_ark_stepper_t) :: my_stepper_t
     integer ::     nx
     ! fft object and differentiaton matrices
     type(pf_fft_t), pointer :: fft_tool
     type(pf_fft_ops_t), pointer :: fft_ops
     complex(pfdp), allocatable :: tmp(:)                ! Temp space for feval
     complex(pfdp), pointer :: p_tmp(:)   ! a useful pointer

     !  Useful for making routines  generic by dimension
     class(pf_zndarray_t), pointer :: f_encap,rhs_encap,y_encap
     complex(pfdp),        pointer :: yvec(:), rhsvec(:), fvec(:)

   contains

     procedure :: f_eval       !  Computes the advection and diffusion terms
     procedure :: f_comp       !  Computes the advection and diffusion terms
     procedure :: initialize
     procedure :: initialize_tmp
     procedure :: destroy

  end type my_stepper_t

contains
  subroutine initialize_tmp(this)
    class(my_stepper_t), intent(inout) :: this

    integer :: istat
    allocate(this%tmp(this%nx),STAT=istat)
    if (istat .ne. 0)  call pf_stop(__FILE__,__LINE__,'Allocate failed ',istat)

  end subroutine initialize_tmp

  ! The rest of the stuff is dimension independent
  include 'stepper_include.f90'

end module pf_my_stepper


text:
Subroutine zpca1a(i)
  Implicit Double Precision (A-H, O-Z)
  Parameter (maxptn=400001)
  Common /para2/xmp, xmu, alpha, rscut2, cutof2
  Common /para5/iconfg, iordsc
  Common /prec2/gx(maxptn), gy(maxptn), gz(maxptn), ft(maxptn), px(maxptn), py(maxptn), pz(maxptn), e(maxptn), xmass(maxptn), ityp(maxptn)
  Common /prec3/gxs(maxptn), gys(maxptn), gzs(maxptn), fts(maxptn), pxs(maxptn), pys(maxptn), pzs(maxptn), es(maxptn), xmasss(maxptn), ityps(maxptn)
  Common /prec5/eta(maxptn), rap(maxptn), tau(maxptn)
  Common /prec6/etas(maxptn), raps(maxptn), taus(maxptn)
  Common /ana1/ts(12)
  Save
  If (iconfg==1) Then
     t1 = fts(i)
     t2 = ft(i)
     ipic = 11
  Else If (iconfg==2 .Or. iconfg==3) Then
     t1 = taus(i)
     t2 = tau(i)
     ipic = 12
  Else If (iconfg==4 .Or. iconfg==5) Then
     t1 = fts(i)
     t2 = ft(i)
     ipic = 12
  End If
  If (iconfg<=3) Then
     Do ian = 1, ipic
        If (t1<=ts(ian) .And. t2>ts(ian)) Then
           rapi = raps(i)
           et = dsqrt(pxs(i)**2+pys(i)**2+xmp**2)
           Call zpca1b(rapi, et, ian)
        End If
     End Do
  Else
     Do ian = 1, ipic
        If (t1<=ts(ian) .And. t2>ts(ian)) Then
           p0 = es(i)
           p1 = pxs(i)
           p2 = pys(i)
           p3 = pzs(i)
           Call zpca1c(p0, p1, p2, p3, ian)
        End If
     End Do
  End If
  Return
End Subroutine zpca1a

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


!***********************************/
!	Name:外部境界のうち流入流速との内積が0以下になる境界を無反射境界に設定しなおす
!	Alias:UPrepareBoundary
!	Description:
!	Type:CellCenter
!	Input:Geometry%Dimension,Geom%CellNumber, CellCenter%ConservedQuantity
!	Output:CellCenter%ConservedQuantity
!	Note:
!	Author:Akitaka Toyota
!	Date:2017.11.11
!	Update:2017.11.11
!	Other:
!***********************************/
subroutine UPrepareBoundary(UG,UCC)
use StructVar_Mod
use LoopVar_Mod
use ConstantVar_Mod, Gamma => SpecificOfHeatRatio
implicit none
    type(UnstructuredGrid), intent(inout) :: UG
    type(CellCenter), intent(inout) :: UCC

    do iCell = UG%GI%RealCells+1, UG%GI%AllCells
        if(UG%GM%CellType(iCell,1,1) == 1) then !Inflow and Outflow
            iEdge = UG%VC%Edge(iCell)
            if(dot_product(-UG%GM%BC%InFlowVariable(2:4), -UG%GM%Normal(iEdge, :)) <= 0.0d0) then    !　流入流速との内積が負になるとき流入境界→無反射境界に変更
                UG%GM%CellType(iCell, 1, 1) = 10   ! non reflect boundary
            end if
        end if
    end do

    return
end subroutine UPrepareBoundary


module variables_mod

  implicit none

  integer, parameter :: &
       NF  = 2, & ! Number of fluids
       NFF = 3, & ! Number of fluid pairs, = NF (NF + 1) / 2
       ND  = 1, & ! Number of spatial dimensions
       NV  = 2, & ! Number of variables, = ND + 1
       NEVOLVE = 4, & ! Number of evolved variables, = NF * NV
       LWA = 33   ! Size of the work array for minpack; depends on NFF as
                  ! LWA >= (NFF (3 NFF + 13)) / 2

  integer, parameter :: &
       N_UP_FILE = 10, &
       N_DOWN_FILE = 11, &
       MU_UP_FILE = 12, &
       MU_DOWN_FILE = 13, &
       N_SQD_FILE = 14, &
       V_UP_FILE = 15, &
       C_CC_FILE = 16, &
       C_ENT_FILE = 17, &
       CS2_FILE = 18, &
       CONSTRAINTS_FILE = 19, &
       FFT_FILE = 20

  integer, parameter :: &
       M_FN_MULTI_POLY = 1, &
       M_FN_TWO_STREAM = 2, &
       M_FN_COSMOLOGICAL_POLY = 3, &
       M_FN_GENERIC_PAPER = 4

end module variables_mod


text:
C-----------------------------------------------------------------------
      SUBROUTINE SPEPS(I,M,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    SPEPS       COMPUTE UTILITY SPECTRAL FIELDS
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 92-10-31
C
C ABSTRACT: COMPUTES CONSTANT FIELDS INDEXED IN THE SPECTRAL DOMAIN
C           IN "IBM ORDER" (ZONAL WAVENUMBER IS THE SLOWER INDEX).
C           IF L IS THE ZONAL WAVENUMBER AND N IS THE TOTAL WAVENUMBER
C           AND A IS THE EARTH RADIUS, THEN THE FIELDS RETURNED ARE:
C           (1) NORMALIZING FACTOR EPSILON=SQRT((N**2-L**2)/(4*N**2-1))
C           (2) LAPLACIAN FACTOR N*(N+1)/A**2
C           (3) ZONAL DERIVATIVE/LAPLACIAN FACTOR L/(N*(N+1))*A
C           (4) MERIDIONAL DERIVATIVE/LAPLACIAN FACTOR EPSILON/N*A
C
C PROGRAM HISTORY LOG:
C   91-10-31  MARK IREDELL
C
C USAGE:    CALL SPEPS(I,M,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
C
C   INPUT ARGUMENT LIST:
C     I        - INTEGER SPECTRAL DOMAIN SHAPE
C                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C     M        - INTEGER SPECTRAL TRUNCATION
C
C   OUTPUT ARGUMENT LIST:
C     EPS      - REAL ((M+1)*((I+1)*M+2)/2) SQRT((N**2-L**2)/(4*N**2-1))
C     EPSTOP   - REAL (M+1) SQRT((N**2-L**2)/(4*N**2-1)) OVER TOP
C     ENN1     - REAL ((M+1)*((I+1)*M+2)/2) N*(N+1)/A**2
C     ELONN1   - REAL ((M+1)*((I+1)*M+2)/2) L/(N*(N+1))*A
C     EON      - REAL ((M+1)*((I+1)*M+2)/2) EPSILON/N*A
C     EONTOP   - REAL (M+1) EPSILON/N*A OVER TOP
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
      REAL EPS((M+1)*((I+1)*M+2)/2),EPSTOP(M+1)
      REAL ENN1((M+1)*((I+1)*M+2)/2),ELONN1((M+1)*((I+1)*M+2)/2)
      REAL EON((M+1)*((I+1)*M+2)/2),EONTOP(M+1)
      PARAMETER(RERTH=6.3712E6,RA2=1./RERTH**2)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DO L=0,M
        K=L*(2*M+(I-1)*(L-1))/2+L+1
        EPS(K)=0.
        ENN1(K)=RA2*L*(L+1)
        ELONN1(K)=RERTH/(L+1)
        EON(K)=0.
      ENDDO
      DO L=0,M
        DO N=L+1,I*L+M
          K=L*(2*M+(I-1)*(L-1))/2+N+1
          EPS(K)=SQRT(FLOAT(N**2-L**2)/FLOAT(4*N**2-1))
          ENN1(K)=RA2*N*(N+1)
          ELONN1(K)=RERTH*L/(N*(N+1))
          EON(K)=RERTH/N*EPS(K)
        ENDDO
      ENDDO
      DO L=0,M
        N=I*L+M+1
        EPSTOP(L+1)=SQRT(FLOAT(N**2-L**2)/FLOAT(4*N**2-1))
        EONTOP(L+1)=RERTH/N*EPSTOP(L+1)
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack


text:
module spectra

    use main
    use crosssection
    use transition
    use constants, only: kev, mN, kilometerpersecond
    use eventrate, only: deventrate
    use quadrature, only: boole
    use settings
    use orbitals, only: bfm
    use nucresponse, only: nucFormFactor, nucFormFactor_transform
    use densities, only: pndens

    implicit none

    real(kind=8) :: x_start
    real(kind=8) :: x_stop
    real(kind=8) :: x_step
    real(kind=8), allocatable :: x_grid(:), energy_grid(:), momentum_grid(:)
    integer :: energy_grid_size

    contains

    function velocitycurve(vlist, q, option)
        implicit none
        real(kind=8) :: q, v
        real(kind=8), dimension(:) :: vlist
        integer :: option

        real(kind=8), dimension(size(vlist)) :: velocitycurve

        integer :: i

        select case(option)
        case(3)
            do i = 1, size(vlist)
                v = vlist(i)
                velocitycurve(i) = transition_probability(v, q)
            end do
        case(2)
            do i = 1, size(vlist)
                v = vlist(i)
                velocitycurve(i) = diffCrossSection(v, q)
            end do
        case default
            stop "Not a velocity curve option."
        end select
    end function velocitycurve

    subroutine velocity_curve(option)
        implicit none

        real(kind=8) :: Er, Qr, vstart, vstop, vstep
        integer :: sizevlist, i
        real(kind=8), allocatable :: vlist(:), cslist(:)
        integer :: option, iunit

        print*,"Enter recoil E (keV):"
        read*,Er
        print '("Er = ",F8.4," (keV)")',Er
        Qr = sqrt(2d0*nuc_target%mass*mN*er*kev)
        print '("Qr = ",F8.4,"(GeV/c)")',qr
        print '("Vmin = ",F8.4,"(km/s)")',1/kilometerpersecond&
            *qr/(2.0* wimp%mass * nuc_target%mass * mN / (wimp%mass + nuc_target%mass * mN))

        print*,"Enter start v (km/s):"
        read*,vstart
        print '("v start = ",F8.4," (km/s)")',vstart
        vstart = vstart * kilometerpersecond
        print '("v start = ",ES12.4," (c)")',vstart

        print*,"Enter stop v (km/s):"
        read*,vstop
        print '("v stop = ",F12.4," (km/s)")',vstop
        vstop = vstop * kilometerpersecond
        print '("v stop = ",ES12.4," (c)")',vstop

        print*,"Enter step v (km/s):"
        read*,vstep
        print '("v step = ",F8.4," (km/s)")',vstep
        vstep = vstep * kilometerpersecond
        print '("v step = ",ES12.4," (c)")',vstep

        sizevlist = int(abs(vstop - vstart)/vstep)
        allocate(vlist(sizevlist))
        allocate(cslist(sizevlist))
        do i = 1, sizevlist
            vlist(i) = vstart + (i-1) * vstep
        end do

        cslist = velocitycurve(vlist, qr, option)

        select case(option)
        case(3)
            open(newunit=iunit, file="transition_probability.dat")
            do i = 1, sizevlist
                write(iunit,*) vlist(i)/kilometerpersecond, cslist(i)
            end do
            close(iunit)
        case(2)
            open(newunit=iunit, file="crosssection.dat")
            do i = 1, sizevlist
                write(iunit,*) vlist(i)/kilometerpersecond, cslist(i)
            end do
            close(iunit)
        case default
            STOP "Not a velocity curve option."
        end select

    end subroutine

    subroutine eventrate_spectra

        implicit none

        integer :: calc_num
        integer :: iunit
        integer :: i, N
        real(kind=8) :: q, dq

        real(kind=8) :: recoil_energy, momentum_transfer
        real(kind=8), allocatable :: event_rate_spectra(:)
        real(kind=8) :: mtarget, totaleventrate

        print*,"Computing differential event rate spectra"

        ! Get parameters
        mtarget = nuc_target%mass

        ! Get domain
        call get_energy_grid(mtarget)

        ! Setup calculation and compute
        print*,'Number of event rates to compute:',energy_grid_size
        allocate(event_rate_spectra(energy_grid_size))

        N = size(momentum_grid)
        !$OMP parallel do private(q) schedule(dynamic, 1)
        do i = 1, N
            q = momentum_grid(i)
            Event_rate_spectra(i) = dEventRate(q)
        end do
        !$OMP end parallel do
        !$OMP barrier

        if (.not.useenergyfile) then
            dq = abs(momentum_grid(2) - momentum_grid(1))
            call boole(energy_grid_size, momentum_grid*Event_rate_spectra/(mn*mtarget),&
                dq, totaleventrate)
            print*,'Total integrated eventrate (events)',totaleventrate
        end if

        write(6,"(A,T30,A,T56,A)")" E-recoil (kev)","q-transfer (gev/c)","Eventrate (events/gev)"
        do calc_num = 1, energy_grid_size
            momentum_transfer = momentum_grid(calc_num)
            recoil_energy = energy_grid(calc_num)
            print*,recoil_energy,momentum_transfer,event_rate_spectra(calc_num)
        end do

        ! Write results to file
        open(newunit=iunit, file='eventrate_spectra.dat')
        write(iunit,"(A,T30,A)")"# Recoil energy (kev)","Event rate (events/gev)"
        do calc_num = 1, energy_grid_size
            recoil_energy = energy_grid(calc_num)
            write(iunit,*)recoil_energy,event_rate_spectra(calc_num)
        end do

        close(iunit)
        print*,"Event rate spectra written to eventrate_spectra.dat"

    end subroutine eventrate_spectra

    subroutine get_energy_grid(mtarget)
        implicit none
        integer :: calc_num
        character(len=100) :: filename
        real(kind=8) :: momentum_transfer, mtarget, recoil_energy

        ! Get recoil energy grid from user or from file
        if (useenergyfile) then
            print*,'Recoil energies will be read from file. Enter filename:'
            read*,filename
            call read_energy_grid(filename)
        else
            if (usemomentum) then
                print*,'What is the range of transfer momenta?'
                print*,'Enter starting momentum, stopping momentum, setp size:'
            else
                print*,'What is the range of recoil energies in kev?'
                print*,'Enter starting energy, stoping energy, step size:'
            end if

            read*,x_start, x_stop, x_step

            if (usemomentum) then
                print*,"q min  (gec/c)",x_start
                print*,"q max  (gev/c)",x_stop
                print*,"q step (gev/c)",x_step
            else
                print*,"E min  (kev)",x_start
                print*,"E max  (kev)",x_stop
                print*,"E step (kev)",x_step
            end if

            energy_grid_size = int((x_stop - x_start) / x_step) + 1

            allocate(energy_grid(energy_grid_size))
            do calc_num = 1, energy_grid_size
                energy_grid(calc_num) = x_start + (calc_num - 1) * x_step
            end do
        end if

        allocate(momentum_grid(energy_grid_size))
        do calc_num = 1, energy_grid_size
            if (usemomentum) then
                momentum_transfer = x_start + (calc_num - 1) * x_step
                recoil_energy = momentum_transfer**2d0 / (2d0*mtarget*mN*kev)
                momentum_grid(calc_num) = momentum_transfer
                energy_grid(calc_num) = recoil_energy
            else
                recoil_energy = energy_grid(calc_num)
                momentum_grid(calc_num) = sqrt(2d0*mtarget*mN*recoil_energy*kev)
            end if
        end do
    end subroutine

    subroutine read_energy_grid(filename)

        implicit none
        character(len=100) :: filename
        integer :: i, io, iunit

        open(newunit=iunit,file=trim(filename))
        do
            read(iunit,*,iostat=io)
            if (io/=0) exit
            energy_grid_size = energy_grid_size + 1
        end do

        allocate(energy_grid(energy_grid_size))
        rewind(iunit)

        do i = 1, energy_grid_size
            read(iunit,*) energy_grid(i)
        end do
        close(iunit)
    end subroutine read_energy_grid

    !-----------------------------------------------------------------------------80
    subroutine nucresponse_spectra

        implicit none
        integer :: iunit
        integer :: ioperator, tau, tau_prime
        integer :: iqq

        real(kind=8) :: yy, qq, xx
        real(kind=8), allocatable :: Wlist(:,:,:,:)
        real(kind=8) :: mtarget
        real(kind=8) :: tolerance

        tolerance = epsilon(qq)

        mtarget = nuc_target%mass
        call get_energy_grid(mtarget)

        allocate(Wlist(energy_grid_size, 8, 0:1, 0:1))

        if (pndens) then
          !$OMP parallel do private(qq, yy, tau, tau_prime, ioperator) schedule(dynamic, 10)
          do iqq = 1, energy_grid_size
            qq = momentum_grid(iqq)
            yy = (qq * bfm / 2d0) ** 2d0
            do tau = 0, 1
              do tau_prime = 0, 1
                do ioperator = 1, 8
                  Wlist(iqq, ioperator, tau, tau_prime) = 0.25d0 * nucFormFactor_transform(&
                      tau, tau_prime, ioperator, yy, nuc_target%densitymats%rho,&
                      nuc_target%groundstate%Tx2, nuc_target%Mt )
                end do
              end do
            end do
          end do
        else
          !$OMP parallel do private(qq, yy, tau, tau_prime, ioperator) schedule(dynamic, 10)
          do iqq = 1, energy_grid_size
            qq = momentum_grid(iqq)
            yy = (qq * bfm / 2d0) ** 2d0
            do tau = 0, 1
              do tau_prime = 0, 1
                do ioperator = 1, 8
                  Wlist(iqq, ioperator, tau, tau_prime) = 0.25d0 * nucFormFactor(&
                      tau, tau_prime, ioperator, yy, nuc_target%densitymats%rho,&
                      nuc_target%groundstate%Tx2, nuc_target%Mt )
                end do
              end do
            end do
          end do
        end if

        open(newunit=iunit, file='nucresponse_spectra.dat')
        write(iunit,'(a,i3)')'# Z = ',nuc_target%Z
        write(iunit,'(a,i3)')'# N = ',nuc_target%N
        write(iunit,'(a)')'# Nuclear response functions W(operator,tau, tau_prime, q)'
        write(iunit,'(a)')'# operator = 1, ..., 8:'
        write(iunit,'(a)')"# M, \Phi'', \tilde{\Phi}', \Delta, \Sigma', \Sigma'', \Phi''M, \Delta \Sigma'"
        write(iunit,'(a)')'# tau = 0, 1 isospin coupling'
        if (usemomentum) then
            write(iunit,'(a)')'# x = momentum transfer (gev/c)'
        else
            write(iunit,'(a)')'# x = recoil energy (kev)'
        end if
        write(iunit,'(a)')  '# Format:'
        write(iunit,'(a)')  '# x, { W(i, 0, 0, x) W(i, 1, 0, x) W(i, 0, 1, x) W(i, 1, 1, x); i = 1, ..., 8}'

        do iqq = 1, energy_grid_size
          if (usemomentum) then
              xx = momentum_grid(iqq)
          else
              xx = energy_grid(iqq)
          end if
          write(iunit,'(ES15.5E3,32(ES15.5E3))') xx, Wlist(iqq, 1:8, 0:1, 0:1)
        end do
        close(iunit)

        deallocate(Wlist)

    end subroutine nucresponse_spectra

end module spectra

max_stars_repo_name: ogorton/dmfortfactor

max_issues_repo_name: ogorton/dmfortfactor

max_forks_repo_name: ogorton/dmfortfactor


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE GRNCRV (LABSID, ICURVE, NPTS, XPTS, YPTS, DRAWLN)
C=======================================================================

C   --*** GRNCRV *** (GRPLIB) Label Curve (PLT)
C   --   Written by Amy Gilkey - revised 10/07/87
C   --
C   --GRNCRV finds the device coordinates of the first, last or middle
C   --point of a curve within the plot window and labels the curve with
C   --the line number.
C   --
C   --Parameters:
C   --   LABSID - IN - label curve on "FIRST", "MIDDLE" or "LAST" point
C   --   ICURVE - IN - the curve number
C   --   NPTS - IN - the number of points in the curve
C   --   XPTS, YPTS - IN - the X and Y points
C   --   DRAWLN - IN - true if vector, false if points

C   --Routines Called:
C   --   MP2VC - (PLTLIB) Get device coordinates of vector
C   --   MP2PT - (PLTLIB) Get device coordinates of point
C   --   PLTGTT - (PLTLIB) Get text parameter
C   --      2 = (KSCHSZ) software character size
C   --   INTSTR - (STRLIB) Convert integer to string

      PARAMETER (KSCHSZ=2)

      CHARACTER*8 LABSID
      INTEGER ICURVE
      INTEGER NPTS
      REAL XPTS(*), YPTS(*)
      LOGICAL DRAWLN

      LOGICAL LDUM, PLTGTT
      CHARACTER*2 STR2

C   --Find the nearest line or point within the window

      IF (LABSID .EQ. 'FIRST') THEN
         ISTART = 1
         IEND = NPTS
         INC = 1
      ELSE IF (LABSID .EQ. 'MIDDLE') THEN
         ISTART = NPTS / 2
         IEND = NPTS
         INC = 1
      ELSE
         ISTART = NPTS
         IEND = 1
         INC = -1
      END IF

      IF (DRAWLN) THEN
         DO 100 I = ISTART, IEND-INC, INC
            MASK = -1
            CALL MP2VC (1, XPTS(I), YPTS(I), XPTS(I+INC), YPTS(I+INC),
     &         DX, DY, DXX, DYY, MASK)
            IF (MOD (MASK,2) .NE. 0) GOTO 120
  100    CONTINUE

      ELSE
         DO 110 I = ISTART, IEND, INC
            MASK = -1
            CALL MP2PT (1, XPTS(I), YPTS(I), DX, DY, MASK)
            IF (MOD (MASK,2) .NE. 0) GOTO 120
  110    CONTINUE
      END IF

      GOTO 130

C   --Write the curve number (location is dependent on the placement)

  120 CONTINUE
C   --Get curve number
      CALL INTSTR (1, 0, ICURVE, STR2, LNUM)

      LDUM = PLTGTT (KSCHSZ, VCS)
      IF (DRAWLN) THEN
         DXO = .25
      ELSE
         DXO = .50
      END IF
      IF (LABSID .EQ. 'FIRST') THEN
         DX = DX - (LNUM+DXO) * VCS
      ELSE
         DX = DX + DXO*VCS
      END IF
      IF (LABSID .EQ. 'MIDDLE') THEN
         DY = DY + .002
      ELSE
         DY = DY - .5*VCS
      END IF

      CALL GRTEXT (DX, DY, STR2)

  130 CONTINUE
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! { dg-do run }
!
! PR 40164: Fortran 2003: "Arrays of procedure pointers" (using PPCs)
!
! Original test case by Barron Bichon <barron.bichon@swri.org>
! Adapted by Janus Weil <janus@gcc.gnu.org>

PROGRAM test_prog

  ABSTRACT INTERFACE
  FUNCTION fn_template(n,x) RESULT(y)
    INTEGER, INTENT(in) :: n
    REAL, INTENT(in) :: x(n)
    REAL :: y(n)
  END FUNCTION fn_template
  END INTERFACE

  TYPE PPA
    PROCEDURE(fn_template), POINTER, NOPASS :: f
  END TYPE PPA

 TYPE ProcPointerArray
   PROCEDURE(add), POINTER, NOPASS :: f
 END TYPE ProcPointerArray

 TYPE (ProcPointerArray) :: f_array(3)
 PROCEDURE(add), POINTER :: f
 real :: r

 f_array(1)%f => add
 f => f_array(1)%f
 f_array(2)%f => sub
 f_array(3)%f => f_array(1)%f

 r = f(1.,2.)
 if (abs(r-3.)>1E-3) call abort()
 r = f_array(1)%f(4.,2.)
 if (abs(r-6.)>1E-3) call abort()
 r = f_array(2)%f(5.,3.)
 if (abs(r-2.)>1E-3) call abort()
 if (abs(f_array(1)%f(1.,3.)-f_array(3)%f(2.,2.))>1E-3) call abort()

CONTAINS

 FUNCTION add(a,b) RESULT(sum)
   REAL, INTENT(in) :: a, b
   REAL :: sum
   sum = a + b
 END FUNCTION add

 FUNCTION sub(a,b) RESULT(diff)
   REAL, INTENT(in) :: a, b
   REAL :: diff
   diff = a - b
 END FUNCTION sub

END PROGRAM test_prog


text:
program ProgramName
    DO 10 i=1,6

    10 END DO
end program ProgramName

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
PROGRAM driver

  use sp_hnfs, only: so_32
  use num_types

  implicit none

  integer :: n, i, j
  integer, allocatable :: spHNFs(:,:,:)
  real(dp), pointer :: B_vecs(:,:)
  integer :: at(1)
  integer :: Co(3,3), Cu(3,3)
  real(dp) :: No(3,3), Nu(3,3), O(3,3), U(3,3), offsets(1,3)
  integer :: n_irr, nhnfs
  real(dp) :: grid(3,3), best_offset(3)
  logical :: all_hnfs

  n = 100
  allocate(B_vecs(1,3))
  B_vecs(1,:) = (/ 0.0_dp, 0.0_dp, 0.0_dp /)
  offsets(1,:) = (/ 0.0_dp, 0.0_dp, 0.0_dp /)
  at(1) = 0
  all_hnfs = .True.

  Co(1,:) = (/ 1, 0, 0 /)
  Co(2,:) = (/ 0, 1, 0 /)
  Co(3,:) = (/ 0, 0, 1 /)

  Cu(1,:) = (/ 1, 0, 0 /)
  Cu(2,:) = (/ 0, 1, 0 /)
  Cu(3,:) = (/ 0, 0, 1 /)

  No(1,:) = (/1.0_dp, 0.0_dp, 0.0_dp /)
  No(2,:) = (/0.0_dp, 1.0_dp, 0.0_dp /)
  No(3,:) = (/0.0_dp, 0.0_dp, 1.0_dp /)

  Nu(1,:) = (/1.0_dp, 0.0_dp, 0.0_dp /)
  Nu(2,:) = (/0.0_dp, 1.0_dp, 0.0_dp /)
  Nu(3,:) = (/0.0_dp, 0.0_dp, 1.0_dp /)

  O(1,:) = (/1.0_dp, 0.0_dp, 0.0_dp /)
  O(2,:) = (/0.0_dp, 1.0_dp, 0.0_dp /)
  O(3,:) = (/0.0_dp, 0.0_dp, 1.0_dp /)

  U(1,:) = (/1.0_dp, 0.0_dp, 0.0_dp /)
  U(2,:) = (/0.0_dp, 1.0_dp, 0.0_dp /)
  U(3,:) = (/0.0_dp, 0.0_dp, 1.0_dp /)

  call so_32(n, No, Nu, Co, Cu, O, U, B_vecs, at, offsets, best_offset, spHNFs, grid, &
       n_irr, nhnfs, all_hnfs_=all_hnfs)

  open(3, file="HNFs.out")
  do i= 1,nhnfs
     write(3,'(6i6)') spHNFs(1,1,i), spHNFs(2,1,i), spHNFs(2,2,i), spHNFs(3,1,i), spHNFs(3,2,i), spHNFs(3,3,i)
  end do
  close(3)

end PROGRAM driver

max_stars_repo_name: glwhart/autoGR

max_issues_repo_name: glwhart/autoGR

max_forks_repo_name: glwhart/autoGR


text:
 ASUCD Office of City and County Affairs

The ASUCD Office of City/County Affairs is an ASUCD unit that researches and analyzes City Council city and county legislation to make sure that undergraduate students are fairly represented, and work to educate students about city and county issues which affect them. The office helps students with housing issues and organizes Housing Day.

The Office of City/County Affairs works to improve relations between university students and the community. They focus on improving student relations with the city of Davis and ensuring that the student population is represented in city affairs. Through special projects, the office seeks to serve as a liaison between the city and ASUCD. If there are any student concerns or particular projects you would like to see implemented, please feel free to contact the Office of City/County Affairs.

The Director is a member of the ASUCD StudentPolice Relations Committee and the City / UCD Student Liaison Commission CityUCD Student Liaison Commission, as well as the ASUCD Cabinet.

City/County is the only ASUCD advocacy unit to actually bring money into the association. So while the budgets for Lobby Corp and University Affairs come exclusively from students, City/County uses Housing Day as a way to bring in funds for the association and relieve some of the burden on students. In fact, the 2010 Housing Day brought in so much money that the proceeds from the event funded the offices expenses for the entire year, making the unit subsidyfree for the first time.

 Current Office

Director: Jeanna Gindi (jagindi@gmail.com)
Assistant Director of Projects: Daniel Wildman (djwildman@ucdavis.edu)
Assistant Director of Policy: Eric Renslo (emrenslo@ucdavis.edu)
Housing Coordinator: Sabrina Dias (sidias@ucdavis.edu)

 Previous Directors
20112012:
20102011: Dylan Schaefer
20092010: Don Gibson
2009: Brian Kim
20072009: users/StevenLee Steven D. Lee
20062007: Jarvis Burston
20052006: Greg Russell
20042005: Richard Cipian
20032004: Clinton Gutierrez
20022003: Andrew Glass
20012002: Amir Bahkshaie
20002001: Molly Anne Meyn
19951996: Dorothy Yee

Taken from CityUCD Student Liaison Commission Agendas

20080731 20:56:05 nbsp Other than housing day, I am curious to know how the city/county affairs office works to educate students about city and county issues which affect them... Users/PeteWillits

20091001 13:41:31 nbsp I know its been a while since you posted this question, but for Pete if youre still out there and anyone else interested. Like EVERYTHING else in ASUCD everything is theoretically and is highly contingent on the work ethic of the particular student in charge. City/County has in the past held forums, but student participation/interest is often low as would be the case with most municipal level governments. Most of their work consists of representing student interests at various bodies.   Users/StevenLee