text:
	module sub
	contains

c	This subroutine is used to calculate p and grad_p which use the data of
c	latitude and longtitude
c	x means lon,y means lat,N_max means Nmax here
	subroutine point(N_max,x,y,p,grad_p)
	implicit none
	real*8 x
	real*8 y
	real*8 p(N_max*(N_max+3)/2+1)
	real*8 grad_p(N_max*(N_max+3)/2+1)
	integer N_max,n,m,i

	p(1)=1
	p(2)=cos(y)
	p(3)=sin(y)
	grad_p(1)=0
	grad_p(2)=(cos(y)*p(2)-p(1))/sin(y)
	grad_p(3)=cos(y)*p(3)/sin(y)
	if (N_max>=2) then
	i=4
	do n=2,N_max
	do m=0,n
	if(m==n)then
	p(i)=p(i-n-1)*sin(y)*((2.0*m-1)/(2.0*m))**0.5
	grad_p(i)=n*cos(y)*p(i-n-1)*((2.0*m-1)/(2.0*m))**0.5
	else
	p(i)=((2*n-1)*cos(y)*p(i-n)-((n-1)**2-m**2)**
     &	0.5*p(i-n-n+1))/(n**2-m**2)**0.5
	grad_p(i)=(n*cos(y)*p(i)-(n**2-m**2)**
     &	0.5*p(i-n))/sin(y)
	end if
	i=i+1
	end do
	end do
	end if
	end subroutine point

c	得到观测方程的系数矩阵LA，j表示计算LA第几行，即第j个观测点的系数行。
	subroutine getLAofx(lat,lon,r,N_max,LA,data_N,j)
	implicit none
	integer N_max,j,data_N
	real*8 lat,lon,r,LA(data_N,N_max*(N_max+2))
	real*8 p(N_max*(N_max+3)/2+1),grad_p(N_max*(N_max+3)/2+1)
	integer i,n,m,k
	call point(N_max,lon,lat,p,grad_p)
c	ux=ux+(g(i)*cos(m*x)+h(i)*sin(m*x))*grad_p(i+1)
	i=1
	k=1
	do n=1,N_max
	do m=0,n
	LA(j,i)=cos(m*lon)*grad_p(i+1)
	if(m>0) then
	LA(j,k+N_max*(N_max+3)/2)=sin(m*lon)*grad_p(i+1)
	k=k+1
	end if
	i=i+1
	end do
	end do
	end subroutine

	end module

ext: f

max_stars_repo_name: DearDon/Geomagnetics

max_issues_repo_name: DearDon/Geomagnetics

max_forks_repo_name: DearDon/Geomagnetics


! { dg-do run }
! Test the fix for PR25746, in which dependency checking was not being
! done for elemental subroutines and therefore for interface assignments.
!
! This test is based on
! http://home.comcast.net/~kmbtib/Fortran_stuff/elem_assign.f90
! as reported by Harald Anlauf <anlauf@gmx.de> in the PR.
!
module elem_assign
   implicit none
   type mytype
      integer x
   end type mytype
   interface assignment(=)
      module procedure myassign
   end interface assignment(=)
   contains
      elemental subroutine myassign(x,y)
         type(mytype), intent(out) :: x
         type(mytype), intent(in) :: y
! Multiply the components by 2 to verify that this is being called.
         x%x = y%x*2
      end subroutine myassign
end module elem_assign

program test
   use elem_assign
   implicit none
   type(mytype) :: y(6), x(6) = (/mytype(1),mytype(20),mytype(300),&
                                  mytype(4000),mytype(50000),&
                                  mytype(1000000)/)
   type(mytype) :: z(2, 3)
! The original case - dependency between lhs and rhs.
   x = x((/2,3,1,4,5,6/))
   if (any(x%x .ne. (/40, 600, 2, 8000, 100000, 2000000/))) call abort ()
! Slightly more elborate case with non-trivial array ref on lhs.
   x(4:1:-1) = x((/1,3,2,4/))
   if (any(x%x .ne. (/16000, 1200, 4, 80, 100000, 2000000/))) call abort ()
! Check that no-dependence case works....
   y = x
   if (any(y%x .ne. (/32000, 2400, 8, 160, 200000, 4000000/))) call abort ()
! ...and now a case that caused headaches during the preparation of the patch
   x(2:5) = x(1:4)
   if (any(x%x .ne. (/16000, 32000, 2400, 8, 160, 2000000/))) call abort ()
! Check offsets are done correctly in multi-dimensional cases
   z = reshape (x, (/2,3/))
   z(:, 3:2:-1) = z(:, 1:2)
   y = reshape (z, (/6/))
   if (any(y%x .ne. (/ 64000, 128000, 19200, 64, 128000, 256000/))) call abort ()
end program test


text:
420 Relief is a marijuana consultation clinic (NOT a dispensary) located downtown. According to http://www.davisenterprise.com/localnews/davisfirstmedicalmarijuanaclinicopensdowntown/ an article in The Davis Enterprise, 420 Relief has a staff of licensed medical doctors who write medical marijuana prescriptions. The owner, Steve Maki, has three other clinics in the Sacramento area.  Treatable ailments include eczema, cancer, herpes, bipolar disorder, alcoholism and obesity.  Patients must show supporting documentation of symptoms from their primary care doctor.  Patients must also be 18 and older, or between 16 and 18 years old and accompanied by a parent or legal caregiver who has to sign the necessary paperwork.

The Davis location opened in 2011. It was previously occupied by Suite E and before that, Jacksonfly Recycled Fashion.

>BEWARE BEWARE!!  THIS PLACE USES DOCTORS THAT CANT WORK FOR A REAL CLINIC OR HOSPITAL DUE TO THEIR LICENSE STATUS!! BEWARE!! CARDS ARE NO GOOD REALLY ANYWHERE!!! <

20110517 10:58:09 nbsp The get your card today! line on their webpage is so shady.  I dont have any problem with medical marijuana, but clinics like this quite obviously set themselves up as subscription services where you pay a little money and get a cannabis card, done and done.  As long as the law limits marijuana use to medical applications, places that pass out cards like gym memberships rate pretty high on the shadiness factor in my book.

Lest anyone get the wrong impression, Im not making any comment on marijuana law.  Only on the way places like this do business given the current state of the law. Users/TomGarberson
  Sometimes laws shouldnt be followed. Users/hankim
This is the LAST thing Davis needs...whats next? Pawn shops? Check cashing stores? users/JoshLawson

20110517 14:49:57 nbsp If dude wants doctors to gamble with their licenses on giving cards to people who arent deserving that is the choice they make Users/StevenDaubert

20110601 20:34:23 nbsp My Herb card is invalid because the Bastard who runs the clinic is under investigation from the State!  Dont go there.  You are much better off going somewhere else.   Users/Friar

20110623 19:34:50 nbsp I currently work at 420 Relief in Davis. If you have any questions please feel free to actually come into our store.
231 E st. Suite 3.

 Users/RobertShiple

20110623 19:52:09 nbsp here is the article from the Davis Aggie regarding our recent opening in Davis.
http://theaggie.org/article/2011/06/02/newmedicinalmarijuanaclinicopensdowntown Users/RobertShiple

 This place uses doctors that cant work anywhere else. If you call and get their license numbers they all have issues or problems, Most dispensaries in Sacramento WILL NOT ACCEPT THEM! Beware!!

20110907 16:14:58 nbsp BEWARE!!! Recommendations from this place is sketchy!! You might be able to use them or maybe not, mine wasnt accepted but one place in Sacramento  BEWARE!!!! Users/BobbyMcAlister

20110912 03:05:13 nbsp They hand out cards like theyre candy and fail to explain to those qualifying for them that they are for home and personal use and do not give them free reign to smoke wherever and whenever they wish. It is a blight on the city of Davis and I am sick of these barely legal doctors and theyre idiotic scams.

Yes, some people are aided by the use of this drug... Hell half of my friends use it and I don;t care at all about its uses or its legality. What I do care about is when some scam artist clinic of doctors with visible problems longer than my arms try to practice by scamming people out of their hard earned cash for the possibility of obtaining a card.

Have I gotten a card? No. Have I been there? No. Do I know many who have? Yes.

