text:
module ionization_mpi

use phys_consts, only: wp,debug
use grid, only: lx1,lx2,lx3,g1,g2,g3
use mpimod, only: mpi_realprec, mpi_cfg, tag=>gemini_mpi, MPI_COMM_WORLD,MPI_STATUS_IGNORE
use neutral, only: Tnmsis

implicit none (type, external)
external :: mpi_send,mpi_recv

contains
  !> Query workers to get a single value for Tninf and gavg (or makes something up)
  subroutine get_gavg_Tinf(gavg,Tninf)
    real(wp), intent(out) :: gavg,Tninf
    real(wp) :: Tninftmp
    integer :: iid, ierr
    real(wp), dimension(:,:,:), allocatable :: g

    ! use an average value for the gravitational field; FIXME: perhaps should be done via averaging over all workers???
    allocate(g(1:lx1,1:lx2,1:lx3))
    g=sqrt(g1**2+g2**2+g3**2)
    !    gavg=sum(g)/(lx1*lx2*lx3)    !single average value for computing column dens.  Interestingly this is a worker average...  Do we need root grav vars. grid mod to prevent tearing?  Should be okay as long as the grid is only sliced along the x3-dimension; problematic for x2divisions...
    gavg=8._wp

    Tninf=maxval(Tnmsis)   !set exospheric temperature based on the max value of the background MSIS atmosphere; note this is a worker max

    !both g and Tinf need to be computed as average over the entire grid...
    if (mpi_cfg%myid==0) then     !root
      ierr=0
      do iid=1,mpi_cfg%lid-1
          call mpi_recv(Tninftmp,1,mpi_realprec,iid,tag%Tninf,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
          if (Tninf < Tninftmp) Tninf=Tninftmp
      end do
      if (ierr /= 0) error stop 'root failed to mpi_recv Tninf'

      ierr=0
      do iid=1,mpi_cfg%lid-1
        call mpi_send(Tninf,1,mpi_realprec,iid,tag%Tninf,MPI_COMM_WORLD,ierr)
      end do
      if (ierr /= 0) error stop 'root failed to mpi_send Tninf'

      if (debug) print *, 'Exospheric temperature used for photoionization:  ',Tninf
    else                  !workders
      call mpi_send(Tninf,1,mpi_realprec,0,tag%Tninf,MPI_COMM_WORLD,ierr)                        !send what I think Tninf should be
      if (ierr /= 0) error stop 'worker failed to mpi_send Tninf'
      call mpi_recv(Tninf,1,mpi_realprec,0,tag%Tninf,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)      !receive roots decision
      if (ierr /= 0) error stop 'worker failed to mpi_recv Tninf'
    end if
    deallocate(g)
  end subroutine
end module ionization_mpi

max_forks_repo_name: mattzett/GEMINI


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

C=======================================================================
      SUBROUTINE GRTEXT (DX, DY, STRING)
C=======================================================================

C   --*** GRTEXT *** (GRPLIB) Write text (PLT)
C   --   Written by Amy Gilkey - revised 02/14/86
C   --
C   --GRTEXT writes a software or hardware character string at a location
C   --(left-justified).
C   --
C   --Parameters:
C   --   DX, DY - IN - the horizontal and vertical string location
C   --      (in device coordinates)
C   --   STRING - IN - the string to be written, may be truncated
C   --
C   --Common Variables:
C   --   Uses ICURDV, SOFTCH of /GRPCOM/

C   --Routines Called:
C   --   PLTXTH - (PLTLIB) Display a hardware string
C   --   PLTXTS - (PLTLIB) Display a software string
C   --   LENSTR - (STRLIB) Find string length

      COMMON /GRPCOC/ DEVNAM(2), DEVCOD(2)
      CHARACTER*3 DEVNAM
      CHARACTER*8 DEVCOD
      COMMON /GRPCOM/ ICURDV, ISHARD, DEVOK(2), TALKOK(2),
     &   NSNAP(2), IFONT(2), SOFTCH(2), AUTOPL(2),
     &   MAXCOL(2), NUMCOL(0:1,2), MAPALT(2), MAPUSE(2)
      LOGICAL ISHARD, DEVOK, TALKOK, SOFTCH, AUTOPL

      REAL DX, DY
      CHARACTER*(*) STRING

      LSTR = LENSTR(STRING)
      IF (STRING(LSTR:LSTR) .EQ. ' ') RETURN

      IF (SOFTCH(ICURDV)) THEN
         CALL PLTXTS (DX, DY, STRING(:LSTR))
      ELSE
         CALL PLTXTH (DX, DY, STRING(:LSTR))
      END IF

      RETURN
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
Mobile Home parks are the resting areas for mobile homes and trailers, not to be confused with RVs.

  Davis Mobile Estates on Olive Drive
  Rancho Yolo located in East Davis
  Royal Oak located in South Davis
  Slatters Court on Olive Drive

Of the past:
  Barthel Mobile Ranch in South Davis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine wrhead
************************************************************************
*
*   WRHEAD -- Write top portion of the output file.
*
************************************************************************
      include "parameters.h"

      integer nnod,ntri,tvtx
      common /mfe010/ nnod,ntri /mfe011/ tvtx(0:2,mtri)

      integer j,k

      open(7,file='mfeout',form='unformatted')
      rewind(7)

      write(7) npde,nnod,ntri
      write(7) ((tvtx(k,j), k=0,2), j=1,ntri)

      return
      end

max_stars_repo_name: mfeproject/legacy-gwmfe

max_issues_repo_name: mfeproject/legacy-gwmfe

max_forks_repo_name: mfeproject/legacy-gwmfe


text: #redirect ABT 49


text:
program prog2

implicit none
real, dimension (1:300) :: v, t, x, xf
real :: vo, a, cont
integer :: c

vo = 2
a = 10

cont = 0
do c = 1,300
    cont = cont + 0.01
    t(c) = cont
    v(c) = 2*vo*t(c) + (4d0/3d0)*a*t(c)**3
end do

open(1,FILE = "dados3.txt")

    do c = 1,300
	write(1,*) t(c), v(c)
    end do

close(1)

call integral (v, t, xf)

    open(2,FILE = "dados4.txt")
	do c =1,299
	    write(2,*) t(c), xf(c)
	end do

    close(2)

do c = 1,299
    write(*,*) t(c), xf(c), v(c)
end do

end program prog2

subroutine integral (v, t, xf)

    real, dimension(1:300) :: v, t, x, ret, tri, xf
    real :: cont
    integer :: c

    cont = 0

    do c = 1,299
	ret(c) = (t(c+1)-t(c))*v(c)
	tri(c) = ((t(c+1)-t(c))*(v(c+1)-v(c)))/2
	x(c) = tri(c) + ret(c)
	cont = cont + x(c)
	xf(c) = cont
    end do

return
end subroutine

max_stars_repo_name: j-rheinheimer/Fisica-Computacional-I

max_issues_repo_name: j-rheinheimer/Fisica-Computacional-I

max_forks_repo_name: j-rheinheimer/Fisica-Computacional-I


text:
      subroutine hepnmg(ID,CHAUG)
c
c  HEPNMG --
c
C...Purpose: to give the particle/parton name as a character string.
c
c SCCS ID: hepnmg.f 1.1 4/6/92
c
c Copyright (c) 1991 Universities Research Association, Inc.
c All rights reserved.
c
c This material resulted from work developed under a Government Contract and
c is subject to the following license:  The Government r~gh~nins a paid-up,
c nonexclusive, irrevocable worldwide license to reproduce, prepare derivative
c works, perform publicly and display publicly by or for the Government,
c including the right to distribute to other Government contractors.  Neither
c the United States nor the United States Department of Energy, nor any of
c their employees, makes any warrenty, express or implied, or assumes any
c legal liability or responsibility for the accuracy, completeness, or
c usefulness of any information, apparatus, product, or process disclosed, or
c represents that its use would not infringe privately owned rights.
c
c Written by Paul Lebrun & Lynn Garren
c
c	Arguments :
c
c	ID  :Standard particle name See Convention in Particle Data book,
c		Physics LetterB vol 239, p. III.67
c
c	CHAUG : Null terminated Character string, the name of
c		particle, using Greek fonts, as defined in module StringUtils.c
c		Limited to 30 characters
c
      IMPLICIT NONE

      INTEGER ID, IDEX, IDA, IDQ1, IDQ2, IDQ3, LEN, LEM, KQN
      INTEGER I, KQ, IDJS, IDCM
      CHARACTER CHAU*25, CHAUG*25
      BYTE ICHAU(25)
      EQUIVALENCE (ICHAU(1), CHAU(1:1))
      character CHSP(1605)*20

      integer hepchg,hepcmp
      external hepchg,hepcmp

      data (CHSP(I), I=1,119) /'d','u','s','c','b','t','l','h',2*' ',
     1'e','~gn~n~de~n','~gm~n','~gn~n~dm~n','~gt~n','~gn~n~dt~n',
     1'~gc~n','~gn~n~dc~n','reggeon','pomeron','g','~gg~n','Z','W',
     2'H~d1~n',' ','~gh~n~dtech~n','LQ','R','Z~d2~n','W~d2~n','H~d2~n',
     3'H~d2~n','H~d2~n','Z~d3~n','H~d3~n','H~d3~n','H~d3~n','H~d4~n',
     4'H~d5~n','L susy d','L susy u','L susy s','L susy c','L susy b',
     5'L susy t','gluino','photino',2*' ','L susy e',
     6'L susy ~gn~n~de~n','L susy ~gm~n','L susy ~gn~n~dm~n',
     7'L susy ~gt~n','L susy ~gn~n~dt~n','susy HL','susy HH','susy HA',
     8' ','R susy d','R susy u','R susy s','R susy c','R susy b',
     9'R susy t','zino','susy Z2','susy Z3','susy Z4','R susy e',
     *'R susy ~gn~n~de~n','R susy ~gm~n','R susy ~gn~n~dm~n',
     1'R susy ~gt~n','R susy ~gn~n~dt~n','wino','susy W2',2*' ',
     2 19*'gen. code',' ','K~ds~n','K~dl~n','diquark','l-baryon',
     3'h-baryon',' ','~gr~n~ddiffr~n','~gp~n~ddiffr~n','~gw~n~ddiffr~n',
     4'~gf~n~ddiffr~n','~gY~n~ddiffr~n','n~ddiffr~n','p~ddiffr~n',
     5'remnant ~gg~n','remnant nucleon',4*' '/
      data (CHSP(I), I=120,219) /'~gp~n','a~d0~n(980)','~gp~n(1300)',
     17*' ','~gr~n(770)','b~d1~n(1235)','a~d1~n(1260)','~gr~n(1450)',
     1'~gr~n(1700)',5*' ','a~d2~n(1320)','~gp~n~d2~n(1670)',8*' ',
     2'~gr~n~d3~n(1690)',19*' ','~gh~n','f~d0~n(975)','~gh~n(1295)',
     3'f~d0~n(1400)','~gh~n(1440)','f~d0~n(1590)','f~d0~n(1710)',3*' ',
     4'~gw~n(783)','h~d1~n(1170)','f~d1~n(1285)','~gw~n(1390)',
     5'f~d1~n(1420)','f~d1~n(1510)','~gw~n(1600)',3*' ','f~d2~n(1270)',
     6' ','f~d2~n(2010)','f~d2~n(2300)','f~d2~n(2340)',5*' ',
     7'~gw~n~d3~n(1670)',9*' ','f~d4~n(2050)',9*' '/
      data (CHSP(I), I=220,369) /'~gh~n''(958)','f''0',8*' ',
     1 '~gf~n(1020)','h~d1~n(1380)','f~d1~n(1420)','~gf~n(1680)',
     1 'f''~d1~n',5*' ','f''~d2~n(1525)', 9*' ','~gf~n~d3~n(1850)',
     2 19*' ','~gh~n~dc~n(1S)','~gc~n~dc0~n(1P)','~gh~n~dc~n(2S)',7*' ',
     3 'J/~gY~n','hc(1P)','~gc~n~dc1~n(1P)','~gY~n(2S)','~gY~n(3770)',
     4 '~gY~n(4040)','~gY~n(4160)','~gY~n(4415)',2*' ',
     4 '~gc~n~dc2~n(1P)',29*' ','~gh~n~db~n(1S)','~gc~n~db0~n(1P)',
     5 '~gh~n~db~n(2S)','~gc~n~db0~n(2P)','~gh~n~db~n(3S)',
     5 5*' ','~gU~n(1S)','h~db~n(1P)','~gc~n~db1~n(1P)',
     6 '~gU~n(2S)','h~db~n(2P)','~gc~n~db1~n(2P)','~gU~n(3S)',
     7 '~gU~n(4S)','~gU~n(10860)','~gU~n(11020)',
     8 '~gc~n~db2~n(1P)','~gc~n~db2~n(2P)',28*' '/
      data (CHSP(I), I=370,419) /'~gh~n~dt~n','~gc~n~dt0~n',8*' ',
     1'~gq~n','h~d1t~n','~gc~n~dt1~n',7*' ','~gc~n~dt2~n',29*' '/
      data (CHSP(I), I=420,449) /'~gh~n~dl~n','~gc~n~dl0~n',' ',
     2'~gq~n~dl~n','h~d1l~n','~gc~n~dl1~n','~gc~n~dl2~n',8*' ',
     3'~gh~n~dh~n','~gc~n~dh0~n',' ','~gq~n~dh~n','h~d1h~n',
     4'~gc~n~dh1~n','~gc~n~h2~n',8*' '/
      data (CHSP(I), I=450,499) /'~gp~n','a~d0~n(980)','~gp~n(1300)',
     1 7*' ','~gr~n(770)','b~d1~n(1235)','a~d1~n(1260)','~gr~n(1450)',
     2 '~gr~n(1700)',5*' ','a~d2~n(1320)','~gp~n~d2~n(1670)',
     2 8*' ','~gr~n~d3~n(1690)',19*' '/
      data (CHSP(I), I=500,599) /'K','K~u*~n~d0~n(1430)',8*' ',
     1'K~u*~n(892)','K~d1~n(1270)','K~d1~n(1400)','K~u*~n(1410)',
     1'K~u*~n(1680)',5*' ','K~u*~n~d2~n(1430)','K~d2~n(1770)',
     2'K~d2~n(1820)',7*' ','K~u*~n~d3~n(1780)',9*' ',
     3 'K~u*~n~d4~n(2045)',9*' ','K','K~u*~n~d0~n(1430)',8*' ',
     3'K~u*~n(892)','K~d1~n(1270)','K~d1~n(1400)','K~u*~n(1410)',
     3'K~u*~n(1680)', 5*' ','K~u*~n~d2~n(1430)','K~d2~n(1770)',
     4'K~d2~n(1820)',7*' ', 'K~u*~n~d3~n(1780)',9*' ',
     5'K~u*~n~d4~n(2045)',9*' '/
      data (CHSP(I), I=600,749) /'D','D~u*~n~d0~n',8*' ','D~u*~n(2010)',
     1'D~d1~n(2420)','D~u*~n~d1~n',7*' ','D~u*~n~d2~n(2460)',29*' ',
     1'D','D~u*~n~d0~n',8*' ','D~u*~n(2010)','D~d1~n(2420)',
     2'D~u*~n~d1~n',7*' ','D~u*~n~d2~n(2460)',29*' ','D~ds~n',
     3'D~u*~n~ds0~n',8*' ','D~u*~n~ds~n','D~ds1~n(2536)',
     4'D~u*~n~ds1~n',7*' ','D~u*~n~ds2~n',29*' '/
      data (CHSP(I), I=750,949) /'B','B~u*~n~d0~n',8*' ','B~u*~n',
     1'B~d1~n','B~u*~n~d1~n',7*' ','B~u*~n~d2~n',29*' ','B',
     1'B~u*~n~d0~n',8*' ','B~u*~n','B~d1~n','B~u*~n~d1~n',7*' ',
     2'B~u*~n~d2~n',29*' ','B~ds~n','B~u*~n~ds0~n',8*' ',
     3'B~u*~n~ds~n','B~ds1~n','B~u*~n~ds1~n',7*' ','B~u*~n~ds2~n',
     3 29*' ','B~dc~n','B~u*~n~dc0~n',8*' ','B~u*~n~dc~n',
     3'B~dc1~n','B~u*~n~dc1~n',7*' ','B~u*~n~dc2~n',29*' '/
      data (CHSP(I), I=950,1199) /'T','T~u*~n~d0~n',8*' ','T~u*~n',
     1'T~d1~n','T~u*~n~d1~n',7*' ','T~u*~n~d2~n',29*' ','T',
     1'T~u*~n~d0~n',8*' ','T~u*~n','T~d1~n','T~u*~n~d1~n',7*' ',
     2'T~u*~n~d2~n',29*' ','T~ds~n','T~u*~n~ds0~n',8*' ',
     3'T~u*~n~ds~n','T~ds1~n','T~u*~n~ds1~n',7*' ','T~u*~n~ds2~n',
     3 29*' ','T~dc~n','T~u*~n~dc0~n',8*' ','T~u*~n~dc~n','T~dc1~n',
     3'T~u*~n~dc1~n',7*' ','T~u*~n~dc2~n',29*' ','T~db~n',
     3 'T~u*~n~db0~n',8*' ','T~u*~n~db~n','T~db1~n','T~u*~n~db1~n',
     4 7*' ','T~u*~n~db2~n',29*' '/
      data (CHSP(I), I=1200,1289) /'L','L~u*~n~d0~n',' ','L~u*~n',
     1'L~d1~n','L~u*~n~d1~n','L~u*~n~d2~n',8*' ','L','L~u*~n~d0~n',
     1' ','L~u*~n','L~d1~n','L~u*~n~d1~n','L~u*~n~d2~n',8*' ','L~ds~n',
     2'L~u*~n~ds0~n',' ','L~u*~n~ds~n','L~ds1~n','L~u*~n~ds1~n',
     2'L~u*~n~ds2~n',8*' ','L~dc~n','L~u*~n~dc0~n',' ','L~u*~n~dc~n',
     3'L~dc1~n','L~u*~n~dc1~n','L~u*~n~dc2~n',8*' ','L~db~n',
     4'L~u*~n~db0~n',' ','L~u*~n~db~n','L~db1~n','L~u*~n~db1~n',
     4'L~u*~n~db2~n',8*' ','L~dt~n','L~u*~n~dt0~n',' ',
     5'L~u*~n~dt~n','L~dt1~n','L~u*~n~dt1~n','L~u*~n~dt2~n',8*' '/
      data (CHSP(I), I=1290,1394) /'H','H~u*~n~d0~n',' ','H~u*~n',
     1'H~d1~n','H~u*~n~d1~n','H~u*~n~d2~n',8*' ','H','H~u*~n~d0~n',
     1' ','H~u*~n','H~d1~n','H~u*~n~d1~n','H~u*~n~d2~n',8*' ','H~ds~n',
     2'H~u*~n~ds0~n',' ','H~u*~n~ds~n','H~ds1~n','H~u*~n~ds1~n',
     2'H~u*~n~ds2~n',8*' ','H~dc~n','H~u*~n~dc0~n',' ','H~u*~n~dc~n',
     3'H~dc1~n','H~u*~n~dc1~n','H~u*~n~dc2~n',8*' ','H~db~n',
     4'H~u*~n~db0~n',' ','H~u*~n~db~n','H~db1~n','H~u*~n~db1~n',
     4'H~u*~n~db2~n',8*' ','H~dt~n','H~u*~n~dt0~n',' ',
     5'H~u*~n~dt~n','H~dt1~n','H~u*~n~dt1~n','H~u*~n~dt2~n',8*' ',
     6'H~dl~n','H~u*~n~dl0~n',' ','H~u*~n~dl~n','H~dl1~n',
     7'H~u*~n~dl1~n','H~u*~n~dl2~n',8*' '/
      data (CHSP(I), I=1395,1415) /' ','~gL~n','~gL~n~dc~n',
     1'~gX~n~dc~n','~gX~n~dc~n','~gL~n~db~n','~gX~n~db~n',
     1'~gX~n~db~n','~gX~n~dbc~n','~gX~n~dbc~n','~gW~n~dbc~n',
     2'~gL~n~dt~n','~gX~n~dt~n','~gX~n~dt~n','~gX~n~dtc~n',
     3'~gX~n~dtc~n','~gW~n~dtc~n','~gX~n~dtb~n','~gX~n~dtb~n',
     4'~gW~n~dtb~n','~gW~n~dtbc~n'/
      data (CHSP(I), I=1416,1475) /' ','n','p',' ','~gS~n','~gS~n',
     1'~gS~n','~gX~n','~gX~n',' ','~gS~n~dc~n','~gS~n~dc~n',
     1'~gS~n~dc~n','~gX~n~dc~n''','~gX~n~dc~n''','~gW~n~dc~n',
     2'~gX~n~dcc~n','~gX~n~dcc~n','~gW~n~dcc~n',' ','~gS~n~db~n',
     3'~gS~n~db~n','~gS~n~db~n','~gX~n~db~n''','~gX~n~db~n''',
     3'~gW~n~db~n','~gX~n~dbc~n''','~gX~n~dbc~n''','~gW~n~dbc~n''',
     4'~gW~n~dbcc~n','~gX~n~dbb~n','~gX~n~dbb~n','~gW~n~dbb~n',
     5'~gW~n~dbbc~n',' ','~gS~n~dt~n','~gS~n~dt~n','~gS~n~dt~n',
     5'~gX~n~dt~n''','~gX~n~dt~n''','~gW~n~dt~n','~gX~n~dtc~n''',
     6'~gX~n~dtc~n''','~gW~n~dtc~n''','~gW~n~dtcc~n','~gX~n~dtb~n''',
     7'~gX~n~dtb~n''','~gW~n~dtb~n''','~gW~n~dtbc~n''','~gW~n~dtbb~n',
     7'~gX~n~dtt~n','~gX~n~dtt~n','~gW~n~dtt~n','~gW~n~dttc~n',
     8'~gW~n~dttb~n',5*' '/
      data (CHSP(I), I=1476,1535) /'~gD~n','~gD~n','~gD~n','~gD~n',
     1'~gS~n~u*~n','~gS~n~u*~n','~gS~n~u*~n','~gX~n~u*~n',
     1'~gX~n~u*~n','~gW~n','~gS~n~u*~n~dc~n','~gS~n~u*~n~dc~n',
     2'~gS~n~u*~n~dc~n','~gX~n~u*~n~dc~n','~gX~n~u*~n~dc~n',
     2'~gW~n~u*~n~dc~n','~gX~n~u*~n~dcc~n','~gX~n~u*~n~dcc~n',
     3'~gW~n~u*~n~dcc~n','~gW~n~u*~n~dccc~n','~gS~n~u*~n~db~n',
     4'~gS~n~u*~n~db~n','~gS~n~u*~n~db~n','~gX~n~u*~n~db~n',
     4'~gX~n~u*~n~db~n','~gW~n~u*~n~db~n','~gX~n~u*~n~dbc~n',
     5'~gX~n~u*~n~dbc~n','~gW~n~u*~n~dbc~n','~gW~n~u*~n~dbcc~n',
     6'~gX~n~u*~n~dbb~n','~gX~n~u*~n~dbb~n','~gW~n~u*~n~dbb~n',
     6'~gW~n~u*~n~dbbc~n','~gW~n~u*~n~dbbb~n','~gS~n~u*~n~dt~n',
     7'~gS~n~u*~n~dt~n','~gS~n~u*~n~dt~n','~gX~n~u*~n~dt~n',
     7'~gX~n~u*~n~dt~n','~gW~n~u*~n~dt~n','~gX~n~u*~n~dtc~n',
     8'~gX~n~u*~n~dtc~n','~gW~n~u*~n~dtc~n','~gW~n~u*~n~dtcc~n',
     8'~gX~n~u*~n~dtb~n','~gX~n~u*~n~dtb~n','~gW~n~u*~n~dtb~n',
     9'~gW~n~u*~n~dtbc~n','~gW~n~u*~n~dtbb~n','~gX~n~u*~n~dtt~n',
     9'~gX~n~u*~n~dtt~n','~gW~n~u*~n~dtt~n','~gW~n~u*~n~dttc~n',
     *'~gW~n~u*~n~dttb~n','~gW~n~u*~n~dttt~n',4*' '/
      data (CHSP(I), I=1536,1600) /65*' '/
      data (CHSP(I), I=1601,1605) /'Hydrogen','Deuteron','Tritium',
     1'He~n~u3~n','Alpha'/

      save CHSP

C...Initial values. Charge (in units of 1/3). Subdivide code.
      CHAU=' '
      LEN=1
      IDA=IABS(ID)
      if(IDA.EQ.0) go to 900

      IDCM=hepcmp(ID)
      if(IDCM.LE.0) go to 900

      KQ=hepchg(ID)
      KQN=MOD(IDA/1000000000,10)
      IDQ3=MOD(IDA/1000,10)
      IDQ2=MOD(IDA/100,10)
      IDQ1=MOD(IDA/10,10)
      IDJS=MOD(IDA,10)
      IDEX=MOD(IDA/10000,10)

C...Read out root name and spin for simple particles and special cases
      if(IDA.LE.100) then
        CHAU=CHSP(IDCM)
        LEN=1
        do 100 LEM=1,20
  100     if(CHAU(LEM:LEM).NE.' ') LEN=LEM
      elseif(KQN.eq.1) then
        CHAU=CHSP(IDCM)
        LEN=0
        do 110 LEM=1,20
  110     if(CHAU(LEM:LEM).NE.' ') LEN=LEM
      elseif(IDJS.EQ.0) then
        CHAU=CHSP(IDCM)
        LEN=1
        do 150 LEM=1,20
  150     if(CHAU(LEM:LEM).NE.' ') LEN=LEM

C...Construct root name for diquark. Add on spin.
      elseif(IDQ1.EQ.0) then
        CHAU(1:2)=CHSP(IDQ3)(1:1)//CHSP(IDQ2)(1:1)
        if(IDJS.EQ.1) CHAU(3:7)='~d0~n'
        if(IDJS.EQ.3) CHAU(3:7)='~d1~n'
        LEN=7

C...Construct root name for meson.
      elseif(IDQ3.EQ.0) then
        CHAU=CHSP(IDCM)
        LEN=1
        do 200 LEM=1,20
  200     if(CHAU(LEM:LEM).NE.' ') LEN=LEM

C...Construct root name for Herwig remnant particles
      elseif(IDQ1.eq.9 .and. IDQ2.eq.9 .and. IDQ3.eq.9) then
        CHAU=CHSP(IDCM)
        LEN=1
        do 250 LEM=1,20
  250     if(CHAU(LEM:LEM).NE.' ') LEN=LEM

C...Construct root name and spin for baryon.
      else
        if(IDQ3.LE.6) then
          CHAU=CHSP(IDCM)
          LEN=1
          do 300 LEM=1,20
  300       if(CHAU(LEM:LEM).NE.' ') LEN=LEM
        else
C...Construct root name and spin for heavy baryon.
          LEN=1
          if(IDQ2.LE.2 .AND. IDQ1.LE.2)then
            CHAU='~gS~n'
            if(IDJS.EQ.4) CHAU='~gS~n~u*~n'
            if(IDQ1.GT.IDQ2) CHAU='~gL~n'
          elseif(IDQ1.LE.2) then
            CHAU='~gX~n'''
            if(IDJS.EQ.4) CHAU='~gX~n~u*~n'
          elseif(IDQ2.LE.2) then
            CHAU='~gX~n'
          elseif(IDQ1.GT.IDQ2) then
            CHAU='~gW~n'
          else
            CHAU='~gW~n'''
            if(IDJS.EQ.4) CHAU='~gW~n~u*~n'
          endif
          do 320 LEM=1,20
  320       if(CHAU(LEM:LEM).NE.' ') LEN=LEM

C...Add on heavy flavour content for heavy baryon.
          CHAU(LEN+1:LEN+5)='~d'//CHSP(IDQ3)(1:1)//'~n'
          LEN=LEN+5
          if(IDQ2.GE.IDQ1.AND.IDQ1.GE.4) then
            CHAU(LEN+1:LEN+2)=CHSP(IDQ2)(1:1)//CHSP(IDQ1)(1:1)
            LEN=LEN+2
          elseif(IDQ2.GE.IDQ1.AND.IDQ2.GE.4) then
            CHAU(LEN+1:LEN+1)=CHSP(IDQ2)(1:1)
            LEN=LEN+1
          elseif(IDQ1.GT.IDQ2.AND.IDQ2.GE.4) then
            CHAU(LEN+1:LEN+2)=CHSP(IDQ1)(1:1)//CHSP(IDQ2)(1:1)
            LEN=LEN+2
          elseif(IDQ1.GT.IDQ2.AND.IDQ1.GE.4) then
            CHAU(LEN+1:LEN+1)=CHSP(IDQ1)(1:1)
            LEN=LEN+1
          endif
        endif
      endif

C...Add on bar sign for antiparticle (where necessary).
      if(ID.GT.0.OR.LEN.EQ.0) then
      elseif(IDA.GT.10.AND.IDA.LE.40.AND.KQ.NE.0) then
      elseif(IDA.GT.46.AND.IDA.LE.60.AND.KQ.NE.0) then
      elseif(IDA.GT.66.AND.IDA.LE.80.AND.KQ.NE.0) then
      elseif(IDA.GE.81.AND.IDA.LE.99) then
      elseif(IDA.GT.100.AND.IDQ3.EQ.0.AND.KQ.NE.0) then
      else
        CHAU(LEN+1:LEN+2)='~~'
        LEN=LEN+2
      endif

C...Add on charge where applicable (conventional cases skipped).
      if(IDA.GE.81.AND.IDA.LE.100) then
C...generator specific codes
      elseif(IDCM.eq.114 .or. IDCM.eq.115) then
C...Herwig remnant particles
      elseif(KQN.eq.1) then
C...ions
      elseif(LEN.GT.0)then
C...everything else
        if(abs(KQ).EQ.6)then
          if(KQ.EQ.6) CHAU(LEN+1:LEN+6)='~u++~n'
          if(KQ.EQ.-6) CHAU(LEN+1:LEN+6)='~u--~n'
          LEN = LEN + 6
        elseif(abs(KQ).eq.3)then
          if(KQ.EQ.3) CHAU(LEN+1:LEN+4)='~u+~n'
          if(KQ.EQ.-3) CHAU(LEN+1:LEN+4)='~u-~n'
          LEN = LEN + 4
        endif

        if(KQ.EQ.0.AND.IDA.LE.22) then
C...quarks and leptons
        elseif(KQ.EQ.0.AND.(IDA.GE.41.AND.IDA.LE.56)) then
C...left squarks, sleptons, etc.
        elseif(KQ.EQ.0.AND.(IDA.GE.61.AND.IDA.LE.67)) then
C...right squarks
        elseif(KQ.EQ.0.AND.(IDA.GE.71.AND.IDA.LE.76)) then
C...right sleptons
        elseif(IDA.GT.100.AND.IDQ3.EQ.0.AND.IDQ2.EQ.IDQ1.AND.
     &      IDQ2.NE.1) then
C... eta, psi, upsilon, etc.
        elseif(KQ.EQ.0) then
          CHAU(LEN+1:LEN+5)='~u0~n'
          LEN = LEN + 5
        endif
      endif

900   continue
      ICHAU(LEN+1) = 0
      CHAUG = CHAU
      return
      end


      subroutine reservoir_output(j)

      use time_module
      use basin_module
      use reservoir_module
      use hydrograph_module
      use water_body_module

      implicit none

      integer, intent (in) :: j   !                |
      integer :: iob              !                |
      real :: const               !                |

      iob = sp_ob1%res + j - 1

!!!!! daily print
         if (pco%day_print == "y" .and. pco%int_day_cur == pco%int_day) then
          if (pco%res%d == "y") then
            write (2540,100) time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, res_wat_d(j), res(j), &
                res_in_d(j), res_out_d(j)
               if (pco%csvout == "y") then
                 write (2544,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, &
                    res_wat_d(j), res(j), res_in_d(j), res_out_d(j)
             end if
          end if
         end if

        res_in_m(j) = res_in_m(j) + res_in_d(j)
        res_out_m(j) = res_out_m(j) + res_out_d(j)
        res_wat_m(j) = res_wat_m(j) + res_wat_d(j)

!!!!! monthly print
        if (time%end_mo == 1) then
          const = float (ndays(time%mo + 1) - ndays(time%mo))
          !res_in_m(j)%flo = res_in_m(j)%flo / const
          res_in_y(j) = res_in_y(j) + res_in_m(j)
          !res_out_m(j)%flo = res_out_m(j)%flo / const
          res_out_y(j) = res_out_y(j) + res_out_m(j)
          res_wat_m(j) = res_wat_m(j) // const
          res_wat_y(j) = res_wat_y(j) + res_wat_m(j)
          if (pco%res%m == "y") then
            write (2541,100) time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, res_wat_m(j), res(j), &
                res_in_m(j), res_out_m(j)
              if (pco%csvout == "y") then
                write (2545,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, &
                    res_wat_m(j), res(j), res_in_m(j), res_out_m(j)
              end if
          end if
          res_in_m(j) = resmz
          res_out_m(j) = resmz
          res_wat_m(j) = wbodz
        end if

!!!!! yearly print
       if (time%end_yr == 1) then
          res_in_a(j) = res_in_a(j) + res_in_y(j)
          res_out_a(j) = res_out_a(j) + res_out_y(j)
          res_wat_y(j) = res_wat_y(j) // 12.
          res_wat_a(j) = res_wat_a(j) + res_wat_y(j)
          if (pco%res%y == "y") then
            write (2542,100) time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, res_wat_y(j), res(j), &
                res_in_y(j), res_out_y(j)
              if (pco%csvout == "y") then
                write (2546,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, &
                    res_wat_y(j), res(j), res_in_y(j), res_out_y(j)
              end if
          end if
          res_in_y(j) = resmz
          res_out_y(j) = resmz
          res_wat_y(j) = wbodz
       end if

!!!!! average annual print
        if (time%end_sim == 1 .and. pco%res%a == "y") then
          res_in_a(j) = res_in_a(j) / time%yrs_prt
          res_out_a(j) = res_out_a(j) / time%yrs_prt
          res_wat_a(j) = res_wat_a(j) / time%yrs_prt
          write (2543,100) time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, res_wat_a(j), res(j), &
            res_in_a(j), res_out_a(j)
          if (pco%csvout == "y") then
            write (2547,'(*(G0.3,:","))')time%day, time%mo, time%day_mo, time%yrc, j, ob(iob)%gis_id, ob(iob)%name, &
                res_wat_a(j), res(j), res_in_a(j), res_out_a(j)
          end if
          res_in_a(j) = resmz
          res_out_a(j) = resmz
          res_wat_a(j) = wbodz
        end if

      return

100   format (4i6,2i10,2x,a,60e15.4)

      end subroutine reservoir_output


text:
MODULE ED_GREENS_FUNCTIONS
  USE ED_GF_SHARED
  USE ED_GF_NORMAL
  !USE ED_GF_CHISPIN
  !
  implicit none
  private

  public :: buildGf_impurity
  !public :: buildChi_impurity

contains

  !+------------------------------------------------------------------+
  ! GF CALCULATIONS
  !+------------------------------------------------------------------+
  subroutine buildGF_impurity()
    !
    call allocate_grids
    !
    impGmats=zero
    impGreal=zero
    !
    impSmats = zero
    impSreal = zero
    !
    impG0mats=zero
    impG0real=zero
    !
    !
    write(LOGfile,"(A)")"Get impurity Greens functions:"
    call build_gf_normal()
    call build_sigma_normal()
    !
    if(MPIMASTER)then
       if(ed_print_Sigma)call ed_print_impSigma()
       if(ed_print_G)call ed_print_impG()
       if(ed_print_G0)call ed_print_impG0()
    endif
    !
    call deallocate_grids
    !
  end subroutine buildgf_impurity

  !+------------------------------------------------------------------+
  ! SUSCEPTIBILITY CALCULATIONS
  !+------------------------------------------------------------------+
  !subroutine buildChi_impurity()
  !!
  !call allocate_grids
  !!
  !!
  !!BUILD SPIN SUSCEPTIBILITY
  !spinChi_tau=zero
  !spinChi_w=zero
  !spinChi_iv=zero
  !call build_chi_spin()
  !!
  !!
  !! !BUILD CHARGE SUSCEPTIBILITY
  !! densChi_tau=zero
  !! densChi_w=zero
  !! densChi_iv=zero
  !! densChi_mix_tau=zero
  !! densChi_mix_w=zero
  !! densChi_mix_iv=zero
  !! densChi_tot_tau=zero
  !! densChi_tot_w=zero
  !! densChi_tot_iv=zero
  !! call build_chi_dens()
  !!
  !!
  !! !BUILD PAIR SUSCEPTIBILITY
  !! pairChi_tau=zero
  !! pairChi_w=zero
  !! pairChi_iv=zero
  !! call build_chi_pair()
  !!
  !!
  !!PRINTING:
  !if(MPIMASTER)call ed_print_impChi()
  !!
  !!
  !call deallocate_grids
  !!
  !end subroutine buildChi_impurity

end MODULE ED_GREENS_FUNCTIONS

max_stars_repo_name: QcmPlab/CDMFT-LANC-ED

max_issues_repo_name: QcmPlab/CDMFT-LANC-ED

max_forks_repo_name: QcmPlab/CDMFT-LANC-ED


text:
! dynamic liveness analysis
      program dla
      real A(10,10)
      integer i, j, answer
!hpf$ dynamic A, T
!hpf$ processors P(2,2)
!hpf$ processors Q(4)
!hpf$ template T(10,10)
!hpf$ align A with T
!hpf$ distribute T(block,block) onto P
!hpf$ independent(j,i)
      do j=1, 10
         do i=1, 10
            A(i,j) = i+j
         enddo
      enddo
      print *, '1 or 0'
      read *, answer
      if (answer.eq.1) then
!
! initial A is not live here, because A is reset.
!
!hpf$    redistribute T(block,*) onto Q
!hpf$    independent(j,i)
         do j=1, 9
            do i=1, 10
               A(i,j) = A(i,j)+A(i,j+1)
            enddo
         enddo
      else
!
! initial A remains live from here
!
!hpf$    redistribute T(*,block) onto Q
         print *, A(5,5)
      endif
!
! bask to initial distribution
!
!hpf$ redistribute T(block,block) onto P
!hpf$ independent(j,i)
      do j=1, 10
         do i=1, 10
            A(i,j) = A(i,j)+1
         enddo
      enddo
      end


text:
AMD Medical Supply provides medical services durable medical supplies, including sales, rentals, and installation services.  They can do free local delivery on purchases of $100 or more.  A second location can be found at 3108 Arden Way in Sacramento. They are now also selling great, very tasty 100% Pure Local Honey from a beekeeper in Woodland.

They sell stethoscopes, scrubs, accessories for healthcare professionals, blood pressure monitors, otoscopes, wheelchairs, walkers, compression stockings & arm sleeves, lift chairs, bath chairs & safety items, bed rails, hospital beds & mattresses, orthopedic braces, arm slings, toilet risers & rails, underpads, chux, adult underwear & diapers, and much more.

Stuff

 Rental items

 For sale

 Related Links
 Medical Services

20130607 12:51:34 nbsp Im sorry. But this just smacks of a pure advertisement. I dont understand how the DavisWiki differentiates this from commercialization and its community mission. Print this in the Enterprise. I have been criticized for puffery and bloated content before. Whats with this? Comments welcome. Users/fknochenhauer

20130607 13:48:17 nbsp Put this company on Yelp! Thats all I have to say. Users/fknochenhauer

20130613 08:20:40 nbsp fknochenhauer: My name is Dmitriy, I am the manager of AMD Medical Supply here in Davis. Our DavisWiki page is not designed to be a pure advertisement. Its purpose is to be informative for Davis residents. Our page is very community oriented and serves a clear community mission: we are the ONLY medical supply store in Davis and we exist and strive to improve the health and wellbeing of our customers in Davis. Thus far, we have only heard positive feedback regarding our mission. I know that medical supply stores have opened in Davis in the past  but all unfortunately closed their doors very soon. I believe this can be attributed to a number of reasons, one of which is a very high cost of rent. Thankfully, with the support of our Davis community we have remained open for over 1 year and plan to stay in Davis for years to come. We opened our Davis store based on a  very much vocalized demand  we had many Davis residents driving all the way out to our store in Sacramento to purchase medical supplies and personally asked that we open a store in Davis. Because my family and I are originally from Yolo County and have deep roots in Davis over 40 years old, we felt that we needed to help Davis residents in need of medical supplies. In opening a medical supply store in Davis, our goal was that residents could have their own, local, friendly, and high quality medical supply store and enjoy this convenience. Prior to our opening, the nearest medical supply store was in Woodland or in Sacramento. It is terrible that our city did not have its own permanent medical supply store prior to our opening. We have a huge medical community with a top hospital, clinics, providers, care homes, caregivers, and students studying to give back professionally to our community  all residing in Davis and serving a large elderly population. It is unfortunate that you so strongly oppose our information being available to those who need it on DavisWiki, however I hope that my short explanation will help you to understand our sincere purpose in serving our Davis community. Users/dmitriymarkiv

20130905 17:35:15 nbsp Beware: this store does not give refunds.  Yes, it says so on your receipt and there is a discreet little disclaimer on the wall next to the cash register.  It would be good business, however, to let the customer know this verbally.  It isnt like there is ever a long line. It never occurred to me that I could not return something that was totally unused (purchased only two hours previously)and not a medically sensitive device.  I was able to exchange the device for another product  so what is the difference?  they still took the first product back.  It also took the clerk several minutes and a calculator to figure out that there was one dollar difference in the price.  It will definitely be the last time I purchase anything from them.  Users/JamieMBoston


text:
C--------------------------------------------------------------------------------------
      SUBROUTINE PARTPRESS
C
C     June,2011  MB
C     currently works only for atoms, not locked in molecules
C     all PP's and PF's from the model atmosphere are kept fixed
C
C     Note: the atomic number is here mapped to the equivalent array-element index,
C           thus also in a model atmosphere all atm.h.eps (element abundances) must
C           be present, but set to zero if UION and PION are not available for these
C           elements
C
C     global variables from lf_param.inc are SUPERSCRIPT
C
C--------------------------------------------------------------------------------------
C
      INTEGER*4   nion_temp(NATOM)
      REAL*4      idx_el_temp(NATOM), eps_temp(NATOM), ref_el(NATOM)
      CHARACTER*2 cel_temp(NATOM)
      REAL*4      uion_temp(NDMAX,NR_ION,NATOM)
      REAL*4      pion_temp(NDMAX,NR_ION,NATOM)

c      INCLUDE 'physcnst.inc'
      INCLUDE 'lf_decl0.inc'        ! includes lf_param.inc, atm_def.inc
c      INCLUDE 'iondat.inc'
c      INCLUDE 'ppcommon.inc'

C reference index for each atom is its atomic number plus ionization stage

      ref_el = atm.h.idx_el         ! save model atmosphere atoms for further cross-check

      k = 0
      DO i = 1, NATOM               ! 92 (lf_param.inc) ==> use as ATOMIC NUMBER identifier
        idx_el_temp(i) = i          ! atomic mass
        nelpos(idx_el_temp(i)) = i  ! == index number
        cel_temp(i)    = cel(i)     ! Element-Label
        nion_temp(i)   = nion(i)    ! Nr of ionisation stages
        eps_temp(i)    = epss(i)    ! solar abundance from elemdata.inc

        DO m = 1, nion_temp(i)      ! ionization stage
          DO n = 1, atm.h.ndepth
            uion_temp(n,m,i) = 0.
            pion_temp(n,m,i) = 0.
          ENDDO
        ENDDO

        DO j = 1, atm.h.nel              ! 16 in all MARCS models
        IF (k.GT.16) WRITE(*,*) 'Something wrong in partpress.f'
        IF (atm.h.idx_el(j).EQ.idx_el_temp(i)) THEN
           k = k + 1
           eps_temp(i)  = atm.h.eps(i)   ! atm.h.neps = 63
           nion_temp(i) = atm.h.nion(j)
          DO l = 1, nion_temp(i)    ! only one element
           uion_temp(n,l,i) = atm.uion(n,atm.h.nion(j),j)
           pion_temp(n,l,i) = atm.pion(n,atm.h.nion(j),j)
          ENDDO
         ENDIF
        ENDDO

      ENDDO                              ! finish re-asssignment of arrays

      WRITE(*,*)
     *   'For ',NATOM-k,' elements UION and PION will be computed'
      WRITE(*,100) (eps_temp(ii),  ii = 1, NATOM)
 100  FORMAT(/,'New abundances',10(F5.2,2X))

      atm.h.nel  = NATOM            ! update total number of elements
      atm.h.nion = nion_temp        ! update ion stages
      atm.h.eps  = eps_temp         ! update abundance array
      atm.h.cel  = cel_temp         ! update element labels
      atm.h.idx_el = idx_el_temp    ! update atomic number array

C --- load atomic partition functions

      CALL ATOM_PARTFUN(atm.h.ndepth,atm.t,atm.h.nel,atm.h.nion)

      DO i = 1, NATOM
        DO m = 1, nion_temp(i)      ! ionization stage
          DO n = 1, atm.h.ndepth
           IF (uion_temp(n,m,i).EQ.0.0.AND.pion_temp(n,m,i).EQ.0.0) THEN
            atm.uion(n,m,i) = uion(n,m,i)       ! new
           ELSE
            atm.uion(n,m,i) = uion_temp(n,m,i)  ! model atmosphere
          ENDDO
        ENDDO
      ENDDO                         ! finish re-asssignment of arrays

C --- compute partial pressures for atoms

      CALL SAHA

      RETURN
      END

C---------------------------------------------------------
      SUBROUTINE ATOM_PARTFUN(ndepth,t,nel,nion,uion,chi)
C---------------------------------------------------------

C --- Calculation of Atomic(Ionic) Partition functions UION
C --- by linear interpolation of tables (iondat.inc)
C --- from KURUCZ

      PARAMETER (FAC=1.5664e-02 * 11607.5)
      INTEGER*4 ndepth,nion(nel)
      REAL*4	t(ndepth)

      INCLUDE 'lf_param.inc'
      INCLUDE 'ppcommon.inc'

      DO n = 1,ndepth
       DO m = 1,nel
	DO iz = 1,nion(m)
	  t_chi	       = FAC * chi(iz,m)
	  tk	       = t(n) / t_chi - 0.5
	  j	       = MAX(1,MIN(9,INT(tk)))
	  du	       = u_tab(j+1,iz,m) - u_tab(j,iz,m)
	  uion(n,iz,m) = u_tab(j,iz,m) + du * (tk - j)
	END DO
       END DO
      END DO
      RETURN
      END
C
C---------------------------------------------------------
      SUBROUTINE IONIS(te,divte,tp,pe,n)
C---------------------------------------------------------

C --- CALCULATION OF THE IONIZATION EQUILIBRIA
C --- ACCORDING TO THE SAHA - EQUATION

      INCLUDE 'physcnst.inc'
      INCLUDE 'lf_decl0.inc'
      INCLUDE 'iondat.inc'
      INCLUDE 'ppcommon.inc'

      tpe = tp / pe
      nh  = nelpos(1)

C-----IONIZATION EQUILIBRIA FOR EVERY ELEMENT, EXCEPT HYDROGEN

      DO m = 1, atm.h.nel
         id     = atm.h.idx_el(m)    ! atomic-numbers of used elements
         nm     = atm.h.nion(id)     ! Ionisation steps taken into account
         s(m)   = 1.
         zs(m)  = 0.
         zzs(m) = 0.
         fmz    = 1.
         IF ( nh .EQ. m ) nm = 2
         DO i   = 2, nm
            iz  = i - 1
            fmz = fmz * tpe * uion(n,i,id) / uion(n,iz,id) *
     *            EXP( - divte * chi(iz,id) )
            dfmz       = FLOAT(iz) * fmz
            sion(iz,m) = fmz            ! ionization fraction N(i+1)/N(i)
            s(m)       = s(m)   + fmz   ! total N
            zs(m)      = zs(m)  + dfmz
            zzs(m)     = zzs(m) + FLOAT(iz) * dfmz
         ENDDO
      ENDDO

C-----IONIZATION EQUILIBRIUM OF HYDROGEN INCLUDING H-MINUS

      sion2nh = EXP( divte * chi(3,1) ) / ( tpe * uion(n,1,1) )
      sion(2,nh) = sion2nh
      s(nh)   = s(NH)   + sion2nh
      zs(mh)  = zs(NH)  - sion2nh
      zzs(mh) = zzs(NH) + sion2nh

      RETURN
      END

C---------------------------------------------------------
      SUBROUTINE SAHA
C---------------------------------------------------------

      INCLUDE 'physcnst.inc'
      INCLUDE 'lf_decl0.inc'
      INCLUDE 'ppcommon.inc'

      REAL*4  te,tp,divte             ! local kinetic temperature
      REAL*4  pe                      ! electron pressure per depth
      REAL*4  pg                      ! gas pressure
      REAL*4  pat                     ! at pressure

      DO n   = 1, atm.h.ndepth

         te    = atm.t(n)
         pg    = atm.pg(n)
         divte = 1. / te
         tp    = .66667 * te ** 2.5
         pe    = atm.pe(n)
         pat   = pg - pe

        CALL IONIS(te,divte,tp,pe,n)

         DO k = 1, atm.h.nel
         m = nelpos(k)               ! index == atomic number
         a = eps(m) * pat/s(m)       ! total N for an element from IONIS
         pion(1,m,n) = a
         DO i = 2, nion(m)
            ii = i - 1
            pion(i,m,n) = a * sion(ii,m)   !sion - ionization fraction
         ENDDO
         ENDDO

      ENDDO

      RETURN
      END

max_stars_repo_name: andycasey/original-oracle

max_issues_repo_name: andycasey/original-oracle

max_forks_repo_name: andycasey/original-oracle


text:
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
      SUBROUTINE NSPXCP ( STRING, ERROR, SCREEN, LOGFIL )

      IMPLICIT NONE
C
C$ Version
C
C-     Command Loop Configured Version 1.0.0, 3-MAY-1994 (WLT)
C
C         This is the configured version of the Command Loop
C         software as of May 4, 1994
C

      CHARACTER*(*)         STRING
      CHARACTER*(*)         ERROR ( 2 )
      CHARACTER*(*)         SCREEN
      CHARACTER*(*)         LOGFIL

      EXTERNAL               NSPWLN

      CHARACTER*(128)       LSTYLE
      CHARACTER*(128)       SSTYLE
      CHARACTER*(128)       MARGIN

      INTEGER               I
      LOGICAL               SCRSTT ( 3 )
      LOGICAL               SAVSTT ( 3 )

      SAVE

      DATA LSTYLE / 'LEFT 1 RIGHT 78' /
      DATA SSTYLE / 'LEFT 1 RIGHT 78' /

      RETURN

      ENTRY NSPERR ( STRING, ERROR )
C$ Version
C
C-     Command Loop Configured Version 1.0.0, 3-MAY-1994 (WLT)
C
C         This is the configured version of the Command Loop
C         software as of May 4, 1994
C

C
C     This entry point is intended to be called once for
C     brief error diagnostics and a second time for more detailed
C     diagnostics.  We can tell which is which by examining the
C     first entry of the error array.  If it is non-blank this
C     must be the first such call (because we set it to blank
C     after we get done doing something with it).  A second
C     call can only happen if the special command
C     was entered by the user ('?').  In this case the command
C     manager will not reset the error array and not pass the
C     command to any other routines.  Instead it returns immediately
C     so that this routine can process the second part of the
C     error message.
C

      IF ( ERROR(1) .NE. ' ' ) THEN
C
C        We automatically clear the procedure stack whenever
C        an error occurs.
C
         CALL PRCLR ()

C
C        First inhibit writing to the log file.
C
         CALL NSPIOH ( 'LOG' )

C
C        Now write out only the first component of the error message.
C
         CALL NSPMRG   (            MARGIN         )
         CALL SUFFIX   ( SSTYLE, 1, MARGIN         )
         CALL NICEPR_1 ( ERROR(1),  MARGIN, NSPWLN )

C
C        Now inhibit writing to the screen or the save file. But
C        fetch their current state so that we can reset them
C        to exactly their current states.
C
         CALL NSPGST ( 'SCREEN', SCRSTT )
         CALL NSPGST ( 'SAVE',   SAVSTT )

         CALL NSPIOH ( 'SCREEN' )
         CALL NSPIOH ( 'SAVE'   )

C
C        Reactivate the log file.
C
         CALL NSPIOA ( 'LOG'    )

         DO I = 1, 2
            CALL NICEPR_1 ( ERROR(I), LSTYLE, NSPWLN )
         END DO

         CALL NSPPST ( 'SCREEN', SCRSTT )
         CALL NSPPST ( 'SAVE',   SAVSTT )

         ERROR(1) = ' '
         RETURN

      END IF

C
C     The only way to get here is for the user to have processed
C     the first half of an error and typed a question mark or
C     blank command. (This relies on all kinds of side effects.
C     Better talk to Bill if you want to be able to figure this out).
C

      IF ( STRING   .EQ. '?' ) THEN

         IF (ERROR(2) .EQ. ' ' ) THEN
            CALL TRNLAT ( 'NOMOREDIAGNOSTICS', ERROR(2) )
         END IF
C
C        We've already written the second part of the error
C        message to the log file, so we shall inhibit writing
C        there now.
C
         CALL NSPIOH ( 'LOG' )

         CALL NSPMRG   (            MARGIN         )
         CALL SUFFIX   ( SSTYLE, 1, MARGIN         )
         CALL NICEPR_1 ( ERROR(2),  MARGIN, NSPWLN )
C
C        Now re-activate the log file.
C
         CALL NSPIOA ( 'LOG' )

         ERROR(2) = ' '

      END IF

      RETURN

C
C     Set the style string that shall be used for printing
C     errors.
C
      ENTRY NSPSTY ( SCREEN, LOGFIL  )
C
C$ Version
C
C-     Command Loop Configured Version 1.0.0, 3-MAY-1994 (WLT)
C
C         This is the configured version of the Command Loop
C         software as of May 4, 1994
C

      SSTYLE = SCREEN
      LSTYLE = LOGFIL

      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
MODULE module_sticky
 !===================================
 ! Global variables for sticky atoms
 !===================================
 IMPLICIT NONE

 REAL(8)         :: sticky_k
 LOGICAL         :: sticky_atoms_wanted

 INTEGER, ALLOCATABLE :: sticky_labels(:)
 REAL(8), ALLOCATABLE :: sticky_fx(:), sticky_fy(:), sticky_fz(:)

END MODULE module_sticky

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
!! SR103029: fixed in 6.2 build 6228
!!
!! TYPE BOUND CONTIGUOUS POINTER FUNCTION NOT RECOGNIZED AS CONTIGUOUS
!!
!! In this example the type bound function contig_res returns a simply
!! contiguous array pointer. However when the function result is an actual
!! argument for a contiguous dummy argument, the 6.2 compiler does not
!! recognize it as simply contiguous. The example compiles fine with 6.1.
!!
!! For the error to occur: 1) the function needs to be type bound. No error
!! if the function is not type bound; and 2) the dummy argument needs to be
!! a pointer. No error if the dummy is not a pointer.
!!
!! $ nagfor -w nag-20190302.f90
!! NAG Fortran Compiler Release 6.2(Chiyoda) Build 6227
!! Error: nag-20190302.f90, line 33: Argument B (no. 1) of CONTIG_ARG is a CONTIGUOUS pointer, but the actual argument BAR%CONTIG_RES() is not simply contiguous
!!

module mod
  type :: foo
  contains
    procedure :: contig_res
  end type
contains
  function contig_res(this) result(a)
    class(foo), intent(in) :: this
    real, pointer, contiguous :: a(:)
    allocate(a(10))
  end function
end module

program main
  use mod
  type(foo) :: bar
  call contig_arg(bar%contig_res())
contains
  subroutine contig_arg(b)
    real, pointer, contiguous, intent(in) :: b(:)
  end subroutine
end program

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

module hecmw_mat_id
  use hecmw_util

  private

  public:: hecmw_mat_id_set
  public:: hecmw_mat_id_get
  public:: hecmw_mat_id_clear

  type mat_mesh
    logical :: used = .false.
    type(hecmwST_matrix), pointer :: mat
    type(hecmwST_local_mesh), pointer :: mesh
  end type mat_mesh

  integer(kind=kint), parameter :: MAX_MM = 8

  type(mat_mesh), save :: mm(MAX_MM)

contains

  subroutine hecmw_mat_id_set(hecMAT, hecMESH, id)
    implicit none
    type(hecmwST_matrix), intent(in), target :: hecMAT
    type(hecmwST_local_mesh), intent(in), target :: hecMESH
    integer(kind=kint), intent(out) :: id
    integer(kind=kint) :: i
    id = 0
    do i = 1, MAX_MM
      if (.not. mm(i)%used) then
        id = i
        exit
      endif
    end do
    if (id == 0) then
      stop 'ERROR: hecmw_mat_id_set: too many matrices set'
    endif
    mm(id)%mat => hecMAT
    mm(id)%mesh => hecMESH
    mm(id)%used = .true.
  end subroutine hecmw_mat_id_set

  subroutine hecmw_mat_id_get(id, hecMAT, hecMESH)
    implicit none
    integer(kind=kint), intent(in) :: id
    type(hecmwST_matrix), pointer :: hecMAT
    type(hecmwST_local_mesh), pointer :: hecMESH
    if (id <= 0 .or. MAX_MM < id) then
      stop 'ERROR: hecmw_mat_id_get: id out of range'
    endif
    if (.not. mm(id)%used) then
      stop 'ERROR: hecmw_mat_id_get: invalid id'
    endif
    hecMAT => mm(id)%mat
    hecMESH => mm(id)%mesh
  end subroutine hecmw_mat_id_get

  subroutine hecmw_mat_id_clear(id)
    implicit none
    integer(kind=kint), intent(in) :: id
    if (.not. mm(id)%used) then
      stop 'ERROR: hecmw_mat_id_clear: invalid id'
    endif
    mm(id)%mat => null()
    mm(id)%mesh => null()
    mm(id)%used = .false.
  end subroutine hecmw_mat_id_clear

end module hecmw_mat_id


text:
        SUBROUTINE ZCLEXP(READNO,NBINS,DATA,EPS,ET,ETEPS,ERR,EXPO)
*
*  Module number:
*
*  Module name: zclexp
*
*  Keyphrase:
*  ----------
*       convert to count rates
*
*  Description:
*  ------------
*       This routine converts the input data to count rates by
*       dividing by the exposure times.  The exposure time is
*       computed for each bin as 0.05 x NCOADD x intper, where
*       NCOADD is the number of coadds to the bin and intper
*       is the integration period in 0.05 millisecond intervals.
*       Both intper and NCOADD are found in the eng. trailer
*       array ET.  If either value contains fill, no exposure time
*       can be computed and the entire bin is flagged as unusable.
*
*  FORTRAN name: zclexp
*
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*
*  SDAS:
*       ZMSPUT
*  Others:
*
*
*  History:
*  --------
*  Version      Date        Author          Description
*       1       Feb. 89   D. Lindler    Designed and coded
*-------------------------------------------------------------------------------
*
* Input parameters
*
*       READNO - Readout number
*       NBINS - number of substep bins
*
* Input/Output parameters
*
*       DATA - data array for main diodes
*       EPS - data quality array for data
*       ET - eng. trailer and special diodes
*       ETEPS - data quality for ET
*       ERR - error array
*
* Output parameters
*
*       EXPO - vector of exposure times for each bin
*
********************************************************************************
        INTEGER NBINS,READNO
        REAL DATA(500,NBINS),EPS(500,NBINS),ET(24,NBINS),
     *          ERR(500,NBINS),ETEPS(24,NBINS),EXPO(7)
C
C     ZMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C HRS epsilons
C
        INTEGER EPSFIL
        PARAMETER (EPSFIL = 800)
C
C LOCAL VARIABLES
C
        CHARACTER*80 CONTXT
        INTEGER INT,NCOADD,INTEG,I,J,ISTAT
        REAL E
C
C------------------------------------------------------------------------------
C
C LOOP ON BINS
C
        DO 100 I=1,NBINS
C
C compute exposure time for the bin
C
            IF((ETEPS(21,I).GT.0.0).OR.(ETEPS(23,I).GT.0.0))THEN
                WRITE(CONTXT,9)I,READNO
9               FORMAT('ERROR: unable to compute exposure time for',
     *                  ' bin:',I2,'  readout:',I6)
                CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
                CALL ZMSPUT('     Data flagged as fill',
     *                  STDOUT+STDERR,0,ISTAT)
C
C trailer contains fill data (no way to compute exposure)
C
                DO 10 J=1,500
                        EPS(J,I)=EPSFIL
                        DATA(J,I)=0.0
                        ERR(J,I)=0.0
10              CONTINUE
                DO 20 J=1,12
20                      ETEPS(J,I)=EPSFIL
                EXPO(I)=0.0
C
C we can compute an exposure time
C
             ELSE
                INT=ET(21,I)
                CALL ZGETBT(INT,8,15,INTEG)
C                                    --->INTEGRATION TIME
                NCOADD=ET(23,I)
                IF(NCOADD.EQ.0)NCOADD=1
C                                    --->DDLINK HAS ZERO RECORDED
                E=(0.05*INTEG-0.002)*NCOADD
                EXPO(I)=E
C
C divide by the exposure time
C
                DO 30 J=1,500
                    IF(EPS(J,I).NE.EPSFIL)THEN
                        DATA(J,I)=DATA(J,I)/E
                        ERR(J,I)=ERR(J,I)/E
                    ENDIF
30              CONTINUE
C
C special diodes
C
                DO 40 J=1,12
40                  IF(ETEPS(J,I).NE.EPSFIL)ET(J,I)=ET(J,I)/E
            ENDIF
100     CONTINUE
        IF(READNO.EQ.1)CALL ZMSPUT('Data converted to count rates',
     *           STDOUT,0,ISTAT)
        RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
program t
  integer::a=1,b=2
  write (*,'(i0,i0)') a,b
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


program test_datetime
    use fplus_datetime
    use fplus_map
    implicit none

    type(datetime) :: now
    type(map) :: dates
    character(len=:), allocatable :: text
    integer :: i, test1, test2

    now = new_datetime()

    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()

    ! put the date into a map
    dates = new_map(10)
    call dates%add(now, 1, .true.)

    ! wait some time
    call sleep(1)

    now = new_datetime()

    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 2, .true.)
    print*, ""

    ! add some time
    print*, "Add one month..."
    call now%add(month=1)
    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 3, .true.)
    print*, ""

    print*, "Substract 17 months and add 100 days..."
    call now%add(month=-17, day=100)
    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 4, .true.)
    print*, ""

    print*, "Substract 60 minutes and add 1 year..."
    call now%add(year=1, minute=-60)
    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 5, .true.)
    print*, ""

    ! create a new date with a given value
    now = new_datetime(hour=12,second=56.0_8)
    print*, "New datetime initialized with hour=12 and second=56..."
    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 6, .true.)
    print*, ""

    ! change some parts of the date
    call now%set(year=2013, month=5)
    print*, "Set year to 2013 and month to 5"
    print*, "Date and time: ", now%to_string()
    print*, "hashcode:      ", now%hashcode()
    call dates%add(now, 7, .true.)
    print*, ""

    ! get some parts of the date
    call now%get(month=test1, hour=test2)
    print*, "Month:", test1, "Hour:", test2, "Year:", now%get_field(DT_YEAR)
    print*, ""

    call dates%printContent()
