text:
The November 2012 Presidential Election occured on Tuesday the 6th. President Barack Obama of the Democratic Party and challenger Mitt Romney of the Republican Party were the two leading candidates, and President Obama won. As usual, an assortment of controversial propositions also shared the ballot space.

This was the first election in which Californians could http://www.sos.ca.gov/elections/elections_vr.htm register to vote online  you had to register by October 22, 2012 to vote in this election. You could also check your http://www.yoloelections.org/regdb/ voter registration status and http://www.yoloelections.org/voting/polling_place polling place online, which is a good thing to do since polling places may have changed since the last time you voted. If you had chosen to vote by mail, you could also turn your mailin ballot in at any polling place in your county on Election Day or at the elections office in your county.

This is also the first election in which Davis had been assigned to the new districts that came about as a result of the passage of Proposition 11 in 2008. This had had an impact on the eligibility of some candidates. Perhaps most notably, Mike Thompson  could no longer run for reelection in our district, whereas John Garamendi could.

Measure E, Prop. 30, and Garamendi all won handily.  Oh, and Barack Obama.  Props. 32, 34, and 37 were rejected.  Apparently, Californians like schools, unions, GMOs, and governmentsponsored deaths.  78.8 percent of Davis residents cast a ballot this election, the highest percentage for any community in Yolo County.

Candidates and measures in bold were elected and approved, respectively.
 U.S. Senator

  Dianne Feinstein (Democratic Party)
  Elizabeth Emken (Republican Party)

 U.S. House of Representatives, 3rd District

  John Garamendi (Democratic Party).  53.7% of the vote.
  Kim Vann (Republican Party). 46.3% of the vote.

State Senate, 3rd District

 Lois Wolk (Democratic Party)
  Frank Miranda (Republican Party)

 State Assembly, 4th District

  Mariko Yamada (Democratic Party)
  John Munn (Republican Party)

State Propositions

  November 2012 Presidential Election/Proposition 30 Proposition 30  The Schools and Local Public Safety Protection Act
  November 2012 Presidential Election/Proposition 31 Proposition 31  Establishes twoyear state budget.
  November 2012 Presidential Election/Proposition 32 Proposition 32  Prohibits unions from using payrolldeducted funds for political purposes.
  November 2012 Presidential Election/Proposition 33 Proposition 33  Changes current law to allow insurance companies to set prices based on whether the driver previously carried auto insurance with any insurance company.
  November 2012 Presidential Election/Proposition 34 Proposition 34  Repeals death penalty and replaces it with life imprisonment without possibility of parole.
  November 2012 Presidential Election/Proposition 35 Proposition 35  Increases prison sentences and fines for human trafficking convictions. Note, this will unfairly make anyone convicted of Lewd Conduct register as a sex offender, including former prostitutes. http://againstthecaseact.com/ More info here.
  November 2012 Presidential Election/Proposition 36 Proposition 36  Revises law to impose life sentence only when new felony conviction is serious or violent.
  November 2012 Presidential Election/Proposition 37 Proposition 37   Mandatory Labeling of Genetically Engineered Food Initiative
  November 2012 Presidential Election/Proposition 38 Proposition 38  Increases taxes on earnings using sliding scale, for twelve years. Revenues go to K–12 schools and early childhood programs, and for four years to repaying state debt.
  November 2012 Presidential Election/Proposition 39 Proposition 39  Requires multistate businesses to pay income taxes based on percentage of their sales in California.
  November 2012 Presidential Election/Proposition 40 Proposition 40  A Yes vote approves, and a No vote rejects, new State Senate districts drawn by the Citizens Redistricting Commission.

 Local Measures

  November 2012 Presidential Election/Measure E Measure E  Davis Joint Unified School District Special Tax Proposal.  Passed with 68.9% of the vote. (2/3 required for passage)
  November 2012 Presidential Election/Measure H Measure H  A Yolo County Measure.  Creation of a Department of Finance.  Passed with 65.8% of the vote.

 Governing Board Member, Davis Joint Unified School District

Two seats were open.


!==============================================================================!
  subroutine Allocate_Memory(grid)
!------------------------------------------------------------------------------!
!   Allocates memory for the convertor.                                        !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Grid_Mod, only: Grid_Type,                &
                      Grid_Mod_Allocate_Nodes,  &
                      Grid_Mod_Allocate_Cells,  &
                      Grid_Mod_Allocate_Faces,  &
                      Grid_Mod_Allocate_New_Numbers
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer, parameter :: F = 3   ! workaround to allocate more memory for bnds
!==============================================================================!

  ! Allocate memory
  ! =--> carefull: there is no checking!
  call Grid_Mod_Allocate_Nodes(grid, grid % n_nodes)
  call Grid_Mod_Allocate_Cells(grid, grid % n_cells,   grid % n_bnd_cells * F)
  call Grid_Mod_Allocate_Faces(grid, grid % n_cells*5, 0)
  call Grid_Mod_Allocate_New_Numbers(grid,                    &
                                     grid % n_nodes,          &
                                     grid % n_bnd_cells * F,  &
                                     grid % n_cells,          &
                                     grid % n_faces)

  allocate(grid % bnd_cond % color(-grid % n_bnd_cells * F:-1))
  grid % bnd_cond % color = 0

  ! (Dirty) trick to allocate additional memory for cities
  grid % n_bnd_cells = grid % n_bnd_cells / F

  end subroutine


text:
C     LAST UPDATE 16/03/89
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE FFT1D (X,Y,TABLE,M,LL,ISN)
      IMPLICIT NONE
C
C Purpose:  FFT is in place DFT computation using Sande algorithm
C           and Markel pruning modification.
C
      REAL    X(1),Y(1),TABLE(1)
      INTEGER M,LL,ISN
C
C X      : Array of length 2**M used to hold real part of complex input
C Y      : Array of length 2**M used to hold imaginary part of input
C          X&Y are real and complex part of output
C TABLE  : Array of length N/4+1 (N=2**M) contains quarter length
C          cosine table.
C M      : Integer power of two which determines size of fft to be
C          performed. (bit reverse table is set for a max of N=2**12)
C LL     : Integer power of 2 which determines 2**ll actual data pts
C          and thus no. of stages in which no pruning is allowable
C ISN    :-1 for forward fft
C         +1 for reverse fft
C
C Calls   0:
C Called by: FRWFFT , INVFFT
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local variables:
C
      REAL    C,S,R,FI,T1,T2
      INTEGER L(12),L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,IARG
      INTEGER J,J1,J2,J3,J4,J5,J6,J7,J8,J9,J10,J11,JR,JN,K1,K2,K3
      INTEGER N,ND4,ND4P1,ND4P2,ND2P2,LLL,LM,LMM,LO,LI,LMX,LIX,ISCL
C
C L     :
C
      EQUIVALENCE (L12,L(1)) , (L11,L(2)) , (L10,L(3)) , (L9,L(4)) ,
     1            (L8,L(5))  , (L7,L(6))  , (L6,L(7))  , (L5,L(8)) ,
     1            (L4,L(9))  , (L3,L(10)) , (L2,L(11)) , (L1,L(12))
C
C-----------------------------------------------------------------------
      N=2**M
      ND4=N/4
      ND4P1=ND4+1
      ND4P2=ND4P1+1
      ND2P2=ND4+ND4P2
      LLL=2**LL
      DO 50 LO=1,M
         LMX=2**(M-LO)
         LMM=LMX
         LIX=2*LMX
         ISCL=N/LIX
C
C=======TEST FOR PRUNING
C
         IF (LO-M+LL) 10,20,20
  10     LMM=LLL
  20     DO 50 LM=1,LMM
            IARG=(LM-1)*ISCL+1
            IF (IARG.LE.ND4P1) GO TO 30
            K1=ND2P2-IARG
            C=-TABLE(K1)
            K3=IARG-ND4
            S=ISN*TABLE(K3)
            GO TO 40
  30        C=TABLE(IARG)
            K2=ND4P2-IARG
            S=ISN*TABLE(K2)
  40        DO 50 LI=LIX,N,LIX
               J1=LI-LIX+LM
               J2=J1+LMX
               T1=X(J1)-X(J2)
               T2=Y(J1)-Y(J2)
               X(J1)=X(J1)+X(J2)
               Y(J1)=Y(J1)+Y(J2)
               X(J2)=C*T1-S*T2
               Y(J2)=C*T2+S*T1
  50  CONTINUE
C
C========PERFORM BIT REVERSAL
C
      DO 70 J=1,12
         L(J)=1
         IF (J-M) 60,60,70
  60     L(J)=2**(M+1-J)
  70  CONTINUE
      JN=1
      DO 120 J1=1,L1
       DO 120 J2=J1,L2,L1
        DO 120 J3=J2,L3,L2
         DO 120 J4=J3,L4,L3
          DO 120 J5=J4,L5,L4
           DO 120 J6=J5,L6,L5
            DO 120 J7=J6,L7,L6
             DO 120 J8=J7,L8,L7
              DO 120 J9=J8,L9,L8
               DO 120 J10=J9,L10,L9
                DO 120 J11=J10,L11,L10
                 DO 120 JR=J11,L12,L11
                 IF (JN-JR) 80,80,90
  80             R=X(JN)
                 X(JN)=X(JR)
                 X(JR)=R
                 FI=Y(JN)
                 Y(JN)=Y(JR)
                 Y(JR)=FI
  90             IF (ISN) 110,110,100
  100            X(JR)=X(JR)/REAL(N)
                 Y(JR)=Y(JR)/REAL(N)
  110            JN=JN+1
  120 CONTINUE
      RETURN
      END

ext: f


text:
The Model United Nations Club at U.C Davis is an organization driven to simulate the United Nations in order to broaden understanding of how global political issues are discussed and solved. While it may at first seem like a club geared specifically towards political science or international relations students, Davis MUN has members representing a great variety of academic disciplines. Communication, psychology, animal science, and engineering majors all come together through a common passion for politics and international affairs. This organization is offered at many colleges and high schools around the world. The club offers a good way for members to practice speech and debate, as well as the development of research skills and a broader understanding of international relations.  We also hold fundraisers to subsidize the cost of attending these conferences, so itll be easier for more members to go!

Conferences
Davis M.U.N. attends several conferences hosted around the world every year. In previous years, we have sent delegates to:
 PAXMUN American West in Las Vegas, NV
 UC Berkeley Model UN Conference in San Francisco, CA
 McGill Model UN Conference in Montreal, Canada
 Harvard National Model UN in Boston, MA
 National Model United Nations in Manhattan, New York City
  UCLA Model United Nations Conference in LA

Davis Model United Nations Conference

For a M.U.N. conference, each attending member typically acts as a U.N. delegate representing a given U.N. member state within a given U.N. committee (assigned by the organization(s) hosting the conference). In accordance with the rules of parlimentary procedure of the U.N., the delegates then spend the conference committee sessions discussing international political issues that vary based on the committee. Through discussion and debate, the delegates hopefully come up with adequate solutions in the form of U.N. draft resolutions which are put to a vote by the end of the conference. Traditionally, Model U.N. conferences offer awards to the delegates that perform the best, with the Best Delegate award often in the form of a gavel.

Aside from broadening knowledge of international politics, sharpening public speaking and debate skills, or exercising diplomacy, M.U.N. Delegates are also extremely social. M.U.N. conferences offer plenty of opportunities for meeting students from all over the world who share a passion for politics. Whether you socialize between committee sessions or at delegate social events (which are scheduled at nearly every conference), M.U.N. is a wonderful way to forge connections across continents and gather perspective of international issues from very different cultural viewpoints. These people also often make wonderful friends.

