text:
      PROGRAM PARTONSS
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      character*5 flav(3)

      data flav/'qq   ','qqp  ','gq   '/

      COMMON/IOUT/IPRINT
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/RESULT/S1,S2,S3,S4
      COMMON/IFLAVOR/IFLAVOR

      EXTERNAL BORN, SOVI, HARD, SOVISCA, HARDSCA

      CALL RSTART(12,34,56,78)

C --- PARAMETER TO CHANGE  -------------------------------------

C***  1:Q QBAR     2:QP QBAR    3:G Q

      DO 100 IFLAVOR = 3,3

C***  THE MASSES:

      MS = 280.D0
      MG = 200.D0
      MT = 175.D0

C***  START VALUE IN LOG(ETA), NUMBERS OF STEPS, STEP WIDTH IN LOG(ETA)

      ETALOG = -2D0
      INUM = 5
      STEP = 1D0

C***  THE NUMBER OF VEGAS CALLS ( DEFAULT = 500 )

      IBORN = 1000
      ISOVI = 500
      IHARD = 500

C***  PRINT VEGAS STATISTICS (10) OR NOT (0)

      IPRINT = 0

C --- CONSTANTS DON'T CHANGE PLEASE ------------------------------

      PI = 4.D0*DATAN(1.D0)
      CONV = 389379660.D0
      ALPHAS = 0.1D0
      SCA = 4.D0

C -------------------------------------------------------------

c$$$      IF (IFLAVOR.EQ.1) OPEN(10,FILE="partonss_qq.dat")
c$$$      IF (IFLAVOR.EQ.2) OPEN(10,FILE="partonss_qp.dat")
c$$$      IF (IFLAVOR.EQ.3) OPEN(10,FILE="partonss_gq.dat")

      PRINT 20, MS,MG,MT
      PRINT 30
c$$$      WRITE(10, 20) MS,MG,MT
c$$$      WRITE(10, 30)

 20   FORMAT(' MS = ',G10.4,' MG= ',G10.4,' MT = ',G10.4)
 30   FORMAT('  ETA        BORN          SOVI       REL     ',
     +   '   HARD       REL       SCALE')

C --- INTEGRATION BY VEGAS ------------------------------------

      DO 100  I = 1,INUM

        ETA = 10.D0**(ETALOG)
        ETALOG = ETALOG +STEP

        S = 4.D0*MS**2*(ETA + 1.D0)
        DEL = 1.D-4*MS**2 * (1 -4.D0*MS**2/S)

      IF (IFLAVOR.LE.2) THEN
         CALL INTEG(BORN,1,IBORN,3,5*IBORN,5,1.D-4,RES)
         RESBORN = MS**2/ALPHAS**2/CONV*S1
         DEVBORN = MS**2/ALPHAS**2/CONV*S2

         CALL INTEG(SOVI,1,ISOVI,3,5*ISOVI,5,5.D-4,RES)
         RESSOVI = MS**2/ALPHAS**3/4.D0/PI/CONV*S1
         DEVSOVI = MS**2/ALPHAS**3/4.D0/PI/CONV*S2
         RELSOVI = DABS(DEVSOVI/RESSOVI)

         CALL INTEG(SOVISCA,1,ISOVI,3,5*ISOVI,5,5.D-4,RES)
         RESSOVISCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S1
         DEVSOVISCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S2

         CALL INTEG(HARD,3,ISOVI,3,5*ISOVI,5,5.D-4,RES)
         RESHARD = MS**2/ALPHAS**3/4.D0/PI/CONV*S1
         DEVHARD = MS**2/ALPHAS**3/4.D0/PI/CONV*S2
         RELHARD = DABS(DEVHARD/RESHARD)

         CALL INTEG(HARDSCA,2,ISOVI,3,5*ISOVI,5,5.D-4,RES)
         RESHARDSCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S1
         DEVHARDSCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S2
      ELSE
         CALL INTEG(HARD,3,IHARD,5,5*IHARD,10,5.D-4,RES)
         RESHARD = MS**2/ALPHAS**3/4.D0/PI/CONV*S1
         DEVHARD = MS**2/ALPHAS**3/4.D0/PI/CONV*S2
         RELHARD = DABS(DEVHARD/RESHARD)

         CALL INTEG(HARDSCA,2,ISOVI,3,5*ISOVI,5,5.D-4,RES)
         RESHARDSCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S1
         DEVHARDSCA = MS**2/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S2

         RESBORN = 0D0
         RESSOVI = 0D0
         RELSOVI = 0D0
         RESSOVISCA = 0D0

      ENDIF

      RESSCALE = RESSOVISCA +RESHARDSCA

c$$$      PRINT 50,ETA,RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE
c$$$      WRITE(10,50)ETA,RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE

      PRINT 51,flav(iflavor),ETA,
     +     RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE

 50   FORMAT(G10.4,' ',G11.5,'  ',G11.5,' ',G8.2,'  ',G11.5,' ',
     +     G8.2,'  ',G11.5)

 51   FORMAT(A5,G9.3,' ',G11.5,'  ',G11.5,' ',F5.4,'  ',G11.5,' ',
     +     F5.4,'  ',G11.5)

 100  CONTINUE

      END