end program test_datetime


text:
      SUBROUTINE FTPAR(PAR,FLAG)

************************************************************************
*   Subroutine extracting the fine-tuning
*
*   At the SUSY scale FTSUSY(i)=d(ln MZ^2)/d(ln PSi^2)
*
*     FTSUSY(1) --> PS=m_Hu
*     FTSUSY(2) --> PS=m_Hd
*     FTSUSY(3) --> PS=m_S
*     FTSUSY(4) --> PS=A_lambda
*     FTSUSY(5) --> PS=A_kappa
*     FTSUSY(6) --> PS=XiF
*     FTSUSY(7) --> PS=XiS
*     FTSUSY(8) --> PS=MUP
*     FTSUSY(9) --> PS=MSP
*     FTSUSY(10) --> PS=M_3H
*     FTSUSY(11) --> PS=lambda
*     FTSUSY(12) --> PS=kappa
*     FTSUSY(13) --> PS=htop
*     FTSUSY(14) --> PS=g
*     FTSUSY(15) --> MAX(FTSUSY(i))
*     FTSUSY(16) --> i FOR FTMAX
*
*   Based on the 3 minimization equations:
*
*     (1)=mhu+mu**2+mz**2*(lh+(tb**2-1)/2)/(1+tb**2)
*        -(mu*b+m3+mu*mup+xfh)/tb-ct*mz**2*tb**2/(1+tb**2)
*     (2)=mhd+mu**2+mz**2*(lh*tb**2+(1-tb**2)/2)/(1+tb**2)
*        -(mu*b+m3+mu*mup+xfh)*tb
*     (3)=ms+msp+mup**2+2*kh*xfh+kh*mu*(ak+2*kh*mu+3*mup)
*        +mz**2*lh*(1-(2*kh*mu+al+mup)/mu*tb/(1+tb**2))
*        +(xsh+xfh*mup)/mu
*
*     where lh=l**2/g**2, kh=k/l, xfh=l*xif, xsh = l*xis, b=al+kh*mu
*     ct0 = (3*ht**4)/(8*pi**2*g**2)*ln(mt**2/QSTSB)
*
*   At the GUT scale FTGUT(i)=d(ln MZ^2)/d(ln PGi^2)
*
*     FTGUT(1) --> PG=m_Hu
*     FTGUT(2) --> PG=m_Hd
*     FTGUT(3) --> PG=m_S
*     FTGUT(4) --> PG=m0
*     FTGUT(5) --> PG=M1
*     FTGUT(6) --> PG=M2
*     FTGUT(7) --> PG=M3
*     FTGUT(8) --> PG=M12
*     FTGUT(9) --> PG=A_lambda
*     FTGUT(10) --> PG=A_kappa
*     FTGUT(11) --> PG=A0
*     FTGUT(12) --> PG=XiF
*     FTGUT(13) --> PG=XiS
*     FTGUT(14) --> PG=MUP
*     FTGUT(15) --> PG=MSP
*     FTGUT(16) --> PG=M3H
*     FTGUT(17) --> PG=lambda
*     FTGUT(18) --> PG=kappa
*     FTGUT(19) --> PG=htop
*     FTGUT(20) --> PG=g
*     FTGUT(21) --> PG=MGUT
*     FTGUT(22) --> MAX(FTGUT(i))
*     FTGUT(23) --> i FOR FTMAX
*
*   JACG(j,i) is the Jacobian defined as d(ln PSj^2)/d(ln PGi^2)
*
*   At the MESS scale FTMES(i)=d(ln MZ^2)/d(ln PMi^2)
*
*     FTMES(1) --> PM=M_susyeff
*     FTMES(2) --> PM=m_S
*     FTMES(3) --> PM=A_lambda
*     FTMES(4) --> PM=XiF
*     FTMES(5) --> PM=XiS
*     FTMES(6) --> PM=MUP
*     FTMES(7) --> PM=MSP
*     FTMES(8) --> PM=lambda
*     FTMES(9) --> PM=kappa
*     FTMES(10) --> PM=htop
*     FTMES(11) --> PM=g1
*     FTMES(12) --> PM=g2
*     FTMES(13) --> PM=MMESS
*     FTMES(14) --> PM=lambda_PP
*     FTMES(15) --> PM=lambda_TT
*     FTMES(16) --> PM=lambda_SU
*     FTMES(17) --> PM=lambda_SD
*     FTMES(18) --> PM=lambda_T
*     FTMES(19) --> PM=lambda_B
*     FTMES(20) --> PM=lambda_L
*     FTMES(21) --> MAX(FTGUT(i))
*     FTMES(22) --> i FOR FTMAX
*
*   JACM(j,i) is the Jacobian defined as d(ln PSj^2)/d(ln PMi^2)
*
*
************************************************************************

      IMPLICIT NONE

      INTEGER I,J,FLAG,M1FLAG,M2FLAG,M3FLAG,MHDFLAG,MHUFLAG
      INTEGER MSFLAG,AKFLAG,ALFLAG,NSUSY,NGUT,NMES,GMFLAG
      PARAMETER (NSUSY=14,NGUT=21,NMES=20)

      DOUBLE PRECISION PAR(*),FTSUSY(NSUSY+2),FTGUT(NGUT+2)
      DOUBLE PRECISION PS(NSUSY),PG(NGUT),JACG(NSUSY,NGUT)
      DOUBLE PRECISION PM(NMES),JACM(NSUSY,NMES),FTMES(NMES+2)
      DOUBLE PRECISION L,K,TB,TB2,MU,AL,AK,PI,ALS,B
      DOUBLE PRECISION LH,KH,XFH,XSH,CT
      DOUBLE PRECISION A1,A2,A3,D1Z,D2Z,D3Z
      DOUBLE PRECISION D1T,D2T,D3T,D1M,D2M,D3M,D
      DOUBLE PRECISION MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW
      DOUBLE PRECISION G,G1,G2,G3,HT,HB,HL
      DOUBLE PRECISION XIF,XIS,MUP,MSP,M3H
      DOUBLE PRECISION MHU2,MHD2,MS2,M0,M12,A0
      DOUBLE PRECISION QSTSB,HTQ,HBQ,MTOPQ,MBOTQ
      DOUBLE PRECISION G1GUT,G2GUT,G3GUT,LGUT,KGUT
      DOUBLE PRECISION HTGUT,HBGUT,HLGUT,MGUT
      DOUBLE PRECISION M1GUT,M2GUT,M3GUT,ALGUT,AKGUT,ATGUT,ABGUT
      DOUBLE PRECISION ATAUGUT,AMUGUT,MHUGUT,MHDGUT,MSGUT,MQ3GUT
      DOUBLE PRECISION MU3GUT,MD3GUT,MQGUT,MUGUT,MDGUT,ML3GUT
      DOUBLE PRECISION ME3GUT,MLGUT,MEGUT
      DOUBLE PRECISION XIFGUT,XISGUT,MUPGUT,MSPGUT,M3HGUT
      DOUBLE PRECISION MSUSYEFF,MMESS,N5
      DOUBLE PRECISION G1MES,G2MES,G3MES,LMES,KMES,HTMES,HBMES,HLMES
      DOUBLE PRECISION XIFMES,XISMES,MSMES,MUPMES,MSPMES,M3HMES
      DOUBLE PRECISION LPPMES,LTTMES,LUMES,LDMES,LTMES,LBMES,LLMES,XIU
      DOUBLE PRECISION M1MES,M2MES,M3MES,ALMES,AKMES,ATMES,ABMES,
     .      ATAUMES,AMUMES,MH1MES,MH2MES,MQ3MES,MU3MES,MD3MES,
     .      MQMES,MUMES,MDMES,ML3MES,ME3MES,MLMES,MEMES
      DOUBLE PRECISION MSREF,DD,DMIN

      COMMON/SMSPEC/MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW
      COMMON/SUSYCOUP/g1,g2,g3,HT,HB,HL
      COMMON/SUSYEXT/XIF,XIS,MUP,MSP,M3H
      COMMON/SUSYMH/MHU2,MHD2,MS2
      COMMON/SOFTGUT/M0,M12,A0
      COMMON/GUTCOUP/G1GUT,G2GUT,G3GUT,LGUT,KGUT,HTGUT,
     .      HBGUT,HLGUT
      COMMON/GUTPAR/M1GUT,M2GUT,M3GUT,ALGUT,AKGUT,ATGUT,ABGUT,
     .      ATAUGUT,AMUGUT,MHUGUT,MHDGUT,MSGUT,MQ3GUT,MU3GUT,MD3GUT,
     .      MQGUT,MUGUT,MDGUT,ML3GUT,ME3GUT,MLGUT,MEGUT
      COMMON/GUTEXT/XIFGUT,XISGUT,MUPGUT,MSPGUT,M3HGUT
      COMMON/SCANFLAGS/M1FLAG,M2FLAG,M3FLAG,MHDFLAG,MHUFLAG,
     .      MSFLAG,AKFLAG,ALFLAG
      COMMON/QQUARK/HTQ,HBQ,MTOPQ,MBOTQ
      COMMON/ALSHIFT/ALS,B
      COMMON/STSBSCALE/QSTSB
      COMMON/MGUT/MGUT
      COMMON/MESCAL/MSUSYEFF,MMESS,N5
      COMMON/MESCOUP/G1MES,G2MES,G3MES,LMES,KMES,HTMES,HBMES,HLMES
      COMMON/MESEXT/XIFMES,XISMES,MSMES,MUPMES,MSPMES,M3HMES
      COMMON/MESGUT/LPPMES,LTTMES,LUMES,LDMES,LTMES,LBMES,LLMES,XIU
      COMMON/SOFTMES/M1MES,M2MES,M3MES,ALMES,AKMES,ATMES,ABMES,
     .      ATAUMES,AMUMES,MH1MES,MH2MES,MQ3MES,MU3MES,MD3MES,
     .      MQMES,MUMES,MDMES,ML3MES,ME3MES,MLMES,MEMES
      COMMON/FINETUN/FTSUSY,FTGUT,FTMES
      COMMON/GMSCEN/MSREF,DD,DMIN,GMFLAG

      PI=4d0*DATAN(1d0)

      DO I=1,NSUSY+2
       FTSUSY(I)=0d0
      ENDDO

      L=PAR(1)
      K=PAR(2)
      TB=PAR(3)
      MU=PAR(4)
      AL=PAR(5)
      AK=PAR(6)
      G=(G1+G2)/2d0

      LH=L**2/G
      KH=K/L
      XFH=L*XIF
      XSH=L*XIS
      TB2=TB**2
      CT=3d0*HT**4/(8d0*PI**2*G)*DLOG(MTOPQ**2/QSTSB)

