ww3_bounc.F90

Go to the documentation of this file.

 
!
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!
PROGRAM w3bounc
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           F. Ardhuin              |
  !/                  |           M. Accensi              |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         21-Jul-2020 |
  !/                  +-----------------------------------+
  !/
  !/    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 )
  !/    20-Mar-2018 : Improve netcdf file reading         ( version 6.02 )
  !/    15-May-2018 : Add namelist feature                ( version 6.05 )
  !/    04-May-2020 : Update spectral conversion          ( version 7.11 )
  !/    21-Jul-2020 : Support rotated pole grid           ( version 7.11 )
  !/
  !/
  !/    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.
  !      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, fnmpre, ndst, ndse
  USE w3cspcmd, ONLY: w3cspc
 
 
  USE w3gdatmd, ONLY: nk, nth, xfr, fr1, dth, th, fachfe,           &
       gname, w3nmod, w3setg,&
       nsea, mapsta, gtype, xgrd, ygrd, x0, y0, &
       sx, sy, mapsf, ungtype, clgtype, rlgtype, flagll
#ifdef W3_RTD
  USE w3gdatmd, ONLY : polat, polon
#endif
  USE w3odatmd, ONLY: ndso, ndse
  USE w3iobcmd, ONLY: verbptbc, idstrbc
  USE w3iogrmd, ONLY: w3iogr
  USE w3timemd
  USE w3servmd, ONLY: itrace, nextln, extcde, dist_sphere
#ifdef W3_RTD
  USE w3servmd, ONLY: w3eqtoll
#endif
  USE w3nmlbouncmd
  USE netcdf
#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                 :: 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, verbose, ibo,    &
       iret, icode, ndsl
  INTEGER                 :: time(2), time2(2), varid(12),         &
       refdate(8), curdate(8), vartype
#ifdef W3_S
  INTEGER, SAVE           :: ient = 0
#endif
  !
  INTEGER, ALLOCATABLE    :: ipbpi(:,:), ipbpo(:,:), ncid(:),      &
       dimid(:,:), dimln(:,:)
  !
  REAL                    :: fr1i, xfri, th1i, factor, offset, dmin,&
       dist, dmin2, cos1, dlon, dlat, dlo,    &
       fillval
  !
  REAL, ALLOCATABLE       :: spec2d(:,:,:,:), lats(:), lons(:),    &
       freq(:), theta(:),                    &
       xbpi(:), ybpi(:), rdbpi(:,:),         &
       xbpo(:), ybpo(:), rdbpo(:,:),         &
       abpin(:,:), abpin2(:,:,:)
#ifdef W3_RTD
  REAL, ALLOCATABLE     :: xtmp(:), ytmp(:), angtmp(:)
  LOGICAL               :: isrtd
#endif
  !
  REAL, ALLOCATABLE       :: tmpspci(:,:),tmpspco(:,:)
 
  !
  DOUBLE PRECISION        :: refjulday, curjulday
  DOUBLE PRECISION, ALLOCATABLE       :: times(:)
  !
  CHARACTER               :: comstr*1, line*512, filename*512,     &
       inxout*5, file*128
  CHARACTER*50            :: timeunits, calendar
  CHARACTER*10            :: vertest  ! = '2018-03-01'
  CHARACTER*32            :: idtst    != 'WAVEWATCH III BOUNDARY DATA FILE'
  CHARACTER*512, ALLOCATABLE          :: specfiles(:)
  CHARACTER, ALLOCATABLE              :: station(:,:)
  !
  LOGICAL                 :: flgnml, spconv
  !
  !/
  !/ ------------------------------------------------------------------- /
 
 
  !/
  ! 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
  ndsl   = 50
  ndso   = 6
  ndse   = 6
  !
  ndstrc =  6
  ntrace = 10
  CALL itrace ( ndstrc, ntrace )
  !
#ifdef W3_S
  CALL strace (ient, 'W3BOUNC')
#endif
  !
  WRITE (ndso,900)
  !
 
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 2.  Read model definition file.
  !
  CALL w3iogr ( 'READ', ndsm )
  WRITE (ndso,920) gname
#ifdef W3_RTD
  !
  isrtd = polat .LT. 90.0
  !
#endif
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 3.  Read requests from input file.
  !
 
  !
  ! process ww3_bounc namelist
  !
  INQUIRE(file=trim(fnmpre)//"ww3_bounc.nml", exist=flgnml)
  IF (flgnml) THEN
    ! Read namelist
    CALL w3nmlbounc (ndsi, trim(fnmpre)//'ww3_bounc.nml', nml_bound, ierr)
 
    inxout = nml_bound%MODE
    interp = nml_bound%INTERP
    verbose = nml_bound%VERBOSE
    file = nml_bound%FILE
 
    nbo2 = 0
    OPEN(ndsl,file=trim(file),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,'(A512)',END=801,ERR=802) SPECFILES(I)
    END DO
    CLOSE(ndsl)
 
  END IF ! FLGNML
 
  !
  ! process old ww3_bounc.inp format
  !
  IF (.NOT. flgnml) THEN
    OPEN (ndsi,file=trim(fnmpre)//'ww3_bounc.inp',status='OLD',err=805,iostat=ierr)
    rewind(ndsi)
 
    READ (ndsi,'(A)',END=801,ERR=802,IOSTAT=IERR) comstr
    IF (comstr.EQ.' ') comstr = '$'
    WRITE (ndso,901) comstr
 
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) inxout
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) interp
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) 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
      !
      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 ) GOTO 803
    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=803,ERR=804) 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))
