c AGUA_GEM6
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                                                      ######
c     ######            PROGRAM FOR SUPPORT DAILY WEATHER         ######
c     ######        SIMULATION IN THE CONTIGOUS UNITED STATES     ######
c     ######                                                      ######
c     ######                     ARS - 114                        ######
c     ######                   (c) Jan   1995                     ######
c     ######                                                      ######
c     ######        UNITED STATES DEPARTMENT OF AGRICULTUE        ######
c     ######            AGRICULTURAL RESEARCH SERVICE.            ######
c     ######                  All rights reserved.                ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

        PROGRAM AGUA_GEM6
c
c#######################################################################
c
c     PURPOSE:
c    
c     This program calculates the parameter matrix Z(I,K) for each
c     station in "USCLIMTE.SITE" file. J from 1 to 4 are P00, P10,
c     beta and mu, K from 1 to 13 are the annual mean, first 6 harmonic's
c     amplitude and phase angle of above 4 parameters.
c
c#######################################################################
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
c     OUTPUT:
c
c       za       parameter (p00,p10,beta,mu) matrix for 
c                USCLIMATE model (upto 6 harmonics).
c                1. Annual mean 
c                2. Amplitude of first harmonic
c                3. Phase angle of first harmonic
c                4. Amplitude of second harmonic
c                5. Phase angle of second harmonic
c                6. Amplitude of third harmonic
c                7. Phase angle of third harmonic
c                8. Amplitude of fourth harmonic
c                9. Phase angle of fourth harmonic
c               10. Amplitude of fifth harmonic
c               11. Phase angle of fifth harmonic
c               12. Amplitude of sixth harmonic
c               13. Phase angle of sixth harmonic
c       alfa    weighting function
c
c     INPUT:
c
c     pre         M years of daily precipitation
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
	implicit  none
	integer   idyr         ! max. # of day in one year.
	integer   iyr          ! max. # of year.
	integer   mper, mmper
        integer ilen           ! length of input file name in characters
	parameter (idyr=365,iyr=30,mper=26,mmper=14)
        common /const1/pi,pi2
	common /use/pre(iyr,idyr),ncount,thres
	common /mrv/pdd(mper),pwd(mper),pddme,pwdme
	common /ocho/pr(mper,400),mx(mper),xmu(mper),alfa(mper)
     !,beta(mper),theta(mper)
	common /one/nper,nnper
	common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
        common /tres/c(mmper,mper),s(mmper,mper),a(mmper)
     !,a1(mmper),ncn,np
	real     c,s,a,a1
	real     pi,pi2
	real     pre,pr,xmu,thres
	real     pdd,pwd,alfa,beta,theta,pddme,pwdme
	real     amp,pan,almean,bemean,themean,xmumean
	real     za(4,7),albar
	integer  iy,ida,mo
	integer  j,i,nper,nnper
	integer  mx,ncount,maxh,ncn,np
c
        character*40 header
        integer numyears
        real alat
        character*25 filename,site,infile
c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c#######################################################################
c      constant defination.
c#######################################################################

        pi=atan(1.)*4
	pi2=pi*2
        ncount=(30-1)             ! real # of year in data set.
	maxh=6                    ! max # of harmonics.
	thres=0.008               ! threshold of wet or dry day.
        nper=26                   ! # of periods.
        nnper=14                  ! # of days in each period.
c#######################################################################
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


         write (6,*)'Enter Location file name (no extension-no blanks) '
         read(5,*) filename
         ilen = index(filename,' ')              ! string position of 1st ' '
c
         infile=filename(1:ilen-1)// '.txt'      ! example boise.txt
         site=filename(1:ilen-1)// '.site'       ! from AGUA46 - ppt output
c        open(unit=12,file='boise61a.txt', status='old')

        open(unit=12,file=infile, status='unknown')   !example: boise.txt
c=======================================================================
          read(12,*) numyears,alat               ! decimal value for alat
c
          read(12,*) header     ! header for boise61a.txt  !skip the header
c                                ! for data input
c========================================================================
	do i=1,59                                            ! offset to march
ccc 1      read(12,*,end=123)imo,ida,iyr,v6,v1,v2,v3,v4,v5   ! boise61a.txt
c
 1       read(12,*,end=123) iy,mo,ida,pre(1,i)
	if(mo.eq.2.and.ida.eq.29) goto 1
 123     continue
        end do

	do i=1,ncount
	do j=1,365