Davis M.U.N. also staffs three conferences over the course of the year. During Fall Quarter, we send members to Diablo Valley College in Contra Costa County to staff a M.U.N. conference for high school students. This not only generates some revenue for the club, but gives our members an opportunity to act as the moderating dias members for an M.U.N. committee before attending a conference as a delegate. We also love sharing our love for international politics and debate with high school students, some of whom may hopefully become future Davis M.U.N. members.
The other two conferences we host ourselves: DCCMUN (Davis Collegiate Conference for Model United Nations) and DMUNC (Davis Model United Nations Conference) both take place right on the U.C. Davis Campus. DCCMUN is a small conference for college students, allowing our own members to be part of the debate.

The Davis Model United Nations Conference (also known as DMUNC) is an annual high school level United Nations simulation presented by MUNUCD each Spring in a weekend in the middle of May. High school delegates from all over California come to represent different countries in a simulation of world bodies discussing and debating world issues this conference also functions as the clubs primary source of revenue, this money is spent to help subsidize trips to various conferences such as conference fees as well as travel and lodging expenses.

We are looking for UC Davis students to help us run and participate in these conferences! Come to our club meetings and learn more! If you are interested in helping out by staffing, running a committee as a chair or just interested in participating, please contact the DCCMUN 2012 SecretaryGeneral Sana Sareshwala at sana.sareshwala@davismun.org or DMUNC 2013 SecretaryGeneral Megan Nguyen at ||MailTo(secgen AT davismun DOT org)|| for more information!

Any questions about the club contact us at sg@davismun.org or check out our website!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*------------------------------------------------------------------------------
*
* Filename             : TEST_HYDJET.F
*
*==============================================================================
*
* Description : Example program to simulate hadron spectra in AA collisions
*               at LHC with HYDJET-code (should be compiled with object files
*               obtained with hydjet1_1.f, pyquen1_1.f, pythia6401.f (or later
*               pythia versions) and jetset_73.f with extended array size of
*               common block LUJETS)
*
*==============================================================================

      double precision ckin,parp,pari
      real A,bmin,bmax,bfix
      external ludata,pydata
      common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
      common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
      common /hyflow/ ytfl,ylfl,fpart
      common /hyjpar/ nhsel,ptmin,njet
      common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
      common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
      common /pypars/ mstp(200),parp(200),msti(200),pari(200)
c      common /pawc/ hmemor(20000)
      save /lujets/,/hyflow/,/hyjpar/,/hyfpar/
      save /pysubs/,/pypars/,/pydat1/

* prepare hbook memory
c      call hlimit(20000)
* open hrout file to write histograms
c      call HROPEN(1,'HISTO','hydjet.hrout','N',1024,ISTAT)
* prepare hbook histograms
c      call hbook1(1,'dN/dy $',100,-10.,10.,0.)     ! rapidity
c      call hbook1(2,'dN/deta $',100,-10.,10.,0.)   ! pseudorapidity
c      call hbook1(3,'dN/dpt $',100,0.,10.,0.)      ! transverse momentum
c      call hbook1(4,'dN/dphi $',100,-3.15,3.15,0.) ! azimuthal angle
c      call hbarx(0)

* set initial beam parameters
      A=207.                        ! atomic weigth
      nh=20000                      ! mean soft multiplicity for central Pb-Pb
      ifb=0                         ! fixed impact parameter
      bfix=0.                       ! in nucleus radius units
c      ifb=1                         ! distribution over impact parameter
c      bmin=0.                       ! from 'bmin'
c      bmax=1.                      ! to 'bmax'

* set hydro parameters
* nhsel=0 - hydro (no jets), nhsel=1 - hydro + pythia jets, nhsel=2 - hydro +
* pyquen jets, nhsel=3 - pythia jets (no hydro), nhsel=4 - pyquen jets (no hydro)
      nhsel=2                        ! flag to include hard scatterings
      ylfl=5.                        ! maximum longitudinal flow rapidity
      ytfl=1.                        ! maximum transverse flow rapidity
      fpart=1.                       ! fraction of soft multiplicity proportional
                                      ! # of nucleons-participants
* set input PYTHIA parameters
      msel=1                        ! QCD-dijet production
      ptmin=10.
      ckin(3)=dble(ptmin)           ! minimum pt in initial hard sub-process
      mstp(51)=7                    ! CTEQ5M pdf
      mstp(81)=0                    ! pp multiple scattering off
      mstu(21)=1                    ! avoid stopping run
      paru(14)=1.d0                 ! tolerance parameter to adjust fragmentation

* set original (rounded) test values and its rms for current model parameters
      pta0=0.56
      eta0=0.
      dna0=29850.
* set initial test values and its rms
      ptam=0.
      ptrms=0.
      etam=0.
      etrms=0.
      dnam=0.
      dnrms=0.

* initialize PYTHIA for hard parton-parton scattering
      if(nhsel.ne.0) call pyinit('CMS','p','p',5500.D0)

* set number of generated events
      ntot=10

      do ne=1,ntot                        ! cycle on events
       call hydro(A,ifb,bmin,bmax,bfix,nh)! single event generation

       call luedit(2)                     ! remove unstable particles and partons

       do ip=1,n                          ! cycle on particles
        pt=plu(ip,10)                     ! transverse momentum...
        ycm=plu(ip,17)                    ! rapidity...
        eta=plu(ip,19)                    ! pseudorapidity...
        phi=plu(ip,15)                    ! azimuthal angle...
        charge=plu(ip,6)                  ! electric charge...

* add current test values of eta, pt and its rms
        etam=etam+eta
	etrms=etrms+(eta-eta0)**2
	ptam=ptam+pt
	ptrms=ptrms+(pt-pta0)**2

* fill histograms for charged particles
c        if(abs(charge).gt.0.) then
c         call hfill(1,ycm,0.,1.)          ! rapidity
c         call hfill(2,eta,0.,1.)          ! pseudorapidity
c         call hfill(3,pt,0.,1.)           ! transverse momentum
c         call hfill(4,phi,0.,1.)          ! azimuthal angle
c        end if
       end do
       write(6,*) 'Event #',ne
       write(6,*) 'Impact parameter',bgen,'*RA',' Total multiplicity',n
       write(6,*) 'Pt hard min',ptmin,' GeV','   Ndijets',njet
       write(6,*) '***************************************************'

* add current test value of event multiplicity and its rms
       dnam=dnam+n
       dnrms=dnrms+(n-dna0)**2
      end do

* test calculating and printing of original "true" (rounded) numbers
* and generated one's (with statistical errors)
      etam=etam/dnam
      etrms=sqrt(etrms)/dnam
      ptam=ptam/dnam
      ptrms=sqrt(ptrms)/dnam
      dnam=dnam/float(ntot)
      dnrms=sqrt(dnrms)/float(ntot)
      write(6,1) dna0
1     format(2x,'True (rounded) mean multiplicity =',f7.0)
      write(6,2) dnam, dnrms
2     format(2x,'Generated mean multiplicity      =',f7.0,3x,
     > '+-  ',f6.0)
      write(6,3) eta0
3     format(2x,'True (rounded) mean pseudorapidity =',f5.2)
       write(6,4) etam, etrms
4     format(2x,'Generated mean pseudorapidity      =',f5.2,3x,
     > '+-  ',f5.2)
      write(6,5) pta0
5     format(2x,'True (rounded) mean transverse momentum =',f5.2)
      write(6,6) ptam, ptrms
6     format(2x,'Generated mean transverse momentum      =',f5.2,3x,
     > '+-  ',f5.2)

* finish histograms writing procedure
c      call hidopt(0,'ERRO')
c      call histdo
c      CALL HROUT(0,ICYCLE,' ')
c      CALL HREND('HISTO')

      end
*******************************************************************************

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot


! Estimates the refreezing of liquid water content based on Illangasekare et al., 1990
subroutine refreezing(volume_water_frozen, theta_1, theta_2, porosity, Vol, T_1, T_2, density, dry_density)
use declarations, only : rho_i, Cpi, rho_w, Cpw, rho_a, Cpa, dp, Lf
implicit none
real(dp), intent(inout) :: volume_water_frozen, &     ! Variable used to compute mass conservation [m3]
                       &   dry_density, &             ! Dry snow density [kg/m3]
                       &   density, &                 ! Bulk snow density [kg/m3]
                       &   porosity                   ! Porosity
real(dp), intent(in):: Vol, &             ! Volume of a cell [m3]
                   &   theta_1, &         ! Initial water content (before refreezing)
                   &   T_1                ! Initial temperature before refreezing [C]
real(dp), intent(out):: theta_2, &        ! Water content after refreezing
                    &   T_2               ! Temperature after refreezing [C]
real(dp) :: ice_content, & ! Ice content of each cell [m3/m3]
        &   rhoCp, &       ! rho*Cp for each snow cell [J/K/m3]
        &   freeze_max     ! Maximum liquid water content that can refreeze, same units as water content

if (theta_1 .gt. 0._dp) then

   rhoCp = rho_i * Cpi * (1._dp - porosity) + rho_w * Cpw * theta_1 + &
      & rho_a * Cpa * (porosity - theta_1)

   freeze_max = (-rhoCp*T_1)/(Lf*rho_w)   ! same units as water content

   if (freeze_max .lt. theta_1) then  ! not all liquid water content refreezes

     T_2 = 0._dp
     theta_2 = theta_1-freeze_max
     volume_water_frozen = volume_water_frozen + freeze_max * Vol
     ice_content = (1.0_dp-porosity) + freeze_max*rho_w/rho_i
     porosity = 1.0_dp - ice_content
     density = ice_content*rho_i + rho_w*theta_1 + rho_a * (porosity - theta_1)
     dry_density = ice_content*rho_i

   elseif (freeze_max .ge. theta_1) then ! all liquid water content refreezes,  i.e. freeze = water_content
      theta_2 = 0.0_dp
      volume_water_frozen = volume_water_frozen + theta_1 * Vol
      T_2 = theta_1* Lf * rho_w / rhoCp + T_1
      ice_content = (1.0_dp-porosity) + theta_1*rho_w/rho_i
      porosity = 1.0_dp-ice_content
      density = ice_content*rho_i + rho_a * porosity
      dry_density = ice_content*rho_i

   endif
else
   T_2 = T_1
   theta_2 = theta_1
endif

end subroutine


submodule(euler_interface_m) euler_prob_0047_m
    implicit none

contains

    module character(len=20) function euler0047()
        write (euler0047, "(i20)") ans(1000000)
    end function euler0047

    integer function ans(n)
        integer, intent(in) :: n
        integer :: i, j, n_factor(n), goal(4)

        n_factor = 0

        outer: do i = 2, n - 1

            if (n_factor(i) == 0) then
                inner: do j = 2*i, n - 1, i
                    n_factor(j) = n_factor(j) + 1
                end do inner
            end if

        end do outer

        goal = [4, 4, 4, 4]

        do i = 2, n - 1
            if (all(n_factor(i:i + 3) == goal(:))) exit
        end do

        ans = i
    end function ans
end submodule euler_prob_0047_m


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! The ODE Jacobian of Chemical Model File
!
! Generated by KPP-2.1 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : redhc_made_soa_vbs_Jacobian.f90
! Time                 : Wed Oct 28 17:42:24 2020
! Working directory    : /mnt/lfs1/BMC/fim/Jordan/models/rapid-refresh/WRFV3.9_RAPCHEM_v7/chem/KPP/mechanisms/redhc_made_soa_vbs
! Equation file        : redhc_made_soa_vbs.kpp
! Output root filename : redhc_made_soa_vbs
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE redhc_made_soa_vbs_Jacobian

  USE redhc_made_soa_vbs_Parameters
  USE redhc_made_soa_vbs_JacobianSP

  IMPLICIT NONE

CONTAINS

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Jac_SP_Vec - function for sparse multiplication: sparse Jacobian times vector
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      UV        - User vector for variables
!      JUV       - Jacobian times user vector
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Jac_SP_Vec ( JVS, UV, JUV )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! UV - User vector for variables
  REAL(kind=dp) :: UV(NVAR)