*   DiT = d(i)/d(tb)

      D1T = 2d0*MZ**2*TB/(1d0+TB2)**2*(1d0-LH-CT)
     .    + (MU*B+M3H+MU*MUP+XFH)/TB2
      D2T = 2d0*MZ**2*TB/(1d0+TB2)**2*(LH-1d0)
     .    - (MU*B+M3H+MU*MUP+XFH)
      D3T = LH*MZ**2*(TB2-1d0)/(1d0+TB2)**2
     .     *(B+KH*MU+MUP)/MU

*   DiM = d(i)/d(mu)

      D1M = 2d0*MU-(B+KH*MU+MUP)/TB
      D2M = 2d0*MU-(B+KH*MU+MUP)*TB
      D3M = KH*AK+4d0*KH**2*MU+3d0*KH*MUP
     .    + LH*MZ**2*TB/(1d0+TB2)*(AL+MUP)/MU**2
     .    - (XSH+XFH*MUP)/MU**2

*   DiZ = d(i)/d(MZ^2)

      D1Z = (LH+(TB2-1d0)/2d0-CT*TB2)/(1d0+TB2)
      D2Z = (LH*TB2+(1d0-TB2)/2d0)/(1d0+TB2)
      D3Z = LH*(1d0-TB/(1d0+TB2)*(B+KH*MU+MUP)/MU)