#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
    !
    IF (isrtd) THEN
      ! Convert grid boundary cell locations to standard pole
      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(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', convert=file_endian,      &
           status='old',iostat=icode)
      IF (icode.NE.0) THEN
        lons(ip)=-999.
        lats(ip)=-999.
        WRITE (ndse,1010) trim(specfiles(ip))
        CALL extcde ( 70 )
      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*nth,nt1,nbo2))
 
        ! instanciates time
        refdate(:)=0.
        iret=nf90_inq_varid(ncid(ip),"time",varid(1))
        CALL check_err(iret)
        iret=nf90_get_var(ncid(ip), varid(1), times(:))
        CALL check_err(iret)
        iret=nf90_get_att(ncid(ip),varid(1),"calendar",calendar)
        IF ( iret/=nf90_noerr ) THEN
          WRITE(ndse,951)
        ELSE IF ((index(calendar, "standard").EQ.0) .AND. &
             (index(calendar, "gregorian").EQ.0)) THEN
          WRITE(ndse,952)
        END IF
        iret=nf90_get_att(ncid(ip),varid(1),"units",timeunits)
        CALL u2d(timeunits,refdate,ierr)
        CALL d2j(refdate,refjulday,ierr)
 
      ELSE
        IF (nki.NE.nk1.OR.nthi.NE.nth1.OR.nt1.NE.nti &
             ) GOTO 805
      END IF
 
      ! position variables : lon/lat or x/y
      IF ( flagll ) THEN
        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
        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=mod(2.5*pi-(pi/180)*theta,tpi)
 
      ! 2D spectra depending on station name or lat/lon
      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_inquire_variable(ncid(ip),varid(7),xtype=vartype)
      CALL check_err(iret)
      iret=nf90_get_att(ncid(ip),varid(7),"_FillValue",fillval)
      CALL check_err(iret)
      iret=nf90_get_att(ncid(ip),varid(7),"scale_factor",factor)
      IF (iret.NE.0) factor=1.
      iret=nf90_get_att(ncid(ip),varid(7),"add_offset",offset)
      IF (iret.NE.0) offset=0.
      iret = nf90_inq_varid(ncid(ip), 'station_name', varid(6))
      IF (iret.NE.0) THEN
        ! efth(time, frequency, direction, latitude, longitude)
        iret=nf90_get_var(ncid(ip),varid(7),spec2d(:,:,:,ip),       &
             start=(/1,1,1,1/),count=(/1,1,nthi,nki,nti/))
        CALL check_err(iret)
      ELSE
        ! efth(time, station, frequency, direction)
        iret=nf90_get_var(ncid(ip),varid(7),spec2d(:,:,:,ip),       &
             start=(/1,1,1,1/),count=(/nthi,nki,1,nti/))
        CALL check_err(iret)
      END IF
      ! apply scale_factor and add_offset
      IF (vartype.EQ.nf90_short) THEN
        WHERE(spec2d(:,:,:,ip).NE.fillval) spec2d(:,:,:,ip)=(exp(spec2d(:,:,:,ip)*factor*log(10.)))-1e-12
      ELSE
        WHERE(spec2d(:,:,:,ip).NE.fillval) spec2d(:,:,:,ip)=(spec2d(:,:,:,ip)*factor)+offset
      END IF
 
      ! close spectra file
      iret=nf90_close(ncid(ip))
      CALL check_err(iret)
      !
    END DO ! IP=1,NBO2
 
 
 
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 6.  Checks on spectral discretization
    !     reminder: fr(NK)=fr1*XFR**(NK-1)
    !
    fr1i=freq(1)
    xfri=exp(alog(freq(nki)/freq(1))/(nki-1))
    th1i=theta(1)
 
    spconv = nki.NE.nk .OR. nthi.NE.nth .OR.                    &
         abs(xfri/xfr-1.).GT.0.01 .OR.                      &
         abs(fr1i/fr1-1.).GT.0.01 .OR.                      &
         abs(th1i-th(1)).GT.0.01*dth
 
    IF (verbose.GE.1) WRITE(ndso,*) 'SPCONV:', spconv, nki, nk, nthi, nth
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 7. Loops on files and instanciate ABPIN2
    !
    IF ( .NOT. spconv ) THEN
 
      DO ip=1,nbo2
        ! Copies spectrum in frequency and direction ranges
        DO i=1,nk
          DO j=1,nth
            abpin2((i-1)*nth+j,:,ip)=spec2d(j,i,:,ip)*tpiinv
          END DO
        END DO
      END DO ! IP=1,NBO2
      !
    ELSE
      ALLOCATE(tmpspci(nki*nthi,nti))
      ALLOCATE(tmpspco(nk*nth,  nti))
      DO ip=1,nbo2
        DO i=1,nki
          DO j=1,nthi
            tmpspci((i-1)*nthi+j,:)=spec2d(j,i,:,ip)*tpiinv
          END DO
        END DO
        CALL w3cspc ( tmpspci, nki, nthi, xfri, fr1i, th1i, &
             tmpspco, nk,  nth,  xfr,  fr1,  th(1),&
             nti, ndst, ndse, fachfe )
        abpin2(:,:,ip)=tmpspco(:,:)
      END DO
      !
    END IF
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 8. Writes header
    !
    !       Writes header in nest.ww3 file
    WRITE(ndsb) idstrbc, verbptbc, nk, nth, xfr, fr1,       &
         th(1), 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
      curjulday=times(it)
      IF (index(timeunits, "seconds").NE.0)   curjulday=curjulday/86400.
      IF (index(timeunits, "minutes").NE.0)   curjulday=curjulday/1440.
      IF (index(timeunits, "hours").NE.0)     curjulday=curjulday/24.
      curjulday=refjulday+curjulday
 
      ! convert julday to date and time
      CALL j2d(curjulday,curdate,ierr)
      CALL d2t(curdate,time,ierr)
 
      ! write to output file nest.ww3
      WRITE(ndso,'(A,2I9,A,I6,A,G16.5)') 'Writing boundary data for time:', &
           time, ' at ',nbo2,' points. Max.: ', maxval(abpin2(:,it,:))
      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'
 
  GOTO 888
 
  !
  ! 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) idtst, idstrbc
  CALL extcde ( 63 )
  !