2       read(12,*,end=321) mo,iy,ida,pre(i,j)          !boise61a.txt
CCC        pre(i,j)=pre(i,j)/100.
	if(mo.eq.2.and.ida.eq.29) goto 2
321     continue
	end do
	end do

c#######################################################################
c      This CALL is for starting the calculation.
c#######################################################################

           CALL MARKOV            ! solve for P00 and P10
           CALL CALCUL            ! solve for beta, mu and alfa

c#######################################################################
c      Here just pick up annual mean and first three harmonics.
c      Users can change DO loop # 40: j=1,n (the max. of n is 6)
c      and print out the largest harmonic insterested.
c#######################################################################

           za(1,1)=pddme
           za(2,1)=pwdme
           za(3,1)=bemean
           za(4,1)=xmumean
           albar=almean
           DO 40 j=1,3
             za(1,2*j)=amp(j,1)
             za(2,2*j)=amp(j,2)
             za(3,2*j)=amp(j,4)
             za(4,2*j)=amp(j,6)
             za(1,2*j+1)=pan(j,1)
             za(2,2*j+1)=pan(j,2)
             za(3,2*j+1)=pan(j,4)
             za(4,2*j+1)=pan(j,6)
40        CONTINUE
c#######################################################################
c      This loop is for printing the final results.
c#######################################################################
        open(unit=16,file=site, status='new')    ! example boise.site
c
        write (16,61) numyears,filename          ! Location / filename
 61     format (1x,i4,2x,a24)
c
c
        write (16,62) alat                    !ADD Longitude & elevation
 62     format (1x,f9.4)
c
          DO 30 j=1,7
            write(16,'(4f11.7)')(za(i,j),i=1,4)
30        CONTINUE
          write(16,'(f11.7)')albar
10      CONTINUE

        STOP
        END


c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE MARKOV                   ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE MARKOV
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the transition probabilities by using first_order
c     Markov chain, to describethe dependence between wet and dry 
c     day occurrences on sussive days.
c
c#######################################################################
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      per        NCOUNT years daily precipitation.
c
c     OUTPUT:
c
c      pdd        transition probability for dry-dry day.
c      pwd        transition probability for wet-dry day.
c      pddme      annual mean of pdd.
c      pwdme      annual mean of pwd.
c
c     WORK ARRAY:
c           
c     For total NCOUNT years:
c     nww         number of wet-wet days in each 14 day period.
c     nwd         number of wet-dry days in each 14 day period.
c     ndw         number of dry-wet days in each 14 day period.
c     ndd         number of dry-dry days in each 14 day period.
c
c     nwwt        total number of wet-wet days.
c     nwdt        total number of wet-dry days.
c     ndwt        total number of dry-wet days.
c     nddt        total number of dry-dry days.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################

	integer iper, iiper
        parameter(iper=26,iiper=14)
	common /mrv/pdd(iper),pwd(iper),pddme,pwdme
	common /use/pre(30,365),ncount,thres
	common /one/nper,nnper
        common /tres/c(iiper,iper),s(iiper,iper),a(iiper)
     !,a1(iiper),ncn,np
	common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
        common /const1/pi,pi2
       integer ndd(iper),ndw(iper),nwd(iper),nww(iper)
       logical logic1,logic2

c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c#######################################################################
c      clean the work space.
c#######################################################################

       
      DO 21 i=1,nper
        ndd(i)=0
        ndw(i)=0
        nwd(i)=0
        nww(i)=0
21    CONTINUE
      nddt=0 
      ndwt=0 
      nwdt=0 
      nwwt=0

c#######################################################################
c      Start to calculate  P00 and P10 for each 14 day period.
c#######################################################################

      logic1=.FALSE.                            ! start day 1 as wet.
      n=1
      DO 10 j=1,NCOUNT
      DO 10 i=1,365
        logic2=.FALSE.                          ! start day 2 as wet.
        n=i/nnper+1
        n=min(n,nper)
c#######################################################################
c       Program from here is for filter out the missing value.
c       For the first day and the first starting day after missing
c       days, save its wet/dry information in "logic1" and check its
c       successive day.
c#######################################################################
        IF(pre(j,i).lt.0.) goto 10              
        IF(pre(j,i-1).lt.0.) THEN               
          IF(pre(j,i).le.thres)THEN            
            logic1=.true.                       
          ELSE
            logic1=.false.                      
          END IF
          goto 10
        END IF  
        IF((i.EQ.1).AND.(j.EQ.1)) THEN
          IF(pre(1,1).LE.thres)logic1=.TRUE.
          GO TO 10
        END IF
