MODULE C_CONSTS
PUBLIC
   REAL,    PARAMETER :: pi = 3.14159265,tpi=6.2831853,d2r=.01745329
   REAL,    PARAMETER :: e = 2.7182818284590455
   REAL,    PARAMETER :: expmin=-20.
   REAL,    PARAMETER :: tpiinv = 1. / tpi   
   REAL,    PARAMETER :: max15bit=32768.,max8bit=256.,max7bit=128.
   REAL,    PARAMETER :: max14bit=16384.
   INTEGER, PARAMETER :: i15bit=32768,i8bit=256,i7bit=128,i14bit=16384
   REAL,    PARAMETER :: g =  9.81          !gravity       : SI  (m/s^2)
   REAL,    PARAMETER :: rhow =  1026.      !water density : SI  (kg/m^3)
   REAL,    PARAMETER :: rhoa =  1.293      !air density (0 C): SI  (kg/m^3)
   REAL,    PARAMETER :: kappa = 0.4        !Von Karman's constant
   REAL,    PARAMETER :: nu  = 3.E-6        !Salt water kinematic viscosity
   REAL, DIMENSION(:),ALLOCATABLE  ::invxtanhx
   INTEGER ntanh

   INTEGER, PARAMETER :: NTAB=1001
   REAL           , DIMENSION(NTAB) :: ELLIPE,ELLIPK

END MODULE C_CONSTS
!*********************************************************************
!*********************************************************************
MODULE c_disper
public 
   CONTAINS

!*********************************************************************
	SUBROUTINE DISP(depth,freq,k,c,cg)
!*********************************************************************
   USE C_CONSTS
   IMPLICIT NONE
   REAL,                INTENT(in)     ::depth,freq  !depth (m) and frequency (Hz)
   REAL,                INTENT(out)    ::k,c,cg      !wave number (m-1),phase speed (m/s)
   REAL dep,sigma,X,Y,dy
   INTEGER I
   
   dep=depth
   IF (depth.LT.0.1) dep=0.1
   sigma=freq*tpi
   Y=(sigma**2*DEP)/G
   IF (Y.LE.(15./REAL(ntanh+1))) THEN 
      X=SQRT(Y)
   ELSE
      IF (Y.GE.14.5) THEN 
         X=Y 
      ELSE
         I=INT(Y*REAL(ntanh+1)/15.)-1
         dy=mod(Y*REAL(ntanh+1)/15.,1.0)
         X=invxtanhx(I)*(1-dy)+invxtanhx(I+1)*dy
      ENDIF
   ENDIF
   k=X/DEP
   c=sigma/k
   cg=g*((y/x)+x*(1-(y/x)**2))/(2.*sigma)
   
   RETURN
   END SUBROUTINE DISP


!*********************************************************************
   SUBROUTINE TABULATE
!*********************************************************************
   USE C_CONSTS
   IMPLICIT NONE

   REAL X,Y,H,F,FD
   INTEGER I,J
   ntanh=10000
   ALLOCATE(invxtanhx(0:ntanh))
   DO J=0,ntanh
      Y=15.*REAL(J+1)/REAL(ntanh+1)
      X=Y
   	F=1.
      DO 10 I=1,100
	   	H=TANH(X)
		   F=Y-X*H
         IF (ABS(F) .LT. 1.E-6) GO TO 20
	   	FD=-H-X/COSH(X)**2
   		X=X-F/FD
10	   CONTINUE
   	PRINT *,' SUBROUTINE DISP DOES NOT CONVERGE !!! '
20	   invxtanhx(J)=X
      ENDDO
      RETURN
	   END SUBROUTINE TABULATE

!**********************************************************************
SUBROUTINE tfact(depth,f,fact,fact2)
!... calculates the transformation factor (dk/dko)*(cgo/cg)
!... to avoid floating point exception in sinh, max value of sig=75
!... first get wave no. and group velocity at transformation site
   IMPLICIT NONE
   REAL,                INTENT(in)     ::depth,f
   REAL,                INTENT(out)    ::fact,fact2
   REAL d,dk,c,cg,dk0,c0,cg0
   
   CALL disp(depth,f,dk,c,cg)
   d=800.
   CALL disp(d,f,dk0,c0,cg0)
   fact=(dk/dk0)*(cg0/cg)
   fact2=(dk/cg)
   RETURN
END SUBROUTINE tfact

END MODULE c_disper
!*********************************************************************
	SUBROUTINE DISP2(f,WNUM,DEPTH)
! For a given frequency SIGMA and depth DEPTH this routine computes
! the wave number WNUM  : SIGMA**2=G*WNUM*tanh(WNUM*DEPTH) (G is the
! acceleration of gravity) with the Newton-Raphson method.
!*********************************************************************
! Code from Thomas Herbers, 1999.
!*********************************************************************
	INTEGER I
   REAL,    PARAMETER :: pi = 3.14159265,tpi=6.2831853,d2r=.01745329
   REAL SIGMA,WNUM,DEPTH,X,Y,G,H,FR,FD
	SIGMA=tpi*f
   G=9.81
	Y=(SIGMA**2*DEPTH)/G
	X=Y
	DO 10 I=1,100
		H=TANH(X)
		FR=Y-X*H
		IF (ABS(F) .LT. 0.000001) GO TO 20
		FD=-H-X/COSH(X)**2
		X=X-FR/FD
10	CONTINUE
	PRINT *,' SUBROUTINE DISP DOES NOT CONVERGE !!! '
20	WNUM=X/DEPTH
   RETURN
	END


!****************************************************************************    
MODULE c_caldat
CONTAINS
!*********************************************************************
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 FUNCTION JULDAY

!*********************************************************************
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 SUBROUTINE CALDAT
END MODULE c_caldat
!*********************************************************************
 !  -convert big_endian
PROGRAM makebulkanl
! this programs reads *.spc files from WW3 and generates bulk parameters
! Hs, Tp, spread, as well as the frequency dependent E(f), th(f) and sth(f)
!*********************************************************************
USE C_CONSTS
USE C_CALDAT
USE C_DISPER
   IMPLICIT NONE

   INTEGER     I,IB,IE,J,II,i1,i2,Ntime,idp1,idp2,id1,id2,im1,im2,nread,nfiles
   INTEGER     np,npar,PIND,IFILE,found,IF02MAX
   INTEGER     iy,im,im0,id,id0,ih,ih2,imin,isec,icode,icode2,icode3,iostat,finished
   INTEGER     iyo,imo,ido,iho
   INTEGER(KIND=4) Jdays,Jdayend,Jdayfile,Jday,Jdayo
   REAL   , DIMENSION(:),     ALLOCATABLE :: freq,freq1,freq2,df,theta
   REAL   , DIMENSION(:,:),     ALLOCATABLE :: Efth,th1,th2,sth1,sth2
   REAL, DIMENSION(:,:),   ALLOCATABLE :: Ef,a1,a2,b1,b2
   REAL, DIMENSION(:),   ALLOCATABLE :: lat,lon,depth,Hs,fp,f02,f0m1,U10,Udir,Curr,Currdir
   REAL, DIMENSION(:),   ALLOCATABLE :: th1p,th2p,sth1p,sth2p,th1tot,sth1tot
   INTEGER, DIMENSION(:),  ALLOCATABLE ::  Jfp
   CHARACTER(120) path,filename,filei,name,string1,string2
   CHARACTER(120), DIMENSION(:),     ALLOCATABLE :: filenames
   CHARACTER(20) , DIMENSION(:),     ALLOCATABLE ::  buoyname
   CHARACTER(15) buoyname1
   CHARACTER(20) buoyname2
   CHARACTER(5) buoyshort,buoyshort2
   CHARACTER(5), DIMENSION(300) :: buoyshort_ECMWF
   REAL , DIMENSION(300) :: XE,YE
   INTEGER, DIMENSION(0:300):: PERIOD_TYPE

   INTEGER nECMWF, START, ICODELIST
   CHARACTER space
   REAL m1,m2,a1p,a2p,b1p,b2p,Ep,a1tot,b1tot,FOUT
   REAL dtheta
   INTEGER nfband,nabin,ip,nchar
   REAL Eband,Etot,emax
   REAL Timezone_data,Time_zone_output
   
   LOGICAL         :: FILEW = .TRUE.
   CHARACTER*10, PARAMETER:: VERBPT = 'III  1.03 '
   CHARACTER*32, PARAMETER::                             &
                        IDSTR = 'WAVEWATCH III BOUNDARY DATA FILE'

   REAL , DIMENSION(:),   ALLOCATABLE :: XBPO,YBPO
   REAL , DIMENSION(:,:),   ALLOCATABLE :: RDBPO
   INTEGER, DIMENSION(:,:),   ALLOCATABLE :: IPBPO
   
   Timezone_data=0.

   path=''
   !path='/local/tmp/ardhuin/DATA/BOUEES/62163/'
  
   nfiles=4
   ALLOCATE(filenames(nfiles))
   OPEN(101,FILE='yearmonth.txt',status='old',iostat=icode2)
   READ(101,'(I4,I2)') iy,im0
   CLOSE(101)
   II=200
   im2=MOD(im0,10)
   im1=im0/10
   WRITE(filename,'(A5,I4,2I1,A4,2I1,A4)') 'SHOM_',iy,im1,im2,'.dat'
   OPEN(II,FILE=filename,status='old',iostat=icode)
   IF (icode.NE.0) THEN 
      OPEN(II,FILE=filename,status='new',iostat=icode)
   ELSE
      CLOSE(II)
      OPEN(II,FILE=filename,status='old',position='append')
   ENDIF
!
! Reads ECMWF buoy list to check positions and get period type
!
   OPEN(101,FILE='ECMWF_buoylist.txt',status='old')
   ICODE=0
   I=0
   DO WHILE (icode.EQ.0) 
      I=I+1
      READ(101,'(A5,F8.3,X,F8.3,X,I1)',iostat=ICODE) buoyshort_ECMWF(I), XE(I), YE(I), &
      PERIOD_TYPE(I) 
   ENDDO
   CLOSE(101)
   nECMWF=I-1
   WRITE(6,*) 'ECMWF buoy list read, number of buoys:',nECMWF
   PERIOD_TYPE(0)=0
!
! Starts looping on model buoy results
!
   START=1
   OPEN(102,FILE='buoylist.txt',status='old')
   ICODELIST=0
   IB=0
   DO WHILE (ICODELIST.EQ.0) 
      READ(102,'(A5)',iostat=ICODELIST) buoyname1
      IF(ICODELIST.EQ.0) THEN 
      IB=IB+1
      nchar=LEN_TRIM(buoyname1)
      READ(buoyname1,'(A5)') buoyshort
      WRITE(6,*) '#',buoyshort,'#',buoyshort_ECMWF(IB)
      IE=0
      DO I=1,nECMWF
         IF (buoyshort_ECMWF(I).EQ.buoyshort) THEN 
            IE=I
         ENDIF
      ENDDO
      WRITE(6,*) 'Current buoy is ECMWF #',IE,IB
      WRITE(filenames(4),'(A,A,A)') 'SPECTRA/ww3.',buoyname1(1:nchar),'.spec'
      WRITE(filenames(3),'(A,A,A)') 'SPECTRA/aww3.',buoyname1(1:nchar),'.spec'
      WRITE(filenames(2),'(A,A,A)') 'SPECTRA/medw3.',buoyname1(1:nchar),'.spec'
      WRITE(filenames(1),'(A,A,A)') 'SPECTRA/gww3.',buoyname1(1:nchar),'.spec'
      !IF (IB.EQ.3.OR.IB.EQ.4.OR.IB.EQ.177.OR.IB.EQ.193.OR.IB.EQ.174) THEN !special case for 62001 and 62163
      !   WRITE(filenames(1),'(A,A,A)') 'SPECTRA/gww3.',buoyname1(1:nchar),'.spec'
      !ELSE
      !   WRITE(filenames(1),'(A,A,A)') 'SPECTRA/aww3.',buoyname1(1:nchar),'.spec'
      !ENDIF
      found=0
      IFILE=1
      DO WHILE (found.EQ.0.AND.IFILE.LE.4)
         OPEN(101,FILE=filenames(Ifile),status='old',iostat=icode)
         IF (icode.EQ.0) THEN 
            found=1
         ELSE
            IFILE=IFILE+1
         ENDIF
      ENDDO
      IF (found.EQ.1) THEN 
         filei=filenames(IFILE)
         WRITE(6,*) 'READING FILE:',icode,filei(1:40)
      READ(101,'(A24,3I6,A33)') string1,nfband,nabin,np,string2
      
      ALLOCATE(Efth(nfband,nabin),freq(nfband),freq1(nfband))
      ALLOCATE(freq2(nfband),df(nfband),Ef(np,nfband),theta(nabin))
   
      IF (START.EQ.1) THEN 
         ALLOCATE(buoyname(np))
         ALLOCATE(lat(np))
         ALLOCATE(lon(np))
         ALLOCATE(depth(np))
         ALLOCATE(fp(np))
         ALLOCATE(Jfp(np))
         ALLOCATE(Hs(np))
         ALLOCATE(f02(np))
         ALLOCATE(f0m1(np))
         ALLOCATE(U10(np))
         ALLOCATE(Udir(np))
         ALLOCATE(Curr(np))
         ALLOCATE(Currdir(np))
         START=2
      ENDIF
      READ(101,*) freq(1:nfband)
      IF02MAX=1
      DO J=1,NFBAND
         IF(freq(J).LT.0.42) IF02MAX=J
      ENDDO 
      READ(101,*) theta(1:nabin)  !read directions in radians
      WRITE(6,*) 'IF02MAX:',IF02MAX,FREQ(IF02MAX)

      freq2(1:nfband-1)=sqrt(freq(1:nfband-1)*freq(2:nfband))
      freq1(2:nfband)=freq2(1:nfband-1)
      freq2(nfband)=freq(nfband)*sqrt(freq(nfband)/freq(nfband-1))
      freq1(1)=freq(1)*sqrt(freq(1)/freq(2))
      df=freq2-freq1

      dtheta=2*pi/nabin
      Ntime=0
      DO WHILE (ICODE.GE.0)
         READ(101,'(I4,2I2,A1,3I2)',iostat=icode) iy,im,id,space,ih,imin,isec
         IF (Ntime.EQ.0) THEN 
            Jdays=julday(id,im,iy)+1
         ENDIF
         Jday=julday(id,im,iy)
         ih2=0 ! THESE ARE "ANALYSES" , this is ttt: forecast horizon
         IF (ICODE.GE.0) THEN 
            Ntime=Ntime+1
            IP=1
            READ(101,'(A1,A10,A1,2F7.2,F10.1,F7.2,F6.1,F7.2,F6.1)') &
            space,buoyname(IP),space,lat(IP),lon(IP),depth(IP),U10(IP), &
            Udir(IP),Curr(IP),Currdir(IP)
            buoyname2=buoyname(IP)
            df(1)=freq(2)-freq(1)
            DO J=2,nfband-1 
               df(J)=0.5*(freq(J+1)-freq(J-1))
            ENDDO
            df(nfband)=freq(nfband)-freq(nfband-1)
            freq1=freq-0.5*df
            freq2(1:nfband-1)=freq1(2:nfband)
            freq2(nfband)=freq(nfband)+df(nfband)
            df=freq2-freq1

            READ(101,*) Efth(:,:)
      
            IF (mod(ih,6).EQ.0.AND.im.EQ.im0) THEN 
               Ef(Ip,1:nfband)=SUM(Efth(:,:),2)*dtheta
               emax=0.
               Jfp(IP)=nfband
               DO J=1,nfband 
                  IF (Ef(IP,J).GT.emax) THEN 
                     Jfp(IP)=J
                     emax=Ef(IP,J)
                  ENDIF
               ENDDO
               IF (emax.GT.0.) THEN 
                  f02(IP)=SQRT(SUM(freq(1:IF02MAX)**2 &
                     *Ef(IP,1:IF02MAX)*df(1:IF02MAX)) &
                     /SUM(Ef(IP,1:IF02MAX)*df(1:IF02MAX)))
                  f0m1(IP)=1./(SUM(Ef(IP,:)*df(:)/freq(:))/SUM(Ef(IP,:)*df(:)))
               ENDIF
               Hs(IP)=4*SQRT(SUM(Ef(IP,:)*df(:)))
               IF (Jfp(IP).NE.nfband) THEN 
                  J=Jfp(IP)
                  IF (J.GT.1) THEN 
                     Ep=SUM(Ef(IP,J-1:J+1)*df(J-1:J+1))
                     fp(IP)=SUM(Ef(IP,J-1:J+1)*freq(J-1:J+1)*df(J-1:J+1))/ Ep
                  ELSE
                     fp(IP)=freq(J)
                  ENDIF
               ELSE
                  fp(IP)=0.
               ENDIF 
               Etot=SUM(Ef(IP,:))

               IF (PERIOD_TYPE(IE).EQ.0)  THEN 
                  PIND=PERIOD_TYPE(IE)
                  IF (F02(IP).GT.0) THEN 
                     FOUT=1/F02(IP)
                  ELSE
                     FOUT=0.
                  ENDIF
               ELSE
                  PIND=PERIOD_TYPE(IE)
                  IF (FP(IP).GT.0) THEN 
                     FOUT=1/FP(IP)
                  ELSE
                     FOUT=0.
                  ENDIF
               ENDIF
               buoyshort2=buoyname2(1:5)
               IF (buoyshort2.NE.buoyshort) WRITE(6,*) '###### ',buoyshort, &
                  ' ',buoyshort2
               IF (buoyshort2.EQ.'BILBA') buoyshort2='62024'
               IF (buoyshort2.EQ.'BARES') buoyshort2='62082'
               IF (buoyshort2.EQ.'VILLA') buoyshort2='62083'
               WRITE(II,'(a5,X,a5,X,i3,X,i4,X,i2,X,i2,x,i2,X,F4.1,X,I3,X,F4.1,X,F4.1,A,I1)') & 
                   'SHOM ',buoyshort2,ih2,iy,im,id,ih,U10(IP), &
                   int(Udir(IP)+0.5),Hs(IP),FOUT,' **** *** **** **** ',PIND
            ENDIF
         ENDIF
      ENDDO
      DEALLOCATE(Efth,freq,freq1,freq2,df,Ef,theta) !,buoyname,lat,lon,depth,fp,Jfp,Hs,f02, &
       !f0m1,U10,Udir,Curr,Currdir)
      ELSE
         WRITE(6,*) 'ERROR: no file for buoy',buoyname1(1:nchar)
      ENDIF
   ENDIF
   ENDDO 
   CLOSE(II)
   CLOSE(201)
   CLOSE(202)
   CLOSE(203)
END
  