804 CONTINUE
  WRITE (ndse,1004)
  CALL extcde ( 64 )
  !
805 CONTINUE
  WRITE (ndse,1005) trim(specfiles(ip)), nki, nk1, nthi, nth1, nti, nt1
  CALL extcde ( 65 )
  !
809 CONTINUE
  WRITE (ndse,1009) file, ierr
  CALL extcde ( 69 )
  !
888 CONTINUE
  WRITE (ndso,999)
 
 
  !
  ! Formats
  !
900 FORMAT (/15x,'   *** WAVEWATCH III Bounday input prep. ***   '/ &
       15x,'==============================================='/)
  !
901 FORMAT ( '  Comment character is ''',a,''''/)
  !
920 FORMAT ( '  Grid name : ',a/)
  !
940 FORMAT ( '  Format version : ',a/)
  !
941 FORMAT ( '  File type : ',a/)
  !
951 FORMAT (/' *** WAVEWATCH III WARNING IN W3BOUNC : '/            &
       '     CALENDAR ATTRIBUTE NOT DEFINED'/                 &
       '     IT MUST RESPECT STANDARD OR GREGORIAN CALENDAR')
  !
952 FORMAT (/' *** WAVEWATCH III WARNING IN W3BOUNC : '/            &
       '     CALENDAR ATTRIBUTE NOT MATCH'/                   &
       '     IT MUST RESPECT STANDARD OR GREGORIAN CALENDAR')
  !
999 FORMAT (/'  End of program '/                                   &
       ' ========================================='/          &
       '         WAVEWATCH III Boundary input '/)
  !
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: '/               &
       '     INCONSISTENT SPECTRAL DIMENSION FOR FILE ',a/    &
       '     NKI =',i3,' DIFFERS FROM NK1 =',i3/      &
       '     OR NTHI =',i3,' DIFFERS FROM NTH1 =',i3/ &
       '     OR NTI =',i5,' DIFFERS FROM NT1 =',i5 /)
  !
1009 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC : '/              &
       '     ERROR IN OPENING SPEC FILE: ', a/                &
       '     IOSTAT =',i5/)
  !
1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3BOUNC : '/              &
       '     SPEC FILE DOES NOT EXIST : ',a/)
  !
  !
  !/
  !/ End of W3BOUNC ---------------------------------------------------- /
  !/
END PROGRAM w3bounc
!/ ------------------------------------------------------------------- /
 
 
!==============================================================================
SUBROUTINE check_err(IRET)
 
  USE netcdf
  USE w3odatmd, ONLY: ndse
  USE w3servmd, ONLY: extcde
 
  IMPLICIT NONE
 
  INTEGER IRET
 
  IF (iret .NE. nf90_noerr) THEN
    WRITE(ndse,*) ' *** WAVEWATCH III ERROR IN BOUNC :'
    WRITE(ndse,*) ' NETCDF ERROR MESSAGE: '
    WRITE(ndse,*) nf90_strerror(iret)
    CALL extcde ( 59 )
  END IF
  RETURN
 
END SUBROUTINE check_err
 
!==============================================================================