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 makebulk
! 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
   INTEGER     np,npar,PIND,IFILE,found,IF02MAX
   INTEGER     iy,im,im0,im1,im2,id,ih,ih2,imin,isec,icode,icode2,icode3,iostat,finished
   INTEGER     iyd,imd,idd,ihd,imind
   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 :: depth,Hs,fp,f02,f0m1,U10,Udir,Curr,Currdir
   REAL :: depthd,Hsd,fpd,f02d,depthdo,Hsdo,fpdo,f02do
   CHARACTER(120) path,filename,filenameo,filei,name,string1,string2
   CHARACTER(120), DIMENSION(:),     ALLOCATABLE :: filenames
   CHARACTER(20) , DIMENSION(:),     ALLOCATABLE ::  buoyname
   CHARACTER(15) buoyname1
   CHARACTER(20) buoyname2,dummystring
   CHARACTER(5) source,buoyshort,buoyshort2
   CHARACTER(5), DIMENSION(300) :: buoyshort_ECMWF
   REAL , DIMENSION(300) :: XE,YE
   INTEGER, DIMENSION(0:300):: PERIOD_TYPE

   REAL FOUT,dt,dtmin
   REAL Timezone_data
   

   
   Timezone_data=0.

   path=''
   !path='/local/tmp/ardhuin/DATA/BOUEES/62163/'
  
   OPEN(101,FILE='yearmonth.txt',status='old')
   READ(101,'(I4,I2)') iy,im0
   CLOSE(101)
   im2=MOD(im0,10)
   im1=im0/10
   WRITE(filename,'(A5,I4,2I1,A4,2I1,A4)') 'SHOM_',iy,im1,im2,'.dat'
   OPEN(200,FILE=filename,status='old',iostat=icode)
   WRITE(filenameo,'(A10,I4,2I1,A4,2I1,A4)') 'SHOM_data_',iy,im1,im2,'.dat'
   OPEN(201,FILE=filenameo,status='unknown')
   DO WHILE (icode.EQ.0) 
      READ(200,'(a5,X,a5,X,i3,X,i4,X,i2,X,i2,x,i2,X,F4.1,X,I3,X,F4.1,X,F4.1,A20,I1)',iostat=icode) & 
                   source,buoyshort2,ih2,iy,im,id,ih,U10, &
                   Udir,Hs,FOUT,dummystring,PIND
      IF (icode.EQ.0) THEN 
      found=0
      !IF (ih2.EQ.0.AND.buoyshort.EQ.'62064') THEN 
      IF (buoyshort2.EQ.'62064') THEN 
         Jdayo=julday(id,im,iy)
         OPEN(102,FILE='avgbulk.obs',status='old')
         icode2=0
         DO I=1,6
           READ(102,'(A10)') path
         ENDDO
         dtmin=6.
         DO WHILE (icode2.EQ.0) 
           READ(102,*,iostat=icode2)  iyd,imd,idd,ihd,imind,Timezone_data,Hsd,fpd,f02d
!           READ(102,*,iostat=icode2)  iyd,imd,idd,ihd,imind,Timezone_data,Hsd,fpd,f02d
           WRITE(6,*) 'TEST1:',icode2,found,iy,im,id,ih,iyd,imd,idd,ihd, dt,dtmin,Hsd,fpd,f02d
           IF (icode2.EQ.0) THEN 
           Jday=julday(idd,imd,iyd)
           imin=0
           dt=ABS(REAL(Jday-Jdayo)*24.+ihd-ih+(imind-imin)/60.)
           IF (dt.LT.dtmin.AND.dt.LE.2) THEN 
           WRITE(6,*) 'TEST2:',icode2,found,iy,im,id,ih,iyd,imd,idd,ihd, dt,dtmin,Hsd,fpd,f02d
              found=1
              dtmin=dt
              Hsdo=Hsd
              fpdo=fpd
              f02do=f02d
           ENDIF
           ENDIF
        END DO
        CLOSE(102)
      ENDIF     
      IF (found.EQ.0) THEN 
         WRITE(201,'(a5,X,a5,X,i3,X,i4,X,i2,X,i2,x,i2,X,F4.1,X,I3,X,F4.1,X,F4.1,A20,I1)') & 
                   source,buoyshort2,ih2,iy,im,id,ih,U10, &
                   Udir,Hs,FOUT,dummystring,PIND
      ELSE
         WRITE(6,*) 'TEST:',iy,im,id,ih,buoyshort2,Hsdo,1./F02do,'#',dummystring,'#',dummystring(1:10),'#'
         WRITE(201,'(a5,X,a5,X,i3,X,i4,X,i2,X,i2,x,i2,X,F4.1,X,I3,X,F4.1,X,F4.1,A10,F4.1,X,F4.1,X,I1)') & 
                   source,buoyshort2,ih2,iy,im,id,ih,U10, &
                   Udir,Hs,FOUT,dummystring(1:10),Hsdo,1./F02do,PIND
      ENDIF
   ENDIF
   ENDDO
         
   CLOSE(200)
   CLOSE(201)
   CLOSE(202)
   CLOSE(203)
END
  