c#######################################################################
c      From here program is checking the second non-missing data day.
c      and save its dry/wet information in "logic2".
c#######################################################################
        IF(pre(j,i).LE.thres)logic2=.TRUE.
c#######################################################################
c       Do the two day checking and save the checking result in the 
c       correponding place, such as saving wet-wet case in nww array
c       and nwwt.
c#######################################################################
        IF(logic1)GO TO 15
        IF(logic2)GO TO 20
        nwwt=nwwt+1
        nww(n)=nww(n)+1
        GO TO 16
 15     IF(logic2)GO TO 18
        ndwt=ndwt+1
        ndw(n)=ndw(n)+1
        GO TO 16
 20     nwdt=nwdt+1
        nwd(n)=nwd(n)+1
        GO TO 16
 18     nddt=nddt+1
        ndd(n)=ndd(n)+1
 16     logic1=logic2
 10   CONTINUE
c#######################################################################
c     Calculate the P00 and P10 values for each period.
c#######################################################################
      DO 22 i=1,nper
        pdd(i)=FLOAT(ndd(i))/FLOAT(ndd(i)+ndw(i))
        pwd(i)=FLOAT(nwd(i))/FLOAT(nwd(i)+nww(i))
   22 CONTINUE
c#######################################################################
c     Annual mean of P00 and P10.
c#######################################################################
      pddme=FLOAT(nddt)/FLOAT(nddt+ndwt)
      pwdme=FLOAT(nwdt)/FLOAT(nwdt+nwwt)
c#######################################################################
c     Calculate the Fourier Coefficients for 6 harmonics.
c#######################################################################
      CALL FOUTER(1)

      RETURN
      END

        
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE CALCUL                   ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE CALCUL
c
c#######################################################################
c
c     PURPOSE:
c
c     Using mixed exponential distribution to computer estimates of
c     alfa, beta and theta for each 14 day period. The values with
c     maximum likelihood were chosen as final output. The path for
c     15 iterations to approch the final results is:
c     1) beta is slightly less than than the observed mean;
c     2) theta is slightly greater than than the observed mean;
c     3) alfa is equal to mean experssion.
c     Finish to calculate the Fourier Coefficients for 6 harmonics.
c
c#######################################################################
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      pre        the NCOUNT years daily precipitation
c
c     OUTPUT:
c
c      alfa       weighting parameter with values between 0 and 1.
c      beta       the means of the smaller exponential distributions.
c      theta      the means of the smaller exponential distributions.
c
c     WORK ARRAY:
c
c      pr         the daily precipitation amount minus a threshold.
c      mu     the mean precipitaion in each 14 day period.
c      xm         the wet days in each 14 day period.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################

       parameter (mper=26, mmper=14)
       common /use/pre(30,365),ncount,thres
       common /mrv/pdd(mper),pwd(mper),pddme,pwdme
       common /ocho/pr(mper,400),mx(mper),xmu(mper),alfa(mper)
     !,beta(mper),theta(mper)
       common /one/nper,nnper
       common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
       common /tres/c(mmper,mper),s(mmper,mper),a(mmper),a1(mmper)
     1,ncn,np
       common /const1/pi,pi2

c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c#######################################################################
c      Calculate the mean precipitation, number of wet days and
c      precipitation amount in wet days for each 14 day period. 
c#######################################################################

      DO 5 j=1,nper
        xmu(j)=0.
        mx(j)=0
5     CONTINUE                  

      n=1
      DO 10 J=1,NCOUNT
        DO 20 i=1,365
          n=(i-1)/nnper+1
          n=min(n,nper)
          IF(pre(j,i).LT.0.) GO TO 20
          IF(pre(j,i).LE.thres)GO TO 20
          mx(n)=mx(n)+1
          IF(mx(n).LE.400)GO TO 8
8         pr(n,mx(n))=pre(j,i)-thres
20      CONTINUE
10    CONTINUE
      
      DO 30 i=1,nper
        mm=mx(i)
        DO 40 j=1,mm
          xmu(i)=xmu(i)+pr(i,j)
40      CONTINUE
        xmu(i)=xmu(i)/float(mm)
30    CONTINUE                  
c#######################################################################
c     Calculate beta, mu and alfa.
c#######################################################################
      CALL MIXEXP
