ww3_bound.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
PROGRAM w3bound
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           F. Ardhuin              |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         27-May-2021 |
  !/                  +-----------------------------------+
  !/
  !/    28-Aug-2012 : adaptation from SHOM/Ifremer code   ( version 4.08 )
  !/    01-Nov-2012 : Bug correction for NKI != NK        ( version 4.08 )
  !/    20-Oct-2016 : Error statement updates             ( version 5.15 )
  !/    21-Jul-2020 : Support rotated pole grid           ( version 7.11 )
  !/                                  Chris Bunney, UKMO.
  !/    27-May-2021 : Add namelist feature                ( version 7.XX )
  !/
  !/    Copyright 2012-2012 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_bound.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.
  !      W3EQTOLL  Subr  W3SERVMD Convert coordinates from rotated pole.
  !     ----------------------------------------------------------------
  !
  !  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
  !
  !     - Input spectra are assumed to be formulated on a standard
  !       pole. However, the model grid can be on a rotated pole.
  !
  !  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, gtype, xgrd, ygrd, x0, y0, &
       sx, sy, mapsf, ungtype, clgtype, rlgtype
#ifdef W3_RTD
  USE w3gdatmd, ONLY : polat, polon
#endif
  USE w3odatmd, ONLY: ndso, ndse, fnmpre
  USE w3iobcmd, ONLY: verbptbc, idstrbc
  USE w3iogrmd, ONLY: w3iogr
  USE w3timemd
  USE w3nmlboundmd
  USE w3servmd, ONLY: itrace, nextln, extcde
#ifdef W3_RTD
  USE w3servmd, ONLY: w3eqtoll
#endif
#ifdef W3_S
  USE w3servmd, ONLY : strace
#endif
  !/
  IMPLICIT NONE
  !
#ifdef W3_MPI
  include "mpif.h"
#endif
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
  TYPE(nml_bound_t)       :: nml_bound
  !
  INTEGER                 :: time1(2)
  INTEGER                 :: time2(2)
  INTEGER                 :: ix
  INTEGER                 :: iy
  INTEGER                 :: isea
  INTEGER                 :: i
  INTEGER                 :: jj
  INTEGER                 :: ip
  INTEGER                 :: ip1
  INTEGER                 :: j
  INTEGER                 :: k
  INTEGER                 :: itime
  INTEGER                 :: ndsi
  INTEGER                 :: ndsm
  INTEGER                 :: ndsi2
  INTEGER                 :: ndss
  INTEGER                 :: ndsb
  INTEGER                 :: ndstrc
  INTEGER                 :: ntrace
  INTEGER                 :: nk1
  INTEGER                 :: nth1
  INTEGER                 :: nspec1
  INTEGER                 :: nbi
  INTEGER                 :: nbi2
  INTEGER                 :: nki
  INTEGER                 :: nthi
  INTEGER                 :: nbo
  INTEGER                 :: nbo2
  INTEGER                 :: ierr
  INTEGER                 :: interp
  INTEGER                 :: iloop
  INTEGER                 :: ifmin
  INTEGER                 :: ifmin2
  INTEGER                 :: ifmax
  INTEGER                 :: verbose
  INTEGER                 :: ndsl
#ifdef W3_S
  INTEGER, SAVE           :: ient = 0
#endif
  INTEGER                 :: ibo
  INTEGER,  ALLOCATABLE   :: ipbpi(:,:)
  INTEGER,  ALLOCATABLE   :: ipbpo(:,:)
  !
  REAL                    :: fr1i
  REAL                    :: th1i
  REAL                    :: depth
  REAL                    :: u10
  REAL                    :: udir
  REAL                    :: curr
  REAL                    :: currdir
  REAL                    :: dmin
  REAL                    :: dist
  REAL                    :: dmin2
  REAL                    :: cos1
  !
  REAL, ALLOCATABLE       :: lats(:)
  REAL, ALLOCATABLE       :: lons(:)
  REAL, ALLOCATABLE       :: spec2d(:,:)
  REAL, ALLOCATABLE       :: freq(:)
  REAL, ALLOCATABLE       :: theta(:)
  REAL, ALLOCATABLE       :: xbpi(:)
  REAL, ALLOCATABLE       :: ybpi(:)
  REAL, ALLOCATABLE       :: rdbpi(:,:)
  REAL, ALLOCATABLE       :: xbpo(:)
  REAL, ALLOCATABLE       :: ybpo(:)
  REAL, ALLOCATABLE       :: rdbpo(:,:)
  REAL, ALLOCATABLE       :: abpin(:,:)