! JUV - Jacobian times user vector
  REAL(kind=dp) :: JUV(NVAR)

  JUV(1) = JVS(1)*UV(1)+JVS(2)*UV(43)
  JUV(2) = JVS(3)*UV(2)+JVS(4)*UV(43)
  JUV(3) = JVS(5)*UV(3)+JVS(6)*UV(4)+JVS(7)*UV(43)
  JUV(4) = JVS(8)*UV(4)+JVS(9)*UV(5)+JVS(10)*UV(43)
  JUV(5) = JVS(11)*UV(5)+JVS(12)*UV(6)+JVS(13)*UV(43)
  JUV(6) = JVS(14)*UV(6)+JVS(15)*UV(43)
  JUV(7) = JVS(16)*UV(7)+JVS(17)*UV(8)+JVS(18)*UV(15)+JVS(19)*UV(43)
  JUV(8) = JVS(20)*UV(8)+JVS(21)*UV(9)+JVS(22)*UV(14)+JVS(23)*UV(43)
  JUV(9) = JVS(24)*UV(9)+JVS(25)*UV(10)+JVS(26)*UV(13)+JVS(27)*UV(43)
  JUV(10) = JVS(28)*UV(10)+JVS(29)*UV(11)+JVS(30)*UV(12)+JVS(31)*UV(43)
  JUV(11) = JVS(32)*UV(11)+JVS(33)*UV(43)
  JUV(12) = JVS(34)*UV(12)+JVS(35)*UV(43)
  JUV(13) = JVS(36)*UV(13)+JVS(37)*UV(43)
  JUV(14) = JVS(38)*UV(14)+JVS(39)*UV(43)
  JUV(15) = JVS(40)*UV(15)+JVS(41)*UV(43)
  JUV(16) = JVS(42)*UV(16)+JVS(43)*UV(17)+JVS(44)*UV(43)
  JUV(17) = JVS(45)*UV(17)+JVS(46)*UV(43)
  JUV(18) = JVS(47)*UV(18)+JVS(48)*UV(43)
  JUV(19) = JVS(49)*UV(19)+JVS(50)*UV(43)+JVS(51)*UV(44)
  JUV(20) = JVS(52)*UV(20)+JVS(53)*UV(46)+JVS(54)*UV(47)
  JUV(21) = JVS(55)*UV(21)+JVS(56)*UV(46)+JVS(57)*UV(48)
  JUV(22) = JVS(58)*UV(22)+JVS(59)*UV(24)+JVS(60)*UV(34)+JVS(61)*UV(37)+JVS(62)*UV(42)+JVS(63)*UV(43)+JVS(64)*UV(48)
  JUV(23) = JVS(65)*UV(23)+JVS(66)*UV(43)
  JUV(24) = JVS(67)*UV(24)+JVS(68)*UV(34)+JVS(69)*UV(43)
  JUV(25) = JVS(70)*UV(25)+JVS(71)*UV(40)+JVS(72)*UV(41)+JVS(73)*UV(43)+JVS(74)*UV(44)+JVS(75)*UV(47)
  JUV(26) = JVS(76)*UV(26)+JVS(77)*UV(43)+JVS(78)*UV(44)+JVS(79)*UV(46)
  JUV(27) = JVS(80)*UV(27)+JVS(81)*UV(42)+JVS(82)*UV(43)
  JUV(28) = JVS(83)*UV(28)+JVS(84)*UV(38)+JVS(85)*UV(39)+JVS(86)*UV(43)+JVS(87)*UV(46)+JVS(88)*UV(48)
  JUV(29) = JVS(89)*UV(23)+JVS(90)*UV(27)+JVS(91)*UV(29)+JVS(92)*UV(30)+JVS(93)*UV(32)+JVS(94)*UV(35)+JVS(95)*UV(36)&
              &+JVS(96)*UV(37)+JVS(97)*UV(38)+JVS(98)*UV(39)+JVS(99)*UV(42)+JVS(100)*UV(43)+JVS(101)*UV(48)
  JUV(30) = JVS(102)*UV(23)+JVS(103)*UV(30)+JVS(104)*UV(33)+JVS(105)*UV(36)+JVS(106)*UV(42)+JVS(107)*UV(43)+JVS(108)&
              &*UV(48)
  JUV(31) = JVS(109)*UV(31)+JVS(110)*UV(32)+JVS(111)*UV(35)+JVS(112)*UV(36)+JVS(113)*UV(40)+JVS(114)*UV(43)+JVS(115)&
              &*UV(45)+JVS(116)*UV(48)
  JUV(32) = JVS(117)*UV(32)+JVS(118)*UV(42)+JVS(119)*UV(43)+JVS(120)*UV(48)
  JUV(33) = JVS(121)*UV(33)+JVS(122)*UV(42)+JVS(123)*UV(43)+JVS(124)*UV(48)
  JUV(34) = JVS(125)*UV(22)+JVS(127)*UV(32)+JVS(128)*UV(33)+JVS(129)*UV(34)+JVS(130)*UV(35)+JVS(131)*UV(36)+JVS(133)&
              &*UV(42)+JVS(134)*UV(43)
  JUV(35) = JVS(136)*UV(35)+JVS(137)*UV(42)+JVS(138)*UV(43)+JVS(139)*UV(48)
  JUV(36) = JVS(140)*UV(36)+JVS(141)*UV(42)+JVS(142)*UV(43)+JVS(143)*UV(48)
  JUV(37) = JVS(144)*UV(33)+JVS(145)*UV(36)+JVS(146)*UV(37)+JVS(147)*UV(42)+JVS(148)*UV(43)+JVS(149)*UV(48)
  JUV(38) = JVS(150)*UV(24)+JVS(151)*UV(27)+JVS(152)*UV(32)+JVS(153)*UV(33)+JVS(154)*UV(34)+JVS(155)*UV(35)+JVS(156)&
              &*UV(36)+JVS(157)*UV(37)+JVS(158)*UV(38)+JVS(159)*UV(42)+JVS(160)*UV(43)+JVS(161)*UV(48)
  JUV(39) = JVS(162)*UV(27)+JVS(163)*UV(32)+JVS(164)*UV(33)+JVS(165)*UV(35)+JVS(166)*UV(36)+JVS(167)*UV(37)+JVS(168)&
              &*UV(38)+JVS(169)*UV(39)+JVS(170)*UV(42)+JVS(171)*UV(43)+JVS(172)*UV(44)+JVS(173)*UV(45)+JVS(174)*UV(47)&
              &+JVS(175)*UV(48)
  JUV(40) = JVS(176)*UV(23)+JVS(177)*UV(32)+JVS(178)*UV(33)+JVS(179)*UV(34)+JVS(180)*UV(35)+JVS(181)*UV(36)+JVS(182)&
              &*UV(37)+JVS(183)*UV(40)+JVS(184)*UV(42)+JVS(185)*UV(43)+JVS(186)*UV(44)+JVS(187)*UV(45)+JVS(188)*UV(48)
  JUV(41) = JVS(189)*UV(23)+JVS(190)*UV(24)+JVS(191)*UV(25)+JVS(192)*UV(27)+JVS(193)*UV(30)+JVS(194)*UV(31)+JVS(195)&
              &*UV(32)+JVS(196)*UV(33)+JVS(197)*UV(34)+JVS(198)*UV(35)+JVS(199)*UV(36)+JVS(200)*UV(37)+JVS(201)*UV(38)&
              &+JVS(203)*UV(41)+JVS(204)*UV(42)+JVS(205)*UV(43)+JVS(206)*UV(44)+JVS(207)*UV(45)+JVS(208)*UV(47)+JVS(209)&
              &*UV(48)
  JUV(42) = JVS(210)*UV(27)+JVS(211)*UV(32)+JVS(212)*UV(33)+JVS(213)*UV(35)+JVS(214)*UV(36)+JVS(215)*UV(37)+JVS(216)&
              &*UV(42)+JVS(217)*UV(43)+JVS(218)*UV(44)+JVS(219)*UV(45)+JVS(220)*UV(46)+JVS(221)*UV(48)
  JUV(43) = JVS(222)*UV(17)+JVS(223)*UV(18)+JVS(224)*UV(19)+JVS(225)*UV(23)+JVS(226)*UV(25)+JVS(227)*UV(26)+JVS(228)&
              &*UV(27)+JVS(229)*UV(28)+JVS(230)*UV(29)+JVS(231)*UV(30)+JVS(232)*UV(31)+JVS(233)*UV(32)+JVS(234)*UV(33)&
              &+JVS(235)*UV(34)+JVS(236)*UV(35)+JVS(237)*UV(36)+JVS(238)*UV(37)+JVS(239)*UV(38)+JVS(240)*UV(39)+JVS(243)&
              &*UV(42)+JVS(244)*UV(43)+JVS(245)*UV(44)+JVS(246)*UV(45)+JVS(247)*UV(46)+JVS(248)*UV(47)+JVS(249)*UV(48)
  JUV(44) = JVS(250)*UV(19)+JVS(251)*UV(23)+JVS(252)*UV(24)+JVS(253)*UV(26)+JVS(254)*UV(27)+JVS(255)*UV(29)+JVS(256)&
              &*UV(30)+JVS(257)*UV(31)+JVS(258)*UV(32)+JVS(260)*UV(34)+JVS(261)*UV(35)+JVS(262)*UV(36)+JVS(263)*UV(37)&
              &+JVS(264)*UV(38)+JVS(265)*UV(39)+JVS(266)*UV(40)+JVS(267)*UV(41)+JVS(268)*UV(42)+JVS(269)*UV(43)+JVS(270)&
              &*UV(44)+JVS(271)*UV(45)+JVS(272)*UV(46)+JVS(273)*UV(47)+JVS(274)*UV(48)
  JUV(45) = JVS(275)*UV(40)+JVS(276)*UV(41)+JVS(277)*UV(42)+JVS(279)*UV(44)+JVS(280)*UV(45)+JVS(281)*UV(46)+JVS(282)&
              &*UV(47)+JVS(283)*UV(48)
  JUV(46) = JVS(284)*UV(20)+JVS(285)*UV(21)+JVS(286)*UV(26)+JVS(287)*UV(28)+JVS(288)*UV(31)+JVS(290)*UV(35)+JVS(291)&
              &*UV(36)+JVS(292)*UV(37)+JVS(296)*UV(41)+JVS(297)*UV(42)+JVS(298)*UV(43)+JVS(299)*UV(44)+JVS(300)*UV(45)&
              &+JVS(301)*UV(46)+JVS(302)*UV(47)+JVS(303)*UV(48)
  JUV(47) = JVS(304)*UV(20)+JVS(305)*UV(23)+JVS(306)*UV(30)+JVS(308)*UV(35)+JVS(309)*UV(36)+JVS(310)*UV(38)+JVS(311)&
              &*UV(42)+JVS(312)*UV(43)+JVS(313)*UV(44)+JVS(314)*UV(45)+JVS(315)*UV(46)+JVS(316)*UV(47)+JVS(317)*UV(48)
  JUV(48) = JVS(318)*UV(21)+JVS(319)*UV(26)+JVS(320)*UV(28)+JVS(321)*UV(32)+JVS(322)*UV(33)+JVS(323)*UV(35)+JVS(324)&
              &*UV(36)+JVS(325)*UV(37)+JVS(326)*UV(38)+JVS(327)*UV(39)+JVS(328)*UV(42)+JVS(329)*UV(43)+JVS(330)*UV(44)&
              &+JVS(331)*UV(45)+JVS(332)*UV(46)+JVS(334)*UV(48)

END SUBROUTINE Jac_SP_Vec

! End of Jac_SP_Vec function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! JacTR_SP_Vec - sparse multiplication: sparse Jacobian transposed times vector
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      UV        - User vector for variables
!      JTUV      - Jacobian transposed times user vector
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE JacTR_SP_Vec ( JVS, UV, JTUV )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! UV - User vector for variables
  REAL(kind=dp) :: UV(NVAR)
