!*********************************************	
! Program mem_deep_raw: Adaptation of Bill O'Reilly MEM estimator for 
! reading time-series of frequency spectra and directional moments and 
! generating f-dir spectra with a given frequency grid (in input file crest.inp) 
! and direction resolution. These spectra are also back-refracted to deep 
! water using Snel's law (all this is hard-wired in the main routine at the end 
! of the present file). 
!*********************************************	
	subroutine mem_est(nf,nabin,e,ec,es,ec2,es2,ds)
!    mem_est.f
!  ** THIS version really reads in normailzed four. coeff.,s
!   so no division by e
!    subroutine to make MEM estimate of the directional spectrum
!    in 33 .01hz frequency bands
!    Input:   energy and 4 directional moments from a slope array
! (e,ec,es..) in trig coord frame
!    Output: Directional spectrum in 3 degree directional bins (in
!     (ds)   the "compass heading" coordinate frame) for
!            the 33 fbands between .04 and .25hz
!            
!  If nabin=120 directional bands,  1=3 degress....120=360 degrees
!            where the direction is the compass heading from which
!            the wave energy is arriving. e.g. 270=from west, 180=
!            from south  etc.
!*********************************************	
   IMPLICIT NONE
	INTEGER             , intent(in)  :: nf,nabin
   REAL, dimension(:)  , intent(in)  :: e,ec,es,ec2,es2
   REAL, dimension(:,:), intent(inout) :: ds
   
   REAL, dimension(360) :: d 
   REAL a1,b1,a2,b2,chk
   INTEGER I,J,K,Kbin
! loop thru freq bands
	Kbin=360/nabin
   do i=1,nf
! moments are first four fourier coeff. of the directional
!    distribution, MEM uses normalized values
!		a1=ec(i)
!		b1=es(i)
!		a2=ec2(i)
!		b2=es2(i)
		a1=ec(i)/e(i)
		b1=es(i)/e(i)
		a2=ec2(i)/e(i)
		b2=es2(i)/e(i)
! get MEM estimate (d) with 1 deg resolution, d is returned in
!    compass coordinates
		IF (e(i).EQ.0.) THEN 
			ds(:,i)=0.
		ELSE
	   	   call mem(a1,a2,b1,b2,d,chk)
! merge into Kbin degrees directional bins
                   do j=Kbin,360-Kbin,Kbin
	              ds(j/Kbin,i)=e(i)*(SUM(d(j-kbin/2:j+kbin/2)))
                   ENDDO
	           ds(nabin,i)=e(i)*(SUM(d(360-kbin/2:360))+SUM(d(1:kbin/2)))
                ENDIF
   ENDDO
	return
	end

!*********************************************	
	subroutine  mem(a1,a2,b1,b2,s,chk)