C ======================================================================

      DOUBLE PRECISION FUNCTION BORN(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      SB = S * DSQRT(1.D0 -4.D0*MS**2/S)
      T1 = SB*(X -1.D0/2.) -S/2.D0
      IF (IFLAVOR.EQ.1) BORN = SB*DSSQQB(ALPHAS,S,T1,MS,MG,1)
      IF (IFLAVOR.EQ.2) BORN = SB*DSSQQB(ALPHAS,S,T1,MS,MG,0)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION SOVI(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      SB = S * DSQRT(1.D0 -4.D0*MS**2/S)
      T1 = SB*(X -1.D0/2.) -S/2.D0
      IF (IFLAVOR.EQ.1) SOVI = SB*DSSQQV(ALPHAS,S,T1,MS,MG,MT,1)
      IF (IFLAVOR.EQ.2) SOVI = SB*DSSQQV(ALPHAS,S,T1,MS,MG,MT,0)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION HARD(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(3)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      Y = VAR(2)
      Z = VAR(3)

      SB = S * DSQRT(1.D0 -4.D0*MS**2/S)
      T1 = DSQRT((DEL -S)**2 -4*MS**2*S)*(X -1.D0/2.) -(S -DEL)/2.D0
      S4MAX = S +T1 +MS**2*S/T1
      S4 = ( S4MAX -DEL)*Y +DEL
      PREF = S4MAX*SB

C***  EPS = ( GAMMA(GLUINO)/MASS(GLUINO) )**2
      EPS = 1.D-8

      IF (IFLAVOR.EQ.3) THEN
         M2 = MG**2 -MS**2
         S4MAX0 = S - 2.D0*DSQRT(S*MS**2)
         GAM = DSQRT(EPS)*MG**2
         TYMIN = DATAN(-M2/GAM)
         TYMAX = DATAN((S4MAX0-M2)/GAM)
         TY = (TYMAX -TYMIN) * X + TYMIN
         PREFTY = (TYMAX -TYMIN)/GAM
         S40 = GAM*DTAN(TY) + M2
         S41 = M2

         PRET10 = DSQRT( (S-S40)**2 -4*MS**2*S )
         PRET11 = 1
         IF((S-S41)**2-4*MS**2*S.GT.0.D0)THEN
          PRET11 = DSQRT( (S-S41)**2 -4*MS**2*S )
         ENDIF
         T10 = -(S-S40)/2.D0 +(Y -0.5D0)*PRET10
         T11 = -(S-S41)/2.D0 +(Y -0.5D0)*PRET11

         PREF0 = PRET10 * PREFTY
         PREF1 = PRET11 * PREFTY
         S4G02 = (S40 -M2)**2 + EPS*MG**4

         S3MAX0 = S - 2.D0*DSQRT(S*MS**2)
         GAM = DSQRT(EPS)*MG**2
         TZMIN = DATAN(-M2/GAM)
         TZMAX = DATAN((S3MAX0-M2)/GAM)
         TZ = (TZMAX -TZMIN) * X + TZMIN
         PREFTZ = (TZMAX -TZMIN)/GAM
         S32 = GAM*DTAN(TZ) + M2
         S33 = M2

         PRES42 = S32/(S32+MS**2)*DSQRT((S-S32)**2 -4*MS**2*S)
         PRES43 = 1
         IF((S-S33)**2 -4*MS**2*S.GT.0.D0)THEN
          PRES43 = S33/(S33+MS**2)*DSQRT((S-S33)**2 -4*MS**2*S)
         ENDIF
         S42 = S32/2D0/(S32+MS**2)*(S -S32 -2*MS**2) +PRES42 *(Y -0.5D0)
         S43 = S33/2D0/(S33+MS**2)*(S -S33 -2*MS**2) +PRES43 *(Y -0.5D0)
         PRET12 = DSQRT( (S-S42)**2 -4*MS**2*S )
         PRET13 = 1
         IF((S-S43)**2 -4*MS**2*S.GT.0.D0)THEN
          PRET13 = DSQRT( (S-S43)**2 -4*MS**2*S )
         ENDIF
         T12 = -(S-S42)/2.D0 +(Z -0.5D0)*PRET12
         T13 = -(S-S43)/2.D0 +(Z -0.5D0)*PRET13
         PREX2 = S42/2D0/(S42+MS**2)*DSQRT((S-S42)**2-4*MS**2*S)
         PREX3 = 1
         IF((S-S43)**2 -4*MS**2*S.GT.0.D0)THEN
          PREX3 = S43/2D0/(S43+MS**2)*DSQRT((S-S43)**2-4*MS**2*S)
         ENDIF

         PREF2 = PRES42 *PRET12 /PREX2 * PREFTZ
         PREF3 = PRES43 *PRET13 /PREX3 * PREFTZ

         S3G2  = S32 -M2
         S3G3  = 0D0

         SB4 = S * DSQRT(1.D0 -4.D0*MS**2/S)
         T14 = DSQRT(S**2 -4*MS**2*S)*(X -1.D0/2.) -S/2.D0
         S4MAX4 = S +T14 +MS**2*S/T14
         PREF4 = SB4

      END IF

      IF (IFLAVOR.EQ.1) HARD =
     +     + PREF*DSSQQH(ALPHAS,S,T1,S4,MS,MG,1)
     +     + PREF*DSSQQD(ALPHAS,S,T1,S4,MS,MG,DEL,S4MAX,1)/(S4MAX-DEL)
      IF (IFLAVOR.EQ.2) HARD =
     +     + PREF*DSSQQH(ALPHAS,S,T1,S4,MS,MG,0)
     +     + PREF*DSSQQD(ALPHAS,S,T1,S4,MS,MG,DEL,S4MAX,0)/(S4MAX-DEL)

      IF (IFLAVOR.EQ.3) THEN
         HARD1 = PREF0 * DSSGQH(ALPHAS,S,T10,S40,MS,MG,EPS)*S4G02
         IF ((MG.GT.MS).AND.(S.GT.(MS+MG)**2)) THEN
            HARD2=(PREF0 *DSSGQS(ALPHAS,S,T10,S40,MS,MG,EPS)
     +           - PREF1 *DSSGQS(ALPHAS,S,T11,M2 ,MS,MG,EPS) )
            HARD3 =
     +           + PREF2 *DSSGQT(ALPHAS,S,T12,S42,S3G2,MS,MG,EPS)
     +           - PREF3 *DSSGQT(ALPHAS,S,T13,S43,0.D0,MS,MG,EPS)
         ELSE
            HARD2= PREF0 *DSSGQS(ALPHAS,S,T10,S40,MS,MG,EPS)
            HARD3= PREF2 *DSSGQT(ALPHAS,S,T12,S42,S3G2,MS,MG,EPS)
         END IF
         IF ((S4MAX4.GE.M2).AND.(MG.GT.MS).AND.(S.GT.(MS+MG)**2)) THEN
            HARD4 =PREF4*DSSGQU(ALPHAS,S,T14,+M2,MS,MG,EPS)
         ELSE
            HARD4 = 0D0
         END IF

         HARD = HARD1 +HARD2 +HARD3 +HARD4

      END IF

      IF (.NOT.(DABS(HARD).LT.1.D40)) THEN
         HARD = 0.D0
         PRINT *,'NAN'
      END IF

      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION SOVISCA(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      SB = S * DSQRT(1.D0 -4.D0*MS**2/S)
      T1 = SB*(X -1.D0/2.) -S/2.D0
      IF (IFLAVOR.EQ.1) SOVISCA = SB*DSSQQ1(ALPHAS,S,T1,MS,MG,SCA,1)
      IF (IFLAVOR.EQ.2) SOVISCA = SB*DSSQQ1(ALPHAS,S,T1,MS,MG,SCA,0)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION HARDSCA(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(2)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      Y = VAR(2)
      SB = S * DSQRT(1.D0 -4.D0*MS**2/S)
      T1 = DSQRT((DEL -S)**2 -4*MS**2*S)*(X -1.D0/2.) -(S -DEL)/2.D0
      S4MAX = S +T1 +MS**2*S/T1
      S4 = ( S4MAX -DEL)*Y +DEL
      PREF = (S4MAX -DEL)*(SB -S/SB*DEL)
      IF (IFLAVOR.EQ.1) HARDSCA =
     +   + PREF*DSSQQ3(ALPHAS,S,T1,S4,MS,MG,SCA,1)
     +   + PREF*DSSQQ2(ALPHAS,S,T1,S4,MS,MG,DEL,S4MAX,SCA,1)/(S4MAX-DEL)
      IF (IFLAVOR.EQ.2) HARDSCA =
     +   + PREF*DSSQQ3(ALPHAS,S,T1,S4,MS,MG,SCA,0)
     +   + PREF*DSSQQ2(ALPHAS,S,T1,S4,MS,MG,DEL,S4MAX,SCA,0)/(S4MAX-DEL)
      IF (IFLAVOR.EQ.3) HARDSCA =
     +   + PREF*DSSGQ3(ALPHAS,S,T1,S4,MS,MG,SCA)

      RETURN
      END

C ======================================================================

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text: The Davis Tutor is the best place to receive lower division tutoring in Math, Physics, Engineering, Computer Science and High/Middle school course work. Upper division courses in Mechanical Engineering are especially welcomed. Bobby Kian has worked as a physics tutor at his local college, with knowledge in best assisting others. He recently Graduated UC Davis with a Mechanical Engineering degree and works in the Mechatronics lab. He speaks fluent English and is kind of funny too. He makes school fun as it should be, and the learning experience a lot more fulfilling. It is backed up with a guarantee that the tutor is able to fully help you with your class work or your money back. 25/hour for individual and study groups. Please call ahead of time to schedule an appointment.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


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

20090108 18:13:04 nbsp Welcome to the Wiki. Users/JasonAller

20090515 13:09:52 nbsp http://daviswiki.org/UC_Davis_Computer_Labs?actiondiff&version219&version118 Please http://daviswiki.org/UC_Davis_Computer_Labs?actiondiff&version221&version120 stop http://daviswiki.org/UC_Davis_Computer_Labs?actiondiff&version225&version124 removing the picture of the old print interface. Further removals will result in further restorations and youll just be wasting your time as well as that of your fellow editors. Were information junkies here and such epherma holds an attraction for most of us. The picture shows the way the print interface used to look and contains something unexpected if you look at it closer. If you have a really good reason why we should make the wiki less enjoyable and informative, please present it. The fact that some people may find it offensive isnt good enough. Unfortunately, just about EVERYTHING in the world is offensive to somebody. We couldnt exist if we excised everything remotely offensive from the wiki. Users/WilliamLewis


module gs2_gryfx_zonal

  use gs2_main, only: gs2_program_state_type
  use iso_c_binding

  implicit none
  private

  type(gs2_program_state_type) :: state
  public :: state
  public :: init_gs2_gryfx, advance_gs2_gryfx, finish_gs2_gryfx
  public :: allocate_gryfx_zonal_arrays
  public :: deallocate_gryfx_zonal_arrays
  public :: gryfx_parameters_type

  public :: interpolate_theta ! for testing
  public :: gs2_2_gryfx_grid

  public :: test_flag

  logical :: test_flag = .false.

  type, bind(c) :: gryfx_parameters_type
      integer(c_int) :: mpirank
      integer(c_int) :: restart
      integer(c_int) :: nstep

       ! Name of gryfx/gryfx input file
       !character(len=1000) :: input_file
      !Base geometry parameters - not currently set by trinity
      !See geometry.f90
       integer(c_int) :: equilibrium_type
       !character(len=800) :: eqfile
       integer(c_int) :: irho
       real(c_double) :: rhoc
       integer(c_int) :: bishop
       integer(c_int) :: nperiod
       integer(c_int) :: ntheta

      ! Miller parameters
       real(c_double) :: rgeo_lcfs
       real(c_double) :: rgeo_local
       real(c_double) :: akappa
       real(c_double) :: akappri
       real(c_double) :: tri
       real(c_double) :: tripri
       real(c_double) :: shift
       real(c_double) :: qinp
       real(c_double) :: shat
       real(c_double) :: asym
       real(c_double) :: asympri

       ! Circular parameters
       real(c_double) :: eps

       ! Other geometry parameters - Bishop/Greene & Chance
       real(c_double) :: beta_prime_input
       real(c_double) :: s_hat_input

       ! Flow shear
       real(c_double) :: g_exb
       ! Species parameters... I think allowing 20 species should be enough!
       ! Allocating the structs and arrays is tedious and prone to segfaults
       ! and is unnecessary given the tiny memory usage of this data object
       ! Most importantly it is not interoperable with C!
       integer(c_int) :: ntspec
       real(c_double):: dens(20)
       real(c_double):: temp(20)
       real(c_double):: fprim(20)
       real(c_double):: tprim(20)
       real(c_double):: nu(20)

       type(c_ptr) :: everything_struct_address

     end type gryfx_parameters_type

    integer, parameter :: verb = 2

    real, dimension(:,:), allocatable :: gs2_2_gryfx_grid, gryfx_2_gs2_grid

contains

  subroutine create_interp_matrices(gryfx_theta)
    use theta_grid, only: ntgrid, theta
    use splines, only: inter_cspl
    real*8, intent(in), dimension(1:2*ntgrid) :: gryfx_theta
    real, dimension(1:2*ntgrid) :: delta_array, gryfx_theta_real
    integer ::iz

    ! NB Gryfx grid has one less theta point
    gryfx_theta_real = gryfx_theta

    !write (*,*) 'theta', theta, 'gryfx_theta_real', gryfx_theta_real

    do iz = 1,2*ntgrid
      delta_array = 0.
      delta_array(iz) = 1.
      call inter_cspl(2*ntgrid, gryfx_theta_real, delta_array, &
                      2*ntgrid, theta(-ntgrid:ntgrid-1), gryfx_2_gs2_grid(:,iz))
      call inter_cspl(2*ntgrid, theta(-ntgrid:ntgrid-1), delta_array, &
                      2*ntgrid, gryfx_theta_real, gs2_2_gryfx_grid(:,iz))
    end do

  end subroutine create_interp_matrices

  subroutine interpolate_theta(matrix, array, is_field)
    use theta_grid, only: ntgrid
    use kt_grids, only: naky, ntheta0
    use species, only: nspec
    real, dimension(:,:), intent(in) :: matrix
    complex*8, dimension (:), intent (inout) :: &
                    array
    logical, intent(in) :: is_field
    complex*8, dimension (:), allocatable :: array_temp
    integer :: iz_row, iz_col, is, it, ik, index_gryfx_row, index_gryfx_col
    integer :: spec_max

    return

    allocate(array_temp(naky*ntheta0*2*ntgrid*nspec))

    array_temp = array

    spec_max = nspec
    if (is_field) spec_max = 1

    do is = 1,spec_max
      do ik = 1,naky
        do it = 1,ntheta0
          do iz_row = 1,2*ntgrid
            index_gryfx_row = 1 + &
                             (ik-1) + &
                             naky*(it-1) + &
                             naky*ntheta0*(iz_row-1) + &
                             (2*ntgrid)*naky*ntheta0*(is-1)
            array(index_gryfx_row) = 0.0
            do iz_col = 1,2*ntgrid
              index_gryfx_col = 1 + &
                               (ik-1) +  &
                               naky*(it-1) + &
                               naky*ntheta0*(iz_col-1) + &
                               (2*ntgrid)*naky*ntheta0*(is-1)

              array(index_gryfx_row) = array(index_gryfx_row) + &
                 matrix(iz_row, iz_col) * &
                 array_temp(index_gryfx_col)

            end do
          end do
        end do
      end do
    end do

  end subroutine interpolate_theta

  subroutine allocate_gryfx_zonal_arrays
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    use nonlinear_terms, only: gryfx_zonal

    implicit none
    integer :: tot_size

    tot_size = naky*ntheta0*2*ntgrid*nspec

    allocate(gs2_2_gryfx_grid(2*ntgrid, 2*ntgrid))
    allocate(gryfx_2_gs2_grid(2*ntgrid, 2*ntgrid))
    allocate(gryfx_zonal%NLdens_ky0(tot_size))
    allocate(gryfx_zonal%NLupar_ky0(tot_size))
    allocate(gryfx_zonal%NLtpar_ky0(tot_size))
    allocate(gryfx_zonal%NLtprp_ky0(tot_size))
    allocate(gryfx_zonal%NLqpar_ky0(tot_size))
    allocate(gryfx_zonal%NLqprp_ky0(tot_size))
  end subroutine allocate_gryfx_zonal_arrays

  subroutine deallocate_gryfx_zonal_arrays
    use mp, only: mp_abort
    use nonlinear_terms, only: gryfx_zonal
    use unit_tests, only: debug_message

    implicit none

    if(.not. associated(gryfx_zonal%NLqprp_ky0)) &
      call mp_abort("Attempting to deallocated unassociated pointers&
     & in deallocate_gryfx_zonal_arrays", .true.)

    deallocate(gs2_2_gryfx_grid, gryfx_2_gs2_grid)

    deallocate(gryfx_zonal%NLdens_ky0)
    deallocate(gryfx_zonal%NLupar_ky0)
    deallocate(gryfx_zonal%NLtpar_ky0)
    call debug_message(verb, &
      'gs2_gryfx_zonal::deallocate_gryfx_zonal_arrays 50% done')
    deallocate(gryfx_zonal%NLtprp_ky0)
    call debug_message(verb, &
      'gs2_gryfx_zonal::deallocate_gryfx_zonal_arrays 65% done')
    deallocate(gryfx_zonal%NLqpar_ky0)
    call debug_message(verb, &
      'gs2_gryfx_zonal::deallocate_gryfx_zonal_arrays 85% done')
    deallocate(gryfx_zonal%NLqprp_ky0)
    call debug_message(verb, &
      'gs2_gryfx_zonal::deallocate_gryfx_zonal_arrays finished')
  end subroutine deallocate_gryfx_zonal_arrays

  subroutine init_gs2_gryfx_c(strlen, run_name, mp_comm, &
                            nz, &
                            gryfx_theta, &
                            gryfx_parameters) &
                            bind(c, name='init_gs2')
    use iso_c_binding
    use theta_grid, only: ntgrid
    use mp, only: proc0, broadcast
    implicit none
    integer(c_int), intent(in) :: strlen
    character(kind=c_char), intent(in) :: run_name
    integer(c_int), intent(in) :: mp_comm
    real(c_float), intent(in), dimension(nz) :: gryfx_theta
    type(gryfx_parameters_type), intent(in) :: gryfx_parameters
    integer(c_int), intent(in) :: nz
    real*8, dimension(nz) :: gryfx_theta_real

    gryfx_theta_real = gryfx_theta
    !call broadcast(gryfx_theta_real)

    !write (*,*) 'gryfx_theta_real is ', gryfx_theta_real
    call init_gs2_gryfx(strlen, run_name, mp_comm, gryfx_theta_real ,gryfx_parameters)

  end subroutine init_gs2_gryfx_c

  subroutine init_gs2_gryfx(strlen, file_name, mp_comm, &
                            gryfx_theta, gryfx_parameters )
    use gs2_main, only: initialize_gs2
    use gs2_main, only: initialize_equations
    use gs2_main, only: initialize_diagnostics
    use gs2_main, only: prepare_miller_geometry_overrides
    use gs2_main, only: prepare_profiles_overrides
    use gs2_main, only: prepare_initial_values_overrides
    use gs2_main, only: prepare_kt_grids_overrides
    use theta_grid, only: ntgrid
    use nonlinear_terms, only: gryfx_zonal
    use file_utils, only: run_name
    use mp, only: proc0, mp_abort, broadcast
    use unit_tests, only: debug_message
    use geometry, only: equal_arc
    use gs2_save, only: init_dt
    use gs2_time, only: init_delt

    implicit none
    integer, intent(in) :: strlen
    character (len=strlen), intent (in) :: file_name
    integer, intent(in) :: mp_comm
    real*8, intent(inout), dimension(2*ntgrid) :: gryfx_theta
    type(gryfx_parameters_type), intent(in) :: gryfx_parameters
    integer :: istatus
    real :: delt

    !gryfx_zonal%on = .true.
    state%mp_comm_external = .true.
    state%mp_comm = mp_comm
    state%is_trinity_job = .true.

    state%run_name_external = .true.
    state%run_name = file_name

    gryfx_zonal%on = .true.

    call initialize_gs2(state)
    !if(proc0) write(*,*) 'initialize_gs2 complete. gs2 thinks run_name is ', run_name
    !set overrides of naky, x0, y0, dt, vnewk here.
    !currently this is hard-coded in kt_grids.f90 and run_parameters.f90
    !just search for "GRYFX"
    if(.not. test_flag) then
      call prepare_miller_geometry_overrides(state)
      call set_miller_geometry_overrides
    end if
    call prepare_kt_grids_overrides(state)
    call debug_message(verb, 'setting overrides')
    state%init%kt_ov%override_gryfx = .true.
    state%init%kt_ov%gryfx = .true.
    if(.not. test_flag) then
      call prepare_profiles_overrides(state)
      call set_profiles_overrides
    endif
    if (gryfx_parameters%restart==1) then
      ! Load the saved value of dt
      call init_dt(delt, istatus)
      call init_delt(delt)
      ! Set up gs2 state so that it reads the
      ! initial values from the restart file.
      call prepare_initial_values_overrides(state)
      if (state%init%initval_ov%in_memory) then
        call mp_abort('ERROR: Gryfx/GS2 restart does not currently work with &
          & in_memory set to true', .true.)
      end if
      state%init%initval_ov%override = .true.
    else
      state%init%initval_ov%override = .false.
    end if
    call initialize_equations(state)
    call debug_message(verb, 'initialize_equations complete.')
    call initialize_diagnostics(state)
    call debug_message(verb, 'initialize_diagnostics complete.')

    call broadcast(equal_arc)
    if (equal_arc) then
      write (*,*) 'equal_arc true'
      call mp_abort(&
      "gs2_gryfx_zonal doesn't work with equal_arc yet", .true.)
    end if

    call allocate_gryfx_zonal_arrays

    call broadcast(gryfx_theta)
    call create_interp_matrices(gryfx_theta)

contains

      subroutine set_miller_geometry_overrides
        state%init%mgeo_ov%override_rhoc = .true.
        state%init%mgeo_ov%rhoc = gryfx_parameters%rhoc

        state%init%mgeo_ov%override_qinp = .true.
        state%init%mgeo_ov%qinp = gryfx_parameters%qinp ! qval_gs2

        state%init%mgeo_ov%override_shat = .true.
        state%init%mgeo_ov%shat = gryfx_parameters%shat ! shat_gs2

        state%init%mgeo_ov%override_rgeo_lcfs = .true.
        state%init%mgeo_ov%rgeo_lcfs = gryfx_parameters%rgeo_lcfs

        state%init%mgeo_ov%override_rgeo_local = .true.
        state%init%mgeo_ov%rgeo_local = gryfx_parameters%rgeo_local

        ! OK
        state%init%mgeo_ov%override_akappa = .true.
        state%init%mgeo_ov%akappa = gryfx_parameters%akappa

        state%init%mgeo_ov%override_akappri = .true.
        state%init%mgeo_ov%akappri = gryfx_parameters%akappri

        state%init%mgeo_ov%override_tri = .true.
        state%init%mgeo_ov%tri = gryfx_parameters%tri

        state%init%mgeo_ov%override_tripri = .true.
        state%init%mgeo_ov%tripri = gryfx_parameters%tripri

        state%init%mgeo_ov%override_shift = .true.
        state%init%mgeo_ov%shift = gryfx_parameters%shift

        state%init%mgeo_ov%override_betaprim = .true.
        state%init%mgeo_ov%betaprim = gryfx_parameters%beta_prime_input
      end subroutine set_miller_geometry_overrides

      subroutine set_profiles_overrides
        use species, only: nspec
        integer :: is, isg, idens
        !write (*,*) 'SET_PROFILES_OVERRIDES'
        do is = 1,nspec
            isg = is
          state%init%prof_ov%override_dens(isg) = .true.
          state%init%prof_ov%dens(isg) = gryfx_parameters%dens(is) ! dens_ms(1,is)/dens_ms(1,idens)

          state%init%prof_ov%override_temp(isg) = .true.
          state%init%prof_ov%temp(isg) = gryfx_parameters%temp(is) ! temp_ms(1,is)/temp_ms(1,1)

          state%init%prof_ov%override_fprim(isg) = .true.
          state%init%prof_ov%fprim(isg) = gryfx_parameters%fprim(is) ! fprim_gs2(is)

          state%init%prof_ov%override_tprim(isg) = .true.
          state%init%prof_ov%tprim(isg) = gryfx_parameters%tprim(is) ! tprim_gs2(is)

          state%init%prof_ov%override_vnewk(isg) = .true.
          state%init%prof_ov%vnewk(isg) = gryfx_parameters%nu(is) ! nu_gs2(is)
        end do
        !gs2_state%init%prof_ov%override_g_exb = .true.
        !gs2_state%init%prof_ov%g_exb = gexb_gs2
        !gs2_state%init%prof_ov%override_mach = .true.
        !gs2_state%init%prof_ov%mach = mach_gs2

      end subroutine set_profiles_overrides

  end subroutine init_gs2_gryfx

  subroutine advance_gs2_gryfx_c(istep, dens_ky0, upar_ky0, tpar_ky0, tprp_ky0, &
                               qpar_ky0, qprp_ky0, phi_ky0, first_half_step) &
                                    bind(c, name='advance_gs2')
    use iso_c_binding
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    implicit none
    complex (c_float_complex), intent (inout) :: &
                                dens_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                upar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                tpar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                tprp_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                qpar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                qprp_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (out) :: phi_ky0(naky*ntheta0*2*ntgrid)
    integer(c_int), intent (in) :: first_half_step
    integer(c_int), intent (in) :: istep
    logical :: first_half_step_l

    first_half_step_l = .false.
    if (first_half_step.eq.1) first_half_step_l = .true.

    call advance_gs2_gryfx(istep, dens_ky0, upar_ky0, tpar_ky0, tprp_ky0, &
                               qpar_ky0, qprp_ky0, phi_ky0, first_half_step_l)

  end subroutine advance_gs2_gryfx_c

  subroutine advance_gs2_gryfx(istep, dens_ky0, upar_ky0, tpar_ky0, tprp_ky0, &
                               qpar_ky0, qprp_ky0, phi_ky0, first_half_step)

    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    use nonlinear_terms, only: gryfx_zonal
    use gs2_main, only: evolve_equations
    use dist_fn, only: getmoms_gryfx_dist
    use mp, only: proc0, broadcast

    implicit none
    complex*8, dimension (:), intent (inout) :: &
                    dens_ky0, upar_ky0, tpar_ky0, tprp_ky0, qpar_ky0, qprp_ky0
    complex*8, dimension (:), intent (inout) :: phi_ky0
    logical, intent (in) :: first_half_step
    integer, intent (in) :: istep

    !state%istep_end = istep
    state%print_times = .false.

    ! only proc0 knows values of arrays, so broadcast
    call broadcast(dens_ky0)
    call broadcast(upar_ky0)
    call broadcast(tpar_ky0)
    call broadcast(tprp_ky0)
    call broadcast(qpar_ky0)
    call broadcast(qprp_ky0)
    ! all procs know value of first_half_step already

    call interpolate_theta(gryfx_2_gs2_grid, dens_ky0, .false.)
    call interpolate_theta(gryfx_2_gs2_grid, upar_ky0, .false.)
    call interpolate_theta(gryfx_2_gs2_grid, tpar_ky0, .false.)
    call interpolate_theta(gryfx_2_gs2_grid, tprp_ky0, .false.)
    call interpolate_theta(gryfx_2_gs2_grid, qpar_ky0, .false.)
    call interpolate_theta(gryfx_2_gs2_grid, qprp_ky0, .false.)

    ! now that all procs know values, copy into gryfx_zonal structure
    gryfx_zonal%first_half_step = first_half_step
    gryfx_zonal%NLdens_ky0 = dens_ky0
    gryfx_zonal%NLupar_ky0 = upar_ky0
    gryfx_zonal%NLtpar_ky0 = tpar_ky0
    gryfx_zonal%NLtprp_ky0 = tprp_ky0
    gryfx_zonal%NLqpar_ky0 = qpar_ky0
    gryfx_zonal%NLqprp_ky0 = qprp_ky0

    if(first_half_step) then
      state%dont_change_timestep = .true.
      state%skip_diagnostics = .true.
    else
      state%dont_change_timestep = .false.
      state%skip_diagnostics = .false.
    endif
    ! only advance one timestep
    call evolve_equations(state, 1)
    ! getmoms_gryfx takes moments of g and puts results into arguments.
    ! since these are the same arguments that are passed into
    ! advance_gs2_gryfx, we don't need to copy them like on the way in

    !write (*,*) 'Calling getmoms_gryfx'
    call getmoms_gryfx_dist(dens_ky0, upar_ky0, tpar_ky0, &
                 tprp_ky0, qpar_ky0, qprp_ky0, phi_ky0)

    call interpolate_theta(gs2_2_gryfx_grid, dens_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, upar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, tpar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, tprp_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, qpar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, qprp_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, phi_ky0, .true.)
  end subroutine advance_gs2_gryfx

  subroutine finish_gs2_gryfx_c() bind(c, name='finish_gs2')

    use iso_c_binding
    implicit none
    call finish_gs2_gryfx

  end subroutine finish_gs2_gryfx_c

  subroutine finish_gs2_gryfx
    use gs2_main, only: finalize_diagnostics
    use gs2_main, only: finalize_equations
    use gs2_main, only: finalize_gs2
    use nonlinear_terms, only: gryfx_zonal
    use unit_tests, only: debug_message
    implicit none

    call debug_message(verb, 'gs2_gryfx_zonal::finish_gs2_gryfx starting')
    call deallocate_gryfx_zonal_arrays
    call debug_message(verb, 'gs2_gryfx_zonal::finish_gs2_gryfx deallocated arrays')
    call finalize_diagnostics(state)
    call debug_message(verb, 'gs2_gryfx_zonal::finish_gs2_gryfx finalized diagnostics')
    call finalize_equations(state)
    call debug_message(verb, 'gs2_gryfx_zonal::finish_gs2_gryfx finalized equations')
    state%print_times = .true.
    state%print_full_timers = .true.
    call finalize_gs2(state)
    call debug_message(verb, 'gs2_gryfx_zonal::finish_gs2_gryfx finalized gs2')
    gryfx_zonal%on = .false.

  end subroutine finish_gs2_gryfx

  subroutine getmoms_gryfx_c (dens_ky0, upar_ky0, tpar_ky0, tprp_ky0, &
                              qpar_ky0, qprp_ky0, phi_ky0) &
                                   bind(c, name='getmoms_gryfx')
    use dist_fn, only: getmoms_gryfx_dist
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    use iso_c_binding
    implicit none
    complex (c_float_complex), intent (inout) :: &
                                dens_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                upar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                tpar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                tprp_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                qpar_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (inout) :: &
                                qprp_ky0(naky*ntheta0*2*ntgrid*nspec)
    complex (c_float_complex), intent (out) :: phi_ky0(naky*ntheta0*2*ntgrid)

    call getmoms_gryfx_dist(dens_ky0, upar_ky0, tpar_ky0, &
            tprp_ky0, qpar_ky0, qprp_ky0, phi_ky0)

    call interpolate_theta(gs2_2_gryfx_grid, dens_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, upar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, tpar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, tprp_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, qpar_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, qprp_ky0, .false.)
    call interpolate_theta(gs2_2_gryfx_grid, phi_ky0, .true.)
  end subroutine getmoms_gryfx_c

  subroutine broadcast_integer_c(a) bind(c, name='broadcast_integer')
    use iso_c_binding
    use mp, only : broadcast
    implicit none
    integer(c_int), intent(in out) :: a

    call broadcast(a)
  end subroutine broadcast_integer_c

end module gs2_gryfx_zonal


module integration

! This module contains the routines for the time integration and the computation
! of numerical fluxes (Rusanov or HLL in this example).

  use pde
  use global_module

  implicit none

  real(kind=8), dimension(Neq,Nx,Ny) :: U_new ! Working variable

  contains

  ! ========================================================================

  subroutine forward_Euler_step(U, dt)

    ! This function performs one step of the Forward Euler explicit time integrator

    implicit none

    real(kind=8), dimension(Neq,Nx,Ny), intent(inout)  :: U

    real(kind=8), intent(in) :: dt
    integer :: i, j, eqID

    real(kind=8), dimension(Neq) :: F_N, F_S, F_W, F_E, Src
    real(kind=8) :: y_j, y_jmhalf, y_jphalf

    ! First, assign BCs. Wall BCs need to be updated in time, taking the value
    ! from neighboring cells.
    call assign_BCs(U) ! See the pde.f03 module

    ! Update solution using the forward Euler integrator - INTERNAL CELLS only
    do j = 3, Ny-2
      do i = 3, Nx-2

        ! call compute_fluxes_Rusanov(U, i, j, F_N, F_S, F_W, F_E)
        call compute_fluxes_HLL(U, i, j, F_N, F_S, F_W, F_E)

        ! Symmetry BC
        if (j.eq.3) then
          call compute_flux_y_SYM(U(:,i,j), F_S)
        end if

        ! Compute geometrical source term (basically, the pressure)
        call compute_source_term(U,i,j, Src)

        y_j      = y_min + dy/2.0d0 + (real(j)-3.0d0)*dy
        y_jmhalf = y_j - dy/2.0d0
        y_jphalf = y_j + dy/2.0d0

        U_new(:,i,j) = U(:,i,j) - dt/dx*(F_E - F_W) - dt/(dy*y_j)*(F_N*y_jphalf - F_S*y_jmhalf) + dt*Src

        ! Check that the solution did not diverge
        do eqID = 1, Neq
          if (isnan(U_new(eqID,i,j))) then
            print*, 'Solution diverged, try with a smaller time step! Aborting.'
            stop
          end if
        end do

      end do
    end do

    ! Save solution (internal cells ONLY! Do not overwrite ghost cells!)
    U(:,3:Nx-2,3:Ny-2) = U_new(:,3:Nx-2,3:Ny-2)

  end subroutine

  ! ========================================================================

  subroutine compute_fluxes_HLL(U, i, j, F_N, F_S, F_W, F_E)

    ! Computes HLL numerical fluxes for the cell (i,j)

    implicit none

    real(kind=8), dimension(Neq,Nx,Ny), intent(in)  :: U
    real(kind=8), dimension(Neq),       intent(out) :: F_N, F_S, F_W, F_E
    integer, intent(in) :: i, j

    integer :: eqID
    real(kind=8), dimension(Neq) :: U_L, U_R, F_L, F_R

    real(kind=8) :: ws_min_L, ws_max_L, ws_min_R, ws_max_R, ws_min, ws_max ! wave speeds

    ! ------ North interface

    ! Reconstruct the solution, passing the four cells that neighbor the interface
    call reconstruct_sol_interface(U(:,i,j-1), U(:,i,j), U(:,i,j+1), U(:,i,j+2), U_L, U_R)

    call compute_flux_ws_y(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_y(U_R, F_R, ws_max_R, ws_min_R)

    ws_min = MIN(ws_min_L, ws_min_R)
    ws_max = MAX(ws_max_L, ws_max_R)

    ! HLL fluxes
    if (ws_min .ge. 0.0) then
      F_N = F_L
    else if (ws_max .lt. 0.0) then
      F_N = F_R
    else
      F_N = (ws_min*ws_max*(U_R - U_L) + ws_max*F_L - ws_min*F_R)/(ws_max - ws_min);
    end if

    ! Update global maximum wave speed (used for setting the time step)

    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    ! ------ South interface

    ! Reconstruct the solution, passing the four cells that neighbor the interface
    call reconstruct_sol_interface(U(:,i,j-2), U(:,i,j-1), U(:,i,j), U(:,i,j+1), U_L, U_R)

    call compute_flux_ws_y(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_y(U_R, F_R, ws_max_R, ws_min_R)

    ws_min = MIN(ws_min_L, ws_min_R)
    ws_max = MAX(ws_max_L, ws_max_R)

    ! HLL fluxes
    if (ws_min .ge. 0.0) then
      F_S = F_L
    else if (ws_max .lt. 0.0) then
      F_S = F_R
    else
      F_S = (ws_min*ws_max*(U_R - U_L) + ws_max*F_L - ws_min*F_R)/(ws_max - ws_min);
    end if

    ! Update global maximum wave speed (used for setting the time step)
    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    ! ------ East interface

    ! Reconstruct the solution, passing the four cells that neighbor the interface
    call reconstruct_sol_interface(U(:,i-1,j), U(:,i,j), U(:,i+1,j), U(:,i+2,j), U_L, U_R)

    call compute_flux_ws_x(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_x(U_R, F_R, ws_max_R, ws_min_R)

    ws_min = MIN(ws_min_L, ws_min_R)
    ws_max = MAX(ws_max_L, ws_max_R)

    ! HLL fluxes
    if (ws_min .ge. 0.0) then
      F_E = F_L
    else if (ws_max .lt. 0.0) then
      F_E = F_R
    else
      F_E = (ws_min*ws_max*(U_R - U_L) + ws_max*F_L - ws_min*F_R)/(ws_max - ws_min);
    end if

    ! Update global maximum wave speed (used for setting the time step)
    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    ! ------ West interface

    ! Reconstruct the solution, passing the four cells that neighbor the interface
    call reconstruct_sol_interface(U(:,i-2,j), U(:,i-1,j), U(:,i,j), U(:,i+1,j), U_L, U_R)

    call compute_flux_ws_x(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_x(U_R, F_R, ws_max_R, ws_min_R)

    ws_min = MIN(ws_min_L, ws_min_R)
    ws_max = MAX(ws_max_L, ws_max_R)

    ! HLL fluxes
    if (ws_min .ge. 0.0) then
      F_W = F_L
    else if (ws_max .lt. 0.0) then
      F_W = F_R
    else
      F_W = (ws_min*ws_max*(U_R - U_L) + ws_max*F_L - ws_min*F_R)/(ws_max - ws_min);
    end if

    ! Update global maximum wave speed (used for setting the time step)
    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

  end subroutine

  ! ==============================================================

  subroutine compute_fluxes_Rusanov(U, i, j, F_N, F_S, F_W, F_E)

    ! Computes Rusanov numerical fluxes for the cell (i,j)

    implicit none

    real(kind=8), dimension(Neq,Nx,Ny), intent(in)  :: U
    real(kind=8), dimension(Neq),       intent(out) :: F_N, F_S, F_W, F_E
    integer, intent(in) :: i, j

    integer :: eqID
    real(kind=8), dimension(Neq) :: U_L, U_R, F_L, F_R

    real(kind=8) :: ws_min_L, ws_max_L, ws_min_R, ws_max_R, ws_max ! wave speeds

    ! ------ North interface
    U_L = U(:,i,j)
    U_R = U(:,i,j+1)

    call compute_flux_ws_y(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_y(U_R, F_R, ws_max_R, ws_min_R)

    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    F_N = 0.5*(F_R + F_L) - ws_max/2.0*(U_R-U_L) ! Rusanov flux

    ! ------ South interface
    U_L = U(:,i,j-1)
    U_R = U(:,i,j)

    call compute_flux_ws_y(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_y(U_R, F_R, ws_max_R, ws_min_R)

    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    F_S = 0.5*(F_R + F_L) - ws_max/2.0*(U_R-U_L) ! Rusanov flux

    ! ------ East interface
    U_L = U(:,i,j)
    U_R = U(:,i+1,j)

    call compute_flux_ws_x(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_x(U_R, F_R, ws_max_R, ws_min_R)

    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    F_E = 0.5*(F_R + F_L) - ws_max/2.0*(U_R-U_L) ! Rusanov flux

    ! ------ West interface
    U_L = U(:,i-1,j)
    U_R = U(:,i,j)

    call compute_flux_ws_x(U_L, F_L, ws_max_L, ws_min_L)
    call compute_flux_ws_x(U_R, F_R, ws_max_R, ws_min_R)

    ws_max = MAX(ABS(ws_max_L), ABS(ws_min_L), ABS(ws_max_R), ABS(ws_min_R))

    F_W = 0.5*(F_R + F_L) - ws_max/2.0*(U_R-U_L) ! Rusanov flux

  end subroutine

  ! ========================================================================

  subroutine reconstruct_sol_interface(U_m2, U_m1, U_p1, U_p2, U_L, U_R)

  ! Takes the four points around an interface (minus_1, minus_2, plus_1, plus_2) and
  ! computes the left and right states using first order or MUSCL TVD-limited reconstruction

    real(kind=8), dimension(Neq), intent(in)  :: U_m2, U_m1, U_p1, U_p2
    real(kind=8), dimension(Neq), intent(out) :: U_L, U_R

    ! Working variables
    integer      :: eID              ! equation ID
    real(kind=8) :: theta_L, theta_R ! ratio of consecutive gradients
    real(kind=8), dimension(Neq) :: P_m2, P_m1, P_p1, P_p2, P_L, P_R ! primitive variables

    if (bool_MUSCL .EQV. .FALSE.) then
    ! ***** First order solution *****

      U_L = U_m1
      U_R = U_p1

    else
      ! ***** Second order solution *****

      ! Perform reconstruction in primitive variables
      call compute_primitive_from_conserved(U_m2, P_m2)
      call compute_primitive_from_conserved(U_m1, P_m1)
      call compute_primitive_from_conserved(U_p1, P_p1)
      call compute_primitive_from_conserved(U_p2, P_p2)

      do eID = 1, Neq ! Loop on the quantities

        ! Compute ratio of consecutive gradients
        theta_L = (P_p1(eID) - P_m1(eID))/(P_m1(eID) - P_m2(eID) + 1.0d-15) ! Add a small tolerance
        theta_R = (P_p2(eID) - P_p1(eID))/(P_p1(eID) - P_m1(eID) + 1.0d-15) ! Add a small tolerance

        ! Limited linear reconstruction
        P_L(eID) = P_m1(eID) + phi_limiter(theta_L)*(P_m1(eID) - P_m2(eID))/2.0
        P_R(eID) = P_p1(eID) - phi_limiter(1.0d0/(theta_R+1.0d-15))*(P_p2(eID) - P_p1(eID))/2.0

      end do

      ! Compute conservative variables
      call compute_conserved_from_primitive(P_L, U_L)
      call compute_conserved_from_primitive(P_R, U_R)

    end if

  end subroutine

  ! ========================================================================

  real(kind=8) function phi_limiter(theta)

    real(kind=8), intent(in) :: theta ! ratio of consecutive gradients

    phi_limiter = (theta*theta + theta)/(theta*theta + 1.0) ! van Albada symmetric limiter
    ! phi_limiter = (theta + abs(theta))/(1.0 + abs(theta)) ! van Leer limiter

  end function

end module


text:
!==============================================================================!
  subroutine Control_Mod_Open_Root_File(file_name)
!------------------------------------------------------------------------------!
!   It rarely gets simpler than this - just opens a file for reading.          !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  character(len=*) :: file_name
!==============================================================================!

  ! Unused
  ! control_file_name = file_name

  call File_Mod_Open_File_For_Reading(file_name, root_control_file_unit)

  ! Make root default to begin with
  control_file_unit = root_control_file_unit

  ! Set default values for domain 1
  ! dom_control_file_name(1) = control_file_name
  dom_control_file_unit(1) = root_control_file_unit

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
C
C     This file is a part of HAXX
C
C     Copyright (c) 2017 David Williams-Young
C     All rights reserved.
C
C     See LICENSE.txt
C
      Subroutine HGEMVHH(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
     $  Y, INCY)
C
      LOGICAL LSAME
      EXTERNAL LSAME
C
      CHARACTER TRANS
      LOGICAL   NOCONJ, AZERO,AONE,BZERO,BONE
      INTEGER*4 M,N, LDA,INCX,INCY, LENX,LENY, I,J,IY,JX,IX,JY, KX,KY
      COMPLEX*16    A(2,LDA,*), X(2,*), Y(2,*), ALPHA(2), BETA(2)
C
      COMPLEX*16    HTMP1S, HTMP1J, HTMP2S, HTMP2J
      COMPLEX*16    HTMP3S, HTMP3J, HTMP4S, HTMP4J
C
      COMPLEX*16    ONE, ZERO
      PARAMETER (ONE=(1.0D+0,0.0D+0), ZERO=(0.0D+0,0.0D+0))
C
      NOCONJ = LSAME('T',TRANS)
      IF ( LSAME('T',TRANS) ) THEN
        LENX = N
        LENY = M
      ELSE
        LENX = M
        LENY = N
      ENDIF

      AZERO = (ALPHA(1).EQ.ZERO).AND.((ALPHA(2).EQ.ZERO))
      AONE  = (ALPHA(1).EQ.ONE).AND.((ALPHA(2).EQ.ZERO))
      BZERO = (BETA(1).EQ.ZERO).AND.((BETA(2).EQ.ZERO))
      BONE  = (BETA(1).EQ.ONE).AND.((BETA(2).EQ.ZERO))
C
      IF ( (M.EQ.0) .OR. (N.EQ.0) .OR.
     $     ( AZERO .AND. BONE ) ) RETURN
C
C
      IF ( INCX.GT.0 ) THEN
        KX = 1
      ELSE
        KX = 1 - (LENX-1)*INCX
      ENDIF

      IF ( INCY.GT.0 ) THEN
        KY = 1
      ELSE
        KY = 1 - (LENY-1)*INCY
      ENDIF
C

      IF ( .NOT.BONE ) THEN
        IF ( INCY.EQ.1 ) THEN
          IF ( BZERO ) THEN
            DO 10 I = 1,LENY
              Y(1,I) = ZERO
              Y(2,I) = ZERO
  10        CONTINUE
          ELSE
            DO 20 I = 1,LENY
              HTMP1S = BETA(1) * Y(1,I) - BETA(2) * CONJG(Y(2,I))
              Y(2,I) = BETA(1) * Y(2,I) + BETA(2) * CONJG(Y(1,I))
              Y(1,I) = HTMP1S
  20        CONTINUE
          ENDIF
        ELSE
          IY = KY
          IF ( BZERO ) THEN
            DO 30 I = 1,LENY
              Y(1,IY) = ZERO
              Y(2,IY) = ZERO
              IY = IY + INCY
  30        CONTINUE
          ELSE
            DO 40 I = 1,LENY
              HTMP1S = BETA(1) * Y(1,IY) - BETA(2) * CONJG(Y(2,IY))
              Y(2,IY) = BETA(1) * Y(2,IY) + BETA(2) * CONJG(Y(1,IY))
              Y(1,IY) = HTMP1S
              IY = IY + INCY
  40        CONTINUE
          ENDIF
        ENDIF
      ENDIF
C
      IF ( AZERO )  RETURN ! Nothing left to do
C
      IF ( LSAME(TRANS,'N') ) THEN
        JX = KX
        IF ( INCY.EQ.1 ) THEN
          DO 60 J = 1,N
            HTMP1S = ALPHA(1) * X(1,JX)
            HTMP1J = ALPHA(1) * X(2,JX)
            HTMP2S = ALPHA(1) * CONJG(X(1,JX))
            HTMP2J = ALPHA(1) * CONJG(X(2,JX))
            HTMP3S = ALPHA(2) * X(1,JX)
            HTMP3J = ALPHA(2) * X(2,JX)
            HTMP4S = ALPHA(2) * CONJG(X(1,JX))
            HTMP4J = ALPHA(2) * CONJG(X(2,JX))
            Do 50 I = 1,M
C             Hamilton Product
              Y(1,I) = Y(1,I) +
     $          A(1,I,J)*HTMP1S - A(2,I,J)*HTMP2J -
     $          CONJG(A(1,I,J))*HTMP4J - CONJG(A(2,I,J))*HTMP3S
              Y(2,I) = Y(2,I) +
     $          A(1,I,J)*HTMP1J + A(2,I,J)*HTMP2S +
     $          CONJG(A(1,I,J))*HTMP4S - CONJG(A(2,I,J))*HTMP3J
  50        CONTINUE
            JX = JX + INCX
  60      CONTINUE
        ELSE
          DO 80 J = 1,N
            HTMP1S = ALPHA(1) * X(1,JX)
            HTMP1J = ALPHA(1) * X(2,JX)
            HTMP2S = ALPHA(1) * CONJG(X(1,JX))
            HTMP2J = ALPHA(1) * CONJG(X(2,JX))
            HTMP3S = ALPHA(2) * X(1,JX)
            HTMP3J = ALPHA(2) * X(2,JX)
            HTMP4S = ALPHA(2) * CONJG(X(1,JX))
            HTMP4J = ALPHA(2) * CONJG(X(2,JX))
            IY = KY
            Do 70 I = 1,M
C             Hamilton Product
              Y(1,IY) = Y(1,IY) +
     $          A(1,I,J)*HTMP1S - A(2,I,J)*HTMP2J -
     $          CONJG(A(1,I,J))*HTMP4J - CONJG(A(2,I,J))*HTMP3S
              Y(2,IY) = Y(2,IY) +
     $          A(1,I,J)*HTMP1J + A(2,I,J)*CONJG(HTMP2S) +
     $          CONJG(A(1,I,J))*HTMP4S - CONJG(A(2,I,J))*HTMP3J
              IY = IY + INCY
  70        CONTINUE
            JX = JX + INCX
  80      CONTINUE
        ENDIF
      ELSE
        JY = KY
        IF ( INCX.EQ.1 ) THEN
          DO 110 J = 1,N
            HTMP1S = ZERO
            HTMP1J = ZERO
            IF ( NOCONJ ) THEN
              DO 90 I = 1,M
C             Hamilton Product
                HTMP1S = HTMP1S +
     $            A(1,I,J)*X(1,I) - A(2,I,J)*CONJG(X(2,I))
                HTMP1J = HTMP1J +
     $            A(1,I,J)*X(2,I) + A(2,I,J)*CONJG(X(1,I))
  90          CONTINUE
            ELSE
              DO 100 I = 1,M
C               Hamilton Product (IMPLIED CONJ(A))
                HTMP1S = HTMP1S +
     $            CONJG(A(1,I,J))*X(1,I) + A(2,I,J)*CONJG(X(2,I))
                HTMP1J = HTMP1J +
     $            CONJG(A(1,I,J))*X(2,I) - A(2,I,J)*CONJG(X(1,I))
  100         CONTINUE
            ENDIF
            Y(1,JY) = Y(1,JY) + ALPHA(1)*HTMP1S - ALPHA(2)*CONJG(HTMP1J)
            Y(2,JY) = Y(2,JY) + ALPHA(1)*HTMP1J + ALPHA(2)*CONJG(HTMP1S)
            JY = JY + INCY
  110     CONTINUE
        ELSE
          DO 140 J = 1,N
            HTMP1S = ZERO
            HTMP1J = ZERO
            IX = KX
            IF ( NOCONJ ) THEN
              DO 120 I = 1,M
C             Hamilton Product
                HTMP1S = HTMP1S +
     $            A(1,I,J)*X(1,IX) - A(2,I,J)*CONJG(X(2,IX))
                HTMP1J = HTMP1J +
     $            A(1,I,J)*X(2,IX) + A(2,I,J)*CONJG(X(1,IX))
                IX = IX + INCX
  120         CONTINUE
            ELSE
              DO 130 I = 1,M
C               Hamilton Product (IMPLIED CONJ(A))
                HTMP1S = HTMP1S +
     $            CONJG(A(1,I,J))*X(1,IX) + A(2,I,J)*CONJG(X(2,IX))
                HTMP1J = HTMP1J +
     $            CONJG(A(1,I,J))*X(2,IX) - A(2,I,J)*CONJG(X(1,IX))
                IX = IX + INCX
  130         CONTINUE
            ENDIF
            Y(1,JY) = Y(1,JY) + ALPHA(1)*HTMP1S - ALPHA(2)*CONJG(HTMP1J)
            Y(2,JY) = Y(2,JY) + ALPHA(1)*HTMP1J + ALPHA(2)*CONJG(HTMP1S)
            JY = JY + INCY
  140     CONTINUE
        ENDIF
      ENDIF
C
      END


text:
	SUBROUTINE GCLOSP  ( iret )
C************************************************************************
C* GCLOSP								*
C* 									*
C* This subroutine closes any intermediate plot file used by a		*
C* sequential access type device.  When plotting is continued, a new	*
C* file is opened.  No action is taken if the current device does not   *
C* use an intermediate plot file.				        *
C*									*
C* GCLOSP  ( IRET )							*
C*									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. Goodman/RDS	 5/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 5/88	Documentation				*
C* L. Williams /EAi      3/94   Removed blank comments from header      *
C* A. Hardy/GSC          6/98   Cleaned up prolog                       *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C
	INTEGER		isend (2)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FCLOSP
C
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  iret = ierr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak