#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      PROGRAM W3BOUNC
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           F. Ardhuin              |
!/                  |           M. Accensi              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         16-Oct-2016 |
!/                  +-----------------------------------+
!/
!/    24-May-2013 : Adaptation from ww3_bound.ftn    ( version 4.08 )
!/     1-Apr-2015 : Add checks on lat lon xfr        ( version 5.05 )
!/    11-May-2015 : Allow use of cartesian grids     ( version 5.08 )
!/    17-Aug-2016 : Bug correction on RDBPO          ( version 5.10 ) 
!/    20-Oct-2016 : Error statement updates          ( version 5.15 )
!/
!/    Copyright 2012-2013 National Weather Service (NWS),
!/       National Oceanic and Atmospheric Administration.  All rights
!/       reserved.  WAVEWATCH III is a trademark of the NWS. 
!/       No unauthorized use without permission.
!/
!  1. Purpose :
!
!     Combines spectra files into a nest.ww3 file for boundary conditions
!
!  2. Method :
!
!     Finds nearest points and performs linear interpolation
!
!     The initial conditions are written to the restart.ww3 using the
!     subroutine W3IORS. Note that the name of the restart file is set
!     in W3IORS.
!
!  3. Parameters :
!
!     Local parameters.
!     ----------------------------------------------------------------
!       NDSI    Int.  Input unit number ("ww3_assm.inp").
!       ITYPE   Int.  Type of data
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      STRACE    Subr.   Id.    Subroutine tracing.
!      NEXTLN    Subr.   Id.    Get next line from input filw
!      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
!      WAVNU1    Subr. W3DISPMD Solve dispersion relation.
!      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     None, stand-alone program.
!
!  6. Error messages :
!
!  7. Remarks :
!
!     - Can be used also to diagnose contents of nest.ww3 file
!       in read mode
!
!  8. Structure :
!
!     ----------------------------------------------------
!        1.a  Set up data structures.
!                            ( W3NMOD , W3NDAT , W3NOUT 
!                              W3SETG , W3SETW , W3SETO )
!          b  I-O setup.
!        ....
!        9.   Convert energy to action
!       10.   Write restart file.              ( W3IORS )
!     ----------------------------------------------------
!
!  9. Switches :
!
!     !/SHRD  Switch for shared / distributed memory architecture.
!     !/DIST  Id.
!
!     !/SHRD  Switch for message passing method.
!     !/MPI   Id.
!
!     !/S     Enable subroutine tracing.
!
!     !/O4    Output normalized 1-D energy spectrum.
!     !/O5    Output normalized 2-D energy spectrum.
!     !/O6    Output normalized wave heights (not MPP adapted).
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS
      USE W3WDATMD, ONLY: W3NDAT, W3SETW
      USE W3ADATMD, ONLY: W3NAUX, W3SETA
      USE W3ODATMD, ONLY: W3NOUT, W3SETO, FLBPI

      USE W3GDATMD, ONLY: NK, NTH, XFR, FR1, GNAME, W3NMOD, W3SETG,     &
                          NSEA, MAPSTA, XYB, GTYPE, XGRD, YGRD, X0, Y0, &
                          SX, SY, MAPSF, UNGTYPE, CLGTYPE, RLGTYPE, FLAGLL
      USE W3ODATMD, ONLY: NDSO, NDSE
      USE W3IOBCMD, ONLY: VERBPTBC, IDSTRBC
      USE W3IOGRMD, ONLY: W3IOGR
      USE W3TIMEMD
      USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE, DIST_SPHERE
      USE NETCDF
!/S      USE W3SERVMD, ONLY : STRACE

!/
      IMPLICIT NONE
!
!/MPI      INCLUDE "mpif.h"
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/

! character
      CHARACTER                           :: COMSTR*1, LINE*512
      CHARACTER*5                         :: INXOUT
      CHARACTER*10                        :: VERTEST  ! = 'III  1.03 '
      CHARACTER*32                        :: IDTST    != 'WAVEWATCH III BOUNDARY DATA FILE'
      CHARACTER*512, ALLOCATABLE          :: SPECFILES(:)
      CHARACTER*512                       :: FILENAME
      CHARACTER*14                        :: TIME_STR