!	Maximum entropy method (MEM) for the estimation of a directional
!     distribution from pitch-and-roll data.
!     (Lygre and Krogstad, JPO v16 1986: NOTE - there is a typo in the
!     paper...BOTH exponetials in eq. 13 should have negative signs.
!     This is fixed in Krogstad, 1991- THH has the book)
!	variables:
!            a1,b1,b1,b2  fourier coef. of dir. dist.
!            (normalized ala Long [1980]) - TRIG COORDINATES
!  
!            s= MEM estimate od directional spectrum, 1 deg.
!               resolution - COMPASS COORDINATES
!     chk = check factor: MEM likes to make narrow spikes for
!           directional spectra if it can (it fits the a1,a2,
!           b1 abd b2's exactly).  If it does this, then the
!           .1 deg stepsize for making the 1 deg resolution
!           estimate is too coarse.  The check factor should
!           be close to 1.
!*********************************************	
	dimension s(360)
	complex c1,c2,p1,p2,e1,e2,x,y  
	data dr/0.0174533/

	do i=1,360
	 s(i)=0.
	end do
! switch to Lygre & Krogstad notation
	d1=a1
	d2=b1
	d3=a2
	d4=b2
	c1=(1.,0)*d1+(0,1.)*d2
	c2=(1.,0)*d3+(0,1.)*d4
	p1=(c1-c2*conjg(c1))/(1.-cabs(c1)**2)
	p2=c2-c1*p1
	x=1.-p1*conjg(c1)-p2*conjg(c2) 
! sum in .1 deg steps, get distribution with 1 degree resolution
	tot=0
	nn=0
	do 100 n=5,3604
		a=real(n)*.1*dr
		e1=(1.,0)*cos(a)-(0,1.)*sin(a)
		e2=(1.,0)*cos(2*a)-(0,1.)*sin(2*a)
		y=cabs((1.,0)-p1*e1-p2*e2)**2
! put in proper 1 deg directional band
	       ndir=1+nn/10
! switch from trig to compass corrdinates
		ndir=270-ndir
		if(ndir.gt.360) ndir=ndir-360
		if(ndir.lt.1) ndir=ndir+360
! normalize by 360/(step size)
		s(ndir)=s(ndir)+cabs(x/y)/3600.
		tot=tot+cabs(x/y)
	       nn=nn+1
100	continue
! tot should = 3600.  If directional peak is extremely narrow then
!    1 deg resolution may be insufficient and tot .ne. 360
       chk=tot/3600.
	return
	end

!****************************************************************************    
!*********************************************************************
INTEGER FUNCTION julday(id,mm,iyyy)
! See numerical recipes 2nd ed. The order of month and day have been swapped!
!*********************************************************************
   IMPLICIT NONE
	INTEGER(KIND=4),    INTENT(in)  :: id,mm,iyyy
   INTEGER(KIND=4), PARAMETER :: IGREG=15+31*(10+12*1582)
   INTEGER(KIND=4) ja,jm,jy
   jy=iyyy
   IF (jy.EQ.0) WRITE(6,*) 'There is no zero year !!'
   IF (jy.LT.0) jy=jy+1
   IF (mm.GT.2) THEN
      jm=mm+1
   ELSE
      jy=jy-1
      jm=mm+13
   ENDIF
   julday=INT(365.25*jy)+int(30.6001*jm)+id+1720995
   IF (id+31*(mm+12*iyyy).GE.IGREG) THEN
      ja=INT(0.01*jy)
      julday=julday+2-ja+INT(0.25*ja)
   ENDIF
   RETURN
END

!*********************************************************************
SUBROUTINE caldat(julian,id,mm,iyyy)
! See numerical recipes 2nd ed. The order of month and day have been swapped!
!*********************************************************************
   IMPLICIT NONE
	INTEGER(KIND=4),    INTENT(in)  :: julian
   INTEGER(KIND=4),    INTENT(out) :: id,mm,iyyy
   INTEGER(KIND=4), PARAMETER :: IGREG=2299161
   INTEGER(KIND=4) ja,jalpha,jb,jc,jd,je
   if (julian.GE.IGREG) THEN
      jalpha=INT(((julian-1867216)-0.25)/36524.25)
      ja=julian+1+jalpha-INT(0.25*jalpha)
   ELSE
      ja=julian
   ENDIF
   jb=ja+1524
   jc=INT(6680.+((jb-2439870)-122.1)/365.25)
   jd=365*jc+INT(0.25*jc)
   je=INT((jb-jd)/30.6001)
   id=jb-jd-INT(30.6001*je)
   mm=je-1
   IF (mm.GT.12) mm=mm-12
   iyyy=jc-4715
   IF (mm.GT.2) iyyy=iyyy-1
   IF (iyyy.LE.0) iyyy=iyyy-1
   RETURN
   END
!*********************************************************************
   
PROGRAM average
   IMPLICIT NONE

   INTEGER     I,J,k,irec,entier,time
   INTEGER     np,nfband,npar,Jfp,nabin,kfp
   INTEGER     iyc,imc,idc,ihc,iminc,icode,iostat,finished,iminlast,ihlast
   INTEGER   iyend,imend,idend,ihend
   INTEGER(KIND=4) julday,Jdayend
   REAL,    PARAMETER :: pi = 3.14159265,tpi=6.2831853,d2r=.01745329
   REAL   , DIMENSION(:,:),     ALLOCATABLE :: ds,ds2
   REAL   , DIMENSION(:),     ALLOCATABLE :: freq,theta,freq1,freq2,df
   REAL, DIMENSION(:),   ALLOCATABLE :: Ef,a1,a2,b1,b2,m1,m2,th1,sth1,th2,sth2
   REAL, DIMENSION(:),   ALLOCATABLE :: Efa,a1a,a2a,b1a,b2a,th1a,th2a,sth1a,sth2a
   INTEGER  iy,im,id,ih,imin,fileend,iyo,imo,ido,iho,imino
   INTEGER(KIND=4)   Jday,Jdayc,Jdayi,Jdaylast
   
   REAL HS1,HS2,HS3,HS4,HS5,HS6,HS7,HS8  !various estimates of Hs
   REAL fp,fmean,emax,emax0,emax2,dt,tavg               !peak frequency
   REAL Timezone_data,w,w2
   REAL Hs,Hs10,fp10,sth1p,th1p,f02,th2p,sth2p,lat,lon,depth,U10, &
                  Udir,Curr,Currdir
   CHARACTER(80) finp,filename1,filename2,format1,format2
   CHARACTER(120) chaine,string1,string2
   CHARACTER(8), DIMENSION(:),   ALLOCATABLE :: buoyname
   CHARACTER space
   interface
   	subroutine mem_est(nf,nabin,e,ec,es,ec2,es2,ds)
      INTEGER             , intent(in)  :: nf,nabin
      REAL, dimension(:)  , intent(in)  :: e,ec,es,ec2,es2
      REAL, dimension(:,:), intent(inout) :: ds
      end subroutine mem_est
   end interface
   
   Timezone_data=0.
   
   !this is a guess, the true values are read from the mem file
   dt=10.
   tavg=180.
   lat=0.
   lon=0.
   U10=0.
   Udir=0.
   Curr=0.
   Currdir=0.
   depth=50.
   nabin=24
   space=' '
   ALLOCATE(THETA(0:nabin-1))
   THETA= (/ (pi*(0.5-2.*FLOAT(I)/FLOAT(nabin)), i=0,nabin-1) /) 
   
   !!!!!!!!!! WARNING THIS PROGRAM WORKS ONLY WHEN THE AVERAGING TIME IS 1 hour
   !!!!!!!!!! AND THE DATA IS SUPPOSED TO BE SAMPLED AT LEAST ONCE AN HOUR
   !!!!!!!!!! THIS PROGRAM WILL FILL IN WITH ZEROS THE HOURS WITHOUT ANY DATA
   !!!!!!!!!! THIS NEED TO BE FIXED IN ORDER TO PROCESS THE FRF 8M F-theta SPECTRA   
   
   time=1
   filename2='DWIROISE_w3grid.obs'
   WRITE(6,*) 'Name of input file:'
   READ(*,'(a)') filename1
   WRITE(6,*) 'Name of output file (time averaged):'
   READ(*,*) filename2
   WRITE(6,*) 'Average time in minutes:'
   READ(*,*) tavg
   WRITE(6,*) 'Water depth:'
   READ(*,*) depth
   !WRITE(*,*) 'Enter file name of input file:'
   !READ(*,*) filename1
   !WRITE(6,*) '#',filename1
   !WRITE(*,*) 'Enter file name of output file:'
   !READ(*,*) filename2
   OPEN(3,FILE=filename1,status='old')
   OPEN(4,FILE=filename2,status='unknown')
   OPEN(5,FILE='avgbulk.obs',status='unknown')

   READ(3,*) entier
   WRITE(4,*) entier
   WRITE(5,*) '0'
   READ(3,*) nfband,np,npar

   IF (npar.EQ.5) OPEN(8,FILE='avgobs.spec',status='unknown')
   OPEN(9,FILE='average.txt',status='unknown')

   ALLOCATE(buoyname(np))
   ALLOCATE(freq(nfband))
   ALLOCATE(freq1(nfband))
   ALLOCATE(freq2(nfband))

   WRITE(4,'(3I4)') nfband,np,npar
   IF (npar.EQ.5) THEN 
      WRITE(5,'(3I4)') 1,np,9
   ELSE
      WRITE(5,'(3I4)') 1,np,5
   ENDIF
   READ(3,*) freq1(:)
   READ(3,*) freq2(:)
   WRITE(4,'(5E16.8)') freq1(:)
   WRITE(4,'(5E16.8)') freq2(:)
   WRITE(5,'(5E16.8)') freq1(1)
   WRITE(5,'(5E16.8)') freq2(nfband)
   freq(:)=0.5*(freq1(:)+freq2(:))
   IF (np.GT.1) THEN 
      WRITE(format1,'(a,I2,a)') "(a",np*8,")"
   ELSE 
      WRITE(format1,'(a,I1,a)') "(a",np*8,")"
   ENDIF
   IF (npar.GT.1) THEN 
      WRITE(format2,'(a,I2,a)') "(a",npar*8,")"
   ELSE
      WRITE(format2,'(a,I1,a)') "(a",npar*8,")"
   ENDIF
   READ(3,format1) buoyname(1:np)  !point names
   WRITE(4,format1) buoyname(1:np)
   WRITE(5,format1) buoyname(1:np)
   READ(3,format2) chaine  !field names
   WRITE(4,format2) chaine
   
   IF (npar.EQ.5) THEN 
      WRITE(5,*) 'Hs      hs10    fp       fp10      f02     th1p    sth1p   th2p    sth2p   '
   ELSE
      WRITE(5,*) 'Hs      hs10    fp       fp10      f02     '
   ENDIF
   
   ALLOCATE(df(nfband))
   ALLOCATE(Ef(nfband))
   ALLOCATE(ds(nabin,nfband))
   ALLOCATE(ds2(nabin,nfband))
   ALLOCATE(th1(nfband))
   ALLOCATE(sth1(nfband))
   ALLOCATE(th2(nfband))
   ALLOCATE(sth2(nfband))
   ALLOCATE(m1(nfband))
   ALLOCATE(m2(nfband))
   ALLOCATE(a1(nfband))
   ALLOCATE(b1(nfband))
   ALLOCATE(a2(nfband))
   ALLOCATE(b2(nfband))
   df(:)=freq2(:)-freq1(:)
   
   



   finished=0
   iostat=0
   iminc=0
   ALLOCATE(Efa(nfband))
   ALLOCATE(a1a(nfband))
   ALLOCATE(b1a(nfband))
   ALLOCATE(a2a(nfband))
   ALLOCATE(b2a(nfband))

   ALLOCATE(th1a(nfband))
   ALLOCATE(sth1a(nfband))
   ALLOCATE(th2a(nfband))
   ALLOCATE(sth2a(nfband))
   
   irec=0
   WRITE(8,'(A24,3I6,A33)') "'OBSERVATIONS         '",nfband,nabin,np, &
                           "'SHOM                          '"
   WRITE(8,'(8G10.3)') freq(1:nfband)
   WRITE(8,'(7G11.3)') theta(0:nabin-1)  !read directions in radians

   IF (npar.EQ.5) THEN 
   DO WHILE (icode.GE.0)   
      READ(3,*, iostat=icode)  iyc,imc,idc,ihc,iminc,Timezone_data
      Jdayc=julday(idc,imc,iyc)
      CALL caldat(Jdayc,ido,imo,iyo)
      iho=FLOOR(REAL(INT((REAL(ihc*60+iminc)+tavg/120)/tavg))*tavg/60)
      WRITE(6,*) 'Processing:',irec,w,iyc,imc,ihc,iminc,iho,INT((REAL(ihc*60+iminc)+tavg/120)/tavg)
      imino=0
      IF (irec.EQ.0) THEN 
         Jdayi=Jdayc
      ENDIF

      irec=irec+1
      IF (icode.GE.0) THEN 
         READ(3,*)  Ef(:)
         READ(3,*)  th1(:)
         READ(3,*)  sth1(:)
         READ(3,*)  th2(:)
         READ(3,*)  sth2(:)

         a1(:)=COS(th1(:)*d2r)
         b1(:)=SIN(th1(:)*d2r)
         m1(:)=SQRT(ABS((1.-0.5*(sth1(:)*d2r)**2)))
         DO I=1,nfband
            IF (m1(i).NE.0) THEN 
               a1(I)=a1(I)*m1(I)
               b1(I)=b1(I)*m1(I)
            ENDIF
         ENDDO

         a2(:)=COS(th2(:)*d2r)
         b2(:)=SIN(th2(:)*d2r)
         m2(:)=SQRT(ABS((1.-0.5*(sth2(:)*d2r)**2)))
         DO I=1,nfband
            IF (m2(i).NE.0) THEN 
               a2(I)=a2(I)*m2(I)
               b2(I)=b2(I)*m2(I)
            ENDIF
         ENDDO

         IF (irec.EQ.1) THEN 
            w=iminc
            Efa(:)=ef(:)*w
            a1a(:)=a1(:)*ef(:)*w
            b1a(:)=b1(:)*ef(:)*w
            a2a(:)=a2(:)*ef(:)*w
            b2a(:)=b2(:)*ef(:)*w
         ELSE

            w2=60*(24*REAL(Jdayc-Jdaylast)+(ihc-ihlast))+iminc-iminlast
            IF (w+w2.LT.tavg) THEN 
               w=w+w2
               Efa(:)=efa(:)+ef(:)*w2
               a1a(:)=a1a(:)+a1(:)*ef(:)*w2
               b1a(:)=b1a(:)+b1(:)*ef(:)*w2
               a2a(:)=a2a(:)+a2(:)*ef(:)*w2
               b2a(:)=b2a(:)+b2(:)*ef(:)*w2
            ELSE
               Efa(:)=(efa(:)+ef(:)*(tavg-w))/tavg
               a1a(:)=(a1a(:)+a1(:)*ef(:)*(tavg-w))/tavg  
               b1a(:)=(b1a(:)+b1(:)*ef(:)*(tavg-w))/tavg  

               a2a(:)=(a2a(:)+a2(:)*ef(:)*(tavg-w))/tavg
               b2a(:)=(b2a(:)+b2(:)*ef(:)*(tavg-w))/tavg
               th1a(:)=ATAN2(b1a(:),a1a(:))           
               th2a(:)=ATAN2(b2a(:),a2a(:))/2.          
              DO J=1,nfband
               IF (Efa(J) /= 0) THEN
                     sth1a(J)=SQRT(ABS(2*(1-(a1a(J)*cos(TH1a(J)) &
                         +b1a(J)*sin(TH1a(J))))))/(d2r*Efa(J))
                     sth2a(J)=SQRT(ABS(0.5*(1-(a2a(J)*cos(2*TH2a(J)) &
                         +b2a(J)*sin(2*TH2a(J))))))/(d2r*Efa(J))
               ELSE
                     sth1a(J)=0.
                     sth2a(J)=0.
               ENDIF
               ENDDO

               WRITE(4,'(I5,4I3,f6.2)')  &
                     iyo,imo,ido,iho,imino,Timezone_data
               WRITE(4,'(5E16.8)')  Efa(:)
               WRITE(4,'(5E16.8)')  th1a(:)/d2r
               WRITE(4,'(5E16.8)')  sth1a(:)
               WRITE(4,'(5E16.8)')  th2a(:)/d2r
               WRITE(4,'(5E16.8)')  sth2a(:)
               
               emax=0.
               DO J=1,nfband 
                  IF (Efa(J).GT.emax) THEN 
                     Jfp=J
                     emax=Efa(J)
                  ENDIF
               ENDDO
               emax2=0.
               DO J=1,11 
                  IF (Efa(J).GT.emax2) THEN 
                     kfp=J
                     emax2=Efa(J)
                  ENDIF
               ENDDO
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               IF (emax.GT.0.) f02=SQRT((SUM(freq(:)**2*Efa(:)*df(:))/SUM(Efa(:)*df(:))))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               Hs=4*SQRT(abs(SUM(Efa(:)*df(:))))
               Hs10=4*SQRT(abs(SUM(Efa(1:12)*df(1:12))))
	       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               IF (Jfp.NE.nfband) THEN 
! write(6,*) 'acte 9a'
         fp=freq(Jfp)
! write(6,*) 'acte 9b'
         th1p=th1a(Jfp)/d2r
! write(6,*) 'acte 9c'
         IF (th1p.LT.0) th1p=th1p+360. 
         th2p=th2a(Jfp)/d2r
         IF (th2p.LT.0) th2p=th2p+360. 
         sth1p=sth1a(Jfp)
         sth2p=sth2a(Jfp)
               ELSE
 ! write(6,*) 'acte 10'
        fp=0.
         th1p=0.
         th2p=0.
         sth1p=0.
         sth2p=0.
               ENDIF 
	       fp10=freq(kfp)
               WRITE (5,'(I5,4I3,f6.2,13E16.8)') &
                  iyo,imo,ido,iho,imin,Timezone_data, &
                  Hs,Hs10,fp,fp10,f02,th1p,sth1p,th2p,sth2p
               WRITE (9,'(I3,I5,3I3,6E16.8)') &
                  time,iyo,imo,ido,iho,Hs,Hs10,fp,fp10,th1p,sth1p
		time=time+1
!  write(6,*) 'acte 11'
                CALL mem_est(nfband,nabin,Efa,a1a,b1a,a2a,b2a,ds)
!  write(6,*) 'acte 12'
               ! this is corrected for the change of angular reference
               !CALL mem_est(nfband,nabin,Efa,-1.*b1a,-1.*a1a,-1.*a2a,b2a,ds)
               ds2(1,:)=ds(nabin,:)
               ds2(2:nabin,:)=ds(1:nabin-1,:)
               WRITE(8,'(I4,2I2,A1,3I2)',iostat=icode) iyo,imo,ido,space,iho,imino,0
               WRITE(8,'(A1,A10,A1,2F7.2,F10.1,F7.2,F6.1,F7.2,F6.1)') &
                  "'",buoyname(1),"'",lat,lon,depth,U10, &
                  Udir,Curr,Currdir
              WRITE(8,'(7G11.3)') TRANSPOSE(ds2)*nabin/(2*pi)
!WRITE(6,*) 'Ef:',Efa
!WRITE(6,*) 'DD:',SUM(TRANSPOSE(ds2),2)


               w=w+w2-tavg
               Efa(:)=ef(:)*w
               a1a(:)=a1(:)*w
               b1a(:)=b1(:)*w
               a2a(:)=a2(:)*w
               b2a(:)=b2(:)*w
            ENDIF
         ENDIF
      ENDIF
      Jdaylast=Jdayc
      ihlast=ihc
      iminlast=iminc
   ENDDO      
   ELSE
   
      DO WHILE (icode.GE.0)   
      READ(3,*, iostat=icode)  iyc,imc,idc,ihc,iminc,Timezone_data
      Jdayc=julday(idc,imc,iyc)
      CALL caldat(Jdayc,ido,imo,iyo)
      iho=FLOOR(REAL(INT((REAL(ihc*60+iminc)+tavg/120)/tavg))*tavg/60)
      WRITE(6,*) 'Processing:',irec,w,iyc,imc,ihc,iminc,iho,INT((REAL(ihc*60+iminc)+tavg/120)/tavg)
      imino=0
! write(6,*) 'acte 1'
      IF (irec.EQ.0) THEN 
         Jdayi=Jdayc
      ENDIF
! write(6,*) 'acte 2'
      irec=irec+1
      IF (icode.GE.0) THEN 
         READ(3,*)  Ef(:)
! write(6,*) 'acte 5'
         IF (irec.EQ.1) THEN 
            w=iminc
            Efa(:)=ef(:)*w
         ELSE
! write(6,*) 'acte 6'
            w2=60*(24*REAL(Jdayc-Jdaylast)+(ihc-ihlast))+iminc-iminlast
            IF (w+w2.LT.tavg) THEN 
               w=w+w2
               Efa(:)=efa(:)+ef(:)*w2
           ELSE
               Efa(:)=(efa(:)+ef(:)*(tavg-w))/tavg
               WRITE(4,'(I5,4I3,f6.2)')  &
                     iyo,imo,ido,iho,imino,Timezone_data
               WRITE(4,'(5E16.8)')  Efa(:)
               
               emax=0.
               DO J=1,nfband 
                  IF (Efa(J).GT.emax) THEN 
                     Jfp=J
                     emax=Efa(J)
                  ENDIF
               ENDDO
               emax2=0.
               DO J=1,11 
                  IF (Efa(J).GT.emax2) THEN 
                     kfp=J
                     emax2=Efa(J)
                  ENDIF
               ENDDO
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               IF (emax.GT.0.) f02=SQRT((SUM(freq(:)**2*Efa(:)*df(:))/SUM(Efa(:)*df(:))))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               Hs=4*SQRT(abs(SUM(Efa(:)*df(:))))
               Hs10=4*SQRT(abs(SUM(Efa(1:12)*df(1:12))))
	       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               IF (Jfp.NE.nfband) THEN 
         fp=freq(Jfp)
               ELSE
        fp=0.
               ENDIF 
	       fp10=freq(kfp)
               WRITE (5,'(I5,4I3,f6.2,13E16.8)') &
                  iyo,imo,ido,iho,imin,Timezone_data, &
                  Hs,Hs10,fp,fp10,f02
               WRITE (9,'(I3,I5,3I3,6E16.8)') &
                  time,iyo,imo,ido,iho,Hs,Hs10,fp,fp10
		time=time+1
            ENDIF
         ENDIF
      ENDIF
      Jdaylast=Jdayc
      ihlast=ihc
      iminlast=iminc
   ENDDO      
   ENDIF
   CLOSE(3)
   CLOSE(4)
   CLOSE(5)
   IF (npar.EQ.5) CLOSE(8)
   CLOSE(9)
END
  
