module swn_outnc
!
! --|-----------------------------------------------------------|--
! | BMT ARGOSS |
! | Voorsterweg 28, 8316 PT Vollenhove |
! | http://www.bmtargoss.com |
! | |
! | |
! | Programmer: A.Th.C. Hulst |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! Authors
!
! 0.01: A.Th.C.Hulst
!
! Updates
!
! 0.01, Feb 2012: New Module
!
! Purpose
!
! Write SWAN output into netcdf files
!
! Method
!
! MODULE construct
!
! Modules used
!
! NETCDF : Standard netcdf library
! AGIONCMD : Set of utilities for writing (metocean) data to netcdf
! NCTABLEMD : description table for metocean parameters
use agioncmd
use M_PARALL
use NETCDF
use nctablemd, only: nctable, nctable_record, get_nctable_record
use OUTP_DATA, only: ORQDAT, MAX_OUTP_REQ, LCOMPGRD, OUTP_FILES, NREOQ
use OCPCOMM2, only: PROJNR, PROJID, VERTXT
use OCPCOMM4
use SWCOMM1, only: CHTIME, OVEXCV
use SWCOMM2, only: XOFFS, YOFFS, OPTG, EXCFLD
use SWCOMM3, only: NSTATM, ALPC, MDC, MSC, MCGRD, MXC, MYC, DNORTH, &
PI2, ICUR, BNAUT
use SWCOMM4
use TIMECOMM, only: TFINC, TIMCO, TINIC, DT
!
implicit none
!
! Module parameters
!
! oqi(1) : saved alias for nref
! oqi(2) : index of request
! oqi(3) : number of output variables
! oqi(4) : idla of output request
! oqr(1) : tbegin
! oqr(2) : delta
!
! Module variables
character*40 :: STNAMES(71,2) = ''
logical :: stnames_initialized = .false., &
skip_range_error = .true.
type(recordaxe_type), save :: recordaxe(MAX_OUTP_REQ)
integer,save :: ncoffset(MAX_OUTP_REQ) = 0
!PUN logical :: PUNSWAN = .true.
logical :: PUNSWAN = .false.
type spcaux_type
real, dimension(:), allocatable :: hs, depth, ux, uy, wndx, wndy, &
xc, yc, xp, yp, f, theta
integer, dimension(:), allocatable :: ips
real, dimension(:,:), allocatable :: E
logical :: is2d = .false., relative = .false.
integer :: mip = 0, ndir = 0, nfreq = 0
end type spcaux_type
!
! Source text
!
contains
subroutine swn_outnc_spec(RTYPE, OQI, OQR, MIP, VOQR, VOQ, AC2, &
SPCSIG, SPCDIR, DEP2, KGRPNT, CROSS, IONOD)
!
!
! --|-----------------------------------------------------------|--
! | BMT ARGOSS |
! | Voorsterweg 28, 8316 PT Vollenhove |
! | http://www.bmtargoss.com |
! | |
! | |
! | Programmer: A.Th.C. Hulst |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
!
! 1. Updates
!
!
! 2. Purpose
!
! Write action density spectrum to netcdf file when serial OR
! binary file when PARLL or PUNSWAN
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! IO..name....type...size.......description
! i RTYPE char * type of output request: 'SPEC' for 2-D spectral
! 'SPE1' for 1-D freq. spectrum
! i IONOD int MIP array indicating in which subdomain output points
! are located
! i SPCDIR real MDC (*,1); spectral directions (radians)
! (*,2); cosine of spectral directions
! (*,3); sine of spectral directions
! (*,4); cosine^2 of spectral directions
! (*,5); cosine*sine of spectral directions
! (*,6); sine^2 of spectral directions
! i SPCSIG real MSC Relative frequencies in computational domain in sigma-space
! i VOQR int * Adminstrative list tracks which column in VOQ has which
! variable
! i CROSS bool 4,MIP
! i KGRPNT int MXC,MYC
! i MIP int 1 number of output points
! i VOQ real * Array with "integrated" parameters
! i AC2 real MDC,MSC, 2D spectra
! MCGRD
! i DEP2 real MCGRD
!
real, intent( in) :: SPCDIR(MDC,6), SPCSIG(MSC)
logical, intent( in) :: CROSS(1:4,1:MIP)
character (len=*), intent( in) :: RTYPE
integer, intent( in) :: VOQR(*), &
KGRPNT(MXC,MYC), IONOD(*), &
MIP
integer, intent(inout) :: OQI(4)
real*8, intent( in) :: OQR(2)
real, intent( in) :: VOQ(MIP,*), AC2(MDC,MSC,MCGRD), &
DEP2(MCGRD)
!
! 5. Local variables
!
integer :: irq
integer, save :: IENT=0
type(spcaux_type) :: spcaux
! 8. Subroutines used
!
! swn_outnc_spcaux
! swn_outnc_openspecfile
! swn_outnc_appendspc
! nccheck
! nf90_close
! swn_outnc_deallocate_spcaux
!
logical STPNOW
!
! 9. Subroutines calling
!
! SWOUTP (SWAN/OUTP)
!
! 10. Error messages
!
! ---
!
! 11. Remarks
!
!
! 12. Structure.
!
! a. compute spectra and auxilary data (writes binary file if parll or punswan)
! if serial
! b. open / create netcdf file
! c. write spectra to netcdf
! d. close ncfile
!
! 13. Source text
!
if (LTRACE) call STRACE (IENT,'swn_outnc_spec')
irq = OQI(2)
call swn_outnc_spcaux(RTYPE, OQI, MIP, VOQR, VOQ, AC2, &
SPCSIG, SPCDIR, DEP2, KGRPNT, CROSS, IONOD, &
spcaux)
if (STPNOW()) return
!
! b. open / create netcdf file (generic output file)
!
if ( .not. (PARLL .or. PUNSWAN) ) then
! opening / creating the netcdf file. Sets module variable recordaxe(irq)
call swn_outnc_openspecfile(OUTP_FILES(irq), spcaux, OQI, OQR)
if (STPNOW()) return
call swn_outnc_appendspc(OQI, spcaux)
if (STPNOW()) return
call nccheck(nf90_close(oqi(1) + ncoffset(irq) ))
oqi(1) = 0
end if
call swn_outnc_deallocate_spcaux( spcaux )
end subroutine swn_outnc_spec
subroutine swn_outnc_colspc ( RTYPE, OQI, OQR, MIP, KGRPGL )
use SwanGriddata, only: xcugrdgl, ycugrdgl, nverts
!PUN USE OCPCOMM2, ONLY: LENFNM, DIRCH2
!PUN USE SIZES, ONLY: GLOBALDIR, LOCALDIR
!
!
! --|-----------------------------------------------------------|--
! | BMT ARGOSS |
! | Voorsterweg 28, 8316 PT Vollenhove |
! | http://www.bmtargoss.com |
! | |
! | |
! | Programmer: A.Th.C. Hulst |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
!
! 1. Updates
!
!
! 2. Purpose
!
! Collect spectra from binary files and store in netcdf file
!
! 3. Method
!
! ---
!
! 4. Argument variables
character (len=*), intent( in) :: RTYPE * 4
integer, intent(inout) :: OQI(4)
real*8, intent( in) :: OQR(2)
integer, intent( in) :: MIP, KGRPGL(MXCGL,MYCGL)
!
! IO..name....type...size.......description
! i RTYPE char * type of output request: 'SPEC' for 2-D spectral
! 'SPE1' for 1-D freq. spectrum
! i XP, YP real MIP coordinates of all output points
! i MIP int 1 number of output points
!
! 5. Local variables
!
logical :: lopen, read_header, file_exists, &
EQREAL, STPNOW
integer :: ip, ierr, otype, xpctime2, &
i, tmip, binnr, irq, &
tip, ips(MIP), iproc, nref, ilpos
integer, save :: IENT=0
character*80 :: binfile, ncfile, basefile
character*256 :: errmsg
type(spcaux_type) :: spcaux, lspcaux
real :: xref, yref
!
! 9. Subroutines calling
!
! SWCOLOUT
!
! 10. Error messages
!
!
! 11. Remarks
!
!
! 12. Structure
!
!
! 13. Source text
!
if (LTRACE) call STRACE (IENT,'swn_outnc_colspc')
ierr = 0
!
! a. Prepare spectra/auxilary data struct
!
irq = OQI(2)
spcaux%mip = MIP
allocate(spcaux%hs(MIP), spcaux%depth(MIP), spcaux%ux(MIP), spcaux%uy(MIP), &
spcaux%wndx(MIP), spcaux%wndy(MIP), spcaux%xc(MIP), spcaux%yc(MIP))
if ( LCOMPGRD .and. PUNSWAN) then
allocate(spcaux%xp(nverts), spcaux%yp(nverts))
else
allocate(spcaux%xp(MIP), spcaux%yp(MIP))
end if
spcaux%hs = OVEXCV(10)
spcaux%depth = OVEXCV(4)
spcaux%ux = OVEXCV(5)
spcaux%uy = OVEXCV(5)
spcaux%wndx = OVEXCV(26)
spcaux%wndy = OVEXCV(26)
if (RTYPE(3:3) .eq. 'R') spcaux%relative = .true.
if (RTYPE(4:4) .eq. 'C') then
allocate(spcaux%E(MSC * MDC, MIP))
spcaux%is2d = .true.
else
allocate(spcaux%E(3 * MSC, MIP))
spcaux%is2d = .false.
end if
spcaux%E = 0
spcaux%nfreq = MSC
spcaux%ndir = MDC
allocate( spcaux%f(MSC), spcaux%theta(MDC))
!
! open and read binary files
!
binfile = OUTP_FILES(irq)
! nref is unitnr set while writing the binary file.
nref = HIOPEN + IRQ
basefile = outp_files(irq)
ilpos = index(basefile, '-0')
if (ilpos.ne.0) basefile = basefile(1:ilpos-1)
!PUN ilpos = len(trim(LOCALDIR))+2
!PUN basefile = basefile(ilpos:LENFNM)
!PUN ilpos = index(LOCALDIR, '0')-1
do iproc = 1, NPROC
binnr = HIOPEN+NREOQ+(NREF-HIOPEN-1)*NPROC+iproc
!PUN write(binfile, '(A, I4.4, A, A)') LOCALDIR(1:ilpos), iproc - 1, DIRCH2, trim(basefile)
write(binfile, '(A,"-",I3.3)') trim(basefile), iproc
inquire(unit=binnr, opened=lopen)
if ( .not. lopen) then
open(file=binfile, unit=binnr, form='unformatted', IOSTAT=ierr)
if (ierr /= 0) then
write(errmsg, '("Opening ",A, " on ",I3, " failed with IOSTAT ",I3)') trim(binfile), binnr, ierr
call MSGERR(4, errmsg)
return
end if
read_header = .true.
else
read_header = .false.
end if
inquire(unit=binnr, opened=lopen)
if ( .not. lopen ) then
write(errmsg,'("unit ",I3," is NOT open")') binnr
call MSGERR(4, errmsg)
return
end if
! In the final collocation the spatial, spectral and directional grids are unknown.
! They are exchanged via the binary file of the master process.
if ( read_header .and. iproc == 1 ) then
read(binnr, IOSTAT=ierr) spcaux%xp
if ( ierr /= 0 ) then
write(errmsg,'("Reading grid data from ", A, " returned IOSTAT ", I3.3)') trim(binfile), ierr
call MSGERR(4, errmsg)
return
end if
read(binnr) spcaux%yp
if ( PUNSWAN .and. LCOMPGRD ) then
deallocate(spcaux%xp, spcaux%yp)
allocate(spcaux%xp(MIP), spcaux%yp(MIP))
spcaux%xp = xcugrdgl
spcaux%yp = ycugrdgl
end if
read(binnr) spcaux%f
read(binnr) spcaux%theta
end if
read(binnr, IOSTAT=ierr) xpctime2, tmip
if ( ierr /= 0 ) then
write(errmsg,'("Reading time and tmip from ", A, "(",I3,") , time ", A, ", returned IOSTAT ", I3.3)') &
trim(binfile), binnr, CHTIME, ierr
call MSGERR(4, errmsg)
return
end if
if ( tmip == 0 ) cycle
read(binnr) ips(1:tmip)
do tip = 1, tmip
ip = ips(tip)
read(binnr, iostat = ierr) xref, yref, spcaux%E(:, ip), spcaux%depth(ip), &
spcaux%ux(ip), spcaux%uy(ip), spcaux%hs(ip), spcaux%wndx(ip), spcaux%wndy(ip)
if ( ierr /= 0 ) then
write(errmsg,'("Reading IP ",I3.3, ", TIP ", I3.3, " in ", A, " returned IOSTAT ", I3.3)') &
ip, tip, binfile, ierr
call MSGERR(4, errmsg)
return
end if
end do ! points
end do !nproc
!
! MPI with a cold start shows some interesting initial values
! remove them
!
where(spcaux%E < epsilon(1.)) spcaux%E = 0.
!
! b. netcdf IO
!
ncfile = outp_files(irq)
ilpos = index(ncfile, '-0')
if (ilpos.ne.0) ncfile = ncfile(1:ilpos-1)
!PUN ilpos = len(trim(LOCALDIR))+2
!PUN ncfile = ncfile(ilpos:LENFNM)
!PUN ncfile = trim(GLOBALDIR)//DIRCH2//trim(ncfile)
call swn_outnc_openspecfile(ncfile, spcaux, OQI, OQR)
if (STPNOW()) return
call swn_outnc_appendspc(OQI, spcaux, xpctime2)
if (STPNOW()) return
call swn_outnc_deallocate_spcaux( spcaux )
if ( NSTATM == 0 .or. TIMCO >= TFINC ) then
call close_ncfile(OQI(1) + ncoffset(irq) )
OQI(1) = 0
! The intermediate binary files are closed and deleted in
! swanparll.ftn -> SWCOLOUT
else
call nccheck( NF90_SYNC( OQI(1) + ncoffset(irq) ))
end if
end subroutine swn_outnc_colspc
subroutine swn_outnc_spcaux(RTYPE, OQI, MIP, VOQR, VOQ, AC2, &
SPCSIG, SPCDIR, DEP2, KGRPNT, CROSS, IONOD, &
lspcaux)
USE OCPCOMM2
USE TIMECOMM, only: TINIC, TFINC, TIMCO
!PUN USE SIZES, ONLY: MYPROC
!
!
! --|-----------------------------------------------------------|--
! | BMT ARGOSS |
! | Voorsterweg 28, 8316 PT Vollenhove |
! | http://www.bmtargoss.com |
! | |
! | |
! | Programmer: A.Th.C. Hulst |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
!
! 1. Updates
!
!
! 2. Purpose
!
! Write action density spectrum to netcdf file
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! IO..name....type...size.......description
! i RTYPE char * type of output request: 'SPEC' for 2-D spectral
! 'SPE1' for 1-D freq. spectrum
! i IONOD int MIP array indicating in which subdomain output points
! are located
! i SPCDIR real MDC (*,1); spectral directions (radians)
! (*,2); cosine of spectral directions
! (*,3); sine of spectral directions
! (*,4); cosine^2 of spectral directions
! (*,5); cosine*sine of spectral directions
! (*,6); sine^2 of spectral directions
! i SPCSIG real MSC Relative frequencies in computational domain in sigma-space
! i VOQR int * Adminstrative list tracks which column in VOQ has which
! variable
! i CROSS bool 4,MIP
! i KGRPNT int MXC,MYC
! i MIP int 1 number of output points
! i VOQ real * Array with "integrated" parameters
! i AC2 real MDC,MSC, 2D spectra
! MCGRD
! i DEP2 real MCGRD
!
real, intent( in) :: SPCDIR(MDC,6), SPCSIG(MSC)
logical, intent( in) :: CROSS(1:4,1:MIP)
character (len=*), intent( in) :: RTYPE
integer, intent( in) :: OQI(4), VOQR(*), &
KGRPNT(MXC,MYC), IONOD(*), &
MIP
real, intent( in) :: VOQ(MIP,*), AC2(MDC,MSC,MCGRD), &
DEP2(MCGRD)
type(spcaux_type), intent( out) :: lspcaux
!
! 5. Local variables
!
logical :: EQREAL, lopen, do_open_files, write_header
integer :: ip, ierr, otype, xpctmp(2), xpctime, &
pnr, ri, i, tmip, irq, ips(MIP), iproc, &
binnr, xi, yi, npnts, kgrpnt1d(MXC*MYC)
integer, save :: IENT=0
character*80 :: outfile
character*256 :: errmsg
! 8. Subroutines used
!
! SWCMSP
!
logical STPNOW
!
! 9. Subroutines calling
!
! SWOUTP (SWAN/OUTP)
!
! 10. Error messages
!
! ---
!
! 11. Remarks
!
! - PUNSWAN support was added later on. In punswan, PARLL is always false but the
! IAMMASTER is set. For most output purposes, PARLL should be considered
! true. Hence a logical PUNSWAN is set at the top if this module. In the structure
! below PARLL means if PARLL .or. PUNSWAN.
!
! 12. Structure.
! a. decide wether relative and 2D spectra are required
! b. collect depth, current, xc, yc, f, theta
! c if PARLL open unformatted binary
! d. set output type (1D or 2D, relative or absolute) and allocate Energy array
! e. if PARLL, write time and number of output points in this subgrid to
! unformatted binary file
! for each point...
! f.1. compute relative / absolute 1D or 2D energy
! f.2. if PARLL spectra and auxillary data to unformatted binary file
!
! 13. Source text
!
iproc = INODE
!PUN iproc = MYPROC
if (LTRACE) call STRACE (IENT,'swn_outnc_spcaux')
lspcaux%mip = MIP
allocate(lspcaux%hs(MIP), lspcaux%depth(MIP), lspcaux%ux(MIP), lspcaux%uy(MIP), &
lspcaux%wndx(MIP), lspcaux%wndy(MIP), lspcaux%xc(MIP), lspcaux%yc(MIP), &
lspcaux%xp(MIP), lspcaux%yp(MIP))
!
! a. decide wether relative and 2D spectra are required
!
irq = OQI(2)
ierr = 0
if (RTYPE(3:3) .eq. 'R') lspcaux%relative = .true.
if (RTYPE(4:4) .eq. 'C') lspcaux%is2d = .true.
!
! b. collect auxillary data / set default (task e. set hs and wind values)
!
lspcaux%depth = VOQ(:,VOQR(4))
lspcaux%hs = 0
lspcaux%wndx = 0
lspcaux%wndy = 0
if (ICUR > 0) then
lspcaux%ux = VOQ(:,VOQR(5))
lspcaux%uy = VOQ(:,VOQR(5)+1)
else
lspcaux%ux = 0.
lspcaux%uy = 0.
end if
lspcaux%xc = OVEXCV(1)
lspcaux%yc = OVEXCV(2)
lspcaux%xp = OVEXCV(1)
lspcaux%yp = OVEXCV(2)
do ip = 1, MIP
if (OPTG.ne.5) then
if (.not.EQREAL(VOQ(ip,VOQR(24)),OVEXCV(1) )) then
lspcaux%xc(ip) = VOQ(ip,VOQR(24))
lspcaux%yc(ip) = VOQ(ip,VOQR(25))
end if
else
if (.not.EQREAL(VOQ(ip,1),OVEXCV(1))) then
lspcaux%xc(ip) = VOQ(ip,1) - XOFFS
lspcaux%yc(ip) = VOQ(ip,2) - YOFFS
end if
end if
lspcaux%xp(ip) = VOQ(ip,1)
lspcaux%yp(ip) = VOQ(ip,2)
end do
lspcaux%nfreq = MSC
lspcaux%ndir = MDC
allocate( lspcaux%f(MSC), lspcaux%theta(MDC))
lspcaux%f = SPCSIG/PI2
if ( BNAUT ) then
! meteorological convention
lspcaux%theta = PI + DNORTH*PI/180 - SPCDIR(:,1)
else
lspcaux%theta = SPCDIR(:,1)
end if
!
! c if PARLL or PUNSWAN open unformatted binary if not already opened
!
if ( PARLL .or. PUNSWAN) then
if ( oqi(1) == 0 ) then
call FOR(oqi(1), OUTP_FILES(irq), 'UU', ierr)
if ( ierr > 0 ) then
write(errmsg,'("File ", A, " open returned IOSTAT ", I3)') trim(OUTP_FILES(irq)), ierr
call MSGERR(4, errmsg)
return
end if
write_header = .true.
else
write_header = .false.
end if
end if
! d. set output type (1D or 2D, relative or absolute) and allocate Energy array
if ( lspcaux%is2d ) then
allocate(lspcaux%E(MSC * MDC, MIP))
if ( lspcaux%relative ) then
otype = 2
else
otype = -2
end if
else
allocate(lspcaux%E(3 * MSC, MIP))
if ( lspcaux%relative ) then
otype = 1
else
otype = -1
end if
end if
lspcaux%E = 0
! e. write time and number of output points in this subgrid to unformatted binary file
if ( PARLL .or. PUNSWAN) then
tmip = 0
! if part of this subgrid and an active point
do ip = 1, MIP
if ( IONOD(ip).eq.IPROC ) then
tmip = tmip + 1
ips(tmip) = ip
end if
end do
if ( NSTATM > 0 ) then
READ (chtime, '(I8,1X,I6)') (xpctmp(i), i=1,2)
xpctime = seconds_since_epoch(xpctmp(1), xpctmp(2))
else
xpctime = 0
end if
! In the final collocation the spatial, spectral and directional grids are unknown.
! They are exchanged via the binary file of the master on the first timestep.
!
! note the statement below and its reading counterpart imply that the master is
! always IPROC 1
if ( write_header .and. IAMMASTER ) then
write(OQI(1)) lspcaux%xp
write(OQI(1)) lspcaux%yp
write(OQI(1)) lspcaux%f
write(OQI(1)) lspcaux%theta
end if
write(OQI(1), IOSTAT=ierr) xpctime, tmip
if (ierr /= 0) then
write(errmsg, '(A, "-> IOSTAT=", I3, " time=", I10, " TMIP=", I4, " OQI(1)=", I3)') &
trim(outfile), ierr, xpctime, tmip, oqi(1)
call MSGERR(4, errmsg)
return
end if
if ( tmip /= 0 ) write(OQI(1)) ips
end if
! for each point, compute spectra and collect remaining auxilary
do ip = 1, MIP
if ( (.not. PARLL .and. .not. PUNSWAN) .or. &
IONOD(ip).eq.INODE ) then
!PUN IONOD(ip).eq.MYPROC ) then
! if depth > 0 and not a dummy value
if ( lspcaux%depth(ip) > epsilon(1.) .and. abs(lspcaux%depth(ip) - OVEXCV(4)) > epsilon(1.) ) then
call SWCMSP (otype, lspcaux%xc(ip), lspcaux%yc(ip), AC2, lspcaux%E(:, ip), SPCSIG, &
lspcaux%depth(ip) , dep2, lspcaux%ux(ip), lspcaux%uy(ip), &
SPCDIR(1,2), SPCDIR(1,3), 1., KGRPNT, &
CROSS(1,IP), ierr)
! when ierr > 0 then energy == 0. That's fine.
lspcaux%hs(ip) = VOQ(ip,VOQR(10))
lspcaux%wndx(ip) = VOQ(ip,VOQR(26))
lspcaux%wndy(ip) = VOQ(ip,VOQR(26)+1)
end if
if ( PARLL .or. PUNSWAN ) then
! write to binary
write(OQI(1)) lspcaux%xp(ip), lspcaux%yp(ip), lspcaux%E(:, ip), lspcaux%depth(ip), &
lspcaux%ux(ip), lspcaux%uy(ip), lspcaux%hs(ip), lspcaux%wndx(ip), &
lspcaux%wndy(ip)
end if
end if
end do
if ( PARLL .or. PUNSWAN ) then
! in parallel mode, close binary file when required
if (NSTATM == 0 .or. TIMCO >= TFINC) then
inquire(unit=OQI(1), opened=lopen)
if ( lopen ) then
close(unit=OQI(1))
else
write(errmsg, '("File ", A, " is already closed. This is a code error")') trim(OUTP_FILES(irq))
call MSGERR(4, errmsg)
return
end if
end if
end if
end subroutine swn_outnc_spcaux
subroutine swn_outnc_spcaux_on_wetnodes( lkgrpnt, lspcaux, spcaux )
! lkgrpnt is KGRPNT when called from shared memory mode or KGRPGL
! from colspc in the case of distributed memory. In both cases, it covers
! the full grid
integer, intent( in) :: lkgrpnt(MXCGL,MYCGL)
type(spcaux_type), intent( in) :: lspcaux
type(spcaux_type), intent( out) :: spcaux
integer :: ip, indx, npnts, ix, iy
integer, save :: IENT=0
character*256 :: errmsg
if (LTRACE) call STRACE(IENT,'swn_outnc_spcaux_on_wetnodes')
npnts = MCGRDGL - 1
allocate(spcaux%hs(npnts), spcaux%depth(npnts), spcaux%ux(npnts), spcaux%uy(npnts), &
spcaux%wndx(npnts), spcaux%wndy(npnts), spcaux%xp(npnts), spcaux%yp(npnts), &
spcaux%xc(npnts), spcaux%yc(npnts), spcaux%ips(npnts), &
spcaux%E(size(lspcaux%E,1), npnts) )
do ix = 1, MXCGL
do iy = 1, MYCGL
ip = (iy-1)*MXCGL + ix
if ( lkgrpnt(ix,iy) /= 1 ) then
indx = lkgrpnt(ix,iy) - 1
spcaux%ips(indx) = ip
spcaux%hs(indx) = lspcaux%hs(ip)
spcaux%depth(indx) = lspcaux%depth(ip)
spcaux%ux(indx) = lspcaux%ux(ip)
spcaux%uy(indx) = lspcaux%uy(ip)
spcaux%wndx(indx) = lspcaux%wndx(ip)
spcaux%wndy(indx) = lspcaux%wndy(ip)
spcaux%xc(indx) = lspcaux%xc(ip)
spcaux%yc(indx) = lspcaux%yc(ip)
spcaux%xp(indx) = lspcaux%xp(ip)
spcaux%yp(indx) = lspcaux%yp(ip)
spcaux%E(:,indx) = lspcaux%E(:,ip)
end if
end do
end do
if ( indx /= npnts ) then
write(errmsg, '("Inconsistent number of wet nodes:", i5, i5)') indx, npnts
call MSGERR(4, errmsg)
return
end if
spcaux%is2d = lspcaux%is2d
spcaux%relative = lspcaux%relative
spcaux%nfreq = lspcaux%nfreq
spcaux%ndir = lspcaux%ndir
allocate( spcaux%f(spcaux%nfreq), spcaux%theta(spcaux%ndir))
spcaux%f = lspcaux%f
spcaux%theta = lspcaux%theta
spcaux%mip = npnts
end subroutine swn_outnc_spcaux_on_wetnodes
subroutine swn_outnc_deallocate_spcaux( lspcaux )
type(spcaux_type), intent(inout) :: lspcaux
integer, save :: IENT=0
if (LTRACE) call STRACE (IENT,'swn_outnc_deallocate_spcaux')
! yesyes, derived type fields are deallocated when out of scope
! but sometimes variables stay in scope
if (allocated(lspcaux%hs)) deallocate(lspcaux%hs)
if (allocated(lspcaux%depth)) deallocate(lspcaux%depth)
if (allocated(lspcaux%ux)) deallocate(lspcaux%ux)
if (allocated(lspcaux%uy)) deallocate(lspcaux%uy)
if (allocated(lspcaux%wndx)) deallocate(lspcaux%wndx)
if (allocated(lspcaux%wndy)) deallocate(lspcaux%wndy)
if (allocated(lspcaux%xc)) deallocate(lspcaux%xc)
if (allocated(lspcaux%yc)) deallocate(lspcaux%yc)
if (allocated(lspcaux%xp)) deallocate(lspcaux%xp)
if (allocated(lspcaux%yp)) deallocate(lspcaux%yp)
if (allocated(lspcaux%f)) deallocate(lspcaux%f)
if (allocated(lspcaux%theta)) deallocate(lspcaux%theta)
if (allocated(lspcaux%E)) deallocate(lspcaux%E)
if (allocated(lspcaux%ips)) deallocate(lspcaux%ips)
lspcaux%mip = 0
lspcaux%ndir = 0
lspcaux%nfreq = 0
end subroutine swn_outnc_deallocate_spcaux
subroutine swn_outnc_appendspc(oqi, spcaux, xpctime2)
USE OCPCOMM2
!
!
! --|-----------------------------------------------------------|--
! | BMT ARGOSS |
! | Voorsterweg 28, 8316 PT Vollenhove |
! | http://www.bmtargoss.com |
! | |
! | |
! | Programmer: A.Th.C. Hulst |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
!
! 1. Updates
!
!
! 2. Purpose
!
! Append spectra and auxilary variables to netcdf file.
!
! 3. Method
!
! ---
!
! 4. Argument variables
!
! IO..name....type...size.......description
!
integer, intent( in) :: oqi(4)
type(spcaux_type), target, intent(inout) :: spcaux
integer, optional, intent( in) :: xpctime2
!
! 5. Local variables
!
integer :: ierr, xpctmp(2), xpctime, pnr, ri, irq, i, &
ncid
character*256 :: errmsg
integer, save :: IENT=0
logical :: spc_as_map, wetnode_list, noaux
type(spcaux_type), target :: lspcaux
type(spcaux_type), pointer :: pspcaux
! 8. Subroutines used
!
! SWCMSP
!
logical STPNOW
!
! 9. Subroutines calling
!
! SWOUTP (SWAN/OUTP)
!
! 10. Error messages
!
! ---
!
! 11. Remarks
! Auxillary data is retrieved from the VOQ array. Please modify SWORDC when
! VOQ need to be extended
!
! 12. Structure
! a. Determine the index in the record axe of ncfile
! b. add spectra to ncfile
! c. add auxilary variables (depth, current and Hs) to ncfile
!
! 13. Source text
!
if (LTRACE) call STRACE (IENT,'swn_outnc_appendspc')
irq = oqi(2)
ncid = oqi(1) + ncoffset(irq)
call swn_outnc_spcflags(oqi(4), spc_as_map, wetnode_list, spcaux, noaux=noaux)
!
! Determine the index in the record axe
!
if ( recordaxe(irq)%nstatm ) then
READ (chtime, '(I8,1X,I6)') (xpctmp(i), i=1,2)
xpctime = seconds_since_epoch(xpctmp(1), xpctmp(2))
! If provided, check the time against a reference time
if ( present(xpctime2) ) then
if ( xpctime2 /= xpctime ) then
write(errmsg, '(I12, " does not equal ", I12)') xpctime2, xpctime
call MSGERR(4, errmsg)
return
end if
end if
call record_index(ncid, recordaxe(irq), xpctime, ri)
else
if ( string_is_int(PROJNR) ) then
read( PROJNR, '(i16)' ) pnr
if ( pnr == 0 ) then
call MSGERR(3, 'The projnr number may not be zero')
return
end if
else
pnr = 1
end if
call record_index(ncid, recordaxe(irq), pnr, ri )
end if
if ( wetnode_list ) then
! remove dry points. In memory there's a copy of spcaux with dry points and
! lspcaux without them. If this remains problematic than we need to do something clever
! with KGRPGL at write time
call swn_outnc_spcaux_on_wetnodes(KGRPGL, spcaux, lspcaux )
if (STPNOW()) return
! pointer to wetnode list
pspcaux => lspcaux
! deallocate list with dry points
call swn_outnc_deallocate_spcaux( spcaux )
else
pspcaux => spcaux
end if
!
! add spectra
!
where(pspcaux%E < epsilon(1.)) pspcaux%E = 0.
if ( pspcaux%is2d ) then
! factor 2*PI to account for transition from rad/s to Hz
pspcaux%E = pspcaux%E * 2 * pi
call agnc_add_spcdata_density(ncid, ri, pspcaux%E, spc_as_map)
else
call agnc_add_spcdata(ncid, ri, pspcaux%E(1:MSC,:), &
pspcaux%E(MSC+1:2*MSC,:), &
pspcaux%E(2*MSC+1:3*MSC,:), &
spc_as_map)
end if
!
! add auxilary variables (depth, current and Hs).
!
! hs is provided so that user can check his integration routines. His / hers
! hs should match the one computed by SWAN.
if ( .not.noaux ) then
if ( spc_as_map ) then
call agnc_add_mapdata(ncid, 'depth', ri, pspcaux%depth, OVEXCV( 4), skip_range_error)
call agnc_add_mapdata(ncid, 'xcur', ri, pspcaux%ux, OVEXCV( 5), skip_range_error)
call agnc_add_mapdata(ncid, 'ycur', ri, pspcaux%uy, OVEXCV( 5), skip_range_error)
call agnc_add_mapdata(ncid, 'hs', ri, pspcaux%hs, OVEXCV(10), skip_range_error)
call agnc_add_mapdata(ncid, 'xwnd', ri, pspcaux%wndx, OVEXCV(26), skip_range_error)
call agnc_add_mapdata(ncid, 'ywnd', ri, pspcaux%wndy, OVEXCV(26), skip_range_error)
else
call agnc_add_pntdata(ncid, 'depth', ri, pspcaux%depth, OVEXCV( 4), skip_range_error)
call agnc_add_pntdata(ncid, 'xcur', ri, pspcaux%ux, OVEXCV( 5), skip_range_error)
call agnc_add_pntdata(ncid, 'ycur', ri, pspcaux%uy, OVEXCV( 5), skip_range_error)
call agnc_add_pntdata(ncid, 'hs', ri, pspcaux%hs, OVEXCV(10), skip_range_error)
call agnc_add_pntdata(ncid, 'xwnd', ri, pspcaux%wndx, OVEXCV(26), skip_range_error)
call agnc_add_pntdata(ncid, 'ywnd', ri, pspcaux%wndy, OVEXCV(26), skip_range_error)
end if
end if
! garbage collect
nullify(pspcaux)
if ( wetnode_list ) then
call swn_outnc_deallocate_spcaux( lspcaux )
else
call swn_outnc_deallocate_spcaux( spcaux )
end if
end subroutine swn_outnc_appendspc
subroutine swn_outnc_openspecfile(ncfile, spcaux, oqi, oqr)
!
! Structure:
!
! 1) input arguments
! 2) local variables
! if file does not exist
! 3) Define new record axe (run or time)
! If spectra requested on COMPGRID
! if not rotated
! 4.a) define mspcgrid with 1D coordinate axes
! else
! 4.b) define mspcgrid with 2D coordinate axes
! end
! 4.c) create ncfile
! else
! 5.a) define spcgrid
! if WETCGRID
! 5.b) filter pointlist
! end if
! 5.b) create ncfile
! end if
! 6.a) add global attributes
! 6.b) create variables
! 6.c) close define mode and reopen write
! 6.d) fill dimension variables and close
! end if file does not exist
! 7) open file in write mode
!
! 1) input arguments
!
character*80, intent( in) :: ncfile
type(spcaux_type), target, intent( in) :: spcaux
integer, intent( in) :: oqi(4)
real*8, intent( in) :: oqr(2)
!
! 2) local variables
!
logical :: file_exists, nc_debug, &
Escale, monthly, spc_as_map, &
wetnode_list, STPNOW, noaux
integer :: ncid, xpctmp(2), i, irq
type(spcgrid_type) :: spcgrid
type(mapgrid_type) :: mapgrid
type(pntgrid_type) :: pntgrid
type(spcaux_type), target :: lspcaux
type(spcaux_type), pointer :: pspcaux
integer, save :: IENT=0
file_exists = .false.
nc_debug = .false.
if (LTRACE) call STRACE (IENT,'swn_outnc_openspecfile')
call swn_outnc_spcflags(oqi(4), spc_as_map, wetnode_list, spcaux, Escale, noaux, monthly)
irq = OQI(2)
if ( .not. stnames_initialized ) call stnames_init()
inquire( FILE=ncfile, EXIST=file_exists )
!
! if file does not exists
!
if ( .not. file_exists) then
!
! 3) Define new record axe (run or time) and spectral grid
!
allocate(recordaxe(irq)%content(1))
recordaxe(irq)%ncontent = 1
if ( NSTATM > 0 ) then
READ (chtime, '(I8,1X,I6)') (xpctmp(i), i=1,2)
! should be seconds since epoch.
recordaxe(irq)%content(1) = seconds_since_epoch(xpctmp(1), xpctmp(2))
recordaxe(irq)%delta = maxval((/oqr(2), DT/))
recordaxe(irq)%nstatm = .true.
else
! stationary, record axe is "run"
if ( string_is_int(PROJNR) ) then
read( PROJNR, '(i16)' ) recordaxe(irq)%content(1)
else
recordaxe(irq)%content(1) = 1
end if
recordaxe(irq)%delta = 1
recordaxe(irq)%nstatm = .false.
end if
allocate (spcgrid%frequency(spcaux%nfreq))
spcgrid%nfreq = spcaux%nfreq
spcgrid%frequency = spcaux%f
if ( spcaux%is2d ) then
spcgrid%ndir = spcaux%ndir
allocate (spcgrid%direction(spcaux%ndir))
spcgrid%direction = spcaux%theta
else
spcgrid%ndir = 0
end if
spcgrid%relative = spcaux%relative
!
! If spectra requested on COMPGRID
!
if ( spc_as_map .and. OPTG /= 5 ) then
! COMPGRID is a special case. By specifying this spectra are 4D
! written to the netcdf instead of only the wet points
!
! if not rotated
!
! XGRDGL, YGRDGL, MXCGL and MYCGL from M_PARALL
if ( abs(ALPC) < epsilon(1.) .and. OPTG == 1 ) then
!
! 4.a) define mapgrid with 1D coordinate axes
!
allocate (mapgrid%longitude(MXCGL,1))
allocate (mapgrid%latitude (1,MYCGL))
if ( PARLL ) then
mapgrid%longitude(:,1) = xgrdgl(:,1) + XOFFS
mapgrid%latitude(1,:) = ygrdgl(1,:) + YOFFS
else
mapgrid%longitude(:,1) = xgrdgl(:,1)
mapgrid%latitude(1,:) = ygrdgl(1,:)
end if
mapgrid%mdc = .false.
mapgrid%alpc = 0.
else
!
! 4.b) define mapgrid with 2D coordinate axes
!
allocate (mapgrid%longitude(MXCGL, MYCGL))
allocate (mapgrid%latitude (MXCGL, MYCGL))
if ( PARLL ) then
mapgrid%longitude = XGRDGL + XOFFS
mapgrid%latitude = YGRDGL + YOFFS
else
mapgrid%longitude = XGRDGL
mapgrid%latitude = YGRDGL
end if
mapgrid%mdc = .true.
mapgrid%alpc = ALPC
end if
! curvilinear grids may contain the EXCEPTION value of the bottom
! Replace with the hardcoded float _FillValue
where ( abs(mapgrid%longitude - EXCFLD(1)) < epsilon(1.) ) &
mapgrid%longitude = NF90_FILL_FLOAT
where ( abs(mapgrid%latitude - EXCFLD(1)) < epsilon(1.) ) &
mapgrid%latitude = NF90_FILL_FLOAT
if ( KSPHER == 0 ) mapgrid%lunit = 'meter'
mapgrid%nx = MXCGL
mapgrid%ny = MYCGL
!
! 4.c) create ncfile
!
! Definition mode
call create_ncfile(ncfile, ncid, recordaxe(irq), spcgrid=spcgrid, mapgrid=mapgrid, Escale=Escale, nautical=BNAUT)
else
!
! 5.a) define pntgrid, eg, spectra along a point list
!
if ( wetnode_list ) then
! remove dry points
call swn_outnc_spcaux_on_wetnodes(KGRPGL, spcaux, lspcaux )
if (STPNOW()) return
pspcaux => lspcaux
else
pspcaux => spcaux
end if
allocate (pntgrid%longitude(pspcaux%mip))
allocate (pntgrid%latitude(pspcaux%mip))
pntgrid%npoints = pspcaux%mip
pntgrid%longitude = pspcaux%xp
pntgrid%latitude = pspcaux%yp
if ( wetnode_list ) then
allocate(pntgrid%ips(pspcaux%mip))
pntgrid%ips = pspcaux%ips
pntgrid%xdimlen = MXCGL
pntgrid%ydimlen = MYCGL
end if
nullify(pspcaux)
if ( wetnode_list ) then
call swn_outnc_deallocate_spcaux( lspcaux )
end if
if ( KSPHER == 0 ) pntgrid%lunit = 'meter'
!
! 5.b) create point list ncfile
!
! Definition mode
call create_ncfile(ncfile, ncid, recordaxe(irq), spcgrid=spcgrid, pntgrid=pntgrid, Escale=Escale, nautical=BNAUT)
end if
!
! 6.a) add global attributes
!
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'project', PROJID) )
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'model', VERTXT) )
if ( NSTATM > 0 ) then
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'run', PROJNR) )
end if
!
! 6.b) create variables
!
if ( .not.noaux ) call create_auxvariables(ncid, spc_as_map)
! 6.c) close define mode and reopen write
!
! write mode, e.g, fill dimension variables. Note that record dimension
! variable is updated by open_ncfile
call close_ncfile(ncid)
call open_ncfile( ncfile, "write", ncid, recordaxe(irq), monthly)
!
! 6.d) fill dimension variables and close
!
call agnc_set_spcgrid (ncid, spcgrid)
if ( spc_as_map .and. OPTG /= 5 ) then
call agnc_set_mapgrid (ncid, mapgrid)
else
call agnc_set_pntgrid (ncid, pntgrid)
end if
! close the freshly created file so all dimension variable data is written
call close_ncfile(ncid)
end if
!
! 7) open file in write mode
!
call open_ncfile( ncfile, "write", ncid, recordaxe(irq))
! if only only record found, set delta to delta specified in command file (or 1 on stationary)
if ( recordaxe(irq)%ncontent == 1 ) then
if ( NSTATM > 0 ) then
recordaxe(irq)%delta = maxval((/oqr(2), DT/))
else
recordaxe(irq)%delta = 1
end if
end if
ncoffset(irq) = ncid - oqi(1)
end subroutine swn_outnc_openspecfile
subroutine swn_outnc_spcflags(oqi, spc_as_map, wetnode_list, spcaux, Escale, noaux, monthly)
integer, intent( in) :: oqi
type(spcaux_type), intent( in) :: spcaux
logical, intent( out) :: spc_as_map, wetnode_list
logical, optional, intent( out) :: Escale, monthly, noaux
integer :: iocode
iocode = oqi
wetnode_list = .false.
spc_as_map = .false.
noaux = .false.
if ( iocode >= 16 ) then
noaux = .true.
iocode = iocode - 16
end if
! MDGRID is the special keyword to store spectra (including dry points) on the
! full geographic grid instead of along a pointlist
! On the other hand, when COMPGRID is requested on any grid but unstructured,
! store the point indices in the computational grid in the netcdf file
if ( spcaux%mip == MCGRDGL-1 .or. spcaux%mip == (MXCGL * MYCGL) ) then
if (iocode >= 8 ) then
spc_as_map = .true.
iocode = iocode - 8
elseif ( iocode >= 4 .and. OPTG /= 5 ) then
! request to remove dry points from COMPGRID pointlist, e.g., compress
wetnode_list = .true.
end if
end if
! catch conflicting MDGRID + COMPRESS
iocode = mod(iocode, 4)
! Enery scaling is of creating a new variable that stores a scale per timestep and point
! the energy itself is then stored as int16 instead of float
if ( present( EScale) ) then
Escale = .false.
if ( iocode > 1) Escale = .true.
end if
! By setting the keyword MONTH a netcdf file is created that has a time axis
! starting at the first of the month.
if ( present( monthly ) ) then
monthly = .false.
if ( mod(iocode,2) == 1 .and. NSTATM == 1) monthly = .true.
end if
end subroutine swn_outnc_spcflags
subroutine swn_outnc_appendblock(myk, mxk, ivtype, nref, irq, data, excv, col)
! note that ivtyp is an array and ivtype is ivtyp(JVAR)
! col is used for writing vector data.
integer, intent( in) :: irq, ivtype, nref, col, myk, mxk
real, intent( in) :: data(mxk * myk), excv
integer :: ri, pnr, xpctmp(2), i, &
xpctime, ilpos, ncid
integer, save :: IENT=0
if (LTRACE) call STRACE (IENT,'swn_outnc_appendblock')
if ( ivtype == 40 ) return
ncid = nref + ncoffset(irq)
if ( recordaxe(irq)%nstatm ) then
READ (chtime, '(I8,1X,I6)') (xpctmp(i), i=1,2)
xpctime = seconds_since_epoch(xpctmp(1), xpctmp(2))
call record_index(ncid, recordaxe(irq), xpctime, ri)
else
if ( string_is_int(PROJNR) ) then
read( PROJNR, '(i16)' ) pnr
if ( pnr == 0 ) then
call MSGERR(3, 'The projnr number may not be zero')
return
end if
else
pnr = 1
end if
call record_index(ncid, recordaxe(irq), pnr, ri )
end if
if (myk == 1) then
call agnc_add_pntdata(ncid, STNAMES(ivtype,col), ri, data, excv, skip_range_error)
else
call agnc_add_mapdata(ncid, STNAMES(ivtype,col), ri, data, excv, skip_range_error)
end if
end subroutine swn_outnc_appendblock
subroutine swn_outnc_close_on_end(nref, irq)
integer, intent(inout) :: nref
integer, optional, intent( in) :: irq
integer :: ncid
integer, save :: IENT=0
ncid = nref
if ( present(irq) ) ncid = nref + ncoffset(irq)
! Decided to sync data to file even in non stationary mode because
! otherwise data is buffered during the entire computation,
! increasing the chance on file corruption
if (LTRACE) call STRACE (IENT,'swn_outnc_close_on_end')
if ( nstatm == 0 .or. TIMCO >= TFINC ) then
call close_ncfile(ncid)
nref = 0
else
call nccheck( NF90_SYNC(ncid) )
end if
end subroutine swn_outnc_close_on_end
subroutine record_index(ncid, recordaxe, xpctime, ri)
integer, intent( in) :: ncid, xpctime
type(recordaxe_type), intent(inout) :: recordaxe
integer, intent( out) :: ri
integer, allocatable, dimension(:) :: tmp_time
integer, save :: IENT=0
if (LTRACE) call STRACE (IENT,'record_index')
ri = timeindex(recordaxe%content, xpctime)
if ( ri < 0 ) then
call MSGERR(4, 'agioncmd cannot prepend data, only overwrite '// &
'and append. Consider setting option in order '// &
'to generate monthly files. You''ll need to rerun '// &
'the first period too though')
call close_ncfile(ncid)
end if
if ( ri == 0) then
! record needs to be appended
allocate(tmp_time(recordaxe%ncontent))
tmp_time = recordaxe%content
deallocate(recordaxe%content);
allocate(recordaxe%content(recordaxe%ncontent+1))
recordaxe%content(1:recordaxe%ncontent) = tmp_time
recordaxe%content(recordaxe%ncontent+1) = xpctime
recordaxe%ncontent = recordaxe%ncontent+1
!
! Update netcdf file with new time axe
!
call agnc_set_recordaxe (ncid, recordaxe)
ri = recordaxe%ncontent
end if
end subroutine record_index
subroutine swn_outnc_openblockfile(ncfile, myk, mxk, ovlnam, xgrdgl, &
ygrdgl, oqi, oqr, ivtyp , irq )
character*80, intent( in) :: ncfile
integer, intent( in) :: myk, mxk, irq
integer, dimension(:), intent( in) :: oqi, ivtyp
real*8, dimension(:), intent( in) :: oqr
character*40, intent( in) :: ovlnam(:)
real, intent( in) :: xgrdgl(mxk, myk), &
ygrdgl(mxk, myk)
! local variables
logical :: file_exists = .false., &
nc_debug = .false., &
monthly = .false.
integer :: ncid, xpctmp(2), i
type(mapgrid_type) :: mapgrid
type(pntgrid_type) :: pntgrid
integer, save :: IENT=0
if (LTRACE) call STRACE (IENT,'swn_outnc_openblockfile')
if ( .not. stnames_initialized ) call stnames_init()
! By setting the idla to 5 a netcdf file is created that has a time axis
! that starts at the first of the month. Monthly is only valid on non-stationary
! run that has time axis. Declaring a monthly file is required when computing
! hindcasts in several chunks in time.
if ( oqi(4) == 5 .and. nstatm == 1 ) monthly = .true.
inquire( FILE=ncfile, EXIST=file_exists )
if ( .not. file_exists) then
allocate(recordaxe(irq)%content(1))
recordaxe(irq)%ncontent = 1
if ( nstatm > 0 ) then
READ (chtime, '(I8,1X,I6)') (xpctmp(i), i=1,2)
! should be seconds since epoch. This is probably not it.
recordaxe(irq)%content(1) = seconds_since_epoch(xpctmp(1), xpctmp(2))
recordaxe(irq)%delta = maxval((/oqr(2), DT/))
recordaxe(irq)%nstatm = .true.
else
! stationary, record axe is "run"
if ( string_is_int(PROJNR) ) then
read( PROJNR, '(i16)' ) recordaxe(irq)%content(1)
else
recordaxe(irq)%content(1) = 1
end if
recordaxe(irq)%delta = 1
recordaxe(irq)%nstatm = .false.
end if
if (myk == 1) then
! unstructured mesh is a pointlist. TABLE output also uses this section
pntgrid%npoints = mxk
allocate (pntgrid%longitude(mxk))
allocate (pntgrid%latitude(mxk))
pntgrid%latitude = ygrdgl(:,1)
pntgrid%longitude = xgrdgl(:,1)
if ( KSPHER == 0 ) pntgrid%lunit = 'meter'
! Definition mode
call create_ncfile(ncfile, ncid, recordaxe(irq), pntgrid=pntgrid, nautical=BNAUT)
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'project', PROJID) )
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'run', PROJNR) )
call create_variables(ncid, ivtyp, ovlnam, .false.)
! write mode, e.g, fill dimension variables. Note that record dimension
! variable is updated by open_ncfile
call close_ncfile(ncid)
call open_ncfile( ncfile, "write", ncid, recordaxe(irq), monthly)
call agnc_set_pntgrid (ncid, pntgrid)
else
! Regular mesh. Depending on the angle of the computational grid,
! longitude and latitude can be a vector or a matrix
! curvilinear grids always have multi-dimension coordinate fields
! Definition mode
if ( abs(ALPC) < epsilon(1.) .and. OPTG == 1 ) then
allocate (mapgrid%longitude(mxk,1))
allocate (mapgrid%latitude (1,myk))
if ( PARLL ) then
mapgrid%longitude(:,1) = xgrdgl(:,1) + XOFFS
mapgrid%latitude(1,:) = ygrdgl(1,:) + YOFFS
else
mapgrid%longitude(:,1) = xgrdgl(:,1)
mapgrid%latitude(1,:) = ygrdgl(1,:)
end if
mapgrid%mdc = .false.
else
allocate (mapgrid%longitude(mxk, myk))
allocate (mapgrid%latitude (mxk, myk))
if ( PARLL ) then
mapgrid%longitude = xgrdgl + XOFFS
mapgrid%latitude = ygrdgl + YOFFS
else
mapgrid%longitude = xgrdgl
mapgrid%latitude = ygrdgl
end if
mapgrid%mdc = .true.
end if
! curvilinear grids may contain the EXCEPTION value of the bottom
! Replace with the hardcoded float _FillValue
where ( abs(mapgrid%longitude - EXCFLD(1)) < epsilon(1.) ) &
mapgrid%longitude = NF90_FILL_FLOAT
where ( abs(mapgrid%latitude - EXCFLD(1)) < epsilon(1.) ) &
mapgrid%latitude = NF90_FILL_FLOAT
mapgrid%nx = mxk
mapgrid%ny = myk
mapgrid%alpc = ALPC
if ( KSPHER == 0 ) mapgrid%lunit = 'meter'
call create_ncfile(ncfile, ncid, recordaxe(irq), mapgrid=mapgrid, nautical=BNAUT)
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'project', PROJID) )
call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'run', PROJNR) )
call create_variables(ncid, ivtyp, ovlnam, .true.)
! write mode, e.g, fill dimension variables. Note that record dimension
! variable is updated by open_ncfile
call close_ncfile(ncid)
call open_ncfile( ncfile, "write", ncid, recordaxe(irq), monthly)
call agnc_set_mapgrid (ncid, mapgrid)
end if
! close the freshly created file so all buffers are emptied
call close_ncfile(ncid)
end if
call open_ncfile( ncfile, "write", ncid, recordaxe(irq))
! set the offset to translate between swan's HIOPEN + IRQ and
! the netcdf 65xxx range
ncoffset(irq) = ncid - oqi(1)
! if only only record found, set delta to delta specified in command file (or 1 on stationary)
if ( recordaxe(irq)%ncontent == 1 ) then
if ( NSTATM > 0 ) then
recordaxe(irq)%delta = maxval((/oqr(2), DT/))
else
recordaxe(irq)%delta = 1
end if
end if
end subroutine swn_outnc_openblockfile
subroutine create_auxvariables(ncid, spc_as_map)
integer, intent( in) :: ncid
logical, intent( in) :: spc_as_map
integer :: ivtyp(4)
character*40 :: dummy(4)
ivtyp = (/4, 5, 10, 26/)
call create_variables(ncid, ivtyp, dummy, spc_as_map)
end subroutine create_auxvariables
subroutine create_variables(ncid, ivtyp, ovlnam, spc_as_map)
integer, intent( in) :: ncid
integer,dimension(:), intent( in) :: ivtyp
character*40, intent( in) :: ovlnam(:)
logical, intent( in) :: spc_as_map
integer :: i
character*40 :: stname
integer, save :: IENT=0
if (LTRACE) call STRACE (IENT,'create_variables')
do i=1, size(ivtyp)
stname = STNAMES(ivtyp(i),1)
if ( stname /= '' ) then
if ( spc_as_map ) then
call agnc_define_mapvariable( ncid, stname)
else
call agnc_define_pntvariable( ncid, stname)
end if
! if vector, write second variable
stname = STNAMES(ivtyp(i),2)
if ( stname /= '' ) then
if ( spc_as_map ) then
call agnc_define_mapvariable( ncid, stname)
else
call agnc_define_pntvariable( ncid, stname)
end if
end if
else
! See list at the end of the source file
if ( ivtyp(i) > 3) call MSGERR(1, "Field '"//trim(ovlnam(ivtyp(i)))//"' not supported (yet).")
end if
end do
end subroutine create_variables
subroutine stnames_init()
STNAMES( 4, 1) = 'depth'
STNAMES( 5, 1) = 'xcur'
STNAMES( 5, 2) = 'ycur'
STNAMES( 6, 1) = 'ubot'
STNAMES(10, 1) = 'hs'
STNAMES(11, 1) = 'tm01'
STNAMES(12, 1) = 'tp'
STNAMES(13, 1) = 'theta0'
STNAMES(14, 1) = 'thetap'
STNAMES(16, 1) = 'spread'
STNAMES(17, 1) = 'L'
STNAMES(26, 1) = 'xwnd'
STNAMES(26, 2) = 'ywnd'
STNAMES(28, 1) = 'rtm01'
STNAMES(30, 1) = 'dhs'
STNAMES(31, 1) = 'dtm'
STNAMES(32, 1) = 'tm02'
STNAMES(34, 1) = 'urms'
STNAMES(35, 1) = 'ustar'
STNAMES(38, 1) = 'cdrag'
STNAMES(44, 1) = 'hswe'
STNAMES(47, 1) = 'tmm10'
STNAMES(48, 1) = 'rtmm10'
STNAMES(51, 1) = 'ssh'
STNAMES(52, 1) = 'botl'
STNAMES(53, 1) = 'tps'
stnames_initialized = .true.
end subroutine stnames_init
function string_is_int(str) result (is_int)
character(len=*), intent( in) :: str
logical :: is_int
character(len=10) :: nums = "0123456789"
is_int = .true.
if ( verify(str, nums) > 0 ) is_int = .false.
end function string_is_int
end module swn_outnc
! ( 3) [ 'distance along output curve']
! ( 7) [ 'Energy dissipation']
! ( 8) [ 'Fraction breaking waves']
! ( 9) [ 'Energy leak over spectral boundaries']
! (15) [ 'direction of the energy transport']
! (18) [ 'Wave steepness']
! (19) [ 'Wave energy transport']
! (20) [ 'Wave driven force per unit surface']
! (21) [ 'spectral action density']
! (22) [ 'spectral energy density']
! (23) [ 'auxiliary variable']
! (24) [ 'X computational grid coordinate']
! (25) [ 'Y computational grid coordinate']
! (27) [ 'Bottom friction coefficient']
! (29) [ 'energy density integrated over direction']
! (33) [ 'Frequency spectral width (Kappa)']
! (35) [ 'Friction velocity']
! (40) [ 'Date-time']
! (41) [ 'Time in seconds from reference time']
! (36) ['Zero velocity thickness of boundary layer']
! (37) [ ' ']
! (39) [ 'Setup due to waves']
! (49) [ 'Diffraction parameter']
! (50) [ 'Bottom wave period']
! (54) [ 'Bottom friction dissipation']
! (55) [ 'Surf breaking dissipation']
! (56) [ 'Whitecapping dissipation']
! (57) [ 'Vegetation dissipation']
! (58) [ 'Peakedness']
! (59) [ 'Benjamin-Feir index']
! (60) [ 'Energy generation']
! (61) [ 'Wind source term']
! (62) [ 'Energy redistribution']
! (63) [ 'Total absolute 4-wave interaction']
! (64) [ 'Total absolute 3-wave interaction']
! (65) [ 'Energy propagation']
! (66) [ 'xy-propagation']
! (67) [ 'theta-propagation']
! (68) [ 'sigma-propagation']
! (69) [ 'Radiation stress']
! (70) [ 'Plants per m2']
! (71) [ 'Peak wave length']