This place is nothing more than a circus troupe of imbeciles trying to get a buck for doing nothing. Users/WesP

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
great(j1,ll1).
great(c1,p1).
great(ii1,ff1).
great(u1,cc1).
great(jj1,w1).
great(g1,x1).
great(a1,r1).
great(u1,dd1).
great(a1,dd1).
great(jj1,v1).
great(kk1,bb1).
great(r1,v1).
great(d1,dd1).
great(e1,u1).
great(h1,e1).
great(jj1,q1).
great(i1,s1).
great(hh1,m1).
great(z1,cc1).
great(i1,t1).
great(k1,o1).
great(cc1,x1).
great(kk1,x1).
great(u1,i1).
great(r1,b1).
great(k1,cc1).
great(ii1,t1).
great(cc1,n1).
great(ee1,s1).
great(h1,t1).
great(ff1,b1).
great(ii1,v1).
great(ff1,ee1).
great(hh1,i1).
great(f1,z1).
great(bb1,o1).
great(i1,kk1).
great(hh1,dd1).
great(x1,n1).
great(j1,aa1).
great(ff1,i1).
great(e1,o1).
great(ii1,g1).
great(ii1,i1).
great(v1,ee1).
great(f1,g1).
great(kk1,p1).
great(h1,cc1).
great(l1,cc1).
great(p1,m1).
great(bb1,t1).
great(e1,cc1).
great(j1,r1).
great(f1,d1).
great(c1,q1).
great(bb1,x1).
great(v1,n1).
great(k1,w1).
great(ii1,h1).
great(l1,kk1).
great(a1,b1).
great(u1,p1).
great(ii1,k1).
great(d1,v1).
great(f1,aa1).
great(l1,jj1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


program test_mead
	use mhm

	call meadfit_main()

end program


text:
!=====================================================================!
! Module to handle all MPI tasks
!
! Author: Komahan Boopathy (komahan@gatech.edu)
!=====================================================================!

module mpi_wrapper

  use constants, only: LOG_UNIT
  use variables, only: idproc, numproc, master, mpierr

#if defined MPI
  use mpi
#endif

  implicit none

  private
  public :: MPI_START_ALL, MPI_STOP_ALL

contains

  !-------------------------------------------------------------------!
  ! Perform MPI initialization, identify master and say hello!
  !-------------------------------------------------------------------!

  subroutine MPI_START_ALL()

#if defined MPI
    ! Initialize MPI
    call MPI_Init(mpierr)

    if (mpierr /= MPI_SUCCESS) then
       stop 'MPI Initializatin error'
    endif

    ! Get the ID and number of processors
    call MPI_Comm_rank(MPI_COMM_WORLD, idproc , mpierr)
    call MPI_Comm_size(MPI_COMM_WORLD, numproc, mpierr)

    ! find master
    if (idproc .eq. 0) master = .true.

    ! Say hello
    !if(master) then
    !   write(LOG_UNIT,'(2x,a,i3)') '>> Number of Processors = ', numproc
    !end if
    ! write(LOG_UNIT,'(2x,a,i3)') '>> Processor ', idproc + 1, 'of', numproc,'...[OK]'
#endif

  end subroutine MPI_START_ALL

  !-------------------------------------------------------------------!
  ! Stop all the MPI processors                                       !
  !-------------------------------------------------------------------!

  subroutine MPI_STOP_ALL()

#if defined MPI
    ! Make sure all procs are here
    call MPI_Barrier(MPI_COMM_WORLD, mpierr)
    if (mpierr /= MPI_SUCCESS) then
       stop 'MPI Barrier error'
    endif

    ! Finalize MPI
    call MPI_Finalize(mpierr)
    if (mpierr /= MPI_SUCCESS) then
       stop 'MPI Finalize error'
    endif
#endif

  end subroutine MPI_STOP_ALL

end module mpi_wrapper

max_stars_repo_name: komahanb/time-marching-adjoint

max_issues_repo_name: komahanb/diffusion

max_forks_repo_name: komahanb/diffusion


c    Subroutines for the program sdf_lister
c
c     Contents:
c             logical function open_output_sect
c
c             logical function write_output_sect
c
c
c****************************************************
c
       logical function open_output_sect(unit_sect)
c
c****************************************************
c
c     Open output data file and write the header for sectored matrix rates data
c     returns .true. if successful, .false. otherwise
c
c     12/1/95 by J. Dwyer
c          modified 2/28/96 by J. Dwyer added headerlines output to write_header
c                         4/8/96 by J. Dwyer changed headerlines to 40
c                         4/9/96 by J. Dwyer added year to output
c                         12/31/96 by J. Dwyer changed format to vse,E2,...,E6
c                         1/20/97 by J. Dwyer removed VSE sectors, added E1 tel 1 and 2,
c                                         removed anisotropy data, added VSE_AVE,VSEBAR and START
c                         2/6/97 by J. Dwyerincluded all E1 sectors
c                         2/7/97 by J. Dwyer added E7 - E10 for protons
c                          2/13/97 by J. Dwyer   added new KP list
c                          2/14/97 by J. Dwyer added write_kptitles
c
      include 'sdf_include.inc'   ! include type declarations
c
      integer   unit_sect   ! file unit of output file
      character*64     filename       ! complete name of output file
      character*4       file_prefix   ! prefix to be attached to output file name
      character*400   kpstring     !  string containing KP data names
      integer   headerlines
      integer   stringlength   ! length of KP string
c
      headerlines = 40
c
      file_prefix = 'SECT'
      write(filename, '(a4,a60)') file_prefix, outputfile
      open(unit=unit_sect,name=outputdir//filename,
     1     status='NEW',DISP='KEEP',RECL=32766, err=1000) ! open output data file
      print *, 'Opened output file: ', filename
c
      call write_header(unit_sect,file_prefix,headerlines)
c
      call write_kptitles(kpstring,stringlength)  !  get string containing KP data names
c
      write(unit_sect,200) kpstring(1:stringlength)
200  	format(' Time,',
     1         ' year,',
     1         ' day of the year,',
     1         ' delta_t/2,',
     1         ' saturation flag,',
     1         ' VSE, sigma,',
     1         ' VSEBAR, sigma,',
     1         ' START, sigma,',
     1         ' E1 (sect1;tel1), sigma,',' E1 (sect2;tel1), sigma,',
     1         ' E1 (sect3;tel1), sigma,',' E1 (sect4;tel1), sigma,',
     1         ' E1 (sect5;tel1), sigma,',' E1 (sect6;tel1), sigma,',
     1         ' E1 (sect7;tel1), sigma,',' E1 (sect8;tel1), sigma,',
     1         ' E1 (sect1;tel2), sigma,',' E1 (sect2;tel2), sigma,',
     1         ' E1 (sect3;tel2), sigma,', ' E1 (sect4;tel2), sigma,',
     1         ' E1 (sect5;tel2), sigma,',' E1 (sect6;tel2), sigma,',
     1         ' E1 (sect7;tel2), sigma,',' E1 (sect8;tel2), sigma,',
     1         ' E2 (sect1;tel1), sigma,',' E2 (sect2;tel1), sigma,',
     1         ' E2 (sect3;tel1), sigma,',' E2 (sect4;tel1), sigma,',
     1         ' E2 (sect5;tel1), sigma,',' E2 (sect6;tel1), sigma,',
     1         ' E2 (sect7;tel1), sigma,',' E2 (sect8;tel1), sigma,',
     1         ' E2 (sect1;tel2), sigma,',' E2 (sect2;tel2), sigma,',
     1         ' E2 (sect3;tel2), sigma,',' E2 (sect4;tel2), sigma,',
     1         ' E2 (sect5;tel2), sigma,', ' E2 (sect6;tel2), sigma,',
     1         ' E2 (sect7;tel2), sigma,', ' E2 (sect8;tel2), sigma,',
     1         ' E3 (sect1;tel1), sigma,',' E3 (sect2;tel1), sigma,',
     1         ' E3 (sect3;tel1), sigma,',' E3 (sect4;tel1), sigma,',
     1         ' E3 (sect5;tel1), sigma,',' E3 (sect6;tel1), sigma,',
     1         ' E3 (sect7;tel1), sigma,',' E3 (sect8;tel1), sigma,',
     1         ' E3 (sect1;tel2), sigma,', ' E3 (sect2;tel2), sigma,',
     1         ' E3 (sect3;tel2), sigma,',' E3 (sect4;tel2), sigma,',
     1         ' E3 (sect5;tel2), sigma,', ' E3 (sect6;tel2), sigma,',
     1         ' E3 (sect7;tel2), sigma,',' E3 (sect8;tel2), sigma,',
     1         ' E4 (sect1;tel1), sigma,', ' E4 (sect2;tel1), sigma,',
     1         ' E4 (sect3;tel1), sigma,', ' E4 (sect4;tel1), sigma,',
     1         ' E4 (sect5;tel1), sigma,', ' E4 (sect6;tel1), sigma,',
     1         ' E4 (sect7;tel1), sigma,', ' E4 (sect8;tel1), sigma,',
     1         ' E4 (sect1;tel2), sigma,',' E4 (sect2;tel2), sigma,',
     1         ' E4 (sect3;tel2), sigma,', ' E4 (sect4;tel2), sigma,',
     1         ' E4 (sect5;tel2), sigma,', ' E4 (sect6;tel2), sigma,',
     1         ' E4 (sect7;tel2), sigma,', ' E4 (sect8;tel2), sigma,',
     1         ' E5 (sect1;tel1), sigma,',' E5 (sect2;tel1), sigma,',
     1         ' E5 (sect3;tel1), sigma,', ' E5 (sect4;tel1), sigma,',
     1         ' E5 (sect5;tel1), sigma,', ' E5 (sect6;tel1), sigma,',
     1         ' E5 (sect7;tel1), sigma,', ' E5 (sect8;tel1), sigma,',
     1         ' E5 (sect1;tel2), sigma,', ' E5 (sect2;tel2), sigma,',
     1         ' E5 (sect3;tel2), sigma,', ' E5 (sect4;tel2), sigma,',
     1         ' E5 (sect5;tel2), sigma,', ' E5 (sect6;tel2), sigma,',
     1         ' E5 (sect7;tel2), sigma,', ' E5 (sect8;tel2), sigma,',
     1         ' E6 (sect1;tel1), sigma,', ' E6 (sect2;tel1), sigma,',
     1         ' E6 (sect3;tel1), sigma,', ' E6 (sect4;tel1), sigma,',
     1         ' E6 (sect5;tel1), sigma,',' E6 (sect6;tel1), sigma,',
     1         ' E6 (sect7;tel1), sigma,', ' E6 (sect8;tel1), sigma,',
     1         ' E6 (sect1;tel2), sigma,', ' E6 (sect2;tel2), sigma,',
     1         ' E6 (sect3;tel2), sigma,',' E6 (sect4;tel2), sigma,',
     1         ' E6 (sect5;tel2), sigma,',' E6 (sect6;tel2), sigma,',
     1         ' E6 (sect7;tel2), sigma,', ' E6 (sect8;tel2), sigma,',
     1         ' E7 (sect1;tel1), sigma,', ' E7 (sect2;tel1), sigma,',
     1         ' E7 (sect3;tel1), sigma,', ' E7 (sect4;tel1), sigma,',
     1         ' E7 (sect5;tel1), sigma,', ' E7 (sect6;tel1), sigma,',
     1         ' E7 (sect7;tel1), sigma,', ' E7 (sect8;tel1), sigma,',
     1         ' E7 (sect1;tel2), sigma,',' E7 (sect2;tel2), sigma,',
     1         ' E7 (sect3;tel2), sigma,', ' E7 (sect4;tel2), sigma,',
     1         ' E7 (sect5;tel2), sigma,', ' E7 (sect6;tel2), sigma,',
     1         ' E7 (sect7;tel2), sigma,', ' E7 (sect8;tel2), sigma,',
     1         ' E8 (sect1;tel1), sigma,',' E8 (sect2;tel1), sigma,',
     1         ' E8 (sect3;tel1), sigma,', ' E8 (sect4;tel1), sigma,',
     1         ' E8 (sect5;tel1), sigma,', ' E8 (sect6;tel1), sigma,',
     1         ' E8 (sect7;tel1), sigma,', ' E8 (sect8;tel1), sigma,',
     1         ' E8 (sect1;tel2), sigma,', ' E8 (sect2;tel2), sigma,',
     1         ' E8 (sect3;tel2), sigma,', ' E8 (sect4;tel2), sigma,',
     1         ' E8 (sect5;tel2), sigma,', ' E8 (sect6;tel2), sigma,',
     1         ' E8 (sect7;tel2), sigma,', ' E8 (sect8;tel2), sigma,',
     1         ' E9 (sect1;tel1), sigma,', ' E9 (sect2;tel1), sigma,',
     1         ' E9 (sect3;tel1), sigma,', ' E9 (sect4;tel1), sigma,',
     1         ' E9 (sect5;tel1), sigma,',' E9 (sect6;tel1), sigma,',
     1         ' E9 (sect7;tel1), sigma,', ' E9 (sect8;tel1), sigma,',
     1         ' E9 (sect1;tel2), sigma,', ' E9 (sect2;tel2), sigma,',
     1         ' E9 (sect3;tel2), sigma,', ' E9 (sect4;tel2), sigma,',
     1         ' E9 (sect5;tel2), sigma,', ' E9 (sect6;tel2), sigma,',
     1         ' E9 (sect7;tel2), sigma,', ' E9 (sect8;tel2), sigma,',
     1         ' E10 (sect1;tel1), sigma,', ' E10 (sect2;tel1), sigma,',
     1         ' E10 (sect3;tel1), sigma,', ' E10 (sect4;tel1), sigma,',
     1         ' E10 (sect5;tel1), sigma,', ' E10 (sect6;tel1), sigma,',
     1         ' E10 (sect7;tel1), sigma,', ' E10 (sect8;tel1), sigma,',
     1         ' E10 (sect1;tel2), sigma,', ' E10 (sect2;tel2), sigma,',
     1         ' E10 (sect3;tel2), sigma,', ' E10 (sect4;tel2), sigma,',
     1         ' E10 (sect5;tel2), sigma,',' E10 (sect6;tel2), sigma,',
     1         ' E10 (sect7;tel2), sigma,', ' E10 (sect8;tel2), sigma,',
     1           a)
c
      open_output_sect = .true.
      return  !  succussful
1000   open_output_sect = .false.  ! on error go to here
      print *, 'Error opening output sectored matrix rates file:',
     1      filename
      return
      end  ! end open_output_sect
c
c
c

c
c
c
c
c****************************************************
c
       logical function write_output_sect
     1       (unit_sect, nrow_sect,
     1        vse_ave, sigma_vse_ave, vse_ave_flag,
     1        disc, sigma_disc, disc_flag,
     1        E1_sect_tel1, E1_sect_tel1_sigma,
     1        E1_sect_tel1_flag,
     1        E1_sect_tel2, E1_sect_tel2_sigma,
     1        E1_sect_tel2_flag,
     1        E2_sect_tel1, E2_sect_tel1_sigma,
     1        E2_sect_tel1_flag,
     1        E2_sect_tel2, E2_sect_tel2_sigma,
     1        E2_sect_tel2_flag,
     1        E3_sect_tel1, E3_sect_tel1_sigma,
     1        E3_sect_tel1_flag,
     1        E3_sect_tel2, E3_sect_tel2_sigma,
     1        E3_sect_tel2_flag,
     1        E4_sect_tel1, E4_sect_tel1_sigma,
     1        E4_sect_tel1_flag,
     1        E4_sect_tel2, E4_sect_tel2_sigma,
     1        E4_sect_tel2_flag,
     1        E5_sect_tel1, E5_sect_tel1_sigma,
     1        E5_sect_tel1_flag,
     1        E5_sect_tel2, E5_sect_tel2_sigma,
     1        E5_sect_tel2_flag,
     1        E6_sect_tel1, E6_sect_tel1_sigma,
     1        E6_sect_tel1_flag,
     1        E6_sect_tel2, E6_sect_tel2_sigma,
     1        E6_sect_tel2_flag,
     1        E7_sect_tel1, E7_sect_tel1_sigma,
     1        E7_sect_tel1_flag,
     1        E7_sect_tel2, E7_sect_tel2_sigma,
     1        E7_sect_tel2_flag,
     1        E8_sect_tel1, E8_sect_tel1_sigma,
     1        E8_sect_tel1_flag,
     1        E8_sect_tel2, E8_sect_tel2_sigma,
     1        E8_sect_tel2_flag,
     1        E9_sect_tel1, E9_sect_tel1_sigma,
     1        E9_sect_tel1_flag,
     1        E9_sect_tel2, E9_sect_tel2_sigma,
     1        E9_sect_tel2_flag,
     1        E10_sect_tel1, E10_sect_tel1_sigma,
     1        E10_sect_tel1_flag,
     1        E10_sect_tel2, E10_sect_tel2_sigma,
     1        E10_sect_tel2_flag)
c
c****************************************************
c
c      writes SECT data to output file with file unit=unit_sect
c      returns .false. if number of output lines written exceeds
c      specified value, .true. otherwise.
c
c     11/27/95 by J. Dwyer
c                 modified     4/9/96 by J. Dwyer added year to output
c            5/6/96 changed outputstring to 4096 characters
c            12/31/96 by J. Dwyer changed format to vse,E2,...,E6
c            1/15/97 by J. Dwyer fixed bug with E4_sect_tel2_sigma
c            1/20/97 by J. Dwyer removed VSE sectors, added E1 tel 1 and 2,
c                               removed anisotropy data, added VSE_AVE, VSEBAR and START.
c            2/6/97 by J. Dwyerincluded all E1 sectors
c            2/7/97 by J. Dwyer added E7 - E10 for protons
c            2/13/97 by J. Dwyer   added new KP list
c            2/14/97 by J. Dwyer  added write_kpdata
c
      include 'sdf_include.inc'   ! include type declarations
      include 'hex_include.inc'   ! include type declarations
c
      real     E1_sect_tel1(8), E1_sect_tel1_sigma(8),
     1        E1_sect_tel2(8), E1_sect_tel2_sigma(8),
     1        E2_sect_tel1(8), E2_sect_tel1_sigma(8),
     1        E2_sect_tel2(8), E2_sect_tel2_sigma(8),
     1        E3_sect_tel1(8), E3_sect_tel1_sigma(8),
     1        E3_sect_tel2(8), E3_sect_tel2_sigma(8),
     1        E4_sect_tel1(8), E4_sect_tel1_sigma(8),
     1        E4_sect_tel2(8), E4_sect_tel2_sigma(8),
     1        E5_sect_tel1(8), E5_sect_tel1_sigma(8),
     1        E5_sect_tel2(8), E5_sect_tel2_sigma(8),
     1        E6_sect_tel1(8), E6_sect_tel1_sigma(8),
     1        E6_sect_tel2(8), E6_sect_tel2_sigma(8),
     1        E7_sect_tel1(8), E7_sect_tel1_sigma(8),
     1        E7_sect_tel2(8), E7_sect_tel2_sigma(8),
     1        E8_sect_tel1(8), E8_sect_tel1_sigma(8),
     1        E8_sect_tel2(8), E8_sect_tel2_sigma(8),
     1        E9_sect_tel1(8), E9_sect_tel1_sigma(8),
     1        E9_sect_tel2(8), E9_sect_tel2_sigma(8),
     1        E10_sect_tel1(8), E10_sect_tel1_sigma(8),
     1        E10_sect_tel2(8), E10_sect_tel2_sigma(8)
      real   vse_ave, sigma_vse_ave
      real   disc(4), sigma_disc(4)
      real  doy_real
      real  half_delta_t  ! 1/2 of elapsed time interval (sec)
      integer      E1_sect_tel1_flag(8),
     1             E1_sect_tel2_flag(8),
     1             E2_sect_tel1_flag(8),
     1             E2_sect_tel2_flag(8),
     1             E3_sect_tel1_flag(8),
     1             E3_sect_tel2_flag(8),
     1             E4_sect_tel1_flag(8),
     1             E4_sect_tel2_flag(8),
     1             E5_sect_tel1_flag(8),
     1             E5_sect_tel2_flag(8),
     1             E6_sect_tel1_flag(8),
     1             E6_sect_tel2_flag(8),
     1             E7_sect_tel1_flag(8),
     1             E7_sect_tel2_flag(8),
     1             E8_sect_tel1_flag(8),
     1             E8_sect_tel2_flag(8),
     1             E9_sect_tel1_flag(8),
     1             E9_sect_tel2_flag(8),
     1             E10_sect_tel1_flag(8),
     1             E10_sect_tel2_flag(8)
      integer   vse_ave_flag
      integer   disc_flag(4)
      integer   mid_time
      integer   unit_sect   ! output file unit
      integer    nrow_sect   !  number of rows written to output file
      integer    saturation_flag   ! 1 if any rate is saturated in mf, 0 otherwise
      integer    yr, mo, dy, hr, mn, se, doy
      integer    index     ! last index of outputstring written to
      character*20   isotime   ! sampex_timecon isotime array
      character*4096   outputstring
c
      mid_time = (stop_time_avg+start_time_avg)/2
      half_delta_t = (stop_time_avg-start_time_avg)/2.0
      call sampex_timcon(mid_time,yr,mo,dy,hr,mn,se,doy,isotime)
      doy_real = doy+(3600.0*hr+60.0*mn+se)/86400.0
c
      if (nrow_sect.lt.maximum_rows) then
c
       saturation_flag = 0
       if (vse_ave_flag.eq.-1) saturation_flag = 1
       if (disc_flag(1).eq.-1) saturation_flag = 1   ! VSEBAR
       if (disc_flag(2).eq.-1) saturation_flag = 1   ! START
c
       do 5 k=1,8
           if (E1_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E1_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E2_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E2_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E3_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E3_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E4_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E4_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E5_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E5_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E6_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E6_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E7_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E7_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E8_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E8_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E9_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E9_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E10_sect_tel1_flag(k).eq.-1)
     1            saturation_flag = 1
           if (E10_sect_tel2_flag(k).eq.-1)
     1            saturation_flag = 1
5     continue
c
         index = 1
c
         call write_line(index, outputstring,
     1          vse_ave, (vse_ave_flag.eq.1), .true.,
     1                dataformat_flag)
         call write_line(index, outputstring,
     1          sigma_vse_ave, (vse_ave_flag.eq.1), .true.,
     1                dataformat_flag)
c
            call write_line(index, outputstring,
     1          disc(1), (disc_flag(1).eq.1), .true.,
     1                dataformat_flag)
            call write_line(index, outputstring,
     1          sigma_disc(1), (disc_flag(1).eq.1), .true.,
     1                dataformat_flag)  ! VSEBAR
c
            call write_line(index, outputstring,
     1          disc(2), (disc_flag(2).eq.1), .true.,
     1                dataformat_flag)
            call write_line(index, outputstring,
     1          sigma_disc(2), (disc_flag(2).eq.1), .true.,
     1                dataformat_flag)   ! START
c
      do 10 k=1,8
            call write_line(index, outputstring,
     1         E1_sect_tel1(k), (E1_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E1_sect_tel1_sigma(k),
     1        (E1_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
10     continue
c
      do 20 k=1,8
            call write_line(index, outputstring,
     1         E1_sect_tel2(k), (E1_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E1_sect_tel2_sigma(k),
     1        (E1_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
20     continue
c
      do 70 k=1,8
            call write_line(index, outputstring,
     1         E2_sect_tel1(k), (E2_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E2_sect_tel1_sigma(k),
     1        (E2_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
70     continue
c
      do 80 k=1,8
            call write_line(index, outputstring,
     1         E2_sect_tel2(k), (E2_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E2_sect_tel2_sigma(k),
     1        (E2_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
80     continue

         do 90 k=1,8
            call write_line(index, outputstring,
     1         E3_sect_tel1(k), (E3_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E3_sect_tel1_sigma(k),
     1        (E3_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
90     continue
c
      do 100 k=1,8
             call write_line(index, outputstring,
     1         E3_sect_tel2(k), (E3_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E3_sect_tel2_sigma(k),
     1        (E3_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
100     continue

         do 110 k=1,8
            call write_line(index, outputstring,
     1         E4_sect_tel1(k), (E4_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E4_sect_tel1_sigma(k),
     1        (E4_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
110     continue

      do 120 k=1,8
            call write_line(index, outputstring,
     1         E4_sect_tel2(k), (E4_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E4_sect_tel2_sigma(k),
     1        (E4_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
120     continue
c
      do 130 k=1,8
            call write_line(index, outputstring,
     1         E5_sect_tel1(k), (E5_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E5_sect_tel1_sigma(k),
     1        (E5_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
130     continue
c
      do 140 k=1,8
            call write_line(index, outputstring,
     1         E5_sect_tel2(k), (E5_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E5_sect_tel2_sigma(k),
     1        (E5_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
140     continue
c
      do 150 k=1,8
            call write_line(index, outputstring,
     1         E6_sect_tel1(k), (E6_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E6_sect_tel1_sigma(k),
     1        (E6_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
150     continue
c
      do 160 k=1,8
            call write_line(index, outputstring,
     1         E6_sect_tel2(k), (E6_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E6_sect_tel2_sigma(k),
     1        (E6_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
160     continue
c
      do 170 k=1,8
            call write_line(index, outputstring,
     1         E7_sect_tel1(k), (E7_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E7_sect_tel1_sigma(k),
     1        (E7_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
170     continue
c
      do 180 k=1,8
            call write_line(index, outputstring,
     1         E7_sect_tel2(k), (E7_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E7_sect_tel2_sigma(k),
     1        (E7_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
180     continue
c
      do 190 k=1,8
            call write_line(index, outputstring,
     1         E8_sect_tel1(k), (E8_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E8_sect_tel1_sigma(k),
     1        (E8_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
190     continue
c
      do 200 k=1,8
            call write_line(index, outputstring,
     1         E8_sect_tel2(k), (E8_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E8_sect_tel2_sigma(k),
     1        (E8_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
200     continue
c
      do 210 k=1,8
            call write_line(index, outputstring,
     1         E9_sect_tel1(k), (E9_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E9_sect_tel1_sigma(k),
     1        (E9_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
210     continue
c
      do 220 k=1,8
            call write_line(index, outputstring,
     1         E9_sect_tel2(k), (E9_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E9_sect_tel2_sigma(k),
     1        (E9_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
220     continue
c
      do 230 k=1,8
            call write_line(index, outputstring,
     1         E10_sect_tel1(k), (E10_sect_tel1_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E10_sect_tel1_sigma(k),
     1        (E10_sect_tel1_flag(k).eq.1), .true.,
     1                dataformat_flag)
230     continue
c
      do 240 k=1,8
            call write_line(index, outputstring,
     1         E10_sect_tel2(k), (E10_sect_tel2_flag(k).eq.1),
     1         .true., dataformat_flag)
            call write_line(index, outputstring,
     1        E10_sect_tel2_sigma(k),
     1        (E10_sect_tel2_flag(k).eq.1), .true.,
     1                dataformat_flag)
240     continue
c
      call write_kpdata(index, outputstring)
c
         write(unit_sect,300) mo, dy, yr, hr, mn, se,
     1   yr, doy_real,
     1   half_delta_t, saturation_flag,
     1   outputstring(1:(index-1))
c
300 	   format(1x,i2'/'i2'/'i4,2x,i2':'i2':'i2,
     1       ', 'i4', 'f12.8', 'f12.4', 'i1',',a<index-1>)
              nrow_sect = nrow_sect + 1
      else
         print *, 'Number of rows exceeds maximum specified:'
     1            , maximum_rows
         write_output_sect = .false.
         return
      end if
c
      write_output_sect = .true.
      return
      end   ! end write_output_sect
c


text:
C Adventure - subroutines
C
C 30-Aug-93	RMS	Revised for portability.
C
	SUBROUTINE VOCAB(ID1,INIT,V)
C
C Look up ID1 in the vocabulary (ATAB) and return its "definition" (KTAB),
C or -1 if not found.  If INIT is positive, this is an init call setting
C up a keyword variable, and not finding it constitutes a bug.  It also means
C that only ktab values which taken over 1000 equal INIT may be considered.
C (thus "steps", which is a motion verb as well as an object, may be located
C as an object.)  And it also means the KTAB value is taken mod 1000.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
	CHARACTER*(*) ID1
C
	DO 1 I = 1,TABSIZ
	IF(KTAB(I).EQ.-1) GOTO 2
	IF(INIT.GE.0 .AND. KTAB(I)/1000.NE.INIT) GOTO 1
	IF(ATAB(I).EQ.ID1) GOTO 3
1	CONTINUE
	CALL BUG(21)
C
2	V = -1
	IF(INIT.LT.0) RETURN
	TYPE 100,ID1
100	FORMAT(' Keyword = ',A)
	CALL BUG(5)
C
3	V = KTAB(I)
	IF(INIT.GE.0) V = MOD(V,1000)
	RETURN
	END

	SUBROUTINE DSTROY(OBJECT)
C
C Permanently eliminate OBJECT by moving to a non-existent location.
C
	IMPLICIT INTEGER (A-Z)
C
	CALL MOVE(OBJECT,0)
	RETURN
	END

	SUBROUTINE JUGGLE(OBJECT)
C
C Juggle OBJECT by picking it up and putting it down again, the purpose
C being to get the object to the front of the chain of things at its loc.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	I = PLACE(OBJECT)
	J = FIXED(OBJECT)
	CALL MOVE(OBJECT,I)
	CALL MOVE(OBJECT+OBJMAX,J)
	RETURN
	END

	SUBROUTINE MOVE(OBJECT,WHERE)
C
C Place any OBJECT anyWHERE by picking it up and dropping it.  May already be
C toting, in which case the carry is a no-op.  Mustn't pick up objects which
c are not at any loc, since carry wants to remove objects from ATLOC chains.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	IF(OBJECT.GT.OBJMAX) GOTO 1
	FROM = PLACE(OBJECT)
	GOTO 2
1	FROM = FIXED(OBJECT-OBJMAX)
2	IF(FROM.GT.0. AND. FROM.LE.300) CALL CARRY(OBJECT,FROM)
	CALL DROP(OBJECT,WHERE)
	RETURN
	END

	INTEGER FUNCTION PUT(OBJECT,WHERE,PVAL)
C
C PUT is the same as MOVE, except it returns a value used to set up the
C negated PROP values for the repository objects.
C
	IMPLICIT INTEGER (A-Z)
C
	CALL MOVE(OBJECT,WHERE)
	PUT = (-1)-PVAL
	RETURN
	END

	SUBROUTINE CARRY(OBJECT,WHERE)
C
C Start toting OBJECT, removing it from the list of things at its former
C location.  Incr HOLDNG unless it was already being toted.  If OBJECT>OBJMAX
C (moving "fixed" second loc), don't change PLACE or HOLDNG.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	IF(OBJECT.GT.OBJMAX) GOTO 5
	IF(PLACE(OBJECT).EQ.-1) RETURN
	PLACE(OBJECT) = -1
	HOLDNG = HOLDNG+1
5	IF(ATLOC(WHERE).NE.OBJECT) GOTO 6
	ATLOC(WHERE) = LINK(OBJECT)
	RETURN
6	TEMP = ATLOC(WHERE)
7	IF(LINK(TEMP).EQ.OBJECT) GOTO 8
	TEMP = LINK(TEMP)
	GOTO 7
8	LINK(TEMP) = LINK(OBJECT)
	RETURN
	END

	SUBROUTINE DROP(OBJECT,WHERE)
C
C Place OBJECT at a given loc, prefixing it onto the ATLOC list.  Decrement
C HOLDNG if the object was being toted.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	IF(OBJECT.GT.OBJMAX) GOTO 1
	IF(PLACE(OBJECT).EQ.-1) HOLDNG = HOLDNG-1
	PLACE(OBJECT) = WHERE
	GOTO 2
1	FIXED(OBJECT-OBJMAX) = WHERE
2	IF(WHERE.LE.0) RETURN
	LINK(OBJECT) = ATLOC(WHERE)
	ATLOC(WHERE) = OBJECT
	RETURN
	END

	SUBROUTINE BUG(NUM)
	IMPLICIT INTEGER (A-Z)
C
C The following conditions are currently considered fatal bugs.  Numbers < 20
C are detected while reading the database; the others occur at run time.
C	0	MESSAGE LINE > 70 CHARACTERS
C	1	NULL LINE IN MESSAGE
C	2	TOO MANY WORDS OF MESSAGES
C	3	TOO MANY TRAVEL OPTIONS
C	4	TOO MANY VOCABULARY WORDS
C	5	REQUIRED VOCABULARY WORD NOT FOUND
C	6	TOO MANY RTEXT OR MTEXT MESSAGES
C	7	TOO MANY HINTS
C	8	LOCATION HAS COND BIT BEING SET TWICE
C	9	INVALID SECTION NUMBER IN DATABASE
C	10	INVALID ACTSPK (VERB) NUMBER
C	11	TOO MANY LOCATION (LTEXT OR STEXT) ENTRIES
C	12	TOO MANY PLAYER CLASSES
C	13	INVALID OBJECT NUMBER
C	20	SPECIAL TRAVEL (500>L>300) EXCEEDS GOTO LIST
C	21	RAN OFF END OF VOCABULARY TABLE
C	22	VOCABULARY TYPE (N/1000) NOT BETWEEN 0 AND 3
C	23	INTRANSITIVE ACTION VERB EXCEEDS GOTO LIST
C	24	TRANSITIVE ACTION VERB EXCEEDS GOTO LIST
C	25	CONDITIONAL TRAVEL ENTRY WITH NO ALTERNATIVE
C	26	LOCATION HAS NO TRAVEL ENTRIES
C	27	HINT NUMBER EXCEEDS GOTO LIST
C	28	INVALID MONTH RETURNED BY DATE FUNCTION
C
	TYPE 1, NUM
1	FORMAT (/' Fatal error',I3,'.  Consult your local wizard.'/)
	CALL EXIT
	END
C
	SUBROUTINE INIRND(I,J)
	IMPLICIT INTEGER(A-Z)
	COMMON /SEED/ RNSEED
C
	RNSEED = (ISHFT(I,16) + J) .OR. 1
	RETURN
	END
C
	INTEGER FUNCTION RND(RANGE)
	IMPLICIT INTEGER(A-Z)
	REAL RAN
	COMMON /SEED/ RNSEED
C
	RND = RAN(RNSEED) * FLOAT(RANGE)
	RETURN
	END

	SUBROUTINE SPEAK(N)
C
C Print the message in record n of the random message file.
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	IF(N.EQ.0) RETURN
	IF(TLINES(N)(1:3) .EQ. '>$<') RETURN
	TEMP = N
1	OLDRCN = TRECNO(TEMP)
	TYPE 2,TLINES(TEMP)(1:NBLEN(TLINES(TEMP)))
2	FORMAT(' ',A)
	TEMP = TEMP+1
	IF(TRECNO(TEMP) .EQ. OLDRCN) GO TO 1
	RETURN
	END

	SUBROUTINE PSPEAK(MSG,SKIP)
C
C Find the SKIP+1st message for object msg and print it.
C MSG should be the index of the object.
C (INVEN+N+1 message is PROP = N message).
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	M = PTEXT(MSG)
	IF(SKIP.LT.0) GOTO 9
	OLDRCN = MSG
	DO 3 I = 1,SKIP+1
1	M = M+1
	IF(TRECNO(M-1).EQ.OLDRCN) GO TO 1
	OLDRCN = TRECNO(M-1)
3	CONTINUE
	M = M-1
9	CALL SPEAK(M)
	RETURN
	END

	SUBROUTINE RSPEAK(MSG)
C
C Print the i-th "random" message (section 6 of database).
C
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	IF(MSG.NE.0) CALL SPEAK(RTEXT(MSG))
	RETURN
	END

	SUBROUTINE GETIN(WORD1,WORD2)
C
C Get a command from the adventurer.  Snarf out the first word, pad it with
C Blanks, and return it in WORD1.  If a second word appears, it is returned
C in WORD2, else WORD2 is set to blanks.
C
	IMPLICIT INTEGER (A-Z)
	CHARACTER*72 FRST
	CHARACTER*8 WORD1,WORD2
	CHARACTER*1 J
C
	TYPE 1
1	FORMAT(' >',$)
2	ACCEPT 3,FRST
3	FORMAT(A)
C
	LUCVT = ICHAR('A')-ICHAR('a')
	DO 5 I = 1,72
	IF((FRST(I:I).GE.'a') .AND. (FRST(I:I).LE.'z'))
	1	FRST(I:I) = CHAR(ICHAR(FRST(I:I))+LUCVT)
5	CONTINUE
	INLEN = NBLEN(FRST)
	IF(INLEN.LE.0) GO TO 2
	WORD1 = ' '
	WORD2 = ' '
	I = 1					! input ptr.
	OP = 0					! output ptr.
C
50      OP = OP+1				! adv output ptr.
	CP = 0					! char ptr = 0.
C
200     IF(I.GT.INLEN) GO TO 2000		! end of input?
	J = FRST(I:I)				! no, get character
	I = I+1					! advance ptr.
	IF(J.EQ.' ') GO TO 1000                 ! space?
	CP = CP+1				! adv char ptr.
	IF(CP.GT.8) GO TO 200			! too much input?
	IF(OP.EQ.1) WORD1(CP:CP) = J		! insert char in word1
	IF(OP.EQ.2) WORD2(CP:CP) = J		! insert char in word2
	GO TO 200
C
C Space.
C
1000    IF(CP.EQ.0) GO TO 200                   ! any word yet?
	GO TO 50                                ! yes, adv op.
C
C End of input, see if partial word available.
C
2000	IF(CP.EQ.0) OP = OP-1			! any last word?
	RETURN					! done
C
	END

	LOGICAL FUNCTION YES(X,Y,Z)
C
C Print message X, wait for yes/no answer.  If yes, print Y and return
C TRUE; if no, print Z and return FALSE.
C
	IMPLICIT INTEGER (A-Z)
	CHARACTER*8 REPLY,JUNK1
C
1	IF(X.NE.0) CALL RSPEAK(X)
	CALL GETIN(REPLY,JUNK1)
	IF(REPLY.EQ.'Y' .OR. REPLY.EQ.'YE' .OR. REPLY.EQ.'YES') GOTO 10
	IF(REPLY.EQ.'N' .OR. REPLY.EQ.'NO') GOTO 20
	TYPE 9
9	FORMAT(' Please answer the question.')
	GOTO 1
10	YES = .TRUE.
	IF(Y.NE.0) CALL RSPEAK(Y)
	RETURN
20	YES = .FALSE.
	IF(Z.NE.0) CALL RSPEAK(Z)
	RETURN
C
	END

        INTEGER FUNCTION NBLEN(STRING)
C
C Compute string length without trailing blanks
C
        IMPLICIT INTEGER(A-Z)
        CHARACTER*(*) STRING
C
        NBLEN = LEN(STRING)			! get nominal length
100     IF(NBLEN.LE.0) RETURN                   ! any string left?
        IF(STRING(NBLEN:NBLEN).NE.' ') RETURN   ! found a non-blank?
        NBLEN = NBLEN-1				! no, trim len by 1
        GO TO 100                               ! and continue.
C
        END

	SUBROUTINE SAVEGM(F2)
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	OPEN (UNIT=1,FILE='ASAVE',ACCESS='SEQUENTIAL',
	1	STATUS='UNKNOWN',FORM='UNFORMATTED',ERR=200)
C
	WRITE(1) VMAJ,VMIN,VEDIT
	WRITE(1) RTEXT,KTAB,ATAB
	WRITE(1) ATLOC,LINK,PLACE,FIXED,HOLDNG
	WRITE(1) PTEXT,ABB,LINUSE,TRVS,CLSSES
	WRITE(1) OLDLOC,LOC,CVAL,NEWLOC,KEY
	WRITE(1) PLAC,FIXD,ACTSPK,COND,HINTS
	WRITE(1) HNTMAX,PROP,TALLY,TALLY2,HINTLC
	WRITE(1) CHLOC,CHLOC2,DSEEN,DFLAG,DLOC,DALTLC
	WRITE(1) TURNS,LMWARN,KNFLOC,DETAIL,ABBNUM
	WRITE(1) NUMDIE,MAXDIE,DKILL,FOOBAR,BONUS
	WRITE(1) CLOCK1,CLOCK2,CLOSNG,PANIC,CLOSED
	WRITE(1) GAVEUP,SCORNG,ODLOC,CTEXT,STEXT,LTEXT
	WRITE(1) TRAVEL,TRVCON,TRVLOC,MAXTRS,HINTED
C
	TYPE 100
100	FORMAT(' Saved.')
	CLOSE (UNIT=1)
	F2=0
	RETURN
C
200	TYPE 300
300	FORMAT(' Open error, SAVE fails.')
	F2=-1
	RETURN
	END

	SUBROUTINE RSTRGM(F2)
	IMPLICIT INTEGER (A-Z)
	INCLUDE 'aparam.for'
C
	OPEN (UNIT=1,FILE='ASAVE',ACCESS='SEQUENTIAL',
	1	STATUS='OLD',READONLY,FORM='UNFORMATTED',ERR=200)
C
	READ(1) I1,I2,I3
	IF((I1.NE.VMAJ).OR.(I2.NE.VMIN)) GO TO 50
C
	READ(1) RTEXT,KTAB,ATAB
	READ(1) ATLOC,LINK,PLACE,FIXED,HOLDNG
	READ(1) PTEXT,ABB,LINUSE,TRVS,CLSSES
	READ(1) OLDLOC,LOC,CVAL,NEWLOC,KEY
	READ(1) PLAC,FIXD,ACTSPK,COND,HINTS
	READ(1) HNTMAX,PROP,TALLY,TALLY2,HINTLC
	READ(1) CHLOC,CHLOC2,DSEEN,DFLAG,DLOC,DALTLC
	READ(1) TURNS,LMWARN,KNFLOC,DETAIL,ABBNUM
	READ(1) NUMDIE,MAXDIE,DKILL,FOOBAR,BONUS
	READ(1) CLOCK1,CLOCK2,CLOSNG,PANIC,CLOSED
	READ(1) GAVEUP,SCORNG,ODLOC,CTEXT,STEXT,LTEXT
	READ(1) TRAVEL,TRVCON,TRVLOC,MAXTRS,HINTED
C
	TYPE 100
100	FORMAT(' Restored.')
	CLOSE (UNIT=1)
	F2=0
	RETURN
C
50	TYPE 110
110	FORMAT(' File is obsolete, RESTORE fails.')
	CLOSE (UNIT=1)
	F2=-1
	RETURN
C
200	TYPE 300
300	FORMAT(' Open error, RESTORE fails.')
	F2=-1
	RETURN
	END


module bootstrap

! JLR 25/8/2015
!
! Bootstrap confidence interval for Gaussian Kernel Correlation

use correlation

implicit none

contains

real function cpdf(x)
implicit none
real, intent(in) :: x

! cumulative probability density function
!
! Abramowitz and Stegun 26.2.17

real*8 :: z

if (x.ge.0) then
  z=1d0-0.5d0/(1d0+0.196854*x+0.115194*x**2+0.000344*x**3+0.019527*x**4)**4
else
  z=0.5d0/(1d0-0.196854*x+0.115194*x**2-0.000344*x**3+0.019527*x**4)**4
endif

cpdf=z

end function

real function icpdf(x)
implicit none
real, intent(in) :: x

! inverse cumulative probability density function
!
! Algorithm AS 241 Wichura 1988

real*8 :: q,r,zp
real*8, parameter :: a0=3.3871327179d0, a1=50.434271938d0, a2=159.29113202d0
real*8, parameter :: a3=59.109374720d0, b1=17.895169469d0, b2=78.757757664d0
real*8, parameter :: b3=67.187563600d0, c0=1.4234372777d0, c1=2.7568153900d0
real*8, parameter :: c2=1.3067284816d0, c3=0.17023821103d0, d1=0.73700164250d0
real*8, parameter :: d2=0.12021132975d0, e0=6.6579051150d0, e1=3.0812263860d0
real*8, parameter :: e2=0.42868294337d0, e3=0.017337203997d0
real*8, parameter :: f1=0.24197894225d0, f2=0.012258202635d0

q=x-0.5d0

if (abs(q).le.0.425d0) then
  r=0.425d0**2-q**2
  zp=q*(((a3*r+a2)*r+a1)*r+a0)/(((b3*r+b2)*r+b1)*r+1d0)
else
  if ((x.eq.0d0).or.(x.eq.1d0)) then
    r=27d0
  else
    r=sqrt(-log(min(x,1d0-x)))
    r=min(r,27d0)
  endif
  if (r.le.5.0d0) then
    r=r-1.6d0
    zp=(((c3*r+c2)*r+c1)*r+c0)/((d2*r+d1)*r+1d0)
  else
    r=r-5.0d0
    zp=(((e3*r+e2)*r+e1)*r+e0)/((f2*r+f1)*r+1d0)
  endif
  if (q.lt.0) zp=-zp
endif

icpdf=zp

end function

subroutine add_one_sort(x)

implicit none

! add an element to a sorted list)

real, intent(inout) :: x(:)
real :: temp
integer :: i,j,n

n=size(x,1)

do i=1,n-1
  if (x(n)<x(i)) then
    temp=x(n)
    do j=n-1,i,-1
      x(j+1)=x(j)
    enddo
    x(i)=temp
    exit
  endif
enddo

end subroutine

subroutine joint_sort(a,ia,n)

implicit none

! sort a (union of two sorted lists, joined at n) perform sam operation on ia

real, intent(inout) :: a(:)
integer, intent(inout) :: ia(:)
integer, intent(in) :: n

integer :: i,j,k,m
real :: tmp(size(a,1))
integer :: itmp(size(a,1))

tmp(:)=a(:)
itmp(:)=ia(:)

i=1
j=n+1
k=1
m=size(a,1)

a(:)=-999
ia(:)=-999

do
  if (tmp(i).le.tmp(j)) then
    a(k)=tmp(i)
    ia(k)=itmp(i)
    i=i+1
  else
    a(k)=tmp(j)
    ia(k)=itmp(j)
    j=j+1
  endif
  k=k+1
  if (i>n) then
    a(k:m)=tmp(j:m)
    ia(k:m)=itmp(j:m)
    exit
  endif
  if (j>m) then
    a(k:m)=tmp(i:n)
    ia(k:m)=itmp(i:n)
    exit
  endif
!  if ((i.eq.n).and.(j.eq.m)) exit
enddo

end subroutine

function bootstrap_ci(x,y,tx,ty)
real ::  bootstrap_ci(1:2)

real, intent(in) :: x(:),y(:),tx(:),ty(:)
integer :: n,m,nm,i,j,k,iter,n_low,ierr,ci1,ci2,outer
real :: taux,tauy,txmean,tymean,rnd,p,logp,tmean,dt
real :: orig,pd,z0,a,astar,asum,a_num,a_den,b0,b1
real, allocatable :: temp(:),tstar(:)
real, allocatable :: nx(:),ny(:),ntx(:),nty(:)
real, allocatable :: bse(:)
integer, allocatable :: tmp_indx(:),indx(:,:),seq(:)
integer, parameter :: max_iter=2000
real :: med_min(25),med_max(25)
logical :: init(2)

n=size(x,1)
m=size(y,1)

allocate (temp(m+n),tstar(n+m),tmp_indx(n+m),indx(n+m,2),seq(n+m+1))
allocate (nx(n+m),ny(n+m),ntx(n+m),nty(n+m),bse(max_iter))

orig=correlate_gaussian(x,y,tx,ty)

! estimate persistence times

txmean=(tx(n)-tx(1))/(n-1)
tymean=(ty(m)-ty(1))/(m-1)
tmean=max(txmean,tymean)

do i=1,n/16
  temp(1:n)=tx(1:n)+i*txmean
  taux=correlate_gaussian(x,x,tx,temp(1:n))
  if (abs(taux)<0.368) exit
enddo
taux=i*txmean

do i=1,m/16
  temp(1:m)=ty(1:m)+i*tymean
  tauy=correlate_gaussian(y,y,temp(1:m),ty)
  if (abs(tauy)<0.368) exit
enddo
tauy=i*tymean

p=1.0-tmean/(4*max(taux,tauy))
logp=log(p)

! build integrated time-list
temp(1:n)=tx(1:n)
temp(n+1:n+m)=ty(1:m)
! index elements for series x  are positive integers, and negative for y
do i=1,n
  tmp_indx(i)=i
enddo
do i=1,m
  tmp_indx(n+i)=-i
enddo
call joint_sort(temp,tmp_indx,n)
j=1
indx(:,:)=0 ! set no data value
do i=1,n+m
  tstar(j)=temp(i)
  if (tmp_indx(i)>0) then
    indx(j,1)=tmp_indx(i)
  else
    indx(j,2)=tmp_indx(i)
  endif
  if (i.lt.n+m) then
    if (temp(i).ne.temp(i+1)) j=j+1
  endif
enddo
nm=j

! do the bootstraps
do outer=1,25 ! outer iteration
 n_low=0
 iter=1
 do ! inner iteration
  j=1
  init(:)=.false.
  ! random starting point
  call random_number(rnd)
  i=floor((nm-1)*rnd)+1
  if (indx(i,1).ne.0) then
    seq(j)=indx(i,1)
    j=j+1
    init(1)=.true.
  endif
  if (indx(i,2).ne.0) then
    seq(j)=indx(i,2)
    j=j+1
    init(2)=.true.
  endif
  ! loop for rest of points
  do
    call random_number(rnd)
    if (i.gt.1) then
      dt=tstar(i)-tstar(i-1)
    else
      dt=tstar(2)-tstar(1)
    endif
    if (rnd.gt.exp(dt/tmean*logp)) then ! start a new block
      if (.not.(init(1).and.init(2))) then ! if haven't already sampled from both series, restart iteration
        j=1
        init(:)=.false.
      endif
      call random_number(rnd)
      i=floor((nm-1)*rnd)+1
      if (indx(i,1).ne.0) then
        seq(j)=indx(i,1)
        j=j+1
      endif
      if (indx(i,2).ne.0) then
        seq(j)=indx(i,2)
        j=j+1
      endif
    else ! continue with existing block
      i=i+1
      if (i.gt.nm) i=1
      if (indx(i,1).ne.0) then
        seq(j)=indx(i,1)
        j=j+1
        init(1)=.true.
      endif
      if (indx(i,2).ne.0) then
        seq(j)=indx(i,2)
        j=j+1
        init(2)=.true.
      endif
    endif
    if (j.gt.n+m) exit ! have enough data
  enddo
  ! built up two series and correspondind time bases
  ! find start times for each series
  init(:)=.false.
  j=1
  k=1
  do i=1,n+m
    if (seq(i)>0) then
      if (init(1)) then ! already have time base
        nx(j)=x(seq(i))
        if (seq(i).gt.1) then
          dt=tx(seq(i))-tx(seq(i)-1)
        else
          dt=tx(2)-tx(1)
        endif
        ntx(j)=ntx(j-1)+dt
      else
        nx(j)=x(seq(i))
        init(1)=.true.
        if (init(2)) then
          if (ty(-seq(i-1)).gt.tx(seq(i))) then
            if (seq(i).gt.1) then
              dt=tx(seq(i))-tx(seq(i)-1)
            else
              dt=tx(2)-tx(1)
            endif
            ntx(j)=nty(k-1)+dt
          else
            ntx(j)=nty(k-1)+tx(seq(i))-ty(-seq(i-1))
          endif
        else
          ntx(j)=tx(seq(i))
        endif
      endif
      j=j+1
    else
      if (init(2)) then ! already have time base
        ny(k)=y(-seq(i))
        if (-seq(i).gt.1) then
          dt=ty(-seq(i))-ty(-seq(i)-1)
        else
          dt=ty(2)-ty(1)
        endif
        nty(k)=nty(k-1)+dt
      else
        ny(k)=y(-seq(i))
        init(2)=.true.
        if (init(1)) then
          if (tx(seq(i-1)).gt.ty(-seq(i))) then
            if (-seq(i).gt.1) then
              dt=ty(-seq(i))-ty(-seq(i)-1)
            else
              dt=ty(2)-ty(1)
            endif
            nty(k)=ntx(j-1)+dt
          else
            nty(k)=ntx(j-1)+ty(-seq(i))-tx(seq(i-1))
          endif
        else
          nty(k)=ty(-seq(i))
        endif
      endif
      k=k+1
    endif
  enddo
  j=j-1
  k=k-1
  bse(iter)=correlate_gaussian(nx(1:j),ny(1:k),ntx(1:j),nty(1:k))
  if (bse(iter)<orig) n_low=n_low+1
  if (iter>1) call add_one_sort(bse(1:iter))
  iter=iter+1
  if (iter.gt.max_iter) exit
 enddo

!calcalate bias parameter
 z0=icpdf(real(n_low)/max_iter)

!calculate accelation parameter
 asum=0.0
 do i=1,max_iter
   asum=asum+bse(i)
 enddo
 astar=0.0
 do i=1,max_iter
   astar=astar+(asum-bse(i))/(max_iter-1)
 enddo
 astar=astar/max_iter
 a_num=0.0
 a_den=0.0
 do i=1,max_iter
   a_num=a_num+(astar-(asum-bse(i))/(max_iter-1))**3
   a_den=a_den+(astar-(asum-bse(i))/(max_iter-1))**2
 enddo
 a=a_num/(6.0*exp(1.5*log(a_den)))

! calculate interval points
 b0=z0+(z0-1.960)/(1-a*(z0-1.960))
 b1=z0+(z0+1.960)/(1-a*(z0+1.960))

 ci1=min(max(floor(iter*cpdf(b0)),1),max_iter)
 ci2=max(min(ceiling(iter*cpdf(b1)),max_iter),1)

 med_min(outer)=bse(ci1)
 med_max(outer)=bse(ci2)

 if (outer>1) then
   call add_one_sort(med_min(1:outer))
   call add_one_sort(med_max(1:outer))
 endif

! print *,bse(ci1),bse(ci2)

enddo

bootstrap_ci(1)=med_min(13)
bootstrap_ci(2)=med_max(13)

deallocate (temp,tstar,indx,tmp_indx,seq)
deallocate (nx,ny,ntx,nty,bse)

end function

end module


MODULE funlib
  USE global
  IMPLICIT NONE
CONTAINS
  RECURSIVE FUNCTION factorial(i) RESULT(fac)
    IMPLICIT NONE
    INTEGER::fac
    INTEGER,INTENT(IN)::i
    INTEGER::init=0,j
    LOGICAL,DIMENSION(12)::lfacsave
    INTEGER,DIMENSION(12)::facsave
    SAVE facsave,init,lfacsave
    IF(init.EQ.0)THEN
       lfacsave(1:12)=.FALSE.
       init=1
    ENDIF
    IF(i.LT.0)THEN
       WRITE(*,*)"ERROR: i<0 in factorial with i=",i
       STOP
    ENDIF
    IF(i.GT.12)THEN
       WRITE(*,*)"ERROR: i > 12, please take long type (KIND=LINT) integer"
       STOP
    ENDIF
    IF(i.EQ.0)THEN
       fac=1
    ELSE
       IF(lfacsave(i))THEN
          fac=facsave(i)
       ELSE
          fac=i*factorial(i-1)
          facsave(i)=fac
          lfacsave(i)=.TRUE.
       ENDIF
    ENDIF
  END FUNCTION factorial

  SUBROUTINE PCL2PE(NLOOPLINE,PCL,PE)
    ! transfer the loop momenta q,q+p1,q+p1+p2,q+p1+p2+p3,...
    ! to p1,p2,p3,p4,...
    IMPLICIT NONE
    INTEGER,INTENT(IN)::NLOOPLINE
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PCL
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(OUT)::PE
    INTEGER::i
    DO i=1,NLOOPLINE-1
       PE(i,0:3)=PCL(i+1,0:3)-PCL(i,0:3)
    ENDDO
    PE(NLOOPLINE,0:3)=PCL(1,0:3)-PCL(NLOOPLINE,0:3)
    RETURN
  END SUBROUTINE PCL2PE

  RECURSIVE SUBROUTINE calc_all_integers(n,ntot,i,sol)
    ! finds all the non-negative solutions to x1,...,xn
    ! that x1+x2+...+xn==i
    ! the number of solutions should be ntot=C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
    ! it can be recycled for all phase space point
    IMPLICIT NONE
    INTEGER,INTENT(IN)::n,ntot,i
    INTEGER,DIMENSION(ntot,n),INTENT(OUT)::sol
    INTEGER::j,k,k0
    IF(i.EQ.0)THEN
       sol(1,1:n)=0
       RETURN
    ENDIF
    IF(n.EQ.1)THEN
       sol(1,1)=i
       RETURN
    ENDIF
    k0=0
    DO j=0,i
       k=xiarray_arg1(i-j,n-1)
       sol(k0+1:k0+k,1)=j
       CALL calc_all_integers(n-1,k,i-j,sol(k0+1:k0+k,2:n))
       k0=k0+k
    ENDDO
    RETURN
  END SUBROUTINE calc_all_integers

  SUBROUTINE all_Integers(n,ntot,i,sol,factor)
    ! finds all the non-negative solutions to x1,...,xn
    ! that x1+x2+...+xn==i
    ! the number of solutions should be ntot=C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
    ! it can be recycled for all phase space point
    IMPLICIT NONE
    INTEGER,INTENT(IN)::n,ntot,i
    INTEGER,DIMENSION(ntot,n),INTENT(OUT)::sol
    REAL(KIND(1d0)),DIMENSION(ntot),INTENT(OUT)::factor
    INTEGER::ifirst=0,j,jk,k,ntemptot
    INTEGER::maxxiarray,inum,nnum,maxxiarray1,maxxiarray2
    INTEGER::maxfactor_xiarray
    SAVE ifirst
    ! calculate xiarray_i_n first
    IF(ifirst.EQ.0)THEN
       maxxiarray=(MAXRANK_IREGI+1)*(MAXNLOOP_IREGI-1)
       IF(.NOT.ALLOCATED(xiarray))THEN
          ALLOCATE(xiarray(maxxiarray))
       ENDIF
       maxfactor_xiarray=(MAXRANK_IREGI+1)*(MAXNLOOP_IREGI-2)
       IF(.NOT.ALLOCATED(factor_xiarray))THEN
          ALLOCATE(factor_xiarray(MAXRANK_IREGI+2:maxfactor_xiarray+MAXRANK_IREGI+1))
       ENDIF
       ! x1+x2+...+xn==i
       ! C(i+n-1)^(n-1)=(i+n-1)!/(n-1)!/i!
       DO j=1,maxxiarray
          inum=MOD(j-1,MAXRANK_IREGI+1) ! i
          nnum=(j-1)/(MAXRANK_IREGI+1)+1 ! n
          maxxiarray1=xiarray_arg1(inum,nnum)
          maxxiarray2=nnum
          IF(.NOT.ALLOCATED(xiarray(j)%xiarray_i_n))THEN
             ALLOCATE(xiarray(j)%xiarray_i_n(maxxiarray1,maxxiarray2))
          ENDIF
          IF(nnum.EQ.1)THEN
             xiarray(j)%xiarray_i_n(1,1)=j-1
             CYCLE
          ENDIF
          CALL calc_all_integers(nnum,maxxiarray1,inum,&
               xiarray(j)%xiarray_i_n(1:maxxiarray1,1:maxxiarray2))
          IF(.NOT.ALLOCATED(factor_xiarray(j)%factor_xiarray_i_n))THEN
             ALLOCATE(factor_xiarray(j)%factor_xiarray_i_n(maxxiarray1))
          ENDIF
          DO jk=1,maxxiarray1
             factor_xiarray(j)%factor_xiarray_i_n(jk)=DBLE(factorial(inum))
             DO k=1,nnum
                factor_xiarray(j)%factor_xiarray_i_n(jk)=factor_xiarray(j)%factor_xiarray_i_n(jk)/&
                     DBLE(factorial(xiarray(j)%xiarray_i_n(jk,k)))
             ENDDO
          ENDDO
          ntot_xiarray(inum,nnum)=maxxiarray1
       ENDDO
       ifirst=1
    ENDIF
    IF(n.GT.MAXNLOOP_IREGI-1.OR.n.LT.1)THEN
       WRITE(*,100)"ERROR: n is out of range 1<=n<=",MAXNLOOP_IREGI-1," in all_integers"
       STOP
    ENDIF
    IF(i.GT.MAXRANK_IREGI.OR.i.LT.0)THEN
       WRITE(*,100)"ERROR: r is out of range 0<=r<=",MAXRANK_IREGI," in all_integers"
       STOP
    ENDIF
    IF(n.EQ.1.AND.ntot.NE.1)THEN
       WRITE(*,*)"ERROR:ntot should be 1 when n=1 in all_integers"
       STOP
    ENDIF
    IF(n.GE.2.AND.n.LE.MAXNLOOP_IREGI-1)THEN
       ! Make it work in MadLoop, otherwise it is wrong
       IF(ntot.NE.ntot_xiarray(i,n))THEN
          WRITE(*,*)"ERROR: ntot is not correct in all_integers"
          STOP
       ENDIF
    ENDIF
    j=(n-1)*(MAXRANK_IREGI+1)+i+1
    SELECT CASE(n)
       CASE(1)
          sol(1:ntot,1:n)=xiarray(j)%xiarray_i_n(1:1,1:n)
          factor(1)=1d0
       CASE DEFAULT
          sol(1:ntot,1:n)=xiarray(j)%xiarray_i_n(1:ntot,1:n)
          factor(1:ntot)=factor_xiarray(j)%factor_xiarray_i_n(1:ntot)
    END SELECT
    RETURN
100 FORMAT(2X,A31,I2,A16)
  END SUBROUTINE all_Integers

  SUBROUTINE calc_factorial_pair
    IMPLICIT NONE
    INTEGER::i,j
    factorial_pair(1:MAXINDICES_IREGI,0)=1d0
    DO i=1,MAXINDICES_IREGI
       DO j=1,MAXRANK_IREGI
          factorial_pair(i,j)=DBLE(factorial(i+j-1))&
               /DBLE(factorial(i-1))/DBLE(factorial(j))
       ENDDO
    ENDDO
    RETURN
  END SUBROUTINE calc_factorial_pair

  FUNCTION number_coefs_for_rank(rank) RESULT(res)
    IMPLICIT NONE
    INTEGER,INTENT(IN)::rank
    INTEGER::res,i
    IF(rank.LT.0)res=0
    res=0
    DO i=0,rank
       res=res+((3+i)*(2+i)*(1+i))/6
    ENDDO
    RETURN
  END FUNCTION number_coefs_for_rank

  SUBROUTINE timestamp
    IMPLICIT NONE
    CHARACTER(len=8)::ampm
    INTEGER::d
    INTEGER::h
    INTEGER::m,n
    INTEGER::mm
    CHARACTER( len = 9 ),PARAMETER,DIMENSION(12) :: month = (/ &
         'January  ', 'February ', 'March    ', 'April    ', &
         'May      ', 'June     ', 'July     ', 'August   ', &
         'September', 'October  ', 'November ', 'December ' /)
    INTEGER::dn
    INTEGER::s
    INTEGER,DIMENSION(8)::values
    INTEGER::y
    CALL date_and_time(values=values)
    y = values(1)
    m = values(2)
    d = values(3)
    h = values(5)
    n = values(6)
    s = values(7)
    mm = values(8)
    IF( h.LT.12 )THEN
       ampm = 'AM'
    ELSEIF( h.EQ.12 )THEN
       IF( n.EQ.0 .AND. s.EQ.0 )THEN
          ampm = 'Noon'
       ELSE
          ampm = 'PM'
       ENDIF
    ELSE
       h = h - 12
       IF( h.LT.12 )THEN
          ampm = 'PM'
       ELSEIF( h.EQ.12 )THEN
          IF(n.EQ.0.AND.s.EQ.0)THEN
             ampm = 'Midnight'
          ELSE
             ampm = 'AM'
          ENDIF
       ENDIF
    ENDIF
    WRITE( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
         d, TRIM( month(m) ), y, h, ':', n, ':', s, '.', mm, TRIM( ampm )
    RETURN
  END SUBROUTINE timestamp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                               !
! FOLLOWING SUBROUTINES/FUNCTIONS ARE ONLY FOR THE REORDERING THE TIR COEFS     !
!                                                                               !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE SORT_IREGICOEFS(RANK,NLOOPCOEFS,OLDCOEFS,NEWCOEFS)
    !
    ! CONVERT THE OUTPUT OF IREGI TO THAT OF (NEW) MADLOOP
    !
    ! THE NEW OUTPUT OF COEFS FROM MADLOOP IS
    ! RANK=0: (,)
    ! RANK=1: (0,),(1,),(2,),(3,)
    ! RANK=2: (0,0),(0,1),(1,1),(0,2),(1,2),(2,2),(0,3),(1,3),(2,3),(3,3)
    ! ...
    ! THE OLD OUTPUT OF COEFS FROM MADLOOP IS
    ! RANK=0: (,)
    ! RANK=1: (0,),(1,),(2,),(3,)
    ! RANK=2: (0,0),(0,1),(0,2),(0,3),(1,1),(2,1),(3,1),(2,2),(2,3),(3,3)
    ! ...
    !
    ! ARGUMENTS
    IMPLICIT NONE
    INTEGER,INTENT(IN)::RANK,NLOOPCOEFS
    COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3),INTENT(IN)::OLDCOEFS
    COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3),INTENT(OUT)::NEWCOEFS
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I
    ! MAX RANK SET AS 10
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1)::POS
    SAVE POS
    LOGICAL::INIT=.TRUE.
    SAVE INIT
    ! ----------
    ! BEGIN CODE
    ! ----------

    IF(INIT)THEN
       IF(NLOOPCOEFS.GT.LOOPMAXCOEFS_IREGI)&
            STOP "ERROR:LOOPMAXCOEFS_IREGI IS TOO SMALL!!!"
       INIT=.FALSE.
       ! ASSIGN THE POSITION OF POS FOR SWAP
       CALL ASSIGN_PJPOS(POS)
    ENDIF

    DO I=0,NLOOPCOEFS-1
    !   NEWCOEFS(I,1:3)=OLDCOEFS(POS(I),1:3)
       NEWCOEFS(POS(I),1:3)=OLDCOEFS(I,1:3)
    ENDDO

  END SUBROUTINE SORT_IREGICOEFS

  SUBROUTINE ASSIGN_PJPOS(POS)
    !
    ! ASSIGN THE POSITION OF POS FOR SWAP
    !
    !
    IMPLICIT NONE
    !
    ! CONSTANS
    !
    ! MAX RANK SET AS 10
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001,MAXRANK=10
    !
    ! ARGUMENTS
    !
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1),INTENT(OUT)::POS
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I,J,K,SHIFT,DN
    INTEGER,DIMENSION(MAXRANK)::POSINDEX,PJPOSINDEX
    ! ----------
    ! BEGIN CODE
    ! ----------
    POS(0)=0
    DO I=1,4
       POS(I)=I
    ENDDO
    SHIFT=4
    DO J=2,MAXRANK
       DN=(J+3)*(J+2)*(J+1)/6
       POSINDEX(1:MAXRANK)=0
       PJPOSINDEX(1:MAXRANK)=0
       DO I=1,DN
          IF(I.GT.1)CALL NEXTINDEX(J,POSINDEX)
          CALL CONVERT_PJPOSINDEX(J,POSINDEX,PJPOSINDEX)
          K=DN-QPOLYPOS(J,PJPOSINDEX)+1+SHIFT
          !POS(K)=I+SHIFT
          POS(I+SHIFT)=K
       ENDDO
       SHIFT=SHIFT+DN
    ENDDO

  END SUBROUTINE ASSIGN_PJPOS

  SUBROUTINE NEXTINDEX(RANK,POSINDEX)
    !
    ! CALL FOR THE NEXT INDEX
    !
    IMPLICIT NONE
    !
    ! CONSTANTS
    !
    INTEGER,PARAMETER::MAXRANK=10
    !
    ! ARGUMENTS
    !
    INTEGER,INTENT(IN)::RANK
    INTEGER,DIMENSION(MAXRANK),INTENT(INOUT)::POSINDEX
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I
    ! ----------
    ! BEGIN CODE
    ! ----------
    DO I=1,RANK
       POSINDEX(I)=POSINDEX(I)+1
       IF(POSINDEX(I).GT.3)THEN
          POSINDEX(I)=0
          IF(I.EQ.RANK)THEN
             RETURN
          ENDIF
       ELSE
          IF(I.GT.1)THEN
             POSINDEX(1:I-1)=POSINDEX(I)
          ENDIF
          RETURN
       ENDIF
    ENDDO

  END SUBROUTINE NEXTINDEX

  SUBROUTINE CONVERT_PJPOSINDEX(RANK,POSINDEX,PJPOSINDEX)
    !
    ! CONVERT POSINDEX TO PJPOSINDEX
    !
    IMPLICIT NONE
    !
    ! CONSTANTS
    !
    INTEGER,PARAMETER::MAXRANK=10
    !
    ! ARGUMENTS
    !
    INTEGER,INTENT(IN)::RANK
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
    INTEGER,DIMENSION(MAXRANK),INTENT(OUT)::PJPOSINDEX
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I
    ! ----------
    ! BEGIN CODE
    ! ----------
    DO I=1,RANK
       PJPOSINDEX(RANK+1-I)=3-POSINDEX(I)
    ENDDO
    RETURN
  END SUBROUTINE CONVERT_PJPOSINDEX

  FUNCTION QPOLYPOS(RANK,POSINDEX)
    !
    ! COMPUTATION THE RELATIVE POSITION OF INDEX WITH RANK
    ! IN THE *OLD* MADLOOP CONVENTION
    !
    IMPLICIT NONE
    !
    ! CONSTANTS
    !
    INTEGER,PARAMETER::MAXRANK=10
    !
    ! ARGUMENTS
    !
    INTEGER,INTENT(IN)::RANK
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
    INTEGER::QPOLYPOS
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I,J,IMIN
    ! ----------
    ! BEGIN CODE
    ! ----------

    IF(RANK.EQ.0)THEN
       QPOLYPOS=1
       RETURN
    ENDIF

    IF(RANK.EQ.1)THEN
       QPOLYPOS=POSINDEX(1)+1
       RETURN
    ENDIF

    QPOLYPOS=POSINDEX(1)-POSINDEX(2)+1
    DO I=2,RANK
       IF(I.EQ.RANK)THEN
          IMIN=0
       ELSE
          IMIN=POSINDEX(I+1)
       ENDIF
       DO J=IMIN,POSINDEX(I)-1
          QPOLYPOS=QPOLYPOS+QPOLYNUMBER(J,I-1)
       ENDDO
    ENDDO
    RETURN
  END FUNCTION QPOLYPOS

  FUNCTION NEWQPOLYPOS(RANK,POSINDEX)
    !
    ! COMPUTATION THE RELATIVE POSITION OF INDEX WITH RANK
    ! IN THE *NEW* MADLOOP CONVENTION
    !
    IMPLICIT NONE
    !
    ! CONSTANTS
    !
    INTEGER,PARAMETER::MAXRANK=10
    !
    ! ARGUMENTS
    !
    INTEGER,INTENT(IN)::RANK
    INTEGER,DIMENSION(MAXRANK),INTENT(IN)::POSINDEX
    INTEGER::NEWQPOLYPOS
    !
    ! LOCAL VARIABLES
    !
    INTEGER::I,J,IMIN
    ! ----------
    ! BEGIN CODE
    ! ----------

    IF(RANK.EQ.0)THEN
       NEWQPOLYPOS=1
       RETURN
    ENDIF

    IF(RANK.EQ.1)THEN
       NEWQPOLYPOS=POSINDEX(1)+1
       RETURN
    ENDIF

    NEWQPOLYPOS=1
    DO I=1,RANK
       IF(POSINDEX(RANK-I+1).EQ.0)THEN
          IMIN=0
       ELSE
          ! Eq.(3.4.6) in Valentin's PhD Thesis
          IMIN=factorial(POSINDEX(RANK-I+1)+I-1)/factorial(I)&
               /factorial(POSINDEX(RANK-I+1)-1)
       ENDIF
       NEWQPOLYPOS=NEWQPOLYPOS+IMIN
    ENDDO
    RETURN
  END FUNCTION NEWQPOLYPOS

  FUNCTION QPOLYNUMBER(I,RANK)
    !
    ! THE INDEPENDENT NUMBER OF Q POLY WITH \MU=I,...,3 AND RANK
    !
    IMPLICIT NONE
    !
    ! CONSTANTS
    !
    !
    ! ARGUMENTS
    !
    INTEGER,INTENT(IN)::I,RANK
    INTEGER::QPOLYNUMBER
    !
    ! LOCAL VARIABLES
    !
    ! ----------
    ! BEGIN CODE
    ! ----------
    SELECT CASE(I)
    CASE(0)
       QPOLYNUMBER=(3+RANK)*(2+RANK)*(1+RANK)/6
    CASE(1)
       QPOLYNUMBER=(2+RANK)*(1+RANK)/2
    CASE(2)
       QPOLYNUMBER=(1+RANK)
    CASE(3)
       QPOLYNUMBER=1
    CASE DEFAULT
       STOP 'I must be >= 0 and <=3 in QPOLYNUMBER.'
    END SELECT
    RETURN
  END FUNCTION QPOLYNUMBER

  FUNCTION POS2RANK(ipos)
    IMPLICIT NONE
    INTEGER,INTENT(IN)::ipos
    INTEGER::POS2RANK
    LOGICAL::INIT=.TRUE.
    ! NOW ITS MAX RANK IS 10
    INTEGER,PARAMETER::MAXRANK=10
    INTEGER,DIMENSION(0:MAXRANK)::IRANGE
    INTEGER::I
    SAVE IRANGE,INIT
    IF(INIT)THEN
       IRANGE(0)=0
       DO I=1,MAXRANK
          IRANGE(I)=IRANGE(I-1)+(I+3)*(I+2)*(I+1)/6
       ENDDO
       INIT=.FALSE.
    ENDIF
    IF(ipos.EQ.0)THEN
       POS2RANK=0
       RETURN
    ENDIF
    DO I=1,MAXRANK
       IF(ipos.LE.IRANGE(I))THEN
          POS2RANK=I
          RETURN
       ENDIF
    ENDDO
    WRITE(*,*)"ERROR in POS2RANK,ipos=",ipos
    STOP
  END FUNCTION POS2RANK

  FUNCTION RANK_NUM(MAXRANK)
    IMPLICIT NONE
    INTEGER::RANK_NUM
    INTEGER,INTENT(IN)::MAXRANK
    INTEGER::I
    RANK_NUM=0
    IF(MAXRANK.LT.0)RETURN
    DO I=0,MAXRANK
       RANK_NUM=RANK_NUM+(I+3)*(I+2)*(I+1)/6
    ENDDO
    RETURN
  END FUNCTION RANK_NUM

  SUBROUTINE SHIFT_MOM(MOM,NCOEFS,TICOEFS_SAVE,TICOEFS)
    ! SHIFT MOMENTUM MOM
    IMPLICIT NONE
    !
    ! CONSTANS
    !
    ! MAX RANK SET AS 10
    ! Sum((3+r)*(2+r)*(1+r)/6,{r,0,10})=1001 --> LOOPMAXCOEFS_IREGI
    ! MAX OF (i1+1)*(i2+1)*(i3+1)*(i4+1) with i1+i2+i3+i4<=MAXRANK = 144 --> MAXSPLIT
    INTEGER,PARAMETER::LOOPMAXCOEFS_IREGI=1001,MAXRANK=10,MAXSPLIT=144
    REAL(KIND(1d0)),DIMENSION(0:3),INTENT(IN)::MOM
    INTEGER,INTENT(IN)::NCOEFS
    COMPLEX(KIND(1d0)),DIMENSION(0:NCOEFS-1,1:4),INTENT(IN)::TICOEFS_SAVE
    COMPLEX(KIND(1d0)),DIMENSION(0:NCOEFS-1,1:4),INTENT(OUT)::TICOEFS
    LOGICAL::INIT=.TRUE.
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1)::NTERM
    REAL(KIND(1d0)),DIMENSION(0:LOOPMAXCOEFS_IREGI-1,MAXSPLIT)::SPLIT_FACTOR
    ! i=0-3 -> number of lorentz indices =i
    ! i=4 -> the position of TICOEFS that should be multiplied
    ! e.g. (q+MOM)^mu1(q+MOM)^mu2...(q+MOM)^muk
    ! it can be splitted into q^mu1...q^muk+q^mu1...q^mui Mom^mui+1 ...q^muk ...
    ! for some specific term configurateion mu1=0,mu2=0,...muk=3
    ! NTERM stores the total number of terms
    ! SPLIT_FACTOR is the factor for each term
    ! SPLIT_INFO(J,K,0-3) is the number of 0,1,2,3 with MOM^mu
    ! SPLIT_INFO(J,K,4) is the corresponding q^mu (remaining lorentz indices)
    INTEGER,DIMENSION(0:LOOPMAXCOEFS_IREGI-1,MAXSPLIT,0:4)::SPLIT_INFO
    SAVE INIT,NTERM,SPLIT_FACTOR,SPLIT_INFO
    INTEGER::I,J,K,L,NINDEX,r,incr
    INTEGER::i1,i2,i3,i4,nsplit,itit,rte,ND
    REAL(KIND(1d0))::temp
    INTEGER,DIMENSION(MAXRANK)::POSINDEX
    IF(INIT)THEN
       incr=0
       ND=0
       DO r=0,MAXRANK
          DO I=r,0,-1
             DO J=r-I,0,-1
                DO K=r-I-J,0,-1
                   L=r-I-J-K
                   ! exhaust all of the subsets
                   nsplit=0
                   NTERM(incr)=(I+1)*(J+1)*(K+1)*(L+1)
                   DO i1=0,I
                      !DO itit=1,i1
                      !   POSINDEX(itit)=0
                      !ENDDO
                      DO i2=0,J
                         !DO itit=i1+1,i1+i2
                         !   POSINDEX(itit)=1
                         !ENDDO
                         DO i3=0,K
                            !DO itit=i1+i2+1,i1+i2+i3
                            !   POSINDEX(itit)=2
                            !ENDDO
                            DO i4=0,L
                               !DO itit=i1+i2+i3+1,i1+i2+i3+i4
                               !   POSINDEX(itit)=3
                               !ENDDO
                               nsplit=nsplit+1
                               rte=i1+i2+i3+i4
                               DO itit=1,i1
                                  POSINDEX(rte-itit+1)=0
                               ENDDO
                               DO itit=i1+1,i1+i2
                                  POSINDEX(rte-itit+1)=1
                               ENDDO
                               DO itit=i1+i2+1,i1+i2+i3
                                  POSINDEX(rte-itit+1)=2
                               ENDDO
                               DO itit=i1+i2+i3+1,i1+i2+i3+i4
                                  POSINDEX(rte-itit+1)=3
                               ENDDO
                               SPLIT_INFO(incr,nsplit,4)=QPOLYPOS(rte,POSINDEX)+RANK_NUM(rte-1)-1
                               SPLIT_INFO(incr,nsplit,0)=I-i1
                               SPLIT_INFO(incr,nsplit,1)=J-i2
                               SPLIT_INFO(incr,nsplit,2)=K-i3
                               SPLIT_INFO(incr,nsplit,3)=L-i4
                               SPLIT_FACTOR(incr,nsplit)=DBLE(factorial(I))&
                                    /DBLE(factorial(i1)*factorial(I-i1))
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
                                    DBLE(factorial(J))/DBLE(factorial(i2)*factorial(J-i2))
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
                                    DBLE(factorial(K))/DBLE(factorial(i3)*factorial(K-i3))
                               SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
                                    DBLE(factorial(L))/DBLE(factorial(i4)*factorial(L-i4))
                               IF(.NOT.ML5_CONVENTION)THEN
                                  SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)*&
                                       DBLE(factorial(r))/DBLE(factorial(I)*factorial(J))&
                                       /DBLE(factorial(K)*factorial(L))
                                  SPLIT_FACTOR(incr,nsplit)=SPLIT_FACTOR(incr,nsplit)/&
                                       DBLE(factorial(rte))*DBLE(factorial(i1))*DBLE(factorial(i2))*&
                                       DBLE(factorial(i3))*DBLE(factorial(i4))
                               ENDIF
                            ENDDO
                         ENDDO
                      ENDDO
                   ENDDO
                   incr=incr+1
                ENDDO
             ENDDO
          ENDDO
          ND=ND+(3+r)*(2+r)*(1+r)/6
       ENDDO
       INIT=.FALSE.
    ENDIF
    DO I=0,NCOEFS-1
       TICOEFS(I,1:4)=DCMPLX(0d0)
       DO J=1,NTERM(I)
          temp=SPLIT_FACTOR(I,J)
          DO K=0,3
             NINDEX=SPLIT_INFO(I,J,K)
             IF(NINDEX.GT.0)THEN
                temp=temp*MOM(K)**NINDEX
             ENDIF
          ENDDO
          TICOEFS(I,1:4)=TICOEFS(I,1:4)&
               +TICOEFS_SAVE(SPLIT_INFO(I,J,4),1:4)*temp
       ENDDO
    ENDDO
    RETURN
  END SUBROUTINE SHIFT_MOM

  FUNCTION xiarray_arg1(i,n)
    IMPLICIT NONE
    INTEGER,INTENT(IN)::i,n
    INTEGER::xiarray_arg1
    INTEGER::imax,imin,j
    !INTEGER(KIND=LINT)::itemp
    IF(i+n-1.GT.12)THEN
       imax=MAX(n-1,i)
       imin=MIN(n-1,i)
       xiarray_arg1=1
       DO j=imax+1,i+n-1
          xiarray_arg1=xiarray_arg1*j
       ENDDO
       IF(imin.LE.12)THEN
          xiarray_arg1=xiarray_arg1/factorial(imin)
       ELSE
          DO j=1,imin
             xiarray_arg1=xiarray_arg1/j
          ENDDO
       ENDIF
    ELSE
       xiarray_arg1=factorial(i+n-1)/factorial(n-1)/factorial(i)
    ENDIF
    RETURN
  END FUNCTION xiarray_arg1
END MODULE FUNLIB


text:
C    @(#)pieqiv.f	20.4 2/13/96
      subroutine pieqiv (p, y, error)
      complex * 16  y(2,2)
      integer p, error
c
c     This subroutine computes the 2-port admittance matrix y
c     for brnch(1,aq)
c
c	Modified to handle double precision complex numbers.
c
      include 'ipfinc/parametr.inc'

      include 'ipfinc/branch.inc'
c	Global variables used:
c		brnch, kx, ky, brtype, brnch_ptr
      include 'ipfinc/bus.inc'
c	Global variables used:
c		base
      include 'ipfinc/prt.inc'
c	Global variables used:
c		errbuf

c     set up record type code parameters...
      include 'ipfinc/brtype.inc'

      common /is_batch / is_batch

      complex * 16 cx
      real r, x, z, tk, tm
      integer q, aq, k1, k2, ltype
      character xbuf * 120
c
      error = 0
      k1 = kx(p)
      k2 = ky(p)
      ltype = brtype(p)
      q = brnch_ptr(p)
      aq = iabs(q)

      if (ltype.eq.BRTYP_L .or. ltype.eq.BRTYP_T .or.
     &    ltype.eq.BRTYP_TP .or. ltype.eq.BRTYP_E) then
c
c        compute transfer impedance
c
         r=brnch(5,aq)
         x=brnch(6,aq)
         z=r**2 + x**2
         if (z .eq. 0.0) then
            write (errbuf(1),110)
  110       format(' Branch has zero transfer impedance.')
            errbuf(2) = ' '
            call bcdbrn(p,xbuf)
            write (errbuf(3), 120) xbuf(1:80)
  120       format(13x,'(',a80,')')
            if (is_batch .eq. 0) then
               call prterx ('E',3)
            else
               call prterx ('F',3)
            endif
            error = 1
            z=1.0
         endif

         if (ltype .eq. BRTYP_L) then
c
c           "L" branch
c
            y(1,2)=-dcmplx (r/z, -x/z)
            y(2,1) = y(1,2)
            y(1,1)=dcmplx(dble(brnch(7,aq)),dble(brnch(8,aq))) - y(1,2)
            y(2,2)=y(1,1)

         else if (ltype .eq. BRTYP_T) then
c
c           "T" branch
c
            if (q .gt. 0) then
               tk=brnch(9,aq)/base(k1)
               tm=brnch(10,aq)/base(k2)
            else
               tk=brnch(10,aq)/base(k1)
               tm=brnch(9,aq)/base(k2)
            endif
            if (amax1(tk,tm).ge.1.40 .or. amin1(tk,tm).le.0.60) then
               if (q .gt. 0) then
                  tk=brnch(10,aq)/base(k1)
                  tm=brnch(9,aq)/base(k2)
               else
                  tk=brnch(9,aq)/base(k1)
                  tm=brnch(10,aq)/base(k2)
               endif
               if (amax1(tk,tm).lt.1.40 .and. amin1(tk,tm).gt.0.60)
     &            then
                  write(errbuf(1),130)
  130             format(' Backwards transformer taps are reversed.')
                  errbuf(2) = ' '
                  call bcdbrn(p,xbuf)
                  write (errbuf(3),120) xbuf(1:80)
                  temp = brnch(9,aq)
                  brnch(9,aq) = brnch(10,aq)
                  brnch(10,aq) = temp
                  call bcdbrn(p,xbuf)
                  write (errbuf(4),120) xbuf(1:80)
                  call prterx ('W',4)
               else
                  write (errbuf(1),140)
  140             format(' Transformer taps are beyond 40% of base kv.')
                  errbuf(2) = ' '
                  call bcdbrn(p,xbuf)
                  write (errbuf(3),120) xbuf(1:80)
                  if (is_batch .eq. 0) then
                     call prterx ('E',3)
                  else
                     call prterx ('F',3)
                  endif
                  if (tk.eq.0.0) tk=1.0
                  if (tm.eq.0.0) tm=1.0
               endif
            endif
            tkm=1.0/(tk*tm)
            y(1,1)=dcmplx(0.5*tkm*brnch(7,aq),-0.5*tkm*abs(brnch(8,aq)))
            y(2,2)=y(1,1)
            y(1,2)=-dcmplx(r/(z*tk*tm),-x/(z*tk*tm))
            y(2,1)=y(1,2)
            y(1,1)=y(1,1)+dcmplx (r/(z*tk*tk),-x/(z*tk*tk))
            y(2,2)=y(2,2)+dcmplx (r/(z*tm*tm), -x/(z*tm*tm))

         else if (ltype .eq. BRTYP_TP) then
c
c           "TP" branch
c
            if (q .gt. 0) then
               tk=brnch(9,aq)
               tm=brnch(10,aq)/base(k2)
            else
               tk=-brnch(9,aq)
               tm=brnch(10,aq)/base(k1)
            endif
            if (tm .eq. 0.0) tm=1.0
            if (tm .ge. 1.40 .or. tm .le. 0.60) then
               write (errbuf(1),150)
  150          format(' Phase shifter tap2 is beyond 40% of base kv.')
               errbuf(2) = ' '
               call bcdbrn(p,xbuf)
               write (errbuf(3),120) xbuf(1:80)
               call prterx ('W',3)
               tm=1.0
            endif
            if (tk .ge. 120.1 .or. tk .le. -120.1) then
               write (errbuf(1), 160)
  160          format(' Phase shifter has excessive angle ',
     &                '(> 120 degrees)')
               errbuf(2) = ' '
               call bcdbrn(p,xbuf)
               write (errbuf(3),120) xbuf(1:80)
               call prterx ('W',3)
            endif
            angle=0.0174532925*tk
            y(1,2)=dcmplx(0.5*brnch(7,aq),-0.5*abs(brnch(8,aq)))
            y(1,1)=y(1,2)+dcmplx(r/z,-x/z)
            y(2,2)=y(1,2)+dcmplx(r/z,-x/z)/cmplx(tm**2,0.0)
            y(1,2)=-dcmplx(cos(angle),
     &                    sin(angle))*cmplx(r/(z*tm),-x/(z*tm))
            y(2,1)=-dcmplx(cos(angle),
     &                   -sin(angle))*cmplx(r/(z*tm),-x/(z*tm))
            if (q .lt. 0) then
               cx=y(2,2)
               y(2,2)=y(1,1)
               y(1,1)=cx
            endif

         else
c
c           "E" branch
c
            y(1,2)= -dcmplx (r/z, -x/z)
            y(2,1)=y(1,2)
            if (q .gt. 0) then
               y(1,1)=dcmplx(brnch(7,aq),brnch(8,aq))-y(1,2)
               y(2,2)=dcmplx(brnch(9,aq),brnch(10,aq))-y(2,1)
            else
               y(2,2)=dcmplx(brnch(7,aq),brnch(8,aq))-y(2,1)
               y(1,1)=dcmplx(brnch(9,aq),brnch(10,aq))-y(1,2)
            endif

         endif

      else if (ltype.eq.BRTYP_PEQ) then
c
c        "pi-equivalent" branch
c
         if (q .gt. 0) then
            y(1,1) = dcmplx(brnch(4,aq),brnch(5,aq))
            y(1,2) = dcmplx(brnch(6,aq),brnch(7,aq))
            y(2,1) = dcmplx(brnch(8,aq),brnch(9,aq))
            y(2,2) = dcmplx(brnch(10,aq),brnch(11,aq))
         else
            y(2,2) = dcmplx(brnch(4,aq),brnch(5,aq))
            y(2,1) = dcmplx(brnch(6,aq),brnch(7,aq))
            y(1,2) = dcmplx(brnch(8,aq),brnch(9,aq))
            y(1,1) = dcmplx(brnch(10,aq),brnch(11,aq))
         endif

      else
c
c        error-illegal branch type calling "pieqiv"
c
         write (errbuf(1), 170)
  170    format(' Illegal branch type processed by "PIEQIV"')
         errbuf(2) = ' '
         call bcdbrn(p,xbuf)
         write (errbuf(3),120) xbuf(1:80)
         call prterx ('W',3)
         error = 1
      endif

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
module dualClass
use adMaths
implicit none

	type, public :: dual
		real(kind=8) :: rp
		real(kind=8), dimension(dualSize) :: ip   ! Real and infinitesimal part
	contains
		procedure, pass :: init
	end type dual

! Operator overload
	interface operator(+)
		module procedure additionDualDual, additionDualReal, additionRealDual
	end interface

	interface operator(-)
		module procedure substractionDualDual, substractionDualReal, substractionRealDual, signChangeDual
	end interface

	interface operator(*)
		module procedure productDualDual, productDualReal, productRealDual
	end interface

	interface operator(/)
		module procedure divisionDualDual, divisionDualReal, divisionRealDual
	end interface

	interface operator(**)
		module procedure powerDual, powerIntDual
	end interface

	interface exp
		module procedure expDual
	end interface

	interface sin
		module procedure sinDual
	end interface

	interface cos
		module procedure cosDual
	end interface

	interface log
		module procedure logDual
	end interface

	interface tanh
		module procedure tanhDual
	end interface

	interface sqrt
		module procedure sqrtDual
	end interface

	interface sum
		module procedure sumDual
	end interface

	interface lnGamma
		module procedure lnGammaDual
	end interface

	interface voigt
		module procedure voigtDual
	end interface

	contains

!---------------------------
! Creation of the dual number
! value and whichIndependent are optional arguments
! value is used to initialize the value of the number
! whichIndependent is used to set independent variables. If not present, the variable is dependent
!---------------------------
	elemental subroutine init(this, value, whichIndependent)
	class(dual), intent(inout) :: this
	real(kind=8), intent(in), optional :: value
	integer, intent(in), optional :: whichIndependent

		this%ip = 0.d0

		this%rp = 0.d0
		if (present(value)) then
			this%rp = value
		endif

		this%ip = 0.d0
		if (present(whichIndependent)) then
			this%ip(whichIndependent) = 1.d0
		endif

	end subroutine init

! Functions for defining the basic operations

! ---------------
! MULTIPLICATION
!----------------
	elemental function productDualDual(a,b) result(c)
		type (dual), intent(in) :: a, b
		type (dual) :: c
		c%rp = a%rp * b%rp
		c%ip = a%rp * b%ip + b%rp * a%ip
	end function productDualDual

	elemental function productDualReal(a,b) result(c)
		type (dual), intent(in) :: a
		real(kind=8), intent(in) :: b
		type (dual) :: c
		c%rp = a%rp * b
		c%ip = a%ip * b
	end function productDualReal

	elemental function productRealDual(a,b) result(c)
		real(kind=8), intent(in) :: a
		type (dual), intent(in) :: b
		type (dual) :: c
		c%rp = b%rp * a
		c%ip = b%ip * a
	end function productRealDual

! ---------------
! DIVISION
!----------------
	elemental function divisionDualDual(a,b) result(c)
		type (dual), intent(in) :: a, b
		type (dual) :: c
		c%rp = a%rp / b%rp
		c%ip = (b%rp * a%ip - a%rp * b%ip) / b%rp**2
	end function divisionDualDual

	elemental function divisionDualReal(a,b) result(c)
		type (dual), intent(in) :: a
		real(kind=8), intent(in) :: b
		type (dual) :: c
		c%rp = a%rp / b
		c%ip = a%ip / b
	end function divisionDualReal

	elemental function divisionRealDual(a,b) result(c)
		real(kind=8), intent(in) :: a
		type (dual), intent(in) :: b
		type (dual) :: c
		c%rp = a / b%rp
		c%ip = -a * b%ip / b%rp**2
	end function divisionRealDual

!---------------
! ADDITION
!---------------
	elemental function additionDualDual(a,b) result(c)
		type (dual), intent(in) :: a, b
		type (dual) :: c
		c%rp = a%rp + b%rp
		c%ip = a%ip + b%ip
	end function additionDualDual

	elemental function additionDualReal(a,b) result(c)
		type (dual), intent(in) :: a
		real(kind=8), intent(in) :: b
		type (dual) :: c
		c%rp = a%rp + b
		c%ip = a%ip
	end function additionDualReal

	elemental function additionRealDual(a,b) result(c)
		real(kind=8), intent(in) :: a
		type (dual), intent(in) :: b
		type (dual) :: c
		c%rp = b%rp + a
		c%ip = b%ip
	end function additionRealDual

!---------------
! SUBSTRACTION
!---------------
	elemental function substractionDualDual(a,b) result(c)
		type (dual), intent(in) :: a, b
		type (dual) :: c
		c%rp = a%rp - b%rp
		c%ip = a%ip - b%ip
	end function substractionDualDual

	elemental function substractionDualReal(a,b) result(c)
		type (dual), intent(in) :: a
		real(kind=8), intent(in) :: b
		type (dual) :: c
		c%rp = a%rp - b
		c%ip = a%ip
	end function substractionDualReal

	elemental function substractionRealDual(a,b) result(c)
		real(kind=8), intent(in) :: a
		type (dual), intent(in) :: b
		type (dual) :: c
		c%rp = b%rp - a
		c%ip = b%ip
	end function substractionRealDual

	elemental function signChangeDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		c%rp = -a%rp
		c%ip = -a%ip
	end function signChangeDual

!---------------
! POWER
!---------------
	elemental function powerDual(a,b) result(c)
		type (dual), intent(in) :: a
		real(kind=8), intent(in) :: b
		type (dual) :: c
		c%rp = a%rp**b
		c%ip = b * a%rp**(b-1.d0) * a%ip
	end function powerDual

	elemental function powerIntDual(a,b) result(c)
		type (dual), intent(in) :: a
		integer, intent(in) :: b
		type (dual) :: c
		c%rp = a%rp**b
		c%ip = b * a%rp**(b-1.d0) * a%ip
	end function powerIntDual

!---------------
! EXP
!---------------
	elemental function expDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		c%rp = exp(a%rp)
		c%ip = a%ip * exp(a%rp)
	end function expDual

!---------------
! SIN
!---------------
	elemental function sinDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		c%rp = sin(a%rp)
		c%ip = a%ip * cos(a%rp)
	end function sinDual

!---------------
! COS
!---------------
	elemental function cosDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		c%rp = cos(a%rp)
		c%ip = -a%ip * sin(a%rp)
	end function cosDual

!---------------
! LN
!---------------
	elemental function logDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		c%rp = log(a%rp)
		c%ip = a%ip / a%rp
	end function logDual

!---------------
! TANH
!---------------
	elemental function tanhDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		real(kind=8) :: t
		c%rp = tanh(a%rp)
		c%ip = 1.d0/cosh(a%rp)**2 * a%ip
	end function tanhDual

!---------------
! SQRT
!---------------
	elemental function sqrtDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c

		c%rp = sqrt(a%rp)
		c%ip = 0.5d0/sqrt(a%rp) * a%ip
	end function sqrtDual

!---------------
! ADDITION OF A VECTOR
!---------------
	function sumDual(a) result(c)
		type (dual), intent(in), dimension(:) :: a
		type (dual), dimension(size(a)) :: c
		integer :: i
		c%rp = sum(a%rp)
 		do i = 1, size(a(1)%ip)
 			c%ip(i) = sum(a%ip(i))
 		enddo
	end function sumDual

!---------------
! lnGamma
!---------------
	function lnGammaDual(a) result(c)
		type (dual), intent(in) :: a
		type (dual) :: c
		integer :: ierr

		c%rp = alngam(a%rp, ierr)
		c%ip = digama(a%rp, ierr) * a%ip

	end function lnGammaDual

!---------------
! Voigt profile
!---------------
	function voigtDual(a, v) result(c)
		type (dual), intent(in) :: a, v
		type (dual) :: c
		integer :: ierr
		real(kind=8) :: H(2), ga, gv, PI = 3.14159265359d0

			H = voigtAutomaticDiff(a%rp, v%rp)
			ga = 2.d0 * (-1.d0 / sqrt(PI) + a%rp * H(1) + v%rp * H(2))
			gv = 2.d0 * (-v%rp * H(1) + a%rp * H(2))

			c%rp = H(1)
			c%ip = ga * a%ip + gv * v%ip

	end function voigtDual

end module dualClass

max_stars_repo_name: aasensio/pyiacsun

max_issues_repo_name: aasensio/pyiacsun

max_forks_repo_name: aasensio/pyiacsun


text:
** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
** See https://llvm.org/LICENSE.txt for license information.
** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

*   adjustable POINTER object and ALLOCATE statement

	parameter (N=1)
	integer result(N)
	integer expect(N)

  	maxdat=1000
	call glim3(maxdat, result)

	data expect /1/
	call check(result, expect, N)

	end
	subroutine glim3(imxdat, ires)
	double precision sspmat
	integer imxdat
	pointer (l_sspmat, sspmat(imxdat))

c  ensure ilms are written for sspmat's array descriptor

	data kount/0/, l2dir/0/	 ! cause ilms to be erased

	maxdat=imxdat
	allocate(sspmat, stat=jallo)
	if ( jallo .eq. 0 ) goto 100
	ires = 2
	goto 9999
 100    do 18 i=1,maxdat
	   sspmat(i)=0.d0
  18    continue
	ires = 1
	deallocate(sspmat)
 9999   continue
	return
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


! written by jxzou at 20220220: generate .47 file from Gaussian .fch(k) file

program main
 implicit none
 integer :: i
 character(len=240) :: fchname

 i = iargc()
 if(i /= 1) then
  write(6,'(/,A)') ' ERROR in program fch247: wrong command line argument!'
  write(6,'(A,/)') ' Example: fch247 a.fch'
  stop
 end if

 call getarg(1, fchname)
 call require_file_exist(fchname)
 call fch247(fchname)
 stop
end program main

subroutine fch247(fchname)
 use fch_content
 implicit none
 integer :: i, j, k, m, nhigh, fid
 integer, allocatable :: eff_nuc_charge(:), cen_assign(:)
 real(kind=8), allocatable :: S(:,:),F(:,:),dm_a(:,:),dm_b(:,:), D(:,:,:)
 character(len=1), parameter :: l_str(0:5) = ['S','P','D','F','G','H']
 character(len=240) :: f47
 character(len=240), intent(in) :: fchname
 logical :: uhf

 call check_uhf_in_fch(fchname, uhf)
 call read_fch(fchname, uhf)

 i = index(fchname, '.fch')
 f47 = fchname(1:i-1)//'.47'
 open(newunit=fid,file=TRIM(f47),status='replace')
 write(fid,'(A,I0,A,I0)',advance='no') ' $GENNBO  NATOMS=',natom,'  NBAS=',nbf
 if(uhf) write(fid,'(A)',advance='no') '  OPEN'
 write(fid,'(A)') '  UPPER  BODM  $END'
 write(fid,'(A)') ' $NBO  $END'
 write(fid,'(A)') ' $COORD'
 write(fid,'(A)') ' Generated by utility fch247 in MOKIT'

 ! calculate effective nuclear chages (remember to minus ECP core, if any)
 allocate(eff_nuc_charge(natom), source=0)
 if(allocated(RNFroz)) then
  forall(i = 1:natom) eff_nuc_charge(i) = ielem(i) - INT(RNFroz(i))
 else
  forall(i = 1:natom) eff_nuc_charge(i) = ielem(i)
 end if
 do i = 1, natom, 1
  write(fid,'(1X,2I6,3F15.6)') ielem(i),eff_nuc_charge(i),coor(:,i)
 end do ! for i
 deallocate(eff_nuc_charge)
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $BASIS'

 allocate(cen_assign(nbf))
 k = 0
 do i = 1, ncontr, 1
  j = shltyp2nbf(shell_type(i))
  cen_assign(k+1:k+j) = shell2atom_map(i)
  k = k + j
 end do ! for i
 deallocate(shell2atom_map)
 write(fid,'(A,17I7)') '  CENTER =',cen_assign(1:min(17,nbf))
 if(nbf > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,nbf)

 ! Next, we use the array cen_assign to store LABEL
 cen_assign = 0; k = 0
 do i = 1, ncontr, 1
  j = shell_type(i)
  select case(j)
  case(0) ! S
   cen_assign(k+1) = 1
  case(1) ! P
   forall(m = 1:3) cen_assign(k+m) = 100 + m
  case(2) ! 6D
   cen_assign(k+1:k+6) = [201,204,206,202,203,205]
  case(3) ! 10F
   cen_assign(k+1:k+10) = [301,307,310,304,302,303,306,309,308,305]
  case(4) ! 15G
   forall(m = 1:15) cen_assign(k+m) = 416 - m
  case(5) ! 21H
   forall(m = 1:21) cen_assign(k+m) = 522 - m
  case(-1) ! L/SP
   cen_assign(k+1:k+4) = [1,101,102,103]
  case(-2) ! 5D
   cen_assign(k+1:k+5) = [255,252,253,254,251]
  case(-3) ! 7F
   forall(m = 1:7) cen_assign(k+m) = 350 + m
  case(-4) ! 9G
   forall(m = 1:9) cen_assign(k+m) = 450 + m
  case(-5) ! 11H
   forall(m = 1:11) cen_assign(k+m) = 550 + m
  case default
   stop
  end select
  k = k + shltyp2nbf(j)
 end do ! for i

 write(fid,'(A,17I7)') '   LABEL =',cen_assign(1:min(17,nbf))
 if(nbf > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,nbf)
 deallocate(cen_assign)
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $CONTRACT'
 write(fid,'(A,I7)') '  NSHELL =',ncontr
 write(fid,'(A,I7)') '    NEXP =',nprim

 allocate(cen_assign(ncontr))
 forall(i = 1:ncontr) cen_assign(i) = shltyp2nbf(shell_type(i))
 write(fid,'(A,17I7)') '   NCOMP =',cen_assign(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,ncontr)

 write(fid,'(A,17I7)') '   NPRIM =',prim_per_shell(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (prim_per_shell(i),i=18,ncontr)

 cen_assign(1) = 1
 do i = 2, ncontr, 1
  cen_assign(i) = cen_assign(i-1) + prim_per_shell(i-1)
 end do ! for i

 write(fid,'(A,17I7)') '    NPTR =',cen_assign(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,ncontr)
 deallocate(cen_assign)

 write(fid,'(A,4E16.7)') '     EXP =',prim_exp(1:min(4,nprim))
 if(nprim > 4) write(fid,'(10X,4E16.7)') (prim_exp(i),i=5,nprim)

 ! (1) multiply contraction coeff by normalization factor
 ! (2) use F(:,:) to store CS, CP, CD, CF, CG, CH.
 i = MAXVAL(shell_type)
 if(i < 0) i = -i
 j = MINVAL(shell_type)
 if(j < 0) j = -j
 nhigh = MAX(i,j)
 allocate(F(nprim,0:nhigh),source=0d0)
 m = 0
 do i = 1, ncontr, 1
  j = prim_per_shell(i); k = shell_type(i)
  if(k == -1) then ! L, SP
   call contr_coeff_multiply_norm_fac(j,0,prim_exp(m+1:m+j),contr_coeff(m+1:m+j))
   call contr_coeff_multiply_norm_fac(j,1,prim_exp(m+1:m+j),contr_coeff_sp(m+1:m+j))
   F(m+1:m+j,0) = contr_coeff(m+1:m+j)
   F(m+1:m+j,1) = contr_coeff_sp(m+1:m+j)
  else !  S/P/D/F/G/H
   call contr_coeff_multiply_norm_fac(j,k,prim_exp(m+1:m+j),contr_coeff(m+1:m+j))
   F(m+1:m+j,ABS(k)) = contr_coeff(m+1:m+j)
  end if
  m = m + j
 end do ! for i
 do i = 0, nhigh, 1
  write(fid,'(A,4E16.7)') '      C'//l_str(i)//' =',F(1:min(4,nprim),i)
  if(nprim > 4) write(fid,'(10X,4E16.7)') (F(j,i),j=5,nprim)
 end do ! for i
 deallocate(F)
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $OVERLAP'
 allocate(S(nbf,nbf))
 if(nbf == nif) then
  call solve_ovlp_from_ctsc(nbf, alpha_coeff, S)
 else
  write(6,'(/,A)') 'ERROR in subroutine fch247: nbf/=nif, linear dependence&
                  & detected.'
  write(6,'(A)') 'Computing overlap integrals not supported currently.'
  stop
 end if
 write(fid,'(2X,5E15.7)') ((S(j,i),j=1,i),i=1,nbf)
 deallocate(S)
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $DENSITY'
 if(uhf) then ! unretricted type wave function
  if(allocated(tot_dm) .and. allocated(spin_dm)) then
   dm_a = 0.5d0*(tot_dm + spin_dm) ! alpha density
   dm_b = 0.5d0*(tot_dm - spin_dm) ! beta density
  else
   allocate(dm_a(nbf,nbf), dm_b(nbf,nbf))
   dm_a = MATMUL(alpha_coeff(:,1:na), TRANSPOSE(alpha_coeff(:,1:na)))
   dm_b = MATMUL(beta_coeff(:,1:nb), TRANSPOSE(beta_coeff(:,1:nb)))
  end if
  write(fid,'(2X,5E15.7)') ((dm_a(j,i),j=1,i),i=1,nbf)
  write(fid,'(2X,5E15.7)') ((dm_b(j,i),j=1,i),i=1,nbf)
  deallocate(dm_a, dm_b)
 else ! retricted type wave function
  if(.not. allocated(tot_dm)) then
   allocate(tot_dm(nbf,nbf))
   tot_dm = 2d0*MATMUL(alpha_coeff(:,1:nb), TRANSPOSE(alpha_coeff(:,1:nb)))
   if(na > nb) tot_dm = tot_dm + MATMUL(alpha_coeff(:,nb+1:na), &
                                        TRANSPOSE(alpha_coeff(:,1:nb+1:na)))
  end if
  write(fid,'(2X,5E15.7)') ((tot_dm(j,i),j=1,i),i=1,nbf)
 end if
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $FOCK'
 allocate(F(nbf,nbf))
 call solve_fock_from_ctfc(nbf, nif, alpha_coeff, eigen_e_a, F)
 write(fid,'(2X,5E15.7)') ((F(j,i),j=1,i),i=1,nbf)
 if(uhf) then
  call solve_fock_from_ctfc(nbf, nif, beta_coeff, eigen_e_b, F)
  write(fid,'(2X,5E15.7)') ((F(j,i),j=1,i),i=1,nbf)
 end if
 deallocate(F)
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $LCAOMO'
 write(fid,'(2X,5E15.7)') alpha_coeff
 if(uhf) write(fid,'(2X,5E15.7)') beta_coeff
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $DIPOLE'
 allocate(D(nbf,nbf,3))
 call calc_dipole_using_fch(fchname, nbf, D)
 do i = 1, 3
  write(fid,'(2X,5E15.7)') ((D(k,j,i),k=1,j),j=1,nbf)
 end do ! for i
 deallocate(D)
 write(fid,'(A)') ' $END'
 close(fid)
end subroutine fch247

! multiply each contraction coefficient with corresponding normalization factor
! Note: this subroutine is different with that in file fch2wfn.f90
subroutine contr_coeff_multiply_norm_fac(n, p, prim_exp, contr_coeff)
 use fch_content, only: iout
 implicit none
 integer :: i, p0
 integer, intent(in) :: n, p
 real(kind=8) :: r2a_pi(n)
 ! 2*alpha/PI. Because variable cannot begin with a integer, so 'r' is added
 real(kind=8), parameter :: PI = 4d0*DATAN(1d0)
 real(kind=8), parameter :: sqrt_2pi = DSQRT(8d0*DATAN(1d0))
 real(kind=8), intent(in) :: prim_exp(n)
 real(kind=8), intent(inout) :: contr_coeff(n)

 if(p == -1) then
  write(iout,'(A)') 'ERROR in subroutine contr_coeff_multiply_norm_fac: p=-1.'
  write(iout,'(A)') 'You should divide L/SP into separate S/P before calling &
                    &this subroutine.'
  stop
 end if
 if(p < 0) then
  p0 = -p
 else
  p0 = p
 end if

 forall(i = 1:n)
  r2a_pi(i) = (2d0*prim_exp(i)/PI)**(0.25d0)
  contr_coeff(i) = contr_coeff(i)*(r2a_pi(i)**(2*p0+3))
 end forall
 if(p0 > 0) contr_coeff = contr_coeff*(sqrt_2pi**p0)

 select case(p0)
 case(0,1) ! do nothing
 case(2) ! D
  contr_coeff = contr_coeff/DSQRT(3d0)
 case(3) ! F
  contr_coeff = contr_coeff/DSQRT(15d0)
 case(4) ! G
  contr_coeff = contr_coeff/DSQRT(105d0)
 case(5) ! H
  contr_coeff = contr_coeff/DSQRT(945d0)
 case default
  write(iout,'(A,I0)') 'ERROR in subroutine contr_coeff_multiply_norm_fac: inva&
                       &lid p0=', p0
  stop
 end select
end subroutine contr_coeff_multiply_norm_fac