PROGRAM sp2thpdf
!***********************************************************************
! This program computes: *
! *
! density of S_i,Hi,T_i in a gaussian process i.e. *
! *
! quart wavelength (up-crossing to crest) and crest amplitude *
!
! def = 1, gives half wave period, Tc (default).
! -1, gives half wave period, Tt.
! 2, gives half wave period and wave crest amplitude (Tc,Ac).
! -2, gives half wave period and wave trough amplitude (Tt,At).
! 3, gives crest front period and wave crest amplitude (Tcf,Ac).
! -3, gives trough back period and wave trough amplitude (Ttb,At).
! 4, gives minimum of crest front/back period and wave crest
! amplitude (min(Tcf,Tcb),Ac).
! -4, gives minimum of trough front/back period and wave trough
! amplitude (min(Ttf,Ttb),At).
!***********************************************************************
!History:
! revised Per A. Brodtkorb 04.04.2000
! -
! revised Per A. Brodtkorb 23.11.99
! - fixed a bug in calculating pdf for def = +/- 4
! revised Per A. Brodtkorb 03.11.99
! - added def = +/-4
! revised Per A. Brodtkorb 23.09.99
! - minor changes to covinput
! - removed the calculation of the transformation to spec2thpdf.m
! by Igor Rychlik
use GLOBALDATA, only : rateLHD,SCIS,NSIMmax,COV,ABSEPS
use globalconst
use rind
IMPLICIT NONE
double precision, dimension(:,:),allocatable :: BIG
double precision, dimension(:,:),allocatable :: ansr
double precision, dimension(: ),allocatable :: ex
double precision, dimension(:,:),allocatable :: xc
double precision, dimension(: ),allocatable :: fxind,h
double precision, dimension(: ),allocatable :: R0,R1,R2,R3,R4
double precision ::CC,U,XddInf,XdInf,XtInf
double precision, dimension(2,6) :: a_up=0.d0,a_lo=0.d0
integer, dimension(6) :: INFIN=2
integer , dimension(: ),allocatable :: seed
integer ,dimension(7) :: indI
integer :: Nx,Nt,Nc,Nd,NI,Mb,Ntd, Ntdc
integer :: Nstart,Ntime,tn,ts,speed,ph,def,seed1,seed_size
double precision :: dT, EPSOLD ! lag spacing for covariances
LOGICAL :: init=.TRUE.
! DIGITAL:
! f90 -g2 -C -automatic -o ../wave/alpha/sp2thpdf.exe rind44.f sp2thpdf.f
! SOLARIS:
!f90 -g -O -w3 -Bdynamic -fixed -o ../wave/sol2/sp2thpdf.exe rind44.f sp2thpdf.f
! linux:
! f90 -gline -Nl126 -C -o ../exec/lnx86/sp2thpdf8.exe intmodule.f rind60.f sp2thpdf.f
! f90 -gline -Nl126 -C -o sp2thpdf.exe rind45.f sp2thpdf.f
! f90 -gline -Nl126 -C -o ../exec/lnx86/sp2thpdf3.exe adaptmodule.f krbvrcmod.f krobovmod.f rcrudemod.f rind55.f sp2thpdf.f
! HP700
!f90 -g -C -o ../exec/hp700/sp2thpdf.exe rind45.f sp2thpdf.f
!f90 -g -C +check=all +FPVZID -o ../exec/hp700/sp2thpdf.exe rind45.f sp2thpdf.f
! f90 +gprof +extend_source +Oall +Odataprefetch +Ofastaccess +Oinfo +Oprocelim -C +check=all -o ../exec/hp700/sp2thpdf.exe rind48.f sp2thpdf.f
!print *,'enter sp2thpdf'
CALL INIT_LEVELS(U,def,Ntime,Nstart,speed,SCIS,seed1,
& Nx,dT,rateLHD)
!print *,'U,def,Ntime,Nstart,NIT,speed,SCIS,seed1,Nx,dT'
!print *,U,def,Ntime,Nstart,NIT,speed,SCIS,seed1,Nx,dT
if (SCIS.GT.0) then
!allocate(COV(1:Nx))
call random_seed(SIZE=seed_size)
allocate(seed(seed_size))
call random_seed(GET=seed(1:seed_size)) ! get current seed
seed(1)=seed1 ! change seed
call random_seed(PUT=seed(1:seed_size))
deallocate(seed)
endif
CALL INITDATA(speed)
!print *,ntime,speed,u,NIT
allocate(R0(1:Ntime+1))
allocate(R1(1:Ntime+1))
allocate(R2(1:Ntime+1))
if (abs(def).GT.1) THEN
allocate(h(1:Nx))
allocate(R3(1:Ntime+1))
allocate(R4(1:Ntime+1))
CALL INIT_AMPLITUDES(h,def,Nx)
endif
CALL INIT_COVARIANCES(Ntime,def,R0,R1,R2,R3,R4)
!print *,'Nx',Nx
indI(1)=0
C ***** The bound 'infinity' is set to 10*sigma *****
XdInf=10.d0*SQRT(-R2(1))
XtInf=10.d0*SQRT(R0(1))
!print *,'XdInf,XtInf'
!print *,XdInf,XtInf
! normalizing constant
CC=TWPI*SQRT(-R0(1)/R2(1))*exp(u*u/(2.d0*R0(1)) )
allocate(ansr(1:Ntime,1:Nx))
ansr=0.d0
allocate(fxind(1:Nx))
!fxind=0.d0 this is not needed
if (abs(def).GT.1) GOTO 200
NI=4; Nd=2
Nc=2; Mb=1
Nx=1
allocate(BIG(1:Ntime+Nc,1:Ntime+Nc))
allocate(xc(1:Nc,1:Nx))
allocate(ex(1:Ntime+Nc))
ex=0.d0
xc(1,1)=u
xc(2,1)=u
! INFIN = INTEGER, array of integration limits flags: size 1 x Nb (in)
! if INFIN(I) < 0, Ith limits are (-infinity, infinity);
! if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
! if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
! if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
if (def.GT.0) then
INFIN(1:2) = 1
INFIN(3) = 0
a_up(1,1)= u+XtInf
a_lo(1,1)= u
a_up(1,2)= XdInf
a_lo(1,3)=-XdInf
else
INFIN(1:2) = 0
INFIN(3) = 1
a_up(1,1)=u
a_lo(1,1)=u-XtInf
a_lo(1,2)=-XdInf
a_up(1,3)= XdInf
endif
!print *,'Nstart',Nstart
Nstart=MAX(2,Nstart)
!print *,'Nstart',Nstart
if (ALLOCATED(COV)) then
open (unit=11, file='COV.out', STATUS='unknown')
write(11,*) 0.d0
endif
do Ntd=Nstart,Ntime
!CALL COV_INPUT2(BIG,Ntd, R0,R1,R2)
CALL COV_INPUT(BIG,Ntd,-1,R0,R1,R2,R3,R4) ! positive wave period
Nt=Ntd-Nd;
indI(2)=Nt;
indI(3)=Nt+1;
indI(4)=Ntd;
Ntdc=Ntd+Nc;
! IF (Ntd.GT.5.AND.(INIT)) THEN
! INIT=.FALSE.
! CALL INITDATA(speed)
! ENDIF
!if (SCIS.gt.1) Nj=Nt
!if (SCIS.gt.0) then
! if (SCIS.EQ.2) then
! Nj=max(Nt,0)
! else
! Nj=min(max(Nt-5, 0),0)
! endif
!endif
!Ex=0.d0
!CALL echo(BIG(1:Ntdc,1:min(7,Ntdc)),Ntdc)
CALL RINDD(fxind,Big(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,
& Nt,indI(1:NI),a_lo(1:Mb,1:NI-1),a_up(1:Mb,1:NI-1),
& INFIN(1:NI-1))
ansr(Ntd,1)=fxind(1)*CC
if (ALLOCATED(COV)) then !SCIS.GT.0
write(11,*) COV(1) ! save coefficient of variation
endif
print *,'Ready: ',Ntd,' of ',Ntime
enddo
if (ALLOCATED(COV)) then
close(11)
endif
goto 300
200 continue
XddInf=10.d0*SQRT(R4(1))
NI=7; Nd=3
Nc=4; Mb=2
allocate(BIG(1:Ntime+Nc+1,1:Ntime+Nc+1))
ALLOCATE(xc(1:Nc,1:Nx))
allocate(ex(1:Ntime+Nc+1))
ex=0.d0
xc(1,1:Nx)=h(1:Nx)
xc(2,1:Nx)=u
xc(3,1:Nx)=u
xc(4,1:Nx)=0.d0
! INFIN = INTEGER, array of integration limits flags: size 1 x Nb (in)
! if INFIN(I) < 0, Ith limits are (-infinity, infinity);
! if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
! if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
! if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
if (def.GT.0) then
INFIN(2)=-1
INFIN(4)=0
INFIN(5)=1
INFIN(6)=0
a_up(2,1)=1.d0 !*h
a_lo(1,1)=u
a_up(1,2)=XtInf ! X(ts) is redundant
a_lo(1,2)=-Xtinf
a_up(2,2)=1.d0 ! *h
a_lo(2,2)=1.d0 ! *h
a_up(2,3)=1.d0 !*h
a_lo(1,3)=u
a_lo(1,4)=-XddInf
a_up(1,5)= XdInf
a_lo(1,6)=-XdInf
else !def<0
INFIN(2)=-1
INFIN(4)=1
INFIN(5)=0
INFIN(6)=1
a_up(1,1)=u
a_lo(2,1)=1.d0 !*h
a_up(1,2)=XtInf ! X(ts) is redundant
a_lo(1,2)=-Xtinf
a_up(2,2)=1.d0 ! *h
a_lo(2,2)=1.d0 ! *h
a_up(1,3)=u
a_lo(2,3)=1.d0 !*h
a_up(1,4)=XddInf
a_lo(1,5)=-XdInf
a_up(1,6)=XdInf
endif
EPSOLD=ABSEPS
Nstart=MAX(Nstart,3)
do tn=Nstart,Ntime,1
Ntd=tn+1
Nt=Ntd-Nd
Ntdc=Ntd+Nc
indI(4)=Nt
indI(5)=Nt+1
indI(6)=Nt+2
indI(7)=Ntd
! IF (Ntd.GT.5.AND.INIT) THEN
! INIT=.FALSE.
! CALL INITDATA(speed)
! ENDIF
!if (SCIS.gt.1) Nj=Nt
!if (SCIS.gt.0) then
! if (SCIS.EQ.2) then
! Nj=max(Nt,0)
! else
! Nj=min(max(Nt-5, 0),0)
! endif
!endif
ABSEPS=MIN(SQRT(DBLE(tn))*EPSOLD*0.5D0,0.1D0)
do ts=2,FLOOR(DBLE(tn+1)/2.d0)
!print *,'ts,tn' ,ts,tn
CALL COV_INPUT(Big(1:Ntdc,1:Ntdc),tn,ts,R0,R1,R2,R3,R4) ! positive wave period
indI(2)=ts-2
indI(3)=ts-1
!CALL echo(BIG(1:Ntdc,1:min(7,Ntdc)),Ntdc)
!print *,'sp call rind'
CALL RINDD(fxind,Big(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,
& Nt,indI(1:NI),a_lo(1:Mb,1:NI-1),a_up(1:Mb,1:NI-1),
& INFIN(1:NI-1))
!CALL echo(BIG(1:Ntdc,1:min(7,Ntdc)),Ntdc)
!print *,'sp rind finished',fxind
!goto 900
SELECT CASE (ABS(def))
CASE (:2)
! 2, gives half wave period and wave crest amplitude (Tc,Ac).
! -2, gives half wave period and wave trough amplitude (Tt,At).
if (ts .EQ.tn-ts+1) then
ansr(tn,1:Nx)=ansr(tn,1:Nx)+fxind*CC*dt
else
ansr(tn,1:Nx)=ansr(tn,1:Nx)+fxind*CC*2.d0*dt
endif
CASE (3)
! 3, gives crest front period and wave crest amplitude (Tcf,Ac).
! -3, gives trough back period and wave trough amplitude (Ttb,At).
ansr(ts,1:Nx)=ansr(ts,1:Nx)+fxind*CC*dT
if ((ts.LT.tn-ts+1)) THEN
ansr(tn-ts+1,1:Nx)=ansr(tn-ts+1,1:Nx)+fxind*CC*dT ! exploiting the symmetry
endif
CASE (4:)
! 4, gives minimum of crest front/back period and wave crest amplitude (min(Tcf,Tcb),Ac).
! -4, gives minimum of trough front/back period and wave trough amplitude (min(Ttf,Ttb),At).
if (ts .EQ.tn-ts+1) then
ansr(ts,1:Nx)=ansr(ts,1:Nx)+fxind*CC*dt
else
ansr(ts,1:Nx)=ansr(ts,1:Nx)+fxind*CC*2.0*dt
endif
end select
enddo ! ts
print *,'Ready: ',tn,' of ',Ntime, ' ABSEPS = ', ABSEPS
enddo !tn
!print *,'ansr',ansr
300 open (unit=11, file='dens.out', STATUS='unknown')
!print *, ansr
do ts=1,Ntime
do ph=1,Nx
write(11,*) ansr(ts,ph)
! write(11,111) ansr(ts,ph)
enddo
enddo
!111 FORMAT(2x,F12.8)
close(11)
900 deallocate(big)
deallocate(fxind)
deallocate(ansr)
deallocate(xc)
deallocate(ex)
deallocate(R0)
deallocate(R1)
deallocate(R2)
if (allocated(COV) ) then
deallocate(COV)
endif
if (allocated(R3)) then
deallocate(R3)
deallocate(R4)
deallocate(h)
ENDIF
stop
!return
CONTAINS
SUBROUTINE INIT_LEVELS
& (U,def,Ntime,Nstart,speed,SCIS,seed1,Nx,dT,rateLHD)
IMPLICIT NONE
integer, intent(out):: def,Ntime,Nstart,speed,Nx,SCIS,seed1,
& rateLHD
double precision ,intent(out) :: U,dT
double precision :: XSPLT
integer :: NIT
OPEN(UNIT=14,FILE='reflev.in',STATUS= 'UNKNOWN')
READ (14,*) U
READ (14,*) def
READ (14,*) Ntime
READ (14,*) Nstart
READ (14,*) NIT
READ (14,*) speed
READ (14,*) SCIS
READ (14,*) seed1
READ (14,*) Nx
READ (14,*) dT
READ (14,*) rateLHD
READ (14,*) XSPLT
if (abs(def).GT.1) then
if (Ntime.lt.3) then
print *,'The number of wavelength points is too small, stop'
stop
end if
else
Nx=1
if (Ntime.lt.2) then
print *,'The number of wavelength points is too small, stop'
stop
end if
endif
CLOSE(UNIT=14)
RETURN
END SUBROUTINE INIT_LEVELS
C******************************************************
SUBROUTINE INIT_AMPLITUDES(h,def,Nx)
IMPLICIT NONE
double precision, dimension(:), intent(out) :: h
integer, intent(in) :: def
integer, intent(in) :: Nx
integer :: ix
OPEN(UNIT=4,FILE='h.in',STATUS= 'UNKNOWN')
C
C Reading in amplitudes
C
do ix=1,Nx
READ (4,*) H(ix)
enddo
CLOSE(UNIT=4)
!if (def.LT.0) THEN
! H=-H
!endif
RETURN
END SUBROUTINE INIT_AMPLITUDES
C**************************************************
C***********************************************************************
C***********************************************************************
SUBROUTINE INIT_COVARIANCES(Ntime,def,R0,R1,R2,R3,R4)
IMPLICIT NONE
double precision, dimension(:),intent(out) :: R0,R1,R2
double precision, dimension(:),intent(out) :: R3,R4
integer,intent(in) :: Ntime,def
integer :: i
open (unit=1, file='Cd0.in',STATUS='unknown')
open (unit=2, file='Cd1.in',STATUS='unknown')
open (unit=3, file='Cd2.in',STATUS='unknown')
do i=1,Ntime
read(1,*) R0(i)
read(2,*) R1(i)
read(3,*) R2(i)
enddo
close(1)
close(2)
close(3)
if (abs(def).GT.1) then
open (unit=4, file='Cd3.in',STATUS='unknown')
open (unit=5, file='Cd4.in',STATUS='unknown')
do i=1,Ntime
read(4,*) R3(i)
read(5,*) R4(i)
enddo
close(4)
close(5)
endif
return
END SUBROUTINE INIT_COVARIANCES
C***********************************************************************
C***********************************************************************
C**********************************************************************
SUBROUTINE COV_INPUT(BIG,tn,ts, R0,R1,R2,R3,R4)
IMPLICIT NONE
double precision, dimension(:,:),intent(inout) :: BIG
double precision, dimension(:),intent(in) :: R0,R1,R2
double precision, dimension(:),intent(in) :: R3,R4
integer ,intent(in) :: tn,ts
integer :: i,j,shft,Ntd1,N !=Ntdc
! the order of the variables in the covariance matrix
! are organized as follows:
! For ts>1:
! ||X(t2)..X(ts),..X(tn-1)||X''(ts) X'(t1) X'(tn)||X(ts) X(t1) X(tn) X'(ts)||
! = [Xt Xd Xc]
!
! For ts<=1:
! ||X(t2)..,..X(tn-1)||X'(t1) X'(tn)||X(t1) X(tn)||
! = [Xt Xd Xc]
! where
!
! Xt= time points in the indicator function
! Xd= derivatives
! Xc=variables to condition on
if (ts.LE.1) THEN
Ntd1=tn
N=Ntd1+2;
shft=0 ! def=1 want only crest period Tc
else
Ntd1=tn+1
N=Ntd1+4
shft=1 ! def=2 or 3 want Tc Ac or Tcf, Ac
endif
do i=1,tn-2
!cov(Xt)
do j=i,tn-2
BIG(i,j) = R0(j-i+1) ! cov(X(ti+1),X(tj+1))
enddo
!cov(Xt,Xc)
BIG(i ,Ntd1+1+shft) = R0(i+1) !cov(X(ti+1),X(t1))
BIG(tn-1-i ,Ntd1+2+shft) = R0(i+1) !cov(X(t.. ),X(tn))
!Cov(Xt,Xd)=cov(X(ti+1),x(tj)
BIG(i,Ntd1-1) =-R1(i+1) !cov(X(ti+1),X' (t1))
BIG(tn-1-i,Ntd1)= R1(i+1) !cov(X(ti+1),X' (tn))
enddo
!call echo(big(1:tn,1:tn),tn)
!cov(Xd)
BIG(Ntd1 ,Ntd1 ) = -R2(1)
BIG(Ntd1-1,Ntd1 ) = -R2(tn) !cov(X'(t1),X'(tn))
BIG(Ntd1-1,Ntd1-1) = -R2(1)
!cov(Xc)
BIG(Ntd1+1+shft,Ntd1+1+shft) = R0(1) ! cov(X(t1),X (t1))
BIG(Ntd1+1+shft,Ntd1+2+shft) = R0(tn) ! cov(X(t1),X (tn))
BIG(Ntd1+2+shft,Ntd1+2+shft) = R0(1) ! cov(X(tn),X (tn))
!cov(Xd,Xc)
BIG(Ntd1 ,Ntd1+1+shft) = R1(tn) !cov(X'(tn),X(t1))
BIG(Ntd1 ,Ntd1+2+shft) = 0.d0 !cov(X'(tn),X(tn))
BIG(Ntd1-1,Ntd1+1+shft) = 0.d0 !cov(X'(t1),X(t1))
BIG(Ntd1-1,Ntd1+2+shft) =-R1(tn) !cov(X'(t1),X(tn))
if (ts.GT.1) then
!
!cov(Xc)
BIG(Ntd1+1,Ntd1+1) = R0(1) ! cov(X(ts),X (ts)
BIG(Ntd1+1,Ntd1+2) = R0(ts) ! cov(X(ts),X (t1))
BIG(Ntd1+1,Ntd1+3) = R0(tn+1-ts) ! cov(X(ts),X (tn))
BIG(Ntd1+1,Ntd1+4) = 0.d0 ! cov(X(ts),X'(ts))
BIG(Ntd1+2,Ntd1+4) = R1(ts) ! cov(X(t1),X'(ts))
BIG(Ntd1+3,Ntd1+4) = -R1(tn+1-ts) !cov(X(tn),X'(ts))
BIG(Ntd1+4,Ntd1+4) = -R2(1) ! cov(X'(ts),X'(ts))
!cov(Xd)
BIG(Ntd1-2,Ntd1-1) = -R3(ts) !cov(X''(ts),X'(t1))
BIG(Ntd1-2,Ntd1-2) = R4(1)
BIG(Ntd1-2,Ntd1 ) = R3(tn+1-ts) !cov(X''(ts),X'(tn))
!cov(Xd,Xc)
BIG(Ntd1 ,Ntd1+4) =-R2(tn+1-ts) !cov(X'(tn),X'(ts))
BIG(Ntd1 ,Ntd1+1) = R1(tn+1-ts) !cov(X'(tn),X (ts))
BIG(Ntd1-1,Ntd1+4) =-R2(ts) !cov(X'(t1),X'(ts))
BIG(Ntd1-1,Ntd1+1) =-R1(ts) !cov(X'(t1),X (ts))
BIG(Ntd1-2,Ntd1+1) = R2(1) !cov(X''(ts),X (ts)
BIG(Ntd1-2,Ntd1+2) = R2(ts) !cov(X''(ts),X (t1))
BIG(Ntd1-2,Ntd1+3) = R2(tn+1-ts) !cov(X''(ts),X (tn))
BIG(Ntd1-2,Ntd1+4) = 0.d0 !cov(X''(ts),X'(ts))
!cov(Xt,Xc)
do i=1,tn-2
j=abs(i+1-ts)
BIG(i,Ntd1+1) = R0(j+1) !cov(X(ti+1),X(ts))
BIG(i,Ntd1+4) = sign(R1(j+1),R1(j+1)*dble(ts-i-1)) !cov(X(ti+1),X'(ts)) ! check this
!Cov(Xt,Xd)=cov(X(ti+1),X(ts))
BIG(i,Ntd1-2) = R2(j+1) !cov(X(ti+1),X''(ts))
enddo
endif ! ts>1
!call echo(big(1:N,1:N),N)
! make lower triangular part equal to upper
do j=1,N-1
do i=j+1,N
BIG(i,j) =BIG(j,i)
enddo
!call echo(big(1:N,1:N),N)
enddo
!call echo(big(1:N,1:N),N)
C write (*,10) ((BIG(j,i),i=N+1,N+6),j=N+1,N+6)
C 10 format(6F8.4)
RETURN
END SUBROUTINE COV_INPUT
SUBROUTINE COV_INPUT2(BIG,pt, R0,R1,R2)
IMPLICIT NONE
double precision, dimension(:,:), intent(out) :: BIG
double precision, dimension(:), intent(in) :: R0,R1,R2
integer :: pt,i,j
! the order of the variables in the covariance matrix
! are organized as follows;
! X(t2)...X(tn-1) X'(t1) X'(tn) X(t1) X(tn) = [Xt Xd Xc]
!
! where Xd is the derivatives
!
! Xt= time points in the indicator function
! Xd= derivatives
! Xc=variables to condition on
!cov(Xc)
BIG(pt+2,pt+2) = R0(1)
BIG(pt+1,pt+1) = R0(1)
BIG(pt+1,pt+2) = R0(pt)
!cov(Xd)
BIG(pt,pt) = -R2(1)
BIG(pt-1,pt-1) = -R2(1)
BIG(pt-1,pt) = -R2(pt)
!cov(Xd,Xc)
BIG(pt,pt+2) = 0.d0
BIG(pt,pt+1) = R1(pt)
BIG(pt-1,pt+2) = -R1(pt)
BIG(pt-1,pt+1) = 0.d0
if (pt.GT.2) then
!cov(Xt)
do i=1,pt-2
do j=i,pt-2
BIG(i,j) = R0(j-i+1)
enddo
enddo
!cov(Xt,Xc)
do i=1,pt-2
BIG(i,pt+1) = R0(i+1)
BIG(pt-1-i,pt+2) = R0(i+1)
enddo
!Cov(Xt,Xd)=cov(X(ti+1),x(tj))
do i=1,pt-2
BIG(i,pt-1) = -R1(i+1)
BIG(pt-1-i,pt)= R1(i+1)
enddo
endif
! make lower triangular part equal to upper
do j=1,pt+1
do i=j+1,pt+2
BIG(i,j)=BIG(j,i)
enddo
enddo
C write (*,10) ((BIG(j,i),i=N+1,N+6),j=N+1,N+6)
C 10 format(6F8.4)
RETURN
END SUBROUTINE COV_INPUT2
END PROGRAM sp2thpdf