#ifdef W3_RTD
  REAL, ALLOCATABLE       :: xtmp(:)
  REAL, ALLOCATABLE       :: ytmp(:)
  REAL, ALLOCATABLE       :: angtmp(:)
  LOGICAL                 :: isrtd
#endif
  CHARACTER               :: comstr
  CHARACTER(LEN=80)       :: line
  CHARACTER(LEN=128)      :: bndfile
  CHARACTER(LEN=5)        :: inxout
  CHARACTER(LEN=10)       :: vertest
  CHARACTER(LEN=32)       :: idtst
  CHARACTER(LEN=120)      :: filename
  CHARACTER(LEN=120)      :: string1
  CHARACTER(LEN=120)      :: buoyname
  CHARACTER               :: space
  !
  CHARACTER(LEN=120), ALLOCATABLE :: specfiles(:)
  !
  LOGICAL                  :: flgnml
  !/
  !/ ------------------------------------------------------------------- /
 
 
  !/
  ! 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
  ndsm   = 20
  ndss   = 40
  ndsl   = 50
  !
  ndstrc =  6
  ntrace = 10
  CALL itrace ( ndstrc, ntrace )
  !
#ifdef W3_S
  CALL strace (ient, 'W3BOUND')
#endif
  !
 
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 2.  Read model definition file.
  !
  CALL w3iogr ( 'READ', ndsm )
  WRITE (ndso,920) gname
#ifdef W3_RTD
  !
  isrtd = polat .LT. 90.0
  !
#endif
  !
  ! 3. Read input file
  !
  ! process ww3_bound namelist
  !
  INQUIRE(file=trim(fnmpre)//"ww3_bound.nml", exist=flgnml)
  IF (flgnml) THEN
    ! Read namelist
    CALL w3nmlbound (ndsi, trim(fnmpre)//'ww3_bound.nml', nml_bound, ierr)
 
    inxout = nml_bound%MODE
    interp = nml_bound%INTERP
    verbose = nml_bound%VERBOSE
    bndfile = nml_bound%FILE
 
    nbo2 = 0
    OPEN(ndsl,file=trim(bndfile),status='OLD',err=809,iostat=ierr)
    rewind(ndsl)
    DO
      READ (ndsl,*,END=400,ERR=802)
      nbo2 = nbo2 + 1
    END DO
400 CONTINUE
    ALLOCATE(specfiles(nbo2))
    rewind(ndsl)
    DO i=1,nbo2
      READ (ndsl,'(A120)',END=801,ERR=802) SPECFILES(I)
    END DO
    CLOSE(ndsl)
 
  END IF ! FLGNML
 
  !
  ! process old ww3_bound.inp format
  !
  IF (.NOT. flgnml) THEN
    OPEN (ndsi,file=trim(fnmpre)//'ww3_bound.inp',status='OLD',err=803,iostat=ierr)
    rewind(ndsi)
 
    READ (ndsi,'(A)',END=801,ERR=802) 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_bound.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,'(A120)') 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
 
      IF ( iloop .EQ. 1 ) CLOSE ( ndss)
 
      IF ( iloop .EQ. 2 ) CLOSE ( ndss, status='DELETE' )
    END DO ! ILOOP = 1, 2
    CLOSE(ndsi)
  END IF ! .NOT. FLGNML
 
 
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 4. Tests the reading of the file
  !
  IF ( inxout.EQ.'READ') THEN
    OPEN(ndsb,file='nest.ww3',form='UNFORMATTED', convert=file_endian,status='old')
    READ(ndsb) idtst, vertest, nk1, nth1, xfr, fr1i, th1i, nbi
    nspec1  = nk1 * nth1
    IF ( idtst .NE. idstrbc ) THEN
      WRITE (ndso,901) idtst, idstrbc
    END IF
    WRITE(ndso,*) "FORMAT VERSION: '",vertest,"'"
    WRITE(ndso,*) "FILE TYPE: '",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.EQ.1) WRITE(ndso,*)      'XBPI:',xbpi
    IF (verbose.EQ.1) WRITE(ndso,*)      'YBPI:',ybpi
    IF (verbose.EQ.1) WRITE(ndso,*)      'IPBPI:'
    DO i=1,nbi
      IF (verbose.EQ.1) WRITE(ndso,*) i,' interpolated from:',ipbpi(i,1:4)
      IF (verbose.EQ.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,iostat=ierr) abpin(:,ip)
        END DO
      END IF
    END DO
    CLOSE(ndsb)
  END IF
  !
  !
  IF ( inxout.EQ.'WRITE') THEN
    IF ( flbpi) THEN
      !
      !  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))
#ifdef W3_RTD
      IF (isrtd) ALLOCATE(xtmp(nbo), ytmp(nbo), angtmp(nbo))