! index
      INTEGER                             :: IX, IY, ISEA, I,JJ,IP,IP1,J,IT,              &
                                             NDSI,NDSM, NDSI2,NDSS,NDSB, NDSC, NDSTRC,    &
                                             NTRACE, NK1,NTH1,NT1, NSPEC1, NBI, NBI2,     &
                                             NKI, NTHI, NTI, NBO, NBO2, IERR, INTERP,     &
                                             ILOOP, IFMIN, IFMIN2, IFMAX, VERBOSE, IBO
!/S      INTEGER, SAVE           :: IENT = 0

! arrays
      INTEGER                             :: TIME2(2)
      REAL, ALLOCATABLE                   :: SPEC2D(:,:,:,:), LATS(:), LONS(:)
      REAL, ALLOCATABLE                   :: FREQ(:), THETA(:)
      REAL                                :: FR1I, TH1I
      REAL, ALLOCATABLE                   :: XBPI(:), YBPI(:), RDBPI(:,:),        &
                                             XBPO(:), YBPO(:), RDBPO(:,:),        &
                                             ABPIN(:,:), ABPIN2(:,:,:) 
      INTEGER, ALLOCATABLE                :: IPBPI(:,:), IPBPO(:,:)
      REAL                                :: DMIN, DIST, DMIN2, COS1, DLON, DLAT, DLA, DLO

! netCDF
      INTEGER                             :: iret, ICODE, varid(12)
      CHARACTER, ALLOCATABLE              :: station(:,:)
      INTEGER, ALLOCATABLE                :: ncid(:), dimid(:,:), dimln(:,:)
      DOUBLE PRECISION, ALLOCATABLE       :: times(:)
      CHARACTER(120)                      :: timeunits 
      INTEGER                             :: dateref, yref, mref, dref
      INTEGER                             :: TIME(2)
      DOUBLE PRECISION                    :: datecur, datetmp
      REAL                                :: is, imn, ih
      INTEGER                             :: id, im, iyy

!/
!/ ------------------------------------------------------------------- /


!/
! 1.  IO set-up.
!
      CALL W3NMOD ( 1, 6, 6 )
      CALL W3SETG ( 1, 6, 6 )
      CALL W3NDAT (    6, 6 )
      CALL W3SETW ( 1, 6, 6 )
      CALL W3NAUX (    6, 6 )
      CALL W3SETA ( 1, 6, 6 )
      CALL W3NOUT (    6, 6 )
      CALL W3SETO ( 1, 6, 6 )
!
      NDSI   = 10
      NDSB   = 33
      NDSC   = 44
      NDSM   = 20
      NDSS   = 30
!
      NDSTRC =  6
      NTRACE = 10
      CALL ITRACE ( NDSTRC, NTRACE )
!
!/S      CALL STRACE (IENT, 'W3BOUNC')
!

!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 2.  Read model definition file.
!
      CALL W3IOGR ( 'READ', NDSM )
      WRITE (NDSO,920) GNAME
!
! 3. Read input file
!
      OPEN(NDSI,FILE='ww3_bounc.inp',status='old')
      READ (NDSI,'(A)',END=2001,ERR=2002) COMSTR
      IF (COMSTR.EQ.' ') COMSTR = '$'
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*) INXOUT
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*) INTERP
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*) VERBOSE
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
!
      NBO2 = 0          
!
!     ILOOP = 1 to count NBO2
!     ILOOP = 2 to read the file names
!
      DO ILOOP = 1, 2
        OPEN (NDSS,FILE='ww3_bounc.scratch',FORM='FORMATTED',          &
              status='UNKNOWN')
        IF ( ILOOP .EQ. 1 ) THEN
          NDSI2 = NDSI
        ELSE
          NDSI2 = NDSS
          ALLOCATE(SPECFILES(NBO2))
          NBO2=0
        ENDIF
 
        NBO2=0