c#######################################################################
c     Calculate beta, mu and alfa's first 6 harmonic's Fourier coeffs.
c#######################################################################
      CALL FOUTER(3)

      RETURN
      END


c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE MIXEXP                   ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE MIXEXP
c
c#######################################################################
c
c     PURPOSE:
c
c     Using mixed exponential distribution to computer estimates of
c     alfa, beta and theta for each 14 day period. The values with
c     maximum likelihood were chosen as final output. The path for
c     15 iterations to approch the final results is:
c     1) beta is slightly less than than the observed mean;
c     2) theta is slightly greater than than the observed mean;
c     3) alfa is equal to mean experssion.
c     And mu = alfa*beta+(1-alfa)*theta.
c
c#######################################################################
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      mu     the mean precipitaion in each 14 day period.
c      xm         the wet days in each 14 day period.
c      pr         the daily precipitation amount minus a threshold.
c
c     OUTPUT:
c
c      alfa       weighting parameter with values between 0 and 1.
c      beta       the means of the smaller exponential distributions.
c      theta      the means of the smaller exponential distributions.
c
c     WORK ARRAY:
c
c     rr          
c     ss          
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
       parameter (mper=26,mmper=14)
       common /one/nper,nnper
       common /ocho/pr(mper,400),mx(mper),xmu(mper),alfa(mper)
     !,beta(mper),theta(mper)
       common /tres/c(mmper,mper),s(mmper,mper),a(mmper),a1(mmper)
     !,ncn,np
       common /const1/pi,pi2
       common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean

c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


      DO 15 i=1,nper
c#######################################################################
c      Give the initial values for each period.
c#######################################################################
        iper=i
        mm=mx(i)
        da=0.1
        dtonea=0.1*xmu(i)
     
        DO 16 k=1,15
c#######################################################################
c      Give the step values for each iteration.
c#######################################################################
          da=da+0.05
          dtonea=dtonea+0.05*xmu(i)
          dttwoa=(xmu(i)-da*dtonea)/(1.-da)
          sume=0.
c#######################################################################
c      This loop calculate the mixed exponential distribution.
c       (ref. Smith and Schreiber 1974)
c#######################################################################
          DO 17 J=1,mm
            rr=EXP(-pr(i,j)/dtonea)
            ss=EXP(-pr(i,j)/dttwoa)
            sume=sume+ALOG((da*rr/dtonea)+((1.-da)*ss/dttwoa))
   17     CONTINUE
c#######################################################################
c      Look for the max value.
c#######################################################################
          IF(k.EQ.1.or.sume.GT.xmt) THEN
            xmt=sume
            alfa(i)=da
            IF(dtonea.gt.dttwoa) THEN
	      beta(i)=dttwoa
              theta(i)=dtonea
            ELSE
	      beta(i)=dtonea
              theta(i)=dttwoa
            END IF
c#######################################################################
c      Adjustment for alfa.
c#######################################################################
	    IF((theta(i)-beta(i)).lt.0.01)alfa(i)=1.-da
	    IF(alfa(i).gt.0.51)alfa(i)=1.-da
          END IF
   16   CONTINUE
   15 CONTINUE
      
      RETURN
      END

c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE FOUTER                   ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE FOUTER(nstudy)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the mean, amplitude and phase angle for each 
c     inputed variable (i.e., P00, P10, beta, theta, alfa, mu).
c
c#######################################################################
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      f         The mean values in 26 periods
c
c     OUTPUT:
c
c      c          cosines function array
c      s          sines function array
c      amp        amplitude
c      phi        phase angle
c      ymean      mean
c
c     WORK ARRAY:
c     a           cosine coefficients
c     a1          sine coefficients
c     rm          amplitude
c     phi         phase angle
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      parameter (mper=26,mmper=14)
      common /ocho/pr(mper,400),mx(mper),xmu(mper),alfa(mper)
     !,beta(mper),theta(mper)
      common /one/nper,nnper
      common /tres/c(mmper,mper),s(mmper,mper),a(mmper),a1(mmper)
     !,ncn,np
      common /const1/pi,pi2
      common /mrv/pdd(mper),pwd(mper),pddme,pwdme
      common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
      dimension ymean(4)

c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


      DO 30 i=1,4
	ymean(i)=0.