*   Ai=d(j)/d(tb)*d(k)/d(mu)-d(k)/d(tb)*d(j)/d(mu)

      A1=D2T*D3M-D2M*D3T
      A2=D3T*D1M-D3M*D1T
      A3=D1T*D2M-D1M*D2T

*   Determinant of the linear system

      D=A1*D1Z+A2*D2Z+A3*D3Z

*   Fine-tuning measures at the SUSY scale

      FTSUSY(1)=MHU2/MZ**2*A1/D
      FTSUSY(2)=MHD2/MZ**2*A2/D
      FTSUSY(3)=MS2/MZ**2*A3/D
      FTSUSY(4)=AL/2d0/MZ**2*(-MU/TB*A1-MU*TB*A2
     .          -MZ**2*LH*TB/(1d0+TB2)/MU*A3)/D
      FTSUSY(5)=AK/2d0/MZ**2*KH*MU*A3/D
      FTSUSY(6)=XFH/MZ**2*(-A1/TB-A2*TB+(2d0*KH+MUP/MU)*A3)/D
      FTSUSY(7)=3d0/2d0*XSH/MZ**2*A3/MU/D
      FTSUSY(8)=MUP/2d0/MZ**2*(-MU/TB*A1-MU*TB*A2+(2d0*MUP
     .          +3d0*KH*MU+XFH/MU-LH*MZ**2/MU*TB/(1d0+TB2))*A3)/D
      FTSUSY(9)=MSP/MZ**2*A3/D
      FTSUSY(10)=M3H/MZ**2*(-A1/TB-A2*TB)/D
      FTSUSY(11)=(LH*(A1+A2*TB2)/(1d0+TB2)+D3Z*A3)/D
      FTSUSY(12)=KH/2d0/MZ**2*(-MU**2/TB*A1-MU**2*TB*A2
     .          +(2d0*XFH+AK*MU+4d0*KH*MU**2+3d0*MU*MUP
     .          -2d0*LH*MZ**2*TB/(1d0+TB2))*A3)/D
      FTSUSY(13)=-2d0*CT*TB2/(1d0+TB2)*A1/D
      FTSUSY(14)=-FTSUSY(11)+CT*TB2/(1d0+TB2)*A1/D
      FTSUSY(11)=FTSUSY(11)-FTSUSY(12)

      FTSUSY(NSUSY+1)=DABS(FTSUSY(1))
      FTSUSY(NSUSY+2)=1d0
      DO I=2,NSUSY
       IF(DABS(FTSUSY(I)).GE.FTSUSY(NSUSY+1))THEN
        FTSUSY(NSUSY+1)=DABS(FTSUSY(I))
        FTSUSY(NSUSY+2)=DFLOAT(I)
       ENDIF
      ENDDO