!       Read input file names
        DO 
          CALL NEXTLN ( COMSTR , NDSI2 , NDSE )
          READ (NDSI2,'(A512)') FILENAME
          JJ     = LEN_TRIM(FILENAME)
          IF ( ILOOP .EQ. 1 ) THEN
            BACKSPACE (NDSI)
            READ (NDSI,'(A)') LINE
            WRITE (NDSS,'(A)') LINE
          END IF
          IF (FILENAME(:JJ).EQ."'STOPSTRING'") EXIT 
          NBO2=NBO2+1
          IF (ILOOP.EQ.1) CYCLE 
          SPECFILES(NBO2)=FILENAME
        END DO 
!
! ... End of ILOOP loop
!
        IF ( ILOOP .EQ. 1 ) CLOSE ( NDSS) 

        IF ( ILOOP .EQ. 2 ) CLOSE ( NDSS, STATUS='DELETE' )
      END DO ! ILOOP = 1, 2
      CLOSE(NDSI)

!
! 4. Tests the reading of the file
!
      IF ( INXOUT.EQ.'READ') THEN 
        OPEN(NDSB,FILE='nest.ww3',FORM='UNFORMATTED',status='old')
        READ(NDSB) IDTST, VERTEST, NK1, NTH1, XFR, FR1I, TH1I, NBI
        NSPEC1  = NK1 * NTH1
        IF ( IDTST .NE. IDSTRBC ) GOTO 2003
        WRITE(NDSO,940) VERTEST
        WRITE(NDSO,941) IDTST
        IF (VERBOSE.EQ.1) WRITE(NDSO,'(A,2I5,3F12.6,I5)') 'NK,NTH,XFR, FR1I, TH1I, NBI :', & 
                    NK1,NTH1,XFR, FR1I, TH1I, NBI 
        ALLOCATE (XBPI(NBI),YBPI(NBI))
        ALLOCATE (IPBPI(NBI,4),RDBPI(NBI,4))
        READ(NDSB) (XBPI(I),I=1,NBI),                                 &
                   (YBPI(I),I=1,NBI),                                 &
                   ((IPBPI(I,J),I=1,NBI),J=1,4),                      &
                   ((RDBPI(I,J),I=1,NBI),J=1,4)
        IF (VERBOSE.GE.1) WRITE(NDSO,*)      'XBPI:',XBPI
        IF (VERBOSE.GE.1) WRITE(NDSO,*)      'YBPI:',YBPI
        IF (VERBOSE.GE.1) WRITE(NDSO,*)      'IPBPI:'
        DO I=1,NBI
          IF (VERBOSE.GE.1) WRITE(NDSO,*) I,' interpolated from:',IPBPI(I,1:4)
          IF (VERBOSE.GE.1) WRITE(NDSO,*) I,' with coefficient :',RDBPI(I,1:4)
        END DO
!
        READ (NDSB) TIME2, NBI2
        BACKSPACE (NDSB)
        ALLOCATE (ABPIN(NSPEC1,NBI2))
        IERR=0
        DO WHILE (IERR.EQ.0)
          READ (NDSB,IOSTAT=IERR) TIME2, NBI2
          IF (IERR.EQ.0) THEN 
            IF (VERBOSE.EQ.1) WRITE(NDSO,*)      'TIME2,NBI2:',TIME2, NBI2,IERR      
            DO IP=1, NBI2
              READ (NDSB,END=2003,ERR=2004) ABPIN(:,IP)
            END DO
          END IF
        END DO
        CLOSE(NDSB)
      END IF ! INXOUT.EQ.'READ'
!
!
      IF ( INXOUT.EQ.'WRITE') THEN 