30    CONTINUE
      IF(nstudy.EQ.1)GO TO 10
      CALL SERIES(alfa,1,3,ymean)      
      almean=ymean(1)
      CALL SERIES(beta,2,4,ymean)      
      bemean=ymean(2)
      CALL SERIES(theta,2,5,ymean)      
      themean=ymean(3)
      CALL SERIES(xmu,2,6,ymean)
      xmumean=ymean(4)
      GO TO 20
   10 CALL SERIES(pdd,1,1,ymean)      
      CALL SERIES(pwd,1,2,ymean)      
   20 RETURN
      END

c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE SERIES                   ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE SERIES(f,ipr,nstudy,ymean)      
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the mean, amplitude and phase angle for each 
c     inputed variable (i.e., P00, P10, beta, theta, alfa, mu).
c
c#######################################################################
c
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      f         The value of 14 periods
c
c     OUTPUT:
c
c      c          cosines function array
c      s          sines function array
c      amp        amplitude
c      phi        phase angle
c      ymean      mean
c
c     WORK ARRAY:
c     a           cosine coefficients
c     a1          sine coefficients
c     rm          amplitude
c     phi         phase angle
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      common /one/nper,nnper
      common /tres/c(14,26),s(14,26),a(14),a1(14)
     1,ncn,np
      common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
      common /const1/pi,pi2
      dimension f(26),ymean(4)


c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c#######################################################################
c      Calculate cosine and sine values to save time in
c      COEFF's calculation.
c#######################################################################

      np=nper
      ncn=(nper-1)/2
      t=pi2/FLOAT(np)
      
      DO 20 i=1,ncn
      DO 20 k=1,np
         c(i,k)=COS(FLOAT(i-1)*(k-1)*t)
         IF(i.EQ.ncn) GO TO 20
         s(i,k)=SIN(FLOAT(i)*(k-1)*t)
 20   CONTINUE
      
      CALL COEFF(f,nstudy,ymean)
      
      RETURN
      END
      
      
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                  SUBROUTINE COEFF                    ######
c     ######                                                      ######
c     ######                     (c) 1995                         ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

      SUBROUTINE COEFF(f,nstudy,ymean)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the mean, amplitude and phase angle for each 
c     inputed variable (i.e., P00, P10, beta, theta, alfa, mu).
c
c#######################################################################
c
c
c     AUTHOR:  (1986)
c             David Woolhiser - Tucson, Az
c             Jose Roldan - Spain
c
c     Modify history:
c
c     Yunyun Lu (Jan,1995)
c     transfered program from machine VAX to HP
c     reduced the unnecessary procedure (only for our project research).
c     added the full documentation
c
c#######################################################################
c
c     INPUT:
c
c      f         The value of 14 periods
c      c         cosines function array
c      s         sines function array
c
c     OUTPUT:
c
c      amp        amplitude
c      phi        phase angle
c      ymean      mean
c
c     WORK ARRAY:
c     a           cosine coefficients
c     a1          sine coefficients
c     rm          amplitude
c     phi         phase angle
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      common /tres/c(14,26),s(14,26),a(14),a1(14)
     1,ncn,np
      common /seis/amp(6,6)
     !,pan(6,6),maxh,almean,bemean,themean,xmumean
      common /const1/pi,pi2
      dimension rm(14),phi(14),f(26),ymean(4)

c#######################################################################
c
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C     Beginning of executable code...
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c#######################################################################
c      Calculate cosine and sine coefficients
c#######################################################################

      DO 10 i=1,ncn
        a(i)=0.
        a1(i)=0.
        DO 20 k=1,np
           a(i)=a(i)+f(k)*c(i,k)
           IF(i.EQ.ncn) GO TO 20
           a1(i)=a1(i)+f(k)*s(i,k)
 20     CONTINUE
        a(i)=(a(i)*2.)/FLOAT(np)
        a1(i)=(a1(i)*2.)/FLOAT(np)

c#######################################################################
c      Calculate amplitude and phase angle
c#######################################################################

        IF(i.EQ.1) THEN
           phi(i)=0.2
           rm(i)=0.
        ELSE
           phi(i)=ATAN2(a(i),a1(i-1))
           rm(i)=a(i)/SIN(phi(i))
        END IF
 10   CONTINUE
      DO 30 i=1,MAXH
         amp(i,nstudy)=rm(i+1)
         pan(i,nstudy)=phi(i+1)-0.11189*FLOAT(i)
 30   CONTINUE
c#######################################################################
c      Calculate mean
c#######################################################################
      IF(nstudy.gt.2) ymean(nstudy-2)=a(1)/2.
      RETURN
      END