! JTUV - Jacobian transposed times user vector
  REAL(kind=dp) :: JTUV(NVAR)

  JTUV(1) = JVS(1)*UV(1)
  JTUV(2) = JVS(3)*UV(2)
  JTUV(3) = JVS(5)*UV(3)
  JTUV(4) = JVS(6)*UV(3)+JVS(8)*UV(4)
  JTUV(5) = JVS(9)*UV(4)+JVS(11)*UV(5)
  JTUV(6) = JVS(12)*UV(5)+JVS(14)*UV(6)
  JTUV(7) = JVS(16)*UV(7)
  JTUV(8) = JVS(17)*UV(7)+JVS(20)*UV(8)
  JTUV(9) = JVS(21)*UV(8)+JVS(24)*UV(9)
  JTUV(10) = JVS(25)*UV(9)+JVS(28)*UV(10)
  JTUV(11) = JVS(29)*UV(10)+JVS(32)*UV(11)
  JTUV(12) = JVS(30)*UV(10)+JVS(34)*UV(12)
  JTUV(13) = JVS(26)*UV(9)+JVS(36)*UV(13)
  JTUV(14) = JVS(22)*UV(8)+JVS(38)*UV(14)
  JTUV(15) = JVS(18)*UV(7)+JVS(40)*UV(15)
  JTUV(16) = JVS(42)*UV(16)
  JTUV(17) = JVS(43)*UV(16)+JVS(45)*UV(17)+JVS(222)*UV(43)
  JTUV(18) = JVS(47)*UV(18)+JVS(223)*UV(43)
  JTUV(19) = JVS(49)*UV(19)+JVS(224)*UV(43)+JVS(250)*UV(44)
  JTUV(20) = JVS(52)*UV(20)+JVS(284)*UV(46)+JVS(304)*UV(47)
  JTUV(21) = JVS(55)*UV(21)+JVS(285)*UV(46)+JVS(318)*UV(48)
  JTUV(22) = JVS(58)*UV(22)+JVS(125)*UV(34)
  JTUV(23) = JVS(65)*UV(23)+JVS(89)*UV(29)+JVS(102)*UV(30)+JVS(176)*UV(40)+JVS(189)*UV(41)+JVS(225)*UV(43)+JVS(251)&
               &*UV(44)+JVS(305)*UV(47)
  JTUV(24) = JVS(59)*UV(22)+JVS(67)*UV(24)+JVS(150)*UV(38)+JVS(190)*UV(41)+JVS(252)*UV(44)
  JTUV(25) = JVS(70)*UV(25)+JVS(191)*UV(41)+JVS(226)*UV(43)
  JTUV(26) = JVS(76)*UV(26)+JVS(227)*UV(43)+JVS(253)*UV(44)+JVS(286)*UV(46)+JVS(319)*UV(48)
  JTUV(27) = JVS(80)*UV(27)+JVS(90)*UV(29)+JVS(151)*UV(38)+JVS(162)*UV(39)+JVS(192)*UV(41)+JVS(210)*UV(42)+JVS(228)&
               &*UV(43)+JVS(254)*UV(44)
  JTUV(28) = JVS(83)*UV(28)+JVS(229)*UV(43)+JVS(287)*UV(46)+JVS(320)*UV(48)
  JTUV(29) = JVS(91)*UV(29)+JVS(230)*UV(43)+JVS(255)*UV(44)
  JTUV(30) = JVS(92)*UV(29)+JVS(103)*UV(30)+JVS(193)*UV(41)+JVS(231)*UV(43)+JVS(256)*UV(44)+JVS(306)*UV(47)
  JTUV(31) = JVS(109)*UV(31)+JVS(194)*UV(41)+JVS(232)*UV(43)+JVS(257)*UV(44)+JVS(288)*UV(46)
  JTUV(32) = JVS(93)*UV(29)+JVS(110)*UV(31)+JVS(117)*UV(32)+JVS(127)*UV(34)+JVS(152)*UV(38)+JVS(163)*UV(39)+JVS(177)&
               &*UV(40)+JVS(195)*UV(41)+JVS(211)*UV(42)+JVS(233)*UV(43)+JVS(258)*UV(44)+JVS(321)*UV(48)
  JTUV(33) = JVS(104)*UV(30)+JVS(121)*UV(33)+JVS(128)*UV(34)+JVS(144)*UV(37)+JVS(153)*UV(38)+JVS(164)*UV(39)+JVS(178)&
               &*UV(40)+JVS(196)*UV(41)+JVS(212)*UV(42)+JVS(234)*UV(43)+JVS(322)*UV(48)
  JTUV(34) = JVS(60)*UV(22)+JVS(68)*UV(24)+JVS(129)*UV(34)+JVS(154)*UV(38)+JVS(179)*UV(40)+JVS(197)*UV(41)+JVS(235)&
               &*UV(43)+JVS(260)*UV(44)
  JTUV(35) = JVS(94)*UV(29)+JVS(111)*UV(31)+JVS(130)*UV(34)+JVS(136)*UV(35)+JVS(155)*UV(38)+JVS(165)*UV(39)+JVS(180)&
               &*UV(40)+JVS(198)*UV(41)+JVS(213)*UV(42)+JVS(236)*UV(43)+JVS(261)*UV(44)+JVS(290)*UV(46)+JVS(308)*UV(47)&
               &+JVS(323)*UV(48)
  JTUV(36) = JVS(95)*UV(29)+JVS(105)*UV(30)+JVS(112)*UV(31)+JVS(131)*UV(34)+JVS(140)*UV(36)+JVS(145)*UV(37)+JVS(156)&
               &*UV(38)+JVS(166)*UV(39)+JVS(181)*UV(40)+JVS(199)*UV(41)+JVS(214)*UV(42)+JVS(237)*UV(43)+JVS(262)*UV(44)&
               &+JVS(291)*UV(46)+JVS(309)*UV(47)+JVS(324)*UV(48)
  JTUV(37) = JVS(61)*UV(22)+JVS(96)*UV(29)+JVS(146)*UV(37)+JVS(157)*UV(38)+JVS(167)*UV(39)+JVS(182)*UV(40)+JVS(200)&
               &*UV(41)+JVS(215)*UV(42)+JVS(238)*UV(43)+JVS(263)*UV(44)+JVS(292)*UV(46)+JVS(325)*UV(48)
  JTUV(38) = JVS(84)*UV(28)+JVS(97)*UV(29)+JVS(158)*UV(38)+JVS(168)*UV(39)+JVS(201)*UV(41)+JVS(239)*UV(43)+JVS(264)&
               &*UV(44)+JVS(310)*UV(47)+JVS(326)*UV(48)
  JTUV(39) = JVS(85)*UV(28)+JVS(98)*UV(29)+JVS(169)*UV(39)+JVS(240)*UV(43)+JVS(265)*UV(44)+JVS(327)*UV(48)
  JTUV(40) = JVS(71)*UV(25)+JVS(113)*UV(31)+JVS(183)*UV(40)+JVS(266)*UV(44)+JVS(275)*UV(45)
  JTUV(41) = JVS(72)*UV(25)+JVS(203)*UV(41)+JVS(267)*UV(44)+JVS(276)*UV(45)+JVS(296)*UV(46)
  JTUV(42) = JVS(62)*UV(22)+JVS(81)*UV(27)+JVS(99)*UV(29)+JVS(106)*UV(30)+JVS(118)*UV(32)+JVS(122)*UV(33)+JVS(133)&
               &*UV(34)+JVS(137)*UV(35)+JVS(141)*UV(36)+JVS(147)*UV(37)+JVS(159)*UV(38)+JVS(170)*UV(39)+JVS(184)*UV(40)&
               &+JVS(204)*UV(41)+JVS(216)*UV(42)+JVS(243)*UV(43)+JVS(268)*UV(44)+JVS(277)*UV(45)+JVS(297)*UV(46)+JVS(311)&
               &*UV(47)+JVS(328)*UV(48)
  JTUV(43) = JVS(2)*UV(1)+JVS(4)*UV(2)+JVS(7)*UV(3)+JVS(10)*UV(4)+JVS(13)*UV(5)+JVS(15)*UV(6)+JVS(19)*UV(7)+JVS(23)&
               &*UV(8)+JVS(27)*UV(9)+JVS(31)*UV(10)+JVS(33)*UV(11)+JVS(35)*UV(12)+JVS(37)*UV(13)+JVS(39)*UV(14)+JVS(41)&
               &*UV(15)+JVS(44)*UV(16)+JVS(46)*UV(17)+JVS(48)*UV(18)+JVS(50)*UV(19)+JVS(63)*UV(22)+JVS(66)*UV(23)+JVS(69)&
               &*UV(24)+JVS(73)*UV(25)+JVS(77)*UV(26)+JVS(82)*UV(27)+JVS(86)*UV(28)+JVS(100)*UV(29)+JVS(107)*UV(30)+JVS(114)&
               &*UV(31)+JVS(119)*UV(32)+JVS(123)*UV(33)+JVS(134)*UV(34)+JVS(138)*UV(35)+JVS(142)*UV(36)+JVS(148)*UV(37)&
               &+JVS(160)*UV(38)+JVS(171)*UV(39)+JVS(185)*UV(40)+JVS(205)*UV(41)+JVS(217)*UV(42)+JVS(244)*UV(43)+JVS(269)&
               &*UV(44)+JVS(298)*UV(46)+JVS(312)*UV(47)+JVS(329)*UV(48)
  JTUV(44) = JVS(51)*UV(19)+JVS(74)*UV(25)+JVS(78)*UV(26)+JVS(172)*UV(39)+JVS(186)*UV(40)+JVS(206)*UV(41)+JVS(218)&
               &*UV(42)+JVS(245)*UV(43)+JVS(270)*UV(44)+JVS(279)*UV(45)+JVS(299)*UV(46)+JVS(313)*UV(47)+JVS(330)*UV(48)
  JTUV(45) = JVS(115)*UV(31)+JVS(173)*UV(39)+JVS(187)*UV(40)+JVS(207)*UV(41)+JVS(219)*UV(42)+JVS(246)*UV(43)+JVS(271)&
               &*UV(44)+JVS(280)*UV(45)+JVS(300)*UV(46)+JVS(314)*UV(47)+JVS(331)*UV(48)
  JTUV(46) = JVS(53)*UV(20)+JVS(56)*UV(21)+JVS(79)*UV(26)+JVS(87)*UV(28)+JVS(220)*UV(42)+JVS(247)*UV(43)+JVS(272)*UV(44)&
               &+JVS(281)*UV(45)+JVS(301)*UV(46)+JVS(315)*UV(47)+JVS(332)*UV(48)
  JTUV(47) = JVS(54)*UV(20)+JVS(75)*UV(25)+JVS(174)*UV(39)+JVS(208)*UV(41)+JVS(248)*UV(43)+JVS(273)*UV(44)+JVS(282)&
               &*UV(45)+JVS(302)*UV(46)+JVS(316)*UV(47)
  JTUV(48) = JVS(57)*UV(21)+JVS(64)*UV(22)+JVS(88)*UV(28)+JVS(101)*UV(29)+JVS(108)*UV(30)+JVS(116)*UV(31)+JVS(120)&
               &*UV(32)+JVS(124)*UV(33)+JVS(139)*UV(35)+JVS(143)*UV(36)+JVS(149)*UV(37)+JVS(161)*UV(38)+JVS(175)*UV(39)&
               &+JVS(188)*UV(40)+JVS(209)*UV(41)+JVS(221)*UV(42)+JVS(249)*UV(43)+JVS(274)*UV(44)+JVS(283)*UV(45)+JVS(303)&
               &*UV(46)+JVS(317)*UV(47)+JVS(334)*UV(48)

END SUBROUTINE JacTR_SP_Vec

! End of JacTR_SP_Vec function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

END MODULE redhc_made_soa_vbs_Jacobian