!
! 5. Defines position of active boundary points
!
        NBO = 0
        DO ISEA=1,NSEA
          IX     = MAPSF(ISEA,1)
          IY     = MAPSF(ISEA,2) 
          IF (MAPSTA(IY,IX).EQ.2) THEN 
            NBO=NBO+1
          END IF
        END DO
        ALLOCATE(XBPO(NBO),YBPO(NBO))
        ALLOCATE (IPBPO(NBO,4),RDBPO(NBO,4))
        IBO=0
        DO ISEA=1,NSEA
          IX     = MAPSF(ISEA,1)
          IY     = MAPSF(ISEA,2)
          IF (MAPSTA(IY,IX).EQ.2) THEN 
            IBO=IBO+1
            SELECT CASE ( GTYPE )
            CASE ( RLGTYPE )
              XBPO(IBO)=X0+SX*(IX-1)
              YBPO(IBO)=Y0+SY*(IY-1)
            CASE ( CLGTYPE )
              XBPO(IBO)= XGRD(IY,IX)
              YBPO(IBO)= YGRD(IY,IX)
            CASE (UNGTYPE) 
              XBPO(IBO)= XYB(IX,1)
              YBPO(IBO)= XYB(IX,2)
            END SELECT !GTYPE
          END IF
        END DO

        OPEN(NDSB,FILE='nest.ww3',FORM='UNFORMATTED',status='unknown')
        ALLOCATE(dimid(NBO2,3),dimln(NBO2,3),ncid(NBO2))

        ALLOCATE(lats(NBO2),lons(NBO2),station(16,NBO2))

        DO IP=1,NBO2
          ! open file
          OPEN(NDSC,FILE=trim(SPECFILES(IP)),FORM='UNFORMATTED',      &
               status='old',iostat=ICODE)
          IF (ICODE.NE.0) THEN 
            WRITE(NDSO,950) trim(SPECFILES(IP))
            lons(IP)=-999.
            lats(IP)=-999.
            STOP
          END IF
 
          iret=nf90_open(trim(SPECFILES(IP)),nf90_nowrite,ncid(IP))
          WRITE(6,*) 'Opening file:',trim(SPECFILES(IP))
          call check_err(iret)
          
          ! dimensions
          iret=nf90_inq_dimid(ncid(IP),'time',dimid(IP,1))
          call check_err(iret)
          iret=nf90_inq_dimid(ncid(IP),'frequency',dimid(IP,2))
          call check_err(iret)
          iret=nf90_inq_dimid(ncid(IP),'direction',dimid(IP,3))
          call check_err(iret)
          iret=nf90_inquire_dimension(ncid(IP),dimid(IP,1),len=dimln(IP,1))
          call check_err(iret)
          iret=nf90_inquire_dimension(ncid(IP),dimid(IP,2),len=dimln(IP,2))
          call check_err(iret)
          iret=nf90_inquire_dimension(ncid(IP),dimid(IP,3),len=dimln(IP,3))
          call check_err(iret)

          NTI=dimln(IP,1)
          NKI=dimln(IP,2)
          NTHI=dimln(IP,3)

          IF (IP.EQ.1) THEN
            NT1=NTI
            NK1=NKI
            NTH1=NTHI
            NSPEC1  = NK1 * NTH1
            ALLOCATE(times(NT1))
            ALLOCATE (FREQ(NK1),THETA(NTH1))
            ALLOCATE (SPEC2D(NTH1,NK1,NT1,NBO2))
            ALLOCATE (ABPIN2(NK*NTH1,NT1,NBO2))

            ! read time variable
            iret=nf90_inq_varid(ncid(IP), 'time', varid(1))
            iret=nf90_get_var(ncid(IP), varid(1), times(:))
            iret=nf90_get_att(ncid(IP), varid(1),'units',timeunits)

            ! time reference
            IF (index(timeunits, "seconds").NE.0) THEN
              read(timeunits(15:18),'(I4.4)') yref
              read(timeunits(20:21),'(I2.2)') mref
              read(timeunits(23:24),'(I2.2)') dref
            ELSE IF (index(timeunits, "days").NE.0) THEN
              read(timeunits(12:15),'(I4.4)') yref
              read(timeunits(17:18),'(I2.2)') mref
              read(timeunits(20:21),'(I2.2)') dref
            END IF
            dateref=julday(dref,mref,yref)
          ELSE 
            IF (NKI.NE.NK1.OR.NTHI.NE.NTH1) THEN 
               WRITE(NDSE,*) 'ERROR: inconsistent spectral dimension',NKI,NTHI,'in file ', &
                trim(SPECFILES(IP)),' compared to ',NK1,NT1,' in first file'
               STOP 
               END IF  
            END IF

          ! position variables : lon/lat or x/y
          IF ( FLAGLL ) THEN
            !WRITE (NDSO,951)
            iret=nf90_inq_varid(ncid(IP), 'latitude', varid(2))
            call check_err(iret)
            iret=nf90_get_var(ncid(IP), varid(2), lats(IP))
            call check_err(iret)
            iret=nf90_inq_varid(ncid(IP), 'longitude', varid(3))        
            call check_err(iret)
            iret=nf90_get_var(ncid(IP), varid(3), lons(IP))
            call check_err(iret)
          ELSE
            !WRITE (NDSO,952)
            iret=nf90_inq_varid(ncid(IP), 'y', varid(2))        
            call check_err(iret)
            iret=nf90_get_var(ncid(IP), varid(2), lats(IP))
            call check_err(iret)
            iret=nf90_inq_varid(ncid(IP), 'x', varid(3))        
            call check_err(iret)
            iret=nf90_get_var(ncid(IP), varid(3), lons(IP))
            call check_err(iret)
          END IF

          ! freq and dir variables
          iret=nf90_inq_varid(ncid(IP),"frequency",varid(4))
          call check_err(iret)
          iret=nf90_get_var(ncid(IP),varid(4),FREQ)
          call check_err(iret)
          iret=nf90_inq_varid(ncid(IP),"direction",varid(5))
          call check_err(iret)
          iret=nf90_get_var(ncid(IP),varid(5),THETA)
          call check_err(iret)
          THETA=(PI/180)*THETA

          iret = nf90_inq_varid(ncid(IP), 'station_name', varid(6)) 
          call check_err(iret)
          iret = nf90_get_var(ncid(IP), varid(6), station(IP,:),      &
                              start=(/1,1/),count=(/16,1/))
          call check_err(iret)


          iret=nf90_inq_varid(ncid(IP),"efth",varid(7))
          if (iret.ne.0 ) iret=nf90_inq_varid(ncid(IP),"Efth",varid(7))
          call check_err(iret)
          iret=nf90_get_var(ncid(IP),varid(7),SPEC2D(:,:,:,IP),       &
                            start=(/1,1,1,1/),count=(/NTHI,NKI,1,NTI/))
          call check_err(iret)
          iret=nf90_close(ncid(IP))
          call check_err(iret)