#endif
      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)= xgrd(1,ix)
            ybpo(ibo)= ygrd(1,ix)
          END SELECT !GTYPE
        END IF
      END DO
#ifdef W3_RTD
      ! Convert grid boundary cell locations to standard pole
      IF( isrtd ) THEN
        xtmp = xbpo
        ytmp = ybpo
        CALL w3eqtoll(ytmp, xtmp, ybpo, xbpo, angtmp, polat, polon, nbo)
        DEALLOCATE(xtmp, ytmp, angtmp)
      ENDIF
#endif
      OPEN(ndsb,file='nest.ww3',form='UNFORMATTED', convert=file_endian,status='unknown')
      ALLOCATE(lats(nbo2),lons(nbo2))
      DO ip=1,nbo2
        OPEN(200+ip,file=specfiles(ip),status='old',iostat=ierr)
        IF (verbose.EQ.1) WRITE(ndso,'(A,I5,3A,I5)') &
             'IP, file, I/O stat:',ip,', ', &
             trim(specfiles(ip)), ', ',ierr
        IF (ierr.NE.0) GOTO 810
        READ(200+ip,'(A1,A22,A1,X,2I6)',iostat=ierr)  &
             space,string1,space,nki,nthi
        IF (verbose.EQ.1) WRITE(ndso,'(A,3I5)') 'IP and spectral dimensions:',ip, nki,nthi
        IF (ip.EQ.1) THEN
          nk1=nki
          nth1=nthi
          nspec1  = nk1 * nth1
          ALLOCATE (freq(nk1),theta(nth1))
          ALLOCATE (spec2d(nk1,nth1))
          ALLOCATE (abpin(nk*nth1,nbo2))
        ELSE
          !
          !  To be cleaned up later ...
          !
          IF (nk1.NE.nki.OR.nth1.NE.nthi) THEN
            WRITE(ndse,'(A,A,4I5)') 'ERROR, SPECTRAL GRID IN FILE:',     &
                 trim(specfiles(ip)),nk1,nki,nth1,nthi
            stop
          END IF
        END IF
        !
        READ(200+ip,*) freq(1:nk1)
        READ(200+ip,*) theta(1:nth1)
      END DO
 
      !
      ! Defines frequency range in spectra
      !
 
      ! Checks consistency of NK
      IF (nki.GT.nk) GOTO 808
      !
      ! HERE we define IFMIN IFMIN2 IFMAX and IFMAX2 frequency indices
      ! such that source spec SPEC (read in input) links with output spec
      ! APBIN with APBIN(IFMIN2:IFMAX2) = SPEC(IFMIN:IFMAX)
      ! Then APBIN(1:IFMIN2) = 0 and APBIN(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 806
      !
      ! 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 807
 
      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
      !
      ierr=0
      itime=0
      !
      ! Loop on times
      !
      DO WHILE (ierr.EQ.0)
        DO ip=1,nbo2
          READ(200+ip,*,iostat=ierr) time2
          IF (ierr.EQ.0) THEN
            !
            IF (ip.EQ.1) THEN
              time1=time2
            ELSE
              IF (time1(1).NE.time2(1).OR.time1(2).NE.time2(2)) THEN
                WRITE(ndse,*) 'AT POINT ',ip,', BAD TIMES:',time1, time2
              END IF
            END IF
            !
            READ(200+ip,'(A1,A10,A1,2F7.2,F10.1,F7.2,F6.1,F7.2,F6.1)') &
                 space,buoyname,space,lats(ip),lons(ip),depth,u10,udir,curr,currdir
#ifdef W3_RTD
            IF (isrtd) THEN
              ! Rotated coordinates are scaled in range 0 - 360
              IF(lons(ip) .LT. 0) lons(ip) = lons(ip) + 360.0
              IF(lons(ip) .GT. 360) lons(ip) = lons(ip) - 360.0
            ENDIF
#endif
            READ(200+ip,*,iostat=ierr) spec2d
            IF (ifmin2.GT.1) THEN
              !
              !  Fills in the low frequency end of the spectrum
              !
              abpin(1:(ifmin2-1)*nth,ip)=0.
            END IF
            DO i=ifmin,ifmax
              DO j=1,nth
                abpin((i-ifmin+(ifmin2-1))*nth+j,ip)=spec2d(i,j)*tpiinv
              END DO
            END DO
            IF (ifmax-ifmin+ifmin2.LT.nk1) THEN
              !IF (VERBOSE.EQ.1) WRITE(NDSO,*) 'FILLING TAIL',IFMAX-IFMIN,NK1,IFMAX-IFMIN+(IFMIN2-1)
              abpin((ifmax-ifmin+ifmin2)*nth+1:nk1*nth,ip)=0.
            END IF
          END IF  ! ned of test on IERR
        END DO
        !
        ! Writes header
        !
        IF (ierr.EQ.0) THEN
          IF (itime.EQ.0) THEN
            ! Correction for rounding error in ASCII files ...
            IF (abs(theta(1)-0.5*pi).LT.0.01) theta(1)=0.5*pi
            ! Writes header in nest.ww3 file
            WRITE(ndsb) idstrbc, verbptbc, nk1, nth, xfr, freq(1), &
                 mod(2.5*pi-theta(1),tpi), nbo
            ipbpo(:,:)=1
            rdbpo(:,1)=1.
            rdbpo(:,2:4)=0.
            !
            DO ip1=1,nbo
              dmin=360.+180.
              dmin2=360.+180.
              DO ip=1,nbo2
                !
                ! Searches for the nearest 2 points where spectra are available
                !
                dist=sqrt((lons(ip)-xbpo(ip1))**2+(lats(ip)-ybpo(ip1))**2)
                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
                  ENDIF
                END IF
              END DO
              !IF (VERBOSE.EQ.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)
              !
              !  Computes linear interpolation coefficient between the nearest 2 points
              !
              IF (interp.GT.1.AND.nbo2.GT.1) THEN
                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)
                !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))
              END IF
              !
              IF (verbose.EQ.1) WRITE(ndso,*) 'IPBP:',ip1,(ipbpo(ip1,j),j=1,4)
              IF (verbose.EQ.1) WRITE(ndso,*) 'RDBP:',ip1,(rdbpo(ip1,j),j=1,4)
              !
            END DO
            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)
          END IF
 
          WRITE(ndso,*) 'Writing boundary data for time:', time2, nbo2
          WRITE(ndsb,iostat=ierr) time2, nbo2
          DO ip=1, nbo2
            WRITE (ndsb) abpin(:,ip)
          END DO
 
          itime=itime+1
        END IF
      END DO
      CLOSE(ndsb)
    END IF
  END IF
  stop
  !
  ! Escape locations read errors :
  !
 