text:
subroutine lmn_outtpo(a)
!-----------------------------------------------------------------------
!
! This routine tracks points for the LowMach problem.
!
!-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_LowMach
   implicit none

   class(LowMachProblem) :: a

   integer(ip) :: ndime,npoin
   real(rp) :: tempe(1_ip,a%nptra), veloc(3_ip,a%nptra), press(1_ip,a%nptra)

   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)

   if(a%nptra>0) then
      call a%TrackingInterpolator%Interpolate(1,a%tempe(1:npoin,1),tempe)
      call a%TrackingInterpolator%Interpolate(ndime,a%veloc(:,:,1),veloc(1:ndime,:))
      call a%TrackingInterpolator%Interpolate(1,a%press(1:npoin,1),press)
      if (a%MPIrank == a%MPIroot) then
         write(a%lun_trap,20)(a%ctime),(tempe(1_ip,1:a%nptra)),(press(1_ip,1:a%nptra)),(veloc(1:ndime,1:a%nptra))
         if (a%kfl_flush == 1) call flush(a%lun_trap)
      endif
  endif

   20 format(21(1x,e14.7))
  return

end subroutine lmn_outtpo

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
#include "eiscor.h"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! z_upr1fact_chasedown
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine chases the a single misfit core transformation down
! one row in a factored unitary plus rank one (upr1fact) matrix.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! INPUT VARIABLES:
!
!  VEC             LOGICAL
!                    .TRUE. update schurvectors
!                    .FALSE. no schurvectors
!
!  P               LOGICAL array of dimension (2)
!                    array of position flags for Q
!
!  Q               REAL(8) array of dimension (6)
!                    array of generators for givens rotations
!
!  D               REAL(8) arrays of dimension (4)
!                    array of generators for complex diagonal matrix
!
!  C,B             REAL(8) arrays of dimension (6)
!                    array of generators for upper-triangular part
!
!  M               INTEGER
!                    leading dimesnion of V and W
!
!  V               COMPLEX(8) array of dimension (M,2)
!                    right schur vectors
!                    if VEC = .FALSE. unused
!                    if VEC = .TRUE. update V to store right schurvectors
!
!  MISFIT          REAL(8) array of dimension (3)
!                    array of generators for misfit
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine z_upr1fact_chasedown(VEC,P,Q,D,C,B,M,V,MISFIT)

  implicit none

  ! input variables
  logical, intent(in) :: VEC
  integer, intent(in) :: M
  logical, intent(inout) :: P(2)
  real(8), intent(inout) :: Q(6), D(4), C(6), B(6), MISFIT(3)
  complex(8), intent(inout) :: V(M,2)

  ! compute variables
  real(8) :: G1(3), G2(3), G3(3)
  complex(8) :: A(2,2)

  ! set cores for turnover based on P(1)
  ! hess
  if (.NOT.P(1)) then

    G1 = Q(1:3)
    G2 = Q(4:6)
    G3 = MISFIT

  ! invhess
  else

    G1 = MISFIT
    G2 = Q(4:6)
    G3 = Q(1:3)

  end if

  ! compute turnover
  call z_rot3_turnover(G1,G2,G3)

  ! move misfit to appropriate side based on P(2)
  ! hess
  if (.NOT.P(2)) then

    ! update Q
    Q(1:3) = G1
    Q(4:6) = G2

    ! update V
    if (VEC) then

      A(1,1) = cmplx(G3(1),G3(2),kind=8)
      A(2,1) = cmplx(G3(3),0d0,kind=8)
      A(1,2) = cmplx(-G3(3),0d0,kind=8)
      A(2,2) = cmplx(G3(1),-G3(2),kind=8)

      V = matmul(V,A)

    end if

    ! pass G3 through upper-triangular part
    call z_upr1utri_rot3swap(.FALSE.,D,C,B,G3)

    ! set MISFIT as G3
    MISFIT = G3

  ! invhess
  else

    ! update Q
    Q(1:3) = G1
    Q(4:6) = G3

    ! pass G2 through upper-triangular part
    call z_upr1utri_rot3swap(.TRUE.,D,C,B,G2)

    ! update V using G2inv
    if (VEC) then

      A(1,1) = cmplx(G2(1),-G2(2),kind=8)
      A(2,1) = cmplx(-G2(3),0d0,kind=8)
      A(1,2) = cmplx(G2(3),0d0,kind=8)
      A(2,2) = cmplx(G2(1),G2(2),kind=8)

      V = matmul(V,A)

    end if

    ! copy G2 to MISFIT
    MISFIT = G2

  end if

  ! update position flag
  P(1) = P(2)

end subroutine z_upr1fact_chasedown

max_stars_repo_name: trcameron/FPML

max_issues_repo_name: trcameron/FPML

max_forks_repo_name: trcameron/FPML


text:
subroutine getfield(itim, comfil, grpnam, funam, ierr, fcom, nmaxus, mmax, kmax, lmax, l)
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: getfield.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/getfield.f90 $
!!--description-----------------------------------------------------------------
!
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use sp_buffer
    !
    implicit none
    !
    integer                                                             , intent(in)    :: nmaxus ! First dimension length
    integer                                                             , intent(in)    :: mmax   ! Second dimension length
    integer                                                             , intent(in)    :: kmax   ! Third dimension length
    integer                                                             , intent(in)    :: lmax   ! Fourth dimension length
    integer                                                             , intent(in)    :: l      ! Fourth dimension index
    integer                                                             , intent(in)    :: itim   ! requested index
    real(fp), dimension(nmaxus, mmax, kmax, lmax)                       , intent(inout) :: fcom   ! data array
    character(*)                                                        , intent(in)    :: comfil ! file name
    character(16)                                                       , intent(in)    :: funam  ! element name
    character(16)                                                       , intent(in)    :: grpnam ! group name
!
! Local variables
!
    integer                                       :: fds
    integer                                       :: i
    integer                                       :: ierr   ! Flag for error when writing to Communication file
    integer                                       :: k
    integer                                       :: m
    integer                                       :: n
    integer      , dimension(3,5)                 :: uindex
    integer                        , external     :: getelt
    integer                        , external     :: clsnef
    integer                        , external     :: open_datdef
!
!! executable statements -------------------------------------------------------
!
    uindex (1,1) = itim ! start index
    uindex (2,1) = itim ! end index
    uindex (3,1) = 1 ! increment in time
    !
    ierr = open_datdef(comfil, fds, .true.)
    if (ierr /= 0) return
    !
    call sbuff_checksize(nmaxus*mmax*kmax)
    ierr = getelt(fds, grpnam, funam, uindex, 1, 4*nmaxus*mmax*kmax, sbuff)
    if (ierr/=0) return
    !
    ierr = clsnef(fds)
    if (ierr/=0) return
    !
    i = 0
    do k = 1, kmax
       do m = 1, mmax
          do n = 1, nmaxus
             i = i+1
             fcom(n, m, k, 1) = sbuff(i)
          enddo
       enddo
    enddo
end subroutine getfield

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Services:
  Selling quality new and used sewing machines and sergers
  Repairing all major brands of machines with a major emphasis on older high quality machines
  Offering sewing classes to both groups and individuals of all skill levels
  Performing expert Tailoring alterations and sewing projects
  Selling and repairing all major brands of vacuum cleaners
  Sharpening knives and scissors

Current sewing classes:
Sewing for Teens/Adults (Ages 12 to Adult)  Tuesdays and Saturdays
Sewing for Kids (Ages 811)  Mondays, Tuesdays, and Wednesdays
Sewing for Kids and Moms  Saturdays

Contact Elaine at: elainebusch27@att.net to get more information on class times or to register.

Past classes:
  Elaine taught a beginning sewing class Tuesdays, in April and May, 6:309:00 p.m. for age 15adults @ the Davis Arts Center.  Students were guided through 4 projects: a simple market bag, lounge pants, pillow slip, and a tshirt. She teaches regular classes in Vacaville through the city of Vacaville Parks & Recreation Dept.

20110616 23:21:47 nbsp I recently hit up sew many things to get my zipper of my leather jacket repaired. It was done in a day for dirt cheap and Im sure to recommend them to my friends! Users/StevenDaubert