!
        END DO ! IP=1,NBO2



!
! 6. Defines frequency range in spectra
!

        ! Checks consistency of NK
        IF (NKI.GT.NK) GOTO 2008
!
        ! HERE we define IFMIN IFMIN2 IFMAX and IFMAX2 frequency indices
        ! such that source spec SPEC (read in input) links with output spec
        ! APBIN2 with APBIN2(IFMIN2:IFMAX2) = SPEC(IFMIN:IFMAX)
        ! Then APBIN2(1:IFMIN2) = 0 and APBIN2(IFMAX2:end) = 0
        IFMIN=1  ! index of first freq. in source spectrum
        IFMIN2=1 ! index of first freq. in output spectrum
        IFMAX=NK1 ! index of last freq. in source spectrum
!        IFMAX2=NK ! index of last freq. in output spectrum
!
        ! Checks consistency of XFR
        IF (ABS((FREQ(IFMIN+1)/FREQ(IFMIN))-XFR).GT.0.005) GOTO 2006
!
        ! Checks consistency of NTH
        ! WARNING: check is only done on number of directions, no check 
        ! is done on the relative offset of first direction in terms of 
        ! the directional increment [-0.5,0.5] (last parameter of the 
        ! spectral definition in ww3_grid.inp, on second active line) 
        IF (NTHI.NE.NTH) GOTO 2007
!
        IF ((FR1-FREQ(1))/FR1.GT. 0.03) THEN 
          DO J=1,MIN(NK1,NK)
            IF (ABS(FREQ(J)-FR1) .LT. ABS(FREQ(IFMIN)-FR1)) THEN 
              IFMIN=J
            END IF
          END DO  
        END IF
!     
        IF ((FREQ(1)-FR1)/FR1.GT. 0.03) THEN 
          DO J=1,MIN(NK,NK1)
            IF (ABS(FREQ(J)-FR1*XFR**(J-1)) .LT. ABS(FREQ(IFMIN2)-FR1)) THEN 
              IFMIN2=J
            END IF
          END DO  
        END IF
!
        IF ((FREQ(NK1)-FR1*XFR**(NK-1))/FREQ(NK1) .GT.0.03) THEN 
          DO J=1,NK
            IF (ABS(FREQ(J)-FR1*XFR**(NK1-1)) .LT. &
                ABS(FREQ(IFMAX)-FR1*XFR**(NK1-1))) THEN 
              IFMAX=J
            END IF
          END DO 
        END IF 
!
! 7. Loops on files and instanciate ABPIN2
!
        DO IP=1,NBO2
!         Fills in the low frequency end of the spectrum
          IF (IFMIN2.GT.1) THEN 
            ABPIN2(1:(IFMIN2-1)*NTH,:,IP)=0.
          END IF

!         Copies spectrum in frequency and direction ranges
          DO I=IFMIN,IFMAX
            DO J=1,NTH
              ABPIN2((I-IFMIN+(IFMIN2-1))*NTH+J,:,IP)=SPEC2D(J,I,:,IP)*tpiinv
            END DO
          END DO

!         Fills tail of spectrum            
          IF (IFMAX-IFMIN+IFMIN2.LT.NK1) THEN 
            IF (VERBOSE.GE.2) WRITE(NDSO,*) 'FILLING TAIL',IFMAX-IFMIN,NK1,IFMAX-IFMIN+(IFMIN2-1)
            ABPIN2((IFMAX-IFMIN+IFMIN2)*NTH+1:NK1*NTH,:,IP)=0.
          END IF 
        END DO ! IP=1,NBO2



!
! 8. Writes header
!
!       Writes header in nest.ww3 file 
        WRITE(NDSB) IDSTRBC, VERBPTBC, NK1, NTH1, XFR, FREQ(1),       &
        MOD(2.5*PI-THETA(1),TPI), NBO
        IPBPO(:,:)=1
        RDBPO(:,1)=1.
        RDBPO(:,2:4)=0.

!       Loops on points
        DO IP1=1,NBO
          DMIN=360.+180.
          DMIN2=360.+180.
!         Loops on files
          DO IP=1,NBO2
!           Searches for the nearest 2 points where spectra are available
            IF (FLAGLL)  THEN 
              DIST=DIST_SPHERE ( LONS(IP),LATS(IP),XBPO(IP1),YBPO(IP1) )
            ELSE 
              DIST=SQRT((LONS(IP)-XBPO(IP1))**2+(LATS(IP)-YBPO(IP1))**2)
              END IF
            IF (DMIN.EQ.(360.+180.)) THEN 
              IF(DIST.LT.DMIN) THEN 
                IPBPO(IP1,1)=IP
                DMIN=DIST
              END IF
            ELSE
              IF(DIST.LT.DMIN2) THEN 
                IF(DIST.LT.DMIN) THEN 
                  IPBPO(IP1,2)=IPBPO(IP1,1)
                  DMIN2=DMIN
                  IPBPO(IP1,1)=IP
                  DMIN=DIST
                ELSE
                  IPBPO(IP1,2)=IP
                  DMIN2=DIST
                END IF
              END IF
            END IF
          END DO ! IP1=1,NBO2
          IF (VERBOSE.GE.1) WRITE(NDSO,*) 'DIST:',DMIN,DMIN2,IP1,IPBPO(IP1,1),IPBPO(IP1,2), &
                                          LONS(IPBPO(IP1,1)),LONS(IPBPO(IP1,2)),XBPO(IP1), &
                                          LATS(IPBPO(IP1,1)),LATS(IPBPO(IP1,2)),YBPO(IP1)


!
! 9. Computes linear interpolation coefficient between the nearest 2 points
!
          IF (INTERP.GT.1.AND.NBO2.GT.1) THEN 
            IF (FLAGLL) THEN 
              DLON=LONS(IPBPO(IP1,2))-LONS(IPBPO(IP1,1))
              DLAT=LATS(IPBPO(IP1,2))-LATS(IPBPO(IP1,1))
              DLO=XBPO(IP1)-LONS(IPBPO(IP1,1))
              IF (DLON.GT.180.) DLON=DLON-360
              IF (DLON.LT.-180.) DLON=DLON+360
              IF (DLO.GT.180.) DLO=DLO-360
              IF (DLO.LT.-180.) DLO=DLO+360
              DIST=SQRT(DLON**2+DLAT**2) 
              COS1=( DLO*DLON &
                  + (YBPO(IP1)-LATS(IPBPO(IP1,1)))  &
                  *DLAT )/(DIST**2)
            ELSE 
              DIST=SQRT((LONS(IPBPO(IP1,1))-LONS(IPBPO(IP1,2)))**2   &
               +(LATS(IPBPO(IP1,1))-LATS(IPBPO(IP1,2)))**2)
              COS1=( (XBPO(IP1)-LONS(IPBPO(IP1,1)))  &
                  *(LONS(IPBPO(IP1,2))-LONS(IPBPO(IP1,1))) &
                  + (YBPO(IP1)-LATS(IPBPO(IP1,1)))  &
                  *(LATS(IPBPO(IP1,2))-LATS(IPBPO(IP1,1))) )/(DIST**2)
              END IF
            !COS2=( (XBPO(IP1)-LONS(IPBPO(IP1,2)))  &
            !      *(LONS(IPBPO(IP1,1))-LONS(IPBPO(IP1,2)))
            !      + (YBPO(IP1)-LATS(IPBPO(IP1,2)))  &
            !      *(LATS(IPBPO(IP1,1))-LATS(IPBPO(IP1,2))))/(DIST**2)
            RDBPO(IP1,1)=1-MIN(1.,MAX(0.,COS1))
            RDBPO(IP1,2)=MIN(1.,MAX(0.,COS1))
          ELSE 
            ! in this case: nearest point
            RDBPO(IP1,1)=1.
            RDBPO(IP1,2:4)=0.
            END IF
          IF (VERBOSE.GE.1) WRITE(NDSO,*) 'IPBP:',IP1,(IPBPO(IP1,J),J=1,4)
          IF (VERBOSE.GE.1) WRITE(NDSO,*) 'RDBP:',IP1,(RDBPO(IP1,J),J=1,4)
          !IF (VERBOSE.GE.1) WRITE(NDSO,*) 'RDBP:',COS1,DIST,DLON,DLO,DLAT,XBPO(IP1)-360.,LONS(IPBPO(IP1,1)),LONS(IPBPO(IP1,2))
        END DO ! IP1=1,NBO

        WRITE(NDSB)  (XBPO(I),I=1,NBO),            &
                     (YBPO(I),I=1,NBO),            &
                     ((IPBPO(I,J),I=1,NBO),J=1,4), &
                     ((RDBPO(I,J),I=1,NBO),J=1,4)


!
! 10. Loops on times and files and write to nest.ww3
!          
        DO IT=1,NT1
          datetmp=times(IT)
          IF (index(timeunits, "seconds").NE.0) datetmp=times(IT)/86400.
          datecur=datetmp+dateref
          call caldat(int(datecur),id,im,iyy)
 
          ih=(datecur-INT(datecur))*24.
          IF (ABS(ih-NINT(ih)).LT.1.E-5) ih=NINT(ih)

          imn=(ih-INT(ih))*60.
          IF (ABS(imn-NINT(imn)).LT.1.E-5) imn=NINT(imn)

          is=(imn-INT(imn))*60.
          IF (ABS(is-NINT(is)).LT.1.E-5) is=NINT(is)

          WRITE(TIME_STR,'(I4.4,5I2.2)') INT(iyy), INT(im), INT(id), &
                                        INT(ih), INT(imn), INT(is)

          read(TIME_STR(1:8),'(I8.8)') TIME(1)
          read(TIME_STR(9:14),'(I6.6)') TIME(2)

          WRITE(NDSO,'(A,2I9,A,I6,A)') 'Writing boundary data for time:', TIME, ' at ',NBO2,' points'
          WRITE(NDSB,IOSTAT=IERR) TIME, NBO2
          DO IP=1, NBO2
            WRITE(NDSB) ABPIN2(:,IT,IP)
          END DO
        END DO ! IT=0,NT1
        CLOSE(NDSB)
      END IF ! INXOUT.EQ.'WRITE'
      STOP
!
! Escape locations read errors :
!

 2001 CONTINUE
      WRITE (NDSE,1001)
      CALL EXTCDE ( 61 )
!
 2002 CONTINUE
      WRITE (NDSE,1002) IERR
      CALL EXTCDE ( 62 )
!
 2003 CONTINUE
      WRITE (NDSE,1003) IDTST, IDSTRBC
      CALL EXTCDE ( 63 )
!
 2004 CONTINUE
      WRITE (NDSE,1004)
      CALL EXTCDE ( 64 )
!
 2005 CONTINUE
      WRITE (NDSE,1005) IERR
      CALL EXTCDE ( 65 )
!
 2006 CONTINUE
      WRITE (NDSE,1006) XFR
      CALL EXTCDE ( 66 )
!
 2007 CONTINUE
      WRITE (NDSE,1007) NTH, NTHI
      CALL EXTCDE ( 67 )
!
 2008 CONTINUE
      WRITE (NDSE,1008) NK, NKI
      CALL EXTCDE ( 68 )
!
 2222 CONTINUE
!
  920 FORMAT ( '  Grid name : ',A/)
!
  940 FORMAT ( '  Format version : ',A/)
!
  941 FORMAT ( '  File type : ',A/)
!
  950 FORMAT ( '  Spectra File does not exist : ',A/)
!
  951 FORMAT ( '  Location expected as spherical (lon/lat) in Spectra &
                  file'/)
!
  952 FORMAT ( '  Location expected as cartesian (x/y) in Spectra     &
                  file'/)
!
 1900 FORMAT ( '  Writing boundary data for time : ',I8.8, I6.6/      &
               '  for this number of positions   : ',I4)
!
!
 1001 FORMAT (/' *** WAVEWATCH-III ERROR IN W3BOUNC : '/              &
               '     PREMATURE END OF INPUT FILE'/)
!
 1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC: '/               &
               '     ERROR IN READING ',A,' FROM INPUT FILE'/         &
               '     IOSTAT =',I5/)
!
 1003 FORMAT (/' *** WAVEWATCH-III ERROR IN W3IOBC :'/                &
               '     ILLEGAL IDSTR, READ : ',A/                       &
               '                   CHECK : ',A/)
!
 1004 FORMAT (/' *** WAVEWATCH-III ERROR IN W3BOUNC : '/              &
               '     PREMATURE END OF NEST FILE'/)
!
 1005 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC: '/               &
               '     ERROR IN READING SPECTRUM FROM NEST FILE'/       &
               '     IOSTAT =',I5/)
!
 1006 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC: '/               &
               '     ILLEGAL XFR, XFR =',F12.6/)
!
 1007 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC: '/               &
               '     ILLEGAL NTH, NTH =',I3,' DIFFERS FROM NTHI =',I3/)
!
 1008 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC: '/               &
               '     ILLEGAL NK, NK =',I3,' DIFFERS FROM NKI =',I3/   &
               '     IT WILL BE MANAGED SOON BY SPCONV')
!/
!/ End of W3BOUNC ---------------------------------------------------- /
!/
  END PROGRAM W3BOUNC
!/ ------------------------------------------------------------------- /


!==============================================================================

SUBROUTINE check_err(iret)

  USE NETCDF
  implicit none

  integer iret

  if (iret .ne. NF90_NOERR) then
     print *, nf90_strerror(iret)
     stop
  endif
  return

END SUBROUTINE check_err

!==============================================================================