801 CONTINUE
  WRITE (ndse,1001)
  CALL extcde ( 61 )
  !
802 CONTINUE
  WRITE (ndse,1002) ierr
  CALL extcde ( 62 )
  !
803 CONTINUE
  WRITE (ndse,1003)
  CALL extcde ( 63 )
  !
806 CONTINUE
  WRITE (ndse,1006) xfr
  CALL extcde ( 66 )
  !
807 CONTINUE
  WRITE (ndse,1007) nth, nthi
  CALL extcde ( 67 )
  !
808 CONTINUE
  WRITE (ndse,1008) nk, nki
  CALL extcde ( 68 )
  !
809 CONTINUE
  WRITE (ndse,1009) bndfile, ierr
  CALL extcde ( 69 )
  !
810 CONTINUE
  WRITE (ndse,1010) specfiles(ip)
  CALL extcde ( 70 )
 
  !
  ! Formats
  !
901 FORMAT (/' *** WAVEWATCH-III ERROR IN W3IOBC :'/                &
       '     ILEGAL IDSTR, READ : ',a/                        &
       '                  CHECK : ',a/)
  !
920 FORMAT ( '  Grid name : ',a/)
  !
1001 FORMAT (/' *** WAVEWATCH-III ERROR IN W3BOUND : '/              &
       '     PREMATURE END OF INPUT FILE'/)
  !
1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND: '/               &
       '     ERROR IN READING ',a,' FROM INPUT FILE'/         &
       '     IOSTAT =',i5/)
  !
1003 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC : '/              &
       '     ERROR IN OPENING INPUT FILE: ', a/               &
       '     IOSTAT =',i5/)
  !
1006 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND: '/               &
       '     ILLEGAL XFR, XFR =',f12.6/)
  !
1007 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND: '/               &
       '     ILLEGAL NTH, NTH =',i3,' DIFFERS FROM NTHI =',i3/)
  !
1008 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND: '/               &
       '     ILLEGAL NK, NK =',i3,' DIFFERS FROM NKI =',i3/   &
       '     IT WILL BE MANAGED SOON BY SPCONV')
  !
1009 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND : '/              &
       '     ERROR IN OPENING SPEC FILE: ', a/                &
       '     IOSTAT =',i5/)
  !
1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUND : '/              &
       '     SPEC FILE NOT EXISTING: ', a/)
 
 
  !
  !/
  !/ End of W3BOUND ---------------------------------------------------- /
  !/
END PROGRAM w3bound
!/ ------------------------------------------------------------------- /