20110617 10:06:24 nbsp Aaaaahhhhhh they are able to repair zippers on jackets/hoodies?? I had this hoodie that I really liked and the zipper broke off. Said hoodie is now somewhere I dont know :(
But I will have to check out this shop very soon Users/NikhilDahal
 I would imagine it would be simple, if it can be done on leather having a master do it on fabric shouldnt be a problem! Give her a jingle Daubert
  Thanks Daubert Users/NikhilDahal

20110625 20:16:05 nbsp I took a sewing class from Elaine and John tuned up and repaired 3 of my machines and I am very pleased with the quality of their workThanks Elaine and John! Users/stephaniestoops

20110626 11:05:54 nbsp I got my leather jacket repaired again! It had a few holes and a nice tear and Elaine got a hold of a leather repair kit and it turned out wonderful and it was fast! Users/StevenDaubert

20110723 21:21:02 nbsp I brought Elaine two pairs of pants to be hemmed and they were done that afternoon. Im very happy with the results and the price cant be beat! Users/EmmyMelton

20110823 08:34:34 nbsp My wedding dress came in two weeks short of my wedding date and about two inches too small. Elaine was amazing. She was very professional, thorough, and her services were well priced. She was able to repair my dress on short notice and did great work.
I highly recommend her. Users/nkristis


text: 20110412 17:57:22 nbsp Thanks for sharing your experiences on Aggie Square. Too bad your bosses wont stand for it. Happy editing! Users/WilliamLewis


! (C) Copyright 2018-2019 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

module pressure_vt_mod

use fv3jedi_constants_mod, only: kappa, grav, rad2deg, constoz
use fv3jedi_geom_mod, only: fv3jedi_geom
use fv3jedi_kinds_mod, only: kind_real

implicit none
private

public delp_to_pe_p_logp
public pe_to_pkz
public ps_to_pkz
public pe_to_delp
public delp_to_pe
public pe_to_pk
public ps_to_delp
public ps_to_pe
public ps_to_delp_tl
public ps_to_delp_ad
public tropprs
public tropprs_th

contains

!----------------------------------------------------------------------------
! Pressure thickness to pressure (edge), pressure (mid) and log p (mid) -----
!----------------------------------------------------------------------------

subroutine delp_to_pe_p_logp(geom,delp,pe,p,logp)

 type(fv3jedi_geom)  , intent(in ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in ) :: delp(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness
 real(kind=kind_real), intent(out) ::   pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure edge/interface
 real(kind=kind_real), intent(out) ::    p(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure mid point
 real(kind=kind_real), optional, intent(out) :: logp(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Log of pressure mid point

 !Locals
 integer :: isc,iec,jsc,jec,k

 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec

 !Pressure at layer edge
 pe(isc:iec,jsc:jec,1) = geom%ptop
 do k = 2,geom%npz+1
   pe(isc:iec,jsc:jec,k) = pe(isc:iec,jsc:jec,k-1) + delp(isc:iec,jsc:jec,k-1)
 enddo

 !Midpoint pressure
 p(isc:iec,jsc:jec,:) = 0.5*(pe(isc:iec,jsc:jec,2:geom%npz+1) + pe(isc:iec,jsc:jec,1:geom%npz))

 if (present(logp)) then
   !Log pressure
   logp(isc:iec,jsc:jec,:) = log(p(isc:iec,jsc:jec,:))
 endif

end subroutine delp_to_pe_p_logp

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

subroutine pe_to_pkz(geom,pe,pkz)

type(fv3jedi_geom)  , intent(in ) :: geom !Geometry for the model
real(kind=kind_real), intent(in ) ::  pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure edge/interface
real(kind=kind_real), intent(out) :: pkz(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure to the kappa

integer :: k
real(kind=kind_real) :: peln(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1)
real(kind=kind_real) ::   pk(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1)

peln = log(pe)
pk = exp(kappa*peln)

do k=1,geom%npz
  pkz(:,:,k) = (pk(:,:,k+1)-pk(:,:,k)) / (kappa*(peln(:,:,k+1)-peln(:,:,k)))
enddo

end subroutine pe_to_pkz

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

subroutine pe_to_delp(geom,pe,delp)

 type(fv3jedi_geom)  , intent(in ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in ) ::   pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure edge/interface
 real(kind=kind_real), intent(out) :: delp(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness

 !Locals
 integer :: isc,iec,jsc,jec,k

 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec

 !Pressure at layer edge
 do k = 1,geom%npz
   delp(isc:iec,jsc:jec,k) = pe(isc:iec,jsc:jec,k+1) - pe(isc:iec,jsc:jec,k)
 enddo

end subroutine pe_to_delp

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

subroutine delp_to_pe( geom, delp, pe )

type(fv3jedi_geom)  , intent(in ) :: geom !Geometry for the model
real(kind=kind_real), intent(in ) :: delp(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness
real(kind=kind_real), intent(out) ::   pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure edge/interface

!Locals
integer :: isc,iec,jsc,jec,k

isc = geom%isc
iec = geom%iec
jsc = geom%jsc
jec = geom%jec

!Pressure at layer edge
pe(isc:iec,jsc:jec,1) = geom%ptop
do k = 2,geom%npz+1
  pe(isc:iec,jsc:jec,k) = pe(isc:iec,jsc:jec,k-1) + delp(isc:iec,jsc:jec,k-1)
enddo

end subroutine delp_to_pe

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

subroutine pe_to_pk( geom, pe, pk )

type(fv3jedi_geom)  , intent(in ) :: geom !Geometry for the model
real(kind=kind_real), intent(in ) :: pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure edge/interface
real(kind=kind_real), intent(out) :: pk(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure to the kappa

!Locals
integer :: i, j, k
real(kind=kind_real) :: pel1, pel2
real(kind=kind_real) :: pek1, pek2

do k=1,geom%npz
  do j = geom%jsc,geom%jec
    do i = geom%isc,geom%iec

      pel1 = log(pe(i,j,k+1))
      pel2 = log(pe(i,j,k))

      pek1 = exp(kappa*pel1)
      pek2 = exp(kappa*pel2)

      pk(i,j,k) = (pek1-pek2)/(kappa*(pel1-pel2))

    end do
  end do
end do

end subroutine pe_to_pk

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

subroutine ps_to_pe(geom,ps,pe)

 type(fv3jedi_geom)  , intent(in   ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in   ) :: ps(geom%isc:geom%iec,geom%jsc:geom%jec,1           )   !Surface pressure
 real(kind=kind_real), intent(inout) :: pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1)   !Pressure thickness

 integer :: k

 do k = 1,geom%npz+1
   pe(:,:,k) = geom%ak(k) + geom%bk(k) * ps(:,:,1)
 enddo

endsubroutine ps_to_pe

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

subroutine ps_to_pkz(geom,ps,pkz)

 type(fv3jedi_geom)  , intent(in   ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in   ) ::  ps(geom%isc:geom%iec,geom%jsc:geom%jec,1         )   !Surface pressure
 real(kind=kind_real), intent(inout) :: pkz(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure ^ kappa

 integer :: k
 real(kind=kind_real) :: pe(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz+1) !Pressure

 call ps_to_pe(geom, ps, pe)
 call pe_to_pkz(geom, pe, pkz)

endsubroutine ps_to_pkz

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

subroutine ps_to_delp(geom,ps,delp)

 type(fv3jedi_geom)  , intent(in   ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in   ) ::   ps(geom%isc:geom%iec,geom%jsc:geom%jec           )   !Surface pressure
 real(kind=kind_real), intent(inout) :: delp(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness

 integer :: isc,iec,jsc,jec,i,j,k

 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec

 do k = 1,geom%npz
   do j = jsc,jec
     do i = isc,iec
       delp(i,j,k) = geom%ak(k+1) + geom%bk(k+1)*ps(i,j) - (geom%ak(k) + geom%bk(k)*ps(i,j))
     enddo
   enddo
 enddo

endsubroutine ps_to_delp

subroutine ps_to_delp_tl(geom,ps_tl,delp_tl)

 type(fv3jedi_geom)  , intent(in   ) :: geom !Geometry for the model
 real(kind=kind_real), intent(in   ) ::   ps_tl(geom%isc:geom%iec,geom%jsc:geom%jec           )   !Surface pressure
 real(kind=kind_real), intent(inout) :: delp_tl(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness

 integer :: isc,iec,jsc,jec,i,j,k

 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec

 delp_tl = 0.0_kind_real
 do k = 1,geom%npz
   do j = jsc,jec
     do i = isc,iec
       delp_tl(i,j,k) = geom%bk(k+1)*ps_tl(i,j) - geom%bk(k)*ps_tl(i,j)
     enddo
   enddo
 enddo

endsubroutine ps_to_delp_tl

subroutine ps_to_delp_ad(geom,ps_ad,delp_ad)

 type(fv3jedi_geom)  , intent(in   ) :: geom !Geometry for the model
 real(kind=kind_real), intent(inout) ::   ps_ad(geom%isc:geom%iec,geom%jsc:geom%jec           )   !Surface pressure
 real(kind=kind_real), intent(inout) :: delp_ad(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure thickness

 integer :: isc,iec,jsc,jec,i,j,k

 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec

 ps_ad = 0.0_kind_real

 do k = geom%npz,1,-1
   do j = jec,jsc,-1
     do i = iec,isc,-1
       ps_ad(i,j) = ps_ad(i,j) + (geom%bk(k+1) - geom%bk(k))*delp_ad(i,j,k)
     enddo
   enddo
 enddo

endsubroutine ps_to_delp_ad

! --------------------------------------------------------------------------------------------------
! Locate tropopause using potential vorticity and ozone

subroutine tropprs(geom, ps, prs, tv, o3, vort, tprs)

 type(fv3jedi_geom)  , intent(in   ) :: geom
 real(kind=kind_real), intent(in   ) ::   ps(geom%isc:geom%iec,geom%jsc:geom%jec,1         )   !Surface pressure [Pa]
 real(kind=kind_real), intent(in   ) ::  prs(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure, midpoint [Pa]
 real(kind=kind_real), intent(in   ) ::   tv(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Virtual temperature [K]
 real(kind=kind_real), intent(in   ) ::   o3(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Ozone mixing ratio [ppmv]
 real(kind=kind_real), intent(in   ) :: vort(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Vorticity [1/s]
 real(kind=kind_real), intent(inout) :: tprs(geom%isc:geom%iec,geom%jsc:geom%jec,1         )   !Tropopause pressure [Pa]

 ! Local
 integer :: isc,iec,jsc,jec
 integer :: i,j,k,kk
 integer :: npz
 integer :: itrop_k
 integer :: ifound_pv, ifound_oz
 integer :: itrp_pv, itrp_oz

 real(kind=kind_real), parameter :: r1e5 = 1.0e5
 real(kind=kind_real), parameter :: r2e6 = 2.0e-6
 real(kind=kind_real), parameter :: r3e7 = 3.0e-7
 real(kind=kind_real) :: pm1, pp1, thetam1, thetap1, pv
 real(kind=kind_real) :: psi
 real(kind=kind_real) :: prsl(geom%npz), pvort(geom%npz)
 real(kind=kind_real) :: o3mr(geom%npz)
 real(kind=kind_real) :: lat, wgt

 ! Domain
 isc = geom%isc
 iec = geom%iec
 jsc = geom%jsc
 jec = geom%jec
 npz = geom%npz

 ! Loop through locations
 do j = jsc, jec
   do i = isc, iec

     ! Convert latitude from radian to degree
     lat = geom%grid_lat(i,j) * rad2deg

     do k = 1, npz
       prsl(k) = prs(i,j,k)  ! [Pa]
     enddo
     psi = 1.0 / ps(i,j,1)

     ! Convert O3 unit from ppmv to mixing ratio [kg/kg]
     do k = 1, npz
       o3mr(k) = o3(i,j,k) / constoz
     enddo

     ! Compute potential vortivity (pv) at midpoint pressure
     ! Work from model surface to top
     do k = 2, npz-1   ! index from top to surface:     2,     3, ..., npz-1
       kk = npz-k+1    ! index from surface to top: npz-1, npz-2, ..., 2
       pm1 = prsl(kk-1)
       pp1 = prsl(kk+1)
       thetam1 = tv(i,j,kk-1) * (r1e5 / pm1)**kappa
       thetap1 = tv(i,j,kk-1) * (r1e5 / pp1)**kappa
       pv = grav * vort(i,j,kk) * (thetam1 - thetap1) / (pm1 - pp1)
       pvort(kk) = abs(pv)
     enddo
     pvort(1) = pvort(2)
     pvort(npz) = pvort(npz-1)

     ! Locate tropopause using vorticity and ozone
     ! Search upward (decressing pressure) for tropopause above sigma 0.7
     ifound_pv = 0; itrp_pv = npz
     ifound_oz = 0; itrp_oz = npz
     do k = 2, npz-1   ! index from top to surface:     2,     3, ..., npz-1
       kk = npz-k+1    ! index from surface to top: npz-1, npz-2, ..., 2
       if (prsl(kk) * psi < 0.7) then  ! sigma level < 0.7
         ! Tropopause at level where pv > 2e-6
         if (pvort(kk) > r2e6 .and. ifound_pv == 0) then
           ifound_pv = 1
           itrp_pv = kk
         endif
         !Tropopause at level where o3mr > 3e-7
         if (o3mr(kk) > r3e7 .and. ifound_oz == 0) then
           ifound_oz = 1
           itrp_oz = kk
         endif
       endif
     enddo

     ! Merge pv and o3 tropopause levels between 20 and 40 degree latitudes
     lat =abs(lat)
     if (lat >= 40.0) then
       itrop_k = itrp_pv
     elseif (lat >=20) then
       wgt = (lat - 20.0) / 20.0
       itrop_k = wgt * itrp_pv + (1.0 - wgt) * itrp_oz
     else ! lat
       itrop_k = itrp_oz
     endif
     itrop_k = max(1, min(itrop_k, npz))
     tprs(i,j,1) = prsl(itrop_k)

     ! The tropopause pressures [Pa] are used for the following two things in UFO
     ! (1) deflate the moisture sensitivity vectors for satellite radiance data
     ! (2) cloud detection in IR quality control
     ! Bounds are applied on the tropopause pressure to make sure we are deflating
     ! at the very minimum above the 150 hPa and no wherer below 350 hPa

     tprs(i,j,1) = max(15000.0, min(35000.0, tprs(i,j,1)))

   enddo
 enddo

end subroutine tropprs

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

!..Find tropopause height based on delta-Theta / delta-Z. The 10/1500 ratio
!.. approximates a vertical line on typical SkewT chart near typical (mid-latitude)
!.. tropopause height.  Since messy data could give us a false signal of such a
!.. transition, do the check over at least 35 mb delta-pres, not just a level-to-level
!.. check.  This method has potential failure in arctic-like conditions with extremely
!.. low tropopause height, as would any other diagnostic, so ensure resulting k_tropo
!.. level is above 700hPa.

subroutine tropprs_th(geom, p, z, t, tprs)

type(fv3jedi_geom)  , intent(in   ) :: geom
real(kind=kind_real), intent(in   ) ::    p(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Pressure, midpoint [Pa]
real(kind=kind_real), intent(in   ) ::    t(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Temperature [K]
real(kind=kind_real), intent(in   ) ::    z(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)   !Height/altitude [m]
real(kind=kind_real), intent(inout) :: tprs(geom%isc:geom%iec,geom%jsc:geom%jec,1         )   !Tropopause pressure [Pa]

! Local
integer :: isc,iec,jsc,jec, i, j, k
integer :: k1, k2, k_p50, k_p150, k_tropo, npz
real(kind=kind_real) :: p2(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)
real(kind=kind_real) :: z2(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)
real(kind=kind_real) :: t2(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)
real(kind=kind_real) :: theta(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)

! Domain
isc = geom%isc
iec = geom%iec
jsc = geom%jsc
jec = geom%jec
npz = geom%npz

! Flip to index increasing upwards
do k = 1, npz
  do j = jsc, jec
    do i = isc, iec
      t2(i,j,k) = t(i,j,npz-k+1)
      p2(i,j,k) = p(i,j,npz-k+1)
      z2(i,j,k) = z(i,j,npz-k+1)
    enddo
  enddo
enddo

! Compute tropopause pressure
do j = jsc, jec
  do i = isc, iec
    k_p150 = 0
    k_p50 = 0
    do k = npz, 1, -1
      theta(i,j,k) = t2(i,j,k)*((100000.0/p2(i,j,k))**(287.04/1004.))
      if (p2(i,j,k) .gt. 4999.0  .and. k_p50  .eq. 0) k_p50  = k
      if (p2(i,j,k) .gt. 14999.0 .and. k_p150 .eq. 0) k_p150 = k
    enddo
    if ( (k_p50-k_p150) .lt. 3) k_p150 = k_p50-3
    do k = k_p150-2, 1, -1
      k1 = k
      k2 = k1+2
      do while (k1.gt.1 .and. (p2(i,j,k1)-p2(i,j,k2)).lt.3500.0 .and. p2(i,j,k1).lt.70000.)
        k1 = k1-1
      enddo
      if ( ((theta(i,j,k2)-theta(i,j,k1))/(z2(i,j,k2)-z2(i,j,k1))) .lt. 10./1500.) exit
    enddo
    k_tropo = max(3, min(nint(0.5*(k1+k2+1)), k_p50-1))
    tprs(i,j,1) = p2(i,j,k_tropo)
  enddo
enddo

end subroutine tropprs_th

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

end module pressure_vt_mod


text:
*> \brief \b DLARGE
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLARGE( N, A, LDA, ISEED, WORK, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, N
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( 4 )
*       DOUBLE PRECISION   A( LDA, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLARGE pre- and post-multiplies a real general n by n matrix A
*> with a random orthogonal matrix: A = U*D*U'.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*>          On entry, the original n by n matrix A.
*>          On exit, A is overwritten by U*A*U' for some random
*>          orthogonal matrix U.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= N.
*> \endverbatim
*>
*> \param[in,out] ISEED
*> \verbatim
*>          ISEED is INTEGER array, dimension (4)
*>          On entry, the seed of the random number generator; the array
*>          elements must be between 0 and 4095, and ISEED(4) must be
*>          odd.
*>          On exit, the seed is updated.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup double_matgen
*
*  =====================================================================
      SUBROUTINE DLARGE( N, A, LDA, ISEED, WORK, INFO )
*
*  -- LAPACK auxiliary routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      DOUBLE PRECISION   A( LDA, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      DOUBLE PRECISION   TAU, WA, WB, WN
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEMV, DGER, DLARNV, DSCAL, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, SIGN
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DNRM2
      EXTERNAL           DNRM2
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -3
      END IF
      IF( INFO.LT.0 ) THEN
         CALL XERBLA( 'DLARGE', -INFO )
         RETURN
      END IF
*
*     pre- and post-multiply A by random orthogonal matrix
*
      DO 10 I = N, 1, -1
*
*        generate random reflection
*
         CALL DLARNV( 3, ISEED, N-I+1, WORK )
         WN = DNRM2( N-I+1, WORK, 1 )
         WA = SIGN( WN, WORK( 1 ) )
         IF( WN.EQ.ZERO ) THEN
            TAU = ZERO
         ELSE
            WB = WORK( 1 ) + WA
            CALL DSCAL( N-I, ONE / WB, WORK( 2 ), 1 )
            WORK( 1 ) = ONE
            TAU = WB / WA
         END IF
*
*        multiply A(i:n,1:n) by random reflection from the left
*
         CALL DGEMV( 'Transpose', N-I+1, N, ONE, A( I, 1 ), LDA, WORK,
     $               1, ZERO, WORK( N+1 ), 1 )
         CALL DGER( N-I+1, N, -TAU, WORK, 1, WORK( N+1 ), 1, A( I, 1 ),
     $              LDA )
*
*        multiply A(1:n,i:n) by random reflection from the right
*
         CALL DGEMV( 'No transpose', N, N-I+1, ONE, A( 1, I ), LDA,
     $               WORK, 1, ZERO, WORK( N+1 ), 1 )
         CALL DGER( N, N-I+1, -TAU, WORK( N+1 ), 1, WORK, 1, A( 1, I ),
     $              LDA )
   10 CONTINUE
      RETURN
*
*     End of DLARGE
*
      END

max_stars_repo_name: sydha/Lapack

max_issues_repo_name: sydha/Lapack

max_forks_repo_name: sydha/Lapack


text:
      COMMON/Options/o_el,o_bgd,o_w1,prog
      character*1 prog
      integer o_el,o_bgd,o_w1

max_stars_repo_name: mantidproject/quasielasticbayes

max_issues_repo_name: mantidproject/quasielasticbayes

max_forks_repo_name: mantidproject/quasielasticbayes


text:
!==============================================================================!
  subroutine Matrix_Mod_Topology(grid, matrix)
!------------------------------------------------------------------------------!
!   Determines the topology of the system matrix.                              !
!   It relies only on SideC structure. Try to keep it that way.                !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type), target :: grid
  type(Matrix_Type)       :: matrix
!-----------------------------------[Locals]-----------------------------------!
  integer              :: c, s, pos, pos1, pos2, n
  integer              :: c1, c2  ! cell 1 and 2
  integer              :: n1, n2  ! neighbour 1 and 2
  integer, allocatable :: stencw(:)
!==============================================================================!

  ! Store pointer to the grid
  matrix % pnt_grid => grid

  ! Memory allocation
  allocate(stencw(grid % n_cells)); stencw = 1

  if(this_proc < 2) print *, '# Determining matrix topology.'

  ! Initialize number of nonzeros
  matrix % nonzeros = 0

  ! Compute stencis widths
  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2 > 0) then
      stencw(c1) = stencw(c1) + 1
      stencw(c2) = stencw(c2) + 1
    end if
  end do

  ! Count the nonzero entries and allocate the memory for the array
  n = 0
  do c = 1, grid % n_cells
    n = n + stencw(c)
  end do
  matrix % nonzeros = n + 1
  allocate(matrix % val(n+1)); matrix % val = 0. ! it reffers to matrix % row+1
  allocate(matrix % col(n+1)); matrix % col = 0  ! it reffers to matrix % row+1
  allocate(matrix % mir(n+1)); matrix % mir = 0  ! it reffers to matrix % row+1

  ! Form matrix % row and diagonal only formation of matrix % col
  matrix % row(1) = 1
  do c = 1, grid % n_cells
    matrix % row(c + 1) = matrix % row(c) + stencw(c)
    matrix % col(matrix % row(c)) = c   ! it is first to its own self
    stencw(c) = 1                       ! reset stencw
  end do

  ! Extend matrix % col entries with off-diagonal terms
  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2 > 0) then
      matrix % col(matrix % row(c1) + stencw(c1)) = c2
      matrix % col(matrix % row(c2) + stencw(c2)) = c1
      stencw(c1) = stencw(c1) + 1
      stencw(c2) = stencw(c2) + 1
    end if
  end do

  ! Sort matrix % col to make them nice and neat
  ! and also locate the position of diagonal
  do c = 1, grid % n_cells
    call Sort_Mod_Int(matrix % col(matrix % row(c) :  &
                                   matrix % row(c) + stencw(c) - 1))
    do pos = matrix % row(c),matrix % row(c+1)-1
      if(matrix % col(pos) .eq. c) matrix % dia(c) = pos
    end do
  end do

  ! Find mirror positions
  do c1 = 1, grid % n_cells
    do pos1 = matrix % row(c1), matrix % row(c1 + 1) - 1
      n1 = matrix % col(pos1)  ! at this point you have c1 and n1

      ! Inner loop (it might probably go from 1 to c1-1
      c2 = n1
      do pos2 = matrix % row(c2), matrix % row(c2 + 1) - 1
        n2 = matrix % col(pos2)  ! at this point you have c2 and n2

        if(n2 .eq. c1) then
          matrix % mir(pos1) = pos2
          matrix % mir(pos2) = pos1
          goto 2  ! done with the inner loop, get out
        end if
      end do

2     continue
    end do
  end do

  ! Connect faces with matrix entries
  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2 > 0) then

      ! Where is matrix(c1,c2) and ...
      do c = matrix % row(c1), matrix % row(c1+1)-1
        if(matrix % col(c)  .eq.  c2) matrix % pos(1, s) = c
      end do

      ! Where is matrix(c2,c1) and ...
      do c=matrix % row(c2),matrix % row(c2+1)-1
        if(matrix % col(c)  .eq.  c1) matrix % pos(2, s) = c
      end do
    end if
  end do

  if(this_proc < 2) print *, '# Finished !'

  deallocate(stencw)

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
The juvenile detention facility is the temporary detention and treatment facility in and for Yolo County for minors who have been Juvenile Justice System charged with a violation of the law or who have violated conditions of probation. It is a division of the Probation Department but operated under the authority and regulations of the California Welfare and Institutions Code, California Penal Code and the California Administrative Code. Mandates include that the facility shall not be treated as a jail, prison, or penal institution and shall be operated in all respects as nearly like a home as possible.

In general, juvenile hall staff is responsible for booking decisions, admissions procedures, and the daily care and custody of minors housed in the facility.


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

20090811 17:55:36 nbsp I moved the page you created to S&S Pool Service. Users/WilliamLewis

20090811 18:00:57 nbsp Hey, cut it out. Were friendly folk here. You dont need to go delete all my stuff. Im just trying to help you. Users/WilliamLewis

20090811 19:02:27 nbsp Welcome to the Wiki. You seem to have gotten off to a rough start here on the Wiki. Take some time to browse around and maybe read Welcome to the Wiki/Business Owner and Wiki Community. Calling people names and deleting other peoples contributions is only going to give people around here a negative impression of your customer service skills. Users/JasonAller

  Well put.  Most Davis Wikians are all about forgiveness. Users/SunjeetBaadkar


text:
! Main program

program intersections

    use newton, only: solve
    use functions, only: f_int, fprime_int

    implicit none
    real(kind=8) :: x, x0, fx
    real(kind=8) :: x0vals(4)
    integer :: iters, itest
	logical :: debug

    print *, "Test routine for computing zero of f"
    debug = .true.

    ! values to test as x0:
    x0vals = (/-2.2d0, -1.6d0, -0.8d0, 1.45d0/)

    do itest=1,4
        x0 = x0vals(itest)
		print *, ' '  ! blank line
        call solve(f_int, fprime_int, x0, x, iters, debug)

        print 11, x, iters
11      format('solver returns x = ', es22.15, ' after', i3, ' iterations')

        fx = f_int(x)
        print 12, fx
12      format('the value of f(x) is ', es22.15)

        enddo

end program intersections

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYSPLI
C...Splits a hadron remnant into two (partons or hadron + parton)
C...in case it is more complicated than just a quark or a diquark.

      SUBROUTINE PYSPLI(KF,KFLIN,KFLCH,KFLSP)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)

C...Commonblocks. PYDAT1 temporary
      include 'inc/pypars'
      include 'inc/pyint1'
      include 'inc/pydat1'

C...Local array.
      DIMENSION KFL(3)

C...Preliminaries. Parton composition.
      KFA=ABS(KF)
      KFS=SIGN(1,KF)
      KFL(1)=MOD(KFA/1000,10)
      KFL(2)=MOD(KFA/100,10)
      KFL(3)=MOD(KFA/10,10)
      IF(KFA.EQ.22.AND.MINT(109).EQ.2) THEN
        KFL(2)=INT(1.5D0+PYR(0))
        IF(MINT(105).EQ.333) KFL(2)=3
        IF(MINT(105).EQ.443) KFL(2)=4
        KFL(3)=KFL(2)
      ELSEIF((KFA.EQ.111.OR.KFA.EQ.113).AND.PYR(0).GT.0.5D0) THEN
        KFL(2)=2
        KFL(3)=2
      ELSEIF(KFA.EQ.223.AND.PYR(0).GT.0.5D0) THEN
        KFL(2)=1
        KFL(3)=1
      ELSEIF((KFA.EQ.130.OR.KFA.EQ.310).AND.PYR(0).GT.0.5D0) THEN
        KFL(2)=MOD(KFA/10,10)
        KFL(3)=MOD(KFA/100,10)
      ENDIF
      IF(KFLIN.NE.21.AND.KFLIN.NE.22.AND.KFLIN.NE.23) THEN
        KFLR=KFLIN*KFS
      ELSE
        KFLR=KFLIN
      ENDIF
      KFLCH=0

C...Subdivide lepton.
      IF(KFA.GE.11.AND.KFA.LE.18) THEN
        IF(KFLR.EQ.KFA) THEN
          KFLSP=KFS*22
        ELSEIF(KFLR.EQ.22) THEN
          KFLSP=KFA
        ELSEIF(KFLR.EQ.-24.AND.MOD(KFA,2).EQ.1) THEN
          KFLSP=KFA+1
        ELSEIF(KFLR.EQ.24.AND.MOD(KFA,2).EQ.0) THEN
          KFLSP=KFA-1
        ELSEIF(KFLR.EQ.21) THEN
          KFLSP=KFA
          KFLCH=KFS*21
        ELSE
          KFLSP=KFA
          KFLCH=-KFLR
        ENDIF

C...Subdivide photon.
      ELSEIF(KFA.EQ.22.AND.MINT(109).NE.2) THEN
        IF(KFLR.NE.21) THEN
          KFLSP=-KFLR
        ELSE
          RAGR=0.75D0*PYR(0)
          KFLSP=1
          IF(RAGR.GT.0.125D0) KFLSP=2
          IF(RAGR.GT.0.625D0) KFLSP=3
          IF(PYR(0).GT.0.5D0) KFLSP=-KFLSP
          KFLCH=-KFLSP
        ENDIF

C...Subdivide Reggeon or Pomeron.
      ELSEIF(KFA.EQ.110.OR.KFA.EQ.990) THEN
        IF(KFLIN.EQ.21) THEN
          KFLSP=KFS*21
        ELSE
          KFLSP=-KFLIN
        ENDIF

C...Subdivide meson.
      ELSEIF(KFL(1).EQ.0) THEN
        KFL(2)=KFL(2)*(-1)**KFL(2)
        KFL(3)=-KFL(3)*(-1)**ABS(KFL(2))
        IF(KFLR.EQ.KFL(2)) THEN
          KFLSP=KFL(3)
        ELSEIF(KFLR.EQ.KFL(3)) THEN
          KFLSP=KFL(2)
        ELSEIF(KFLR.EQ.21.AND.PYR(0).GT.0.5D0) THEN
          KFLSP=KFL(2)
          KFLCH=KFL(3)
        ELSEIF(KFLR.EQ.21) THEN
          KFLSP=KFL(3)
          KFLCH=KFL(2)
        ELSEIF(KFLR*KFL(2).GT.0) THEN
          NTRY=0
  100     NTRY=NTRY+1
          CALL PYKFDI(-KFLR,KFL(2),KFDUMP,KFLCH)
          IF(KFLCH.EQ.0.AND.NTRY.LT.100) THEN
            GOTO 100
          ELSEIF(KFLCH.EQ.0) THEN
            CALL PYERRM(14,'(PYSPLI:) caught in infinite loop')
            MINT(51)=1
            RETURN
          ENDIF
          KFLSP=KFL(3)
        ELSE
          NTRY=0
  110     NTRY=NTRY+1
          CALL PYKFDI(-KFLR,KFL(3),KFDUMP,KFLCH)
          IF(KFLCH.EQ.0.AND.NTRY.LT.100) THEN
            GOTO 110
          ELSEIF(KFLCH.EQ.0) THEN
            CALL PYERRM(14,'(PYSPLI:) caught in infinite loop')
            MINT(51)=1
            RETURN
          ENDIF
          KFLSP=KFL(2)
        ENDIF

C...Special case for extracting photon from baryon without splitting
C...the latter. (Currently only used by external programs.)
      ELSEIF(KFLIN.EQ.22.AND.MSTP(98).EQ.1) then
        KFLSP=KFA
        KFLCH=0

C...Subdivide baryon.
      ELSE
        NAGR=0
        DO 120 J=1,3
          IF(KFLR.EQ.KFL(J)) NAGR=NAGR+1
  120   CONTINUE
        IF(NAGR.GE.1) THEN
          RAGR=0.00001D0+(NAGR-0.00002D0)*PYR(0)
          IAGR=0
          DO 130 J=1,3
            IF(KFLR.EQ.KFL(J)) RAGR=RAGR-1D0
            IF(IAGR.EQ.0.AND.RAGR.LE.0D0) IAGR=J
  130     CONTINUE
        ELSE
          IAGR=1.00001D0+2.99998D0*PYR(0)
        ENDIF
        ID1=1
        IF(IAGR.EQ.1) ID1=2
        IF(IAGR.EQ.1.AND.KFL(3).GT.KFL(2)) ID1=3
        ID2=6-IAGR-ID1
        KSP=3
        IF(MOD(KFA,10).EQ.2.AND.KFL(1).EQ.KFL(2)) THEN
          IF(IAGR.NE.3.AND.PYR(0).GT.0.25D0) KSP=1
        ELSEIF(MOD(KFA,10).EQ.2.AND.KFL(2).GE.KFL(3)) THEN
          IF(IAGR.NE.1.AND.PYR(0).GT.0.25D0) KSP=1
        ELSEIF(MOD(KFA,10).EQ.2) THEN
          IF(IAGR.EQ.1) KSP=1
          IF(IAGR.NE.1.AND.PYR(0).GT.0.75D0) KSP=1
        ENDIF
        KFLSP=1000*KFL(ID1)+100*KFL(ID2)+KSP
        IF(KFLR.EQ.21) THEN
          KFLCH=KFL(IAGR)
        ELSEIF(NAGR.EQ.0.AND.KFLR.GT.0) THEN
          NTRY=0
  140     NTRY=NTRY+1
          CALL PYKFDI(-KFLR,KFL(IAGR),KFDUMP,KFLCH)
          IF(KFLCH.EQ.0.AND.NTRY.LT.100) THEN
            GOTO 140
          ELSEIF(KFLCH.EQ.0) THEN
            CALL PYERRM(14,'(PYSPLI:) caught in infinite loop')
            MINT(51)=1
            RETURN
          ENDIF
        ELSEIF(NAGR.EQ.0) THEN
          NTRY=0
  150     NTRY=NTRY+1
          CALL PYKFDI(10000*KFL(ID1)+KFLSP,-KFLR,KFDUMP,KFLCH)
          IF(KFLCH.EQ.0.AND.NTRY.LT.100) THEN
            GOTO 150
          ELSEIF(KFLCH.EQ.0) THEN
            CALL PYERRM(14,'(PYSPLI:) caught in infinite loop')
            MINT(51)=1
            RETURN
          ENDIF
          KFLSP=KFL(IAGR)
        ENDIF
      ENDIF

C...Add on correct sign for result.
      KFLCH=KFLCH*KFS
      KFLSP=KFLSP*KFS

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


SUBROUTINE p43(input_file,output_file)
!-------------------------------------------------------------------------
! Program 4.3 Analysis of elastic beams using 2-node beam elements
!             (elastic foundation optional).
!-------------------------------------------------------------------------
 USE main
 USE geom
 IMPLICIT NONE
 CHARACTER(len=60),INTENT(IN) :: input_file
 CHARACTER(len=60),INTENT(OUT) :: output_file
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndof=4,nels,neq,nod=2,      &
   nodof=2,nn,nprops,np_types,nr,nlen
 REAL(iwp)::penalty=1.0e20_iwp,zero=0.0_iwp
!-----------------------dynamic arrays------------------------------------
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),kdiag(:),nf(:,:),no(:),      &
   node(:),num(:),sense(:)
 REAL(iwp),ALLOCATABLE::action(:),eld(:),ell(:),km(:,:),kv(:),loads(:),   &
   mm(:,:),prop(:,:),value(:)
 CHARACTER(LEN=15)::argv
!-----------------------input and initialisation--------------------------
 OPEN(10,FILE=input_file)
 OPEN(11,FILE=output_file)
 READ(10,*)nels,nprops,np_types
 nn=nels+1
 ALLOCATE(g(ndof),num(nod),nf(nodof,nn),etype(nels),ell(nels),eld(ndof),  &
   km(ndof,ndof),mm(ndof,ndof),action(ndof),g_g(ndof,nels),               &
   prop(nprops,np_types))
 READ(10,*)prop
 etype=1
 IF(np_types>1)READ(10,*)etype
 READ(10,*)ell
 nf=1
 READ(10,*)nr,(k,nf(:,k),i=1,nr)
 CALL formnf(nf)
 neq=MAXVAL(nf)
 ALLOCATE(kdiag(neq),loads(0:neq))
!-----------------------loop the elements to find global array sizes------
 kdiag=0
 elements_1: DO iel=1,nels
   num=(/iel,iel+1/)
   CALL num_to_g(num,nf,g)
   g_g(:,iel)=g
   CALL fkdiag(kdiag,g)
 END DO elements_1
 DO i=2,neq
   kdiag(i)=kdiag(i)+kdiag(i-1)
 END DO
 ALLOCATE(kv(kdiag(neq)))
 WRITE(11,'(2(A,I5))')                                                    &
   " There are",neq," equations and the skyline storage is",kdiag(neq)
!-----------------------global stiffness matrix assembly------------------
 kv=zero
 elements_2: DO iel=1, nels
   CALL beam_km(km,prop(1,etype(iel)),ell(iel))
   mm=zero
   IF(nprops>1)CALL beam_mm(mm,prop(2,etype(iel)),ell(iel))
   g=g_g(:,iel)
   CALL fsparv(kv,km+mm,g,kdiag)
 END DO elements_2
!-----------------------read loads and/or displacements-------------------
 loads=zero
 READ(10,*)loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes)
 READ(10,*)fixed_freedoms
 IF(fixed_freedoms/=0)THEN
   ALLOCATE(node(fixed_freedoms),no(fixed_freedoms),                      &
     sense(fixed_freedoms),value(fixed_freedoms))
   READ(10,*)(node(i),sense(i),value(i),i=1,fixed_freedoms)
   DO i=1,fixed_freedoms
     no(i)=nf(sense(i),node(i))
   END DO
   kv(kdiag(no))=kv(kdiag(no))+penalty
   loads(no)=kv(kdiag(no))*value
 END IF
!-----------------------equation solution --------------------------------
 CALL sparin(kv,kdiag)
 CALL spabac(kv,loads,kdiag)
 loads(0)=zero
 WRITE(11,'(/A)')"  Node Translation Rotation"
 DO k=1,nn
   WRITE(11,'(I5,2E12.4)')k,loads(nf(:,k))
 END DO
!-----------------------retrieve element end actions----------------------
 WRITE(11,'(/A)')" Element Force       Moment      Force       Moment"
 elements_3: DO iel=1,nels
   CALL beam_km(km,prop(1,etype(iel)),ell(iel))
   mm=zero
   IF(nprops>1)CALL beam_mm(mm,prop(2,etype(iel)),ell(iel))
   g=g_g(:,iel)
   eld=loads(g)
   action=MATMUL(km+mm,eld)
   WRITE(11,'(I5,4E12.4)')iel,action
 END DO elements_3
 CALL beamdis(loads,nf,0.40_iwp,10,nels,ell,argv,nlen,12)
END SUBROUTINE p43