*   mSUGRA boundary conditions

      IF(FLAG.EQ.1)THEN

       DO I=1,NGUT+2
        FTGUT(I)=0d0
       ENDDO

*   GUT scale parameters

       IF(MHUFLAG.EQ.1)THEN
        PG(1)=MHUGUT
       ELSE
        PG(1)=0d0
       ENDIF
       IF(MHDFLAG.EQ.1)THEN
        PG(2)=MHDGUT
       ELSE
        PG(2)=0d0
       ENDIF
       IF(MSFLAG.EQ.1)THEN
        PG(3)=MSGUT
       ELSE
        PG(3)=0d0
       ENDIF
       PG(4)=M0**2
       IF(M1FLAG.EQ.1)THEN
        PG(5)=M1GUT
       ELSE
        PG(5)=0d0
       ENDIF
       IF(M2FLAG.EQ.1)THEN
        PG(6)=M2GUT
       ELSE
        PG(6)=0d0
       ENDIF
       IF(M3FLAG.EQ.1)THEN
        PG(7)=M3GUT
       ELSE
        PG(7)=0d0
       ENDIF
       IF(M1FLAG*M2FLAG*M3FLAG.EQ.0)THEN
        PG(8)=M12
       ELSE
        PG(8)=0d0
       ENDIF
       IF(ALFLAG.EQ.1)THEN
        PG(9)=ALGUT
       ELSE
        PG(9)=0d0
       ENDIF
       IF(AKFLAG.EQ.1)THEN
        PG(10)=AKGUT
       ELSE
        PG(10)=0d0
       ENDIF
       PG(11)=A0
       PG(12)=XIFGUT
       PG(13)=XISGUT
       PG(14)=MUPGUT
       PG(15)=MSPGUT
       PG(16)=M3HGUT
       PG(17)=LGUT
       PG(18)=KGUT
       PG(19)=HTGUT
       PG(20)=G2GUT
       PG(21)=MGUT**2

       CALL RGESVAR(NGUT,NSUSY,PG,PS)

       IF(MHU2.NE.0d0)
     .  FTSUSY(1)=PS(1)/MHU2*FTSUSY(1)
       IF(MHD2.NE.0d0)
     .  FTSUSY(2)=PS(2)/MHD2*FTSUSY(2)
       IF(MS2.NE.0d0)
     .  FTSUSY(3)=PS(3)/MS2*FTSUSY(3)
       IF(AL.NE.0d0)
     .  FTSUSY(4)=PS(4)/AL*FTSUSY(4)
       IF(AK.NE.0d0)
     .  FTSUSY(5)=PS(5)/AK*FTSUSY(5)
       IF(XIF.NE.0d0)
     .  FTSUSY(6)=PS(6)/XIF*FTSUSY(6)
       IF(XIS.NE.0d0)
     .  FTSUSY(7)=PS(7)/XIS*FTSUSY(7)
       IF(MUP.NE.0d0)
     .  FTSUSY(8)=PS(8)/MUP*FTSUSY(8)
       IF(MSP.NE.0d0)
     .  FTSUSY(9)=PS(9)/MSP*FTSUSY(9)
       IF(M3H.NE.0d0)
     .  FTSUSY(10)=PS(10)/M3H*FTSUSY(10)
       FTSUSY(11)=PS(11)/L**2*FTSUSY(11)
       IF(K.NE.0d0)
     .  FTSUSY(12)=PS(12)/K*FTSUSY(12)
       FTSUSY(13)=PS(13)/HT**2*FTSUSY(13)
       FTSUSY(14)=PS(14)/G*FTSUSY(14)

       CALL JACOBIAN(NGUT,NSUSY,PG,PS,JACG)

       DO I=1,NGUT
        FTGUT(I)=0d0
        DO J=1,NSUSY
         FTGUT(I)=FTGUT(I)+FTSUSY(J)*JACG(J,I)
        ENDDO
       ENDDO

       FTGUT(NGUT+1)=DABS(FTGUT(1))
       FTGUT(NGUT+2)=1d0
       DO I=2,NGUT-2
        IF(DABS(FTGUT(I)).GE.FTGUT(NGUT+1))THEN
         FTGUT(NGUT+1)=DABS(FTGUT(I))
         FTGUT(NGUT+2)=DFLOAT(I)
        ENDIF
       ENDDO

      ENDIF

*   GMSB boundary conditions

      IF(FLAG.EQ.2)THEN

       DO I=1,NMES+2
        FTMES(I)=0d0
       ENDDO

*   MES scale parameters

       PM(1)=MSUSYEFF
       IF(GMFLAG.EQ.0)THEN
        PM(2)=MSMES
       ELSE
        PM(2)=0d0
       ENDIF
       IF(GMFLAG.EQ.0)THEN
        PM(3)=ALMES
       ELSE
        PM(3)=0d0
       ENDIF
       PM(4)=XIFMES
       PM(5)=XISMES
       PM(6)=MUPMES
       PM(7)=MSPMES
       PM(8)=LMES
       PM(9)=KMES
       PM(10)=HTMES
       PM(11)=G1MES
       PM(12)=G2MES
       PM(13)=MMESS
       PM(14)=LPPMES
       PM(15)=LTTMES
       PM(16)=LUMES
       PM(17)=LDMES
       PM(18)=LTMES
       PM(19)=LBMES
       PM(20)=LLMES

       CALL RGESVARGM(NMES,NSUSY,PM,PS,N5)

       IF(MHU2.NE.0d0)
     .  FTSUSY(1)=PS(1)/MHU2*FTSUSY(1)
       IF(MHD2.NE.0d0)
     .  FTSUSY(2)=PS(2)/MHD2*FTSUSY(2)
       IF(MS2.NE.0d0)
     .  FTSUSY(3)=PS(3)/MS2*FTSUSY(3)
       IF(AL.NE.0d0)
     .  FTSUSY(4)=PS(4)/AL*FTSUSY(4)
       IF(AK.NE.0d0)
     .  FTSUSY(5)=PS(5)/AK*FTSUSY(5)
       IF(XIF.NE.0d0)
     .  FTSUSY(6)=PS(6)/XIF*FTSUSY(6)
       IF(XIS.NE.0d0)
     .  FTSUSY(7)=PS(7)/XIS*FTSUSY(7)
       IF(MUP.NE.0d0)
     .  FTSUSY(8)=PS(8)/MUP*FTSUSY(8)
       IF(MSP.NE.0d0)
     .  FTSUSY(9)=PS(9)/MSP*FTSUSY(9)
       IF(M3H.NE.0d0)
     .  FTSUSY(10)=PS(10)/M3H*FTSUSY(10)
       FTSUSY(11)=PS(11)/L**2*FTSUSY(11)
       IF(K.NE.0d0)
     .  FTSUSY(12)=PS(12)/K*FTSUSY(12)
       FTSUSY(13)=PS(13)/HT**2*FTSUSY(13)
       FTSUSY(14)=PS(14)/G*FTSUSY(14)

       CALL JACOBIANGM(NMES,NSUSY,PM,PS,N5,JACM)

       DO I=1,NMES
        FTMES(I)=0d0
        DO J=1,NSUSY
         FTMES(I)=FTMES(I)+FTSUSY(J)*JACM(J,I)
        ENDDO
       ENDDO

       FTMES(NMES+1)=DABS(FTMES(1))
       FTMES(NMES+2)=1d0
       DO I=2,NMES
        IF(DABS(FTMES(I)).GE.FTMES(NMES+1))THEN
         FTMES(NMES+1)=DABS(FTMES(I))
         FTMES(NMES+2)=DFLOAT(I)
        ENDIF
       ENDDO

      ENDIF

      END

      SUBROUTINE JACOBIAN(M,N,PG,PS,JACG)

      IMPLICIT NONE

      INTEGER N,M,I,J
      DOUBLE PRECISION JACG(N,M),PG(M),PS(N),PP(N),T,H,EPS

      EPS=1d-4

      DO J=1,M
       T=PG(J)
       IF(T.EQ.0d0)THEN
        DO I=1,N
         JACG(I,J)=0d0
        ENDDO
       ELSE
        H=EPS*DABS(T)
        PG(J)=T+H
        CALL RGESVAR(M,N,PG,PP)
        PG(J)=T
        DO I=1,N
         IF(PS(I).EQ.0d0)THEN
          JACG(I,J)=0d0
         ELSE
          JACG(I,J)=PG(J)/PS(I)*PP(I)/(2d0*H)
         ENDIF
        ENDDO
        PG(J)=T-H
        CALL RGESVAR(M,N,PG,PP)
        PG(J)=T
        DO I=1,N
         IF(PS(I).EQ.0d0)THEN
          JACG(I,J)=0d0
         ELSE
          JACG(I,J)=JACG(I,J)-PG(J)/PS(I)*PP(I)/(2d0*H)
         ENDIF
        ENDDO
       ENDIF
      ENDDO

      DO J=1,M
       JACG(4,J)=JACG(4,J)*2d0
       JACG(5,J)=JACG(5,J)*2d0
       JACG(7,J)=JACG(7,J)*2d0/3d0
       JACG(8,J)=JACG(8,J)*2d0
       JACG(12,J)=JACG(12,J)*2d0
      ENDDO
      DO I=1,N
       DO J=5,11
        JACG(I,J)=JACG(I,J)/2d0
       ENDDO
       JACG(I,13)=JACG(I,14)*3d0/2d0
       JACG(I,14)=JACG(I,14)/2d0
       JACG(I,18)=JACG(I,18)/2d0
      ENDDO

      END

      SUBROUTINE RGESVAR(M,N,PG,PS)

      IMPLICIT NONE

      INTEGER I,IFAIL,M,N,M1,M2,M3,MHD,MHU,MS,AK,AL,NN
      PARAMETER (NN=35)

      DOUBLE PRECISION PG(M),PS(N),EPS,X1,X2,Y(NN),PI
      DOUBLE PRECISION G1,G2,G3,L,K,HT,HB,HL,Q2,COEF

      COMMON/GUTCOUP/G1,G2,G3,L,K,HT,HB,HL
      COMMON/RENSCALE/Q2
      COMMON/SCANFLAGS/M1,M2,M3,MHD,MHU,MS,AK,AL

      EXTERNAL DERIVSS,RKQSS

      EPS=1d-8
      PI=4d0*DATAN(1d0)
      COEF=1d0/(16d0*PI**2)

      DO I=1,N
       PS(I)=0d0
      ENDDO

      Y(1)=3d0/5d0*PG(20)
      Y(2)=PG(20)
      Y(3)=G3-G2+PG(20)
      Y(4)=PG(17)
      Y(5)=PG(18)
      Y(6)=PG(19)
      Y(7)=HB
      Y(8)=HL
      IF(M1.EQ.1)THEN
       Y(9)=PG(5)
      ELSE
       Y(9)=PG(8)
      ENDIF
      IF(M2.EQ.1)THEN
       Y(10)=PG(6)
      ELSE
       Y(10)=PG(8)
      ENDIF
      IF(M3.EQ.1)THEN
       Y(11)=PG(7)
      ELSE
       Y(11)=PG(8)
      ENDIF
      IF(AL.EQ.1)THEN
       Y(12)=PG(9)
      ELSE
       Y(12)=PG(11)
      ENDIF
      IF(AK.EQ.1)THEN
       Y(13)=PG(10)
      ELSE
       Y(13)=PG(11)
      ENDIF
      Y(14)=PG(11)
      Y(15)=PG(11)
      Y(16)=PG(11)
      IF(MHU.EQ.1)THEN
       Y(17)=PG(1)
      ELSE
       Y(17)=PG(4)
      ENDIF
      IF(MHD.EQ.1)THEN
       Y(18)=PG(2)
      ELSE
       Y(18)=PG(4)
      ENDIF
      IF(MS.EQ.1)THEN
       Y(19)=PG(3)
      ELSE
       Y(19)=PG(4)
      ENDIF
      DO I=20,29
        Y(I)=PG(4)
      ENDDO
      Y(30)=PG(12)
      Y(31)=PG(13)
      Y(32)=PG(14)
      Y(33)=PG(15)
      Y(34)=PG(16)
      Y(35)=PG(11)

      X1=COEF*DLOG(PG(21)/Q2)
      X2=0d0

      CALL ODEINTS(Y,NN,X1,X2,EPS,DERIVSS,RKQSS,IFAIL)

      IF(IFAIL.NE.0)THEN
       RETURN
      ENDIF

      PS(1)=Y(17)
      PS(2)=Y(18)
      PS(3)=Y(19)
      PS(4)=Y(12)
      PS(5)=Y(13)
      PS(6)=Y(30)
      PS(7)=Y(31)
      PS(8)=Y(32)
      PS(9)=Y(33)
      PS(10)=Y(34)
      PS(11)=Y(4)
      PS(12)=Y(5)
      PS(13)=Y(6)
      PS(14)=(Y(1)+Y(2))/2d0

      END

      SUBROUTINE JACOBIANGM(M,N,PM,PS,N5,JACM)

      IMPLICIT NONE

      INTEGER N,M,I,J
      DOUBLE PRECISION JACM(N,M),PM(M),PS(N),PP(N),N5,T,H,EPS

      EPS=1d-4

      DO J=1,M
       T=PM(J)
       IF(T.EQ.0d0)THEN
        DO I=1,N
         JACM(I,J)=0d0
        ENDDO
       ELSE
        H=EPS*DABS(T)
        PM(J)=T+H
        CALL RGESVARGM(M,N,PM,PP,N5)
        PM(J)=T
        DO I=1,N
         IF(PS(I).EQ.0d0)THEN
          JACM(I,J)=0d0
         ELSE
          JACM(I,J)=PM(J)/PS(I)*PP(I)/(2d0*H)
         ENDIF
        ENDDO
        PM(J)=T-H
        CALL RGESVARGM(M,N,PM,PP,N5)
        PM(J)=T
        DO I=1,N
         IF(PS(I).EQ.0d0)THEN
          JACM(I,J)=0d0
         ELSE
          JACM(I,J)=JACM(I,J)-PM(J)/PS(I)*PP(I)/(2d0*H)
         ENDIF
        ENDDO
       ENDIF
      ENDDO

      DO J=1,M
       JACM(4,J)=JACM(4,J)*2d0
       JACM(5,J)=JACM(5,J)*2d0
       JACM(7,J)=JACM(7,J)*2d0/3d0
       JACM(8,J)=JACM(8,J)*2d0
       JACM(12,J)=JACM(12,J)*2d0
      ENDDO
      DO I=1,N
       JACM(I,1)=JACM(I,1)/2d0
       JACM(I,3)=JACM(I,3)/2d0
       JACM(I,5)=JACM(I,5)*3d0/2d0
       JACM(I,6)=JACM(I,6)/2d0
       JACM(I,9)=JACM(I,9)/2d0
       DO J=13,20
        JACM(I,J)=JACM(I,J)/2d0
       ENDDO
      ENDDO

      END

      SUBROUTINE RGESVARGM(M,N,PM,PS,N5)

      IMPLICIT NONE

      INTEGER GMFLAG,I,IFAIL,M,N,NN
      PARAMETER (NN=35)

      DOUBLE PRECISION PM(M),PS(N),N5,EPS,X1,X2,Y(NN)
      DOUBLE PRECISION G1,G2,G3,L,K,HT,HB,HL,Q2,COEF
      DOUBLE PRECISION X,F1,F2,ALP1,ALP2,ALP3,SP2,PI
      DOUBLE PRECISION MSREF,D,DMIN

      COMMON/MESCOUP/G1,G2,G3,L,K,HT,HB,HL
      COMMON/RENSCALE/Q2
      COMMON/GMSCEN/MSREF,D,DMIN,GMFLAG

      EXTERNAL DERIVSS,RKQSS

      EPS=1d-8
      PI=4d0*DATAN(1d0)
      COEF=1d0/(16d0*PI**2)

      DO I=1,N
       PS(I)=0d0
      ENDDO

      Y(1)=PM(11)
      Y(2)=PM(12)
      Y(3)=G3
      Y(4)=PM(8)
      Y(5)=PM(9)
      Y(6)=PM(10)
      Y(7)=HB
      Y(8)=HL

      X=PM(1)/PM(13)
      IF(X.GE.1d-7) THEN
        IF(X.LT.1d0) THEN
          F1=((1d0+X)*DLOG(1d0+X)+(1d0-X)*DLOG(1d0-X))/X**2
          F2=(1d0+X)/X**2*(DLOG(1d0+X)-2d0*SP2(X/(1d0+X))
     .       +SP2(2d0*X/(1d0+X))/2d0)
     .       +(1d0-X)/X**2*(DLOG(1d0-X)-2d0*SP2(-X/(1d0-X))
     .       +SP2(-2d0*X/(1d0-X))/2d0)
        ELSE
          F1=0d0
          F2=0d0
        ENDIF
      ELSE
        F1=1d0
        F2=1d0
      ENDIF

      ALP1=PM(11)/(4d0*PI)
      ALP2=PM(12)/(4d0*PI)
      ALP3=G3/(4d0*PI)

      Y(9)=5d0/3d0*ALP1*N5*PM(1)*F1/(4d0*PI)
      Y(10)=ALP2*N5*PM(1)*F1/(4d0*PI)
      Y(11)=ALP3*N5*PM(1)*F1/(4d0*PI)
      IF(GMFLAG.NE.0)THEN
       PM(3)=-COEF*PM(1)*
     .       (2d0*PM(14)**2+3d0*PM(15)**2+3d0*PM(16)**2+3d0*PM(17)**2)
      ENDIF
      Y(12)=PM(3)
      IF(PM(9).NE.0d0)THEN
       IF(GMFLAG.NE.0)THEN
        Y(13)=-3d0*COEF*PM(1)*
     .        (2d0*PM(14)**2+3d0*PM(15)**2+2d0*PM(16)**2+2d0*PM(17)**2)
       ELSE
        Y(13)=3d0*PM(3)
       ENDIF
      ELSE
       Y(13)=0d0
      ENDIF
      Y(14)=-COEF*PM(1)*(PM(16)**2+3d0*PM(18)**2+PM(19)**2)
      Y(15)=-COEF*PM(1)*(PM(17)**2+PM(18)**2+3d0*PM(19)**2)
      Y(16)=-COEF*PM(1)*(PM(17)**2+3d0*PM(20)**2)
      Y(17)=COEF*(5d0/6d0*ALP1**2+3d0/2d0*ALP2**2)*N5*F2*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 3d0*PM(10)*PM(19)**2 - 3d0*PM(17)**2*PM(8)
     .     - 2d0*PM(8)*PM(14)**2 - 9d0*PM(10)*PM(18)**2
     .     - 3d0*PM(8)*PM(15)**2 + 3d0*PM(14)**2*PM(16)**2
     .     + 3d0*PM(15)**2*PM(16)**2
     .     + 6d0*DSQRT(HB*PM(8))*PM(19)*PM(16)
     .     + 2d0*DSQRT(HL*PM(8))*PM(20)*PM(16)
     .     + 2d0*PM(17)*DSQRT(PM(8))*PM(14)*PM(16)
     .     - PM(11)*PM(16)**2 - 3d0*PM(12)*PM(16)**2
     .     + 2d0*PM(9)**2*PM(16)**2 + 3d0*PM(19)**2*PM(16)**2
     .     + 2d0*PM(17)**2*PM(16)**2
     .     + PM(20)**2*PM(16)**2 + 2d0*PM(8)*PM(16)**2 + 4d0*PM(16)**4)
     .     - PM(1)**4/PM(13)**2*COEF/6d0*PM(16)**2
      Y(18)=COEF*(5d0/6d0*ALP1**2+3d0/2d0*ALP2**2)*N5*F2*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 9d0*HB*PM(19)**2 - PM(11)*PM(17)**2 - 3d0*PM(12)*PM(17)**2
     .     + 2d0*PM(9)**2*PM(17)**2 + 4d0*PM(17)**4
     .     - 3d0*HL*PM(20)**2 + 2d0*PM(17)**2*PM(8)
     .     + 3d0*PM(17)**2*PM(14)**2 - 2d0*PM(8)*PM(14)**2
     .     + 6d0*DSQRT(PM(10)*PM(8))*PM(17)*PM(18)
     .     - 3d0*HB*PM(18)**2 + 3d0*PM(17)**2*PM(18)**2
     .     + 3d0*PM(17)**2*PM(15)**2 - 3d0*PM(8)*PM(15)**2
     .     + 2d0*PM(17)*DSQRT(PM(8))*PM(14)*PM(16)
     .     + 2d0*PM(17)**2*PM(16)**2 - 3d0*PM(8)*PM(16)**2)
     .     - PM(1)**4/PM(13)**2*COEF/6d0*PM(17)**2
      IF(GMFLAG.NE.0)THEN
       Y(19)=PM(1)**2/(4d0*PI)**4*(
     .      - 2d0*PM(11)*PM(17)**2 - 6d0*PM(12)*PM(17)**2
     .      + 6d0*HB*PM(17)**2 + 2d0*HL*PM(17)**2
     .      - 8d0*PM(9)**2*PM(17)**2 + 8d0*PM(17)**4
     .      + 4d0*PM(17)**2*PM(8)
     .      + 12d0*DSQRT(HB)*PM(19)*PM(17)*PM(14)
     .      + 4d0*DSQRT(HL)*PM(17)*PM(20)*PM(14)
     .      - 2d0*PM(11)*PM(14)**2 - 6d0*PM(12)*PM(14)**2
     .      - 8d0*PM(9)**2*PM(14)**2 + 6d0*PM(19)**2*PM(14)**2
     .      + 16d0*PM(17)**2*PM(14)**2 + 2d0*PM(20)**2*PM(14)**2
     .      + 8d0*PM(14)**4 + 12d0*DSQRT(PM(10)*PM(8))*PM(17)*PM(18)
     .      + 6d0*PM(17)**2*PM(18)**2 + 6d0*PM(14)**2*PM(18)**2
     .      - 4d0/3d0*PM(11)*PM(15)**2 - 16d0*G3*PM(15)**2
     .      - 12d0*PM(9)**2*PM(15)**2 + 12d0*PM(17)**2*PM(15)**2
     .      + 12d0*PM(14)**2*PM(15)**2 + 15d0*PM(15)**4
     .      + 12d0*DSQRT(HB*PM(8))*PM(19)*PM(16)
     .      + 4d0*DSQRT(HL*PM(8))*PM(20)*PM(16)
     .      + 16d0*PM(17)*DSQRT(PM(8))*PM(14)*PM(16)
     .      + 12d0*DSQRT(PM(10))*PM(14)*PM(18)*PM(16)
     .      - 2d0*PM(11)*PM(16)**2 - 6d0*PM(12)*PM(16)**2
     .      + 6d0*PM(10)*PM(16)**2 - 8d0*PM(9)**2*PM(16)**2
     .      + 6d0*PM(19)**2*PM(16)**2 + 8d0*PM(17)**2*PM(16)**2
     .      + 2d0*PM(20)**2*PM(16)**2 + 4d0*PM(8)*PM(16)**2
     .      + 16d0*PM(14)**2*PM(16)**2 + 12d0*PM(15)**2*PM(16)**2
     .      + 8d0*PM(16)**4)
     .      -PM(1)**4/(48d0*PI**2*PM(13)**2)
     .       *(2d0*PM(14)**2+3d0*PM(15)**2+PM(16)**2+PM(17)**2)
      ELSE
       Y(19)=PM(2)
      ENDIF
      Y(20)=COEF*(5d0/54d0*ALP1**2+3d0/2d0*ALP2**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 7d0/9d0*PM(11)*PM(19)**2 - 3d0*PM(12)*PM(19)**2
     .     - 16d0/3d0*G3*PM(19)**2 + 6d0*HB*PM(19)**2
     .     + 6d0*PM(19)**4 - HB*PM(17)**2 + PM(19)**2*PM(20)**2
     .     + 2d0*DSQRT(HB*HL)*PM(19)*PM(20)
     .     + 2d0*DSQRT(HB)*PM(19)*PM(17)*PM(14)
     .     + PM(19)**2*PM(14)**2 + 2d0*DSQRT(PM(10)*PM(8))*PM(17)*PM(18)
     .     - 13d0/9d0*PM(11)*PM(18)**2 - 3d0*PM(12)*PM(18)**2
     .     - 16d0/3d0*G3*PM(18)**2 + 6d0*PM(10)*PM(18)**2
     .     + 2d0*PM(19)**2*PM(18)**2 + PM(17)**2*PM(18)**2
     .     + PM(14)**2*PM(18)**2 + 6d0*PM(18)**4
     .     + 2d0*DSQRT(HB*PM(8))*PM(19)*PM(16)
     .     + 2d0*DSQRT(PM(10))*PM(14)*PM(18)*PM(16)
     .     - PM(10)*PM(16)**2 + PM(19)**2*PM(16)**2)
     .     - PM(1)**4/PM(13)**2*COEF/6d0*(PM(18)**2+PM(19)**2)
      Y(21)=COEF*(40d0/27d0*ALP1**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 2d0*PM(10)*PM(19)**2 - 26d0/9d0*PM(11)*PM(18)**2
     .     + 4d0*DSQRT(PM(10)*PM(8))*PM(17)*PM(18)
     .     - 6d0*PM(12)*PM(18)**2 - 32d0/3d0*G3*PM(18)**2
     .     + 2d0*HB*PM(18)**2 - 2d0*PM(10)*PM(16)**2
     .     + 12d0*PM(10)*PM(18)**2 + 2d0*PM(19)**2*PM(18)**2
     .     + 2d0*PM(17)**2*PM(18)**2 + 2d0*PM(14)**2*PM(18)**2
     .     + 12d0*PM(18)**4 + 4d0*DSQRT(PM(10))*PM(14)*PM(18)*PM(16))
     .     - PM(1)**4/PM(13)**2*COEF/3d0*PM(18)**2
      Y(22)=COEF*(10d0/27d0*ALP1**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 14d0/9d0*PM(11)*PM(19)**2 - 6d0*PM(12)*PM(19)**2
     .     - 32d0/3d0*G3*PM(19)**2 + 12d0*HB*PM(19)**2
     .     + 2d0*PM(10)*PM(19)**2 + 12d0*PM(19)**4 - 2d0*HB*PM(17)**2
     .     + 4d0*DSQRT(HB*HL)*PM(19)*PM(20) + 2d0*PM(19)**2*PM(20)**2
     .     + 4d0*DSQRT(HB)*PM(19)*PM(17)*PM(14) - 2d0*HB*PM(18)**2
     .     + 2d0*PM(19)**2*PM(14)**2 + 2d0*PM(19)**2*PM(18)**2
     .     + 4d0*DSQRT(HB*PM(8))*PM(19)*PM(16)
     .     + 2d0*PM(19)**2*PM(16)**2)
     .     - PM(1)**4/PM(13)**2*COEF/3d0*PM(19)**2
      Y(23)=COEF*(5d0/54d0*ALP1**2+3d0/2d0*ALP2**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
      Y(24)=COEF*(40d0/27d0*ALP1**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
      Y(25)=COEF*(10d0/27d0*ALP1**2
     .     +8d0/3d0*ALP3**2)*F2*N5*PM(1)**2
      Y(26)=COEF*(5d0/6d0*ALP1**2
     .     +3d0/2d0*ALP2**2)*F2*N5*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - HL*PM(17)**2 + 6d0*DSQRT(HB*HL)*PM(19)*PM(20)
     .     - 3d0*PM(11)*PM(20)**2 - 3d0*PM(12)*PM(20)**2
     .     + 2d0*HL*PM(20)**2 + 3d0*PM(19)**2*PM(20)**2 + 4d0*PM(20)**4
     .     + 2d0*DSQRT(HL)*PM(17)*PM(20)*PM(14) + PM(20)**2*PM(14)**2
     .     + 2d0*DSQRT(HL*PM(8))*PM(20)*PM(16) + PM(20)**2*PM(16)**2)
     .     - PM(1)**4/PM(13)**2*COEF/6d0*PM(20)**2
      Y(27)=COEF*(10d0/3d0*ALP1**2)*F2*N5*PM(1)**2
     .     +(COEF*PM(1))**2*(
     .     - 2*HL*PM(17)**2 + 12*DSQRT(HB*HL)*PM(19)*PM(20)
     .     - 6*PM(11)*PM(20)**2 - 6*PM(12)*PM(20)**2
     .     + 4*HL*PM(20)**2 + 6*PM(19)**2*PM(20)**2 + 8*PM(20)**4
     .     + 4*DSQRT(HL)*PM(17)*PM(20)*PM(14) + 2*PM(20)**2*PM(14)**2
     .     + 4*DSQRT(HL*PM(8))*PM(20)*PM(16) + 2*PM(20)**2*PM(16)**2)
     .     - PM(1)**4/PM(13)**2*COEF/3d0*PM(20)**2
      Y(28)=COEF*(5d0/6d0*ALP1**2
     .     +3d0/2d0*ALP2**2)*F2*N5*PM(1)**2
      Y(29)=COEF*(10d0/3d0*ALP1**2)*F2*N5*PM(1)**2
      Y(30)=PM(4)
      Y(31)=PM(5)
      Y(32)=PM(6)
      Y(33)=PM(7)
      Y(34)=0d0
      Y(35)=0d0

      X1=COEF*DLOG(PM(13)**2/Q2)
      X2=0d0

      CALL ODEINTS(Y,NN,X1,X2,EPS,DERIVSS,RKQSS,IFAIL)

      IF(IFAIL.NE.0)THEN
       RETURN
      ENDIF

      PS(1)=Y(17)
      PS(2)=Y(18)
      PS(3)=Y(19)
      PS(4)=Y(12)
      PS(5)=Y(13)
      PS(6)=Y(30)
      PS(7)=Y(31)
      PS(8)=Y(32)
      PS(9)=Y(33)
      PS(10)=Y(34)
      PS(11)=Y(4)
      PS(12)=Y(5)
      PS(13)=Y(6)
      PS(14)=(Y(1)+Y(2))/2d0

      END