!----- AGPL --------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2015.
!
! This file is part of Delft3D (D-Flow Flexible Mesh component).
!
! Delft3D is free software: you can redistribute it and/or modify
! it under the terms of the GNU Affero General Public License as
! published by the Free Software Foundation version 3.
!
! Delft3D 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 Affero General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License
! along with Delft3D. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D",
! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
! Deltares, and remain the property of Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!-------------------------------------------------------------------------------------------------------
! Origin:
! URL: https://svn.oss.deltares.nl/repos/delft3d/trunk/src/engines_gpl/flow2d3d/packages/data/include/fourier.igs
! Revision: 5108
module m_fourier_analysis
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!-------------------------------------------------------------------------------
! TODO:
! * finalizing
! * writing of fourier results
! * call from unstruc --> see d3d calls as examples
! * ucx en ucy on (1:#nods) is altijd de diepte gemiddelde, ook in 2d (dan is het de enige snelheid)
use precision
use string_module, only: str_lower
implicit none
type gd_fourier
!
! integers
!
integer :: fouwrt ! Time to write fourier file (TEKAL/map-NetCDF)
integer :: nofouvar ! Number of parameters to write to NetCDF file
integer :: ibluv
integer :: iblqf
integer :: iblbs
integer :: iblep
integer :: iblwl
integer :: ibleh
integer :: iblcn
integer :: idfile ! Id of the NetCDF file
!
! pointers
!
integer , dimension(:) , pointer :: fconno ! Constituent number for Fourier analysis
integer , dimension(:) , pointer :: flayno ! Layer number for fourier analysis
integer , dimension(:) , pointer :: fnumcy ! Number of cycles for fourier analysis
integer , dimension(:) , pointer :: ftmsto ! Integer time step counter stop time for fourier analysis
integer , dimension(:) , pointer :: ftmstr ! Integer time step counter start time for fourier analysis
integer(pntrsize), dimension(:) , pointer :: ifoupt ! Pointer for FOUNAM(IFOU), points in real array R
integer , dimension(:) , pointer :: iofset ! Offset from pointer address to asked layer and constituent for IFOU
integer , dimension(:) , pointer :: foumask ! 0: no additional mask, 1: initial dry points only
integer , dimension(:) , pointer :: idvar ! Ids of the variables in the NetCDF file
integer , dimension(:,:) , pointer :: fouref ! Reference table: (ifou,1): fouid (line in input file)
! (ifou,2): fouvarstart (first index in fouvarnam/idvar to be used by this ifou
!
real(fp) , dimension(:) , pointer :: fknfac ! Fourier amplitude correction factor
real(fp) , dimension(:,:,:), pointer :: foucomp ! Component in Fourier Analysis
real(fp) , dimension(:) , pointer :: foufas ! Frequency for fourier analysis
real(fp) , dimension(:,:,:), pointer :: fousma ! Suma of fourier analysis
real(fp) , dimension(:,:,:), pointer :: fousmb ! Sumb of fourier analysis
real(fp) , dimension(:,:,:), pointer :: fouvec ! Maximum of vector magnitude for fourier analysis
! For velocity (u,v), discharge (qxk, qyk) and bed shear stress (taubpu,taubpv)
! NB: For discharges the analysis is actually performed on the unit discharges qxk/guu and qyk/gvv
real(fp) , dimension(:) , pointer :: fv0pu ! Fourier phase correction
!
character(1) , dimension(:) , pointer :: fouelp ! Y/N: Yes/No requesting elliptic parameters
! X/I: Max/Min values requested instead of fourier analysis
! E : Max Energy head requested instead of fourier analysis
character(16) , dimension(:) , pointer :: founam ! Names of variables for fourier analysis
character(50) , dimension(:) , pointer :: fouvarnam ! Names of variables for fourier analysis as written to NetCDF file
character(50) , dimension(:) , pointer :: fouvarnamlong ! Part of the long names of variables for fourier analysis as written to NetCDF file
character(50) , dimension(:) , pointer :: fouvarunit ! Unit of variables for fourier analysis as written to NetCDF file
character(1) , dimension(:) , pointer :: foutyp ! Character indicating whether parameter is a scalar (s) or vector (v) quantity
!
end type gd_fourier
!-------------------------------------------------------------------------------------------------------
type(gd_fourier), target :: gdfourier
type(gd_fourier), pointer :: gdfourier_ptr
integer, parameter :: maxMessageLen = 1000
character(maxMessageLen) :: FouMessage = ' '
private
integer :: nofou ! Number of fourier components to be analyzed
real(fp) :: ag_fouana = 9.81d0
public fouini
public fouana
public alloc_fourier_analysis_arrays
public reafou
public postpr_fourier
public gdfourier
public gdfourier_ptr
public nofou
public ag_fouana
contains
subroutine default_fouana()
implicit none
! some default setting
! ....
! Do the rest here
call reset_fouana()
end subroutine default_fouana
subroutine reset_fouana()
implicit none
! ... rest of the initialization of fourier analysis
end subroutine reset_fouana
subroutine alloc_fourier_analysis_arrays(gdfourier,gddimens,nofou)
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
!
!!--declarations----------------------------------------------------------------
use precision
use m_d3ddimens
!
implicit none
!
type(gd_fourier), intent(inout), target :: gdfourier
type(gd_dimens), intent(inout) :: gddimens
integer, intent(in) :: nofou
integer :: istat
integer :: imissval = -1
double precision :: rmissval = -999.999_fp
istat = 0
!
! Arrays for Fourier analysis (fourier.igs)
!
if (istat == 0) allocate (gdfourier%fconno (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%flayno (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fnumcy (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%ftmsto (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%ftmstr (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%ifoupt (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%iofset (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%foumask (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%idvar (1:gdfourier%nofouvar), STAT = istat)
if (istat == 0) allocate (gdfourier%fouref (1:nofou,2), STAT = istat)
!
if (istat == 0) allocate (gdfourier%fknfac ( 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%foucomp(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub, 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%foufas ( 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fousma (gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub, 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fousmb (gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub, 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fouvec (gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub, 1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fv0pu ( 1:nofou), STAT = istat)
!
if (istat == 0) allocate (gdfourier%fouelp (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%founam (1:nofou), STAT = istat)
if (istat == 0) allocate (gdfourier%fouvarnam (1:gdfourier%nofouvar), STAT = istat)
if (istat == 0) allocate (gdfourier%fouvarnamlong (1:gdfourier%nofouvar), STAT = istat)
if (istat == 0) allocate (gdfourier%fouvarunit (1:gdfourier%nofouvar), STAT = istat)
if (istat == 0) allocate (gdfourier%foutyp (1:nofou), STAT = istat)
! Initialise arrays for Fourier analysis
!
gdfourier%fconno = imissval
gdfourier%flayno = imissval
gdfourier%fnumcy = imissval
gdfourier%ftmsto = imissval
gdfourier%ftmstr = imissval
gdfourier%ifoupt = imissval
gdfourier%iofset = imissval
gdfourier%foumask = imissval
gdfourier%idvar = imissval
gdfourier%fouref = imissval
!
gdfourier%fknfac = rmissval
gdfourier%foucomp = rmissval
gdfourier%foufas = rmissval
gdfourier%fousma = rmissval
gdfourier%fousmb = rmissval
gdfourier%fouvec = rmissval
gdfourier%fv0pu = rmissval
!
gdfourier%fouelp = ' '
gdfourier%founam = ' '
gdfourier%fouvarnam = ' '
gdfourier%fouvarnamlong = ' '
gdfourier%fouvarunit = ' '
gdfourier%foutyp = ' '
!
if (istat /= 0) then
! Exception handling for allocation of fourier arrays
endif
gdfourier_ptr => gdfourier
end subroutine alloc_fourier_analysis_arrays
subroutine reafou(lunfou ,filfou ,kmax ,&
& lstsc ,lsal ,ltem ,&
& tstart ,tstop ,dt ,success)
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
!
! Function: - Read fourier input file and stores the
! variables necessary for the analysis in
! arrays.
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use mathconsts
use string_module
!
implicit none
!
! Global variables
!
integer , intent(in) :: lsal ! Description and declaration in dimens.igs
integer , intent(in) :: lstsc ! Description and declaration in dimens.igs
integer , intent(in) :: ltem ! Description and declaration in dimens.igs
integer , intent(in) :: lunfou ! Unit number fourier input file
character(*) , intent(in) :: filfou ! File name for fourier analysis input
integer , intent(in) :: kmax ! number of vertical layers
real(fp) , intent(in) :: tstart ! simulation start time
real(fp) , intent(in) :: tstop ! simulation stop time
real(fp) , intent(in) :: dt ! timestep
logical , intent(out):: success
!
! Local variables
!
integer , dimension(:) , pointer :: fconno
integer , dimension(:) , pointer :: flayno
integer , dimension(:) , pointer :: fnumcy
integer , dimension(:) , pointer :: foumask
integer , dimension(:) , pointer :: idvar
integer , dimension(:,:) , pointer :: fouref
integer , pointer :: fouwrt
integer , pointer :: nofouvar
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
real(fp) , dimension(:) , pointer :: fknfac
real(fp) , dimension(:) , pointer :: foufas
real(fp) , dimension(:) , pointer :: fv0pu
character(1) , dimension(:) , pointer :: fouelp
character(16) , dimension(:) , pointer :: founam
character(50) , dimension(:) , pointer :: fouvarnam
character(50) , dimension(:) , pointer :: fouvarnamlong
character(50) , dimension(:) , pointer :: fouvarunit
character(1) , dimension(:) , pointer :: foutyp
integer, parameter :: maxvld = 40 ! Maximum number of fields/columns that can be read from a record
integer :: fouid ! Counter linenumber-commentlines
integer :: i ! Counter
integer :: ifou ! Counter
integer :: ivar ! Counter
integer :: irelp
integer :: it
integer :: lfile ! Length of file name
integer :: linenumber ! Line number in Fourier input file
integer :: nopos ! Used for format free reading
integer :: nveld ! actual number of fields encountered in a record
integer :: iveld ! loop counter over fields in a record
real(fp) :: rstart ! Start time for fourier analysis
real(fp) :: rstop ! Stop time for fourier analysis
real(fp) :: t
character(4) :: cdummy ! Help string to read FOUELP
character(26) :: errmsg ! Character var. containing the error message to be written to file.
! The message depends on the error.
character(7) :: fmt ! Used for format free reading
character(300) :: message
character(132) :: record ! Used for format free reading
character(30) , dimension(maxvld) :: columns ! each record is split into separate fields (columns)
integer :: iostat
!
!! executable statements -------------------------------------------------------
!
fknfac => gdfourier%fknfac
foufas => gdfourier%foufas
fv0pu => gdfourier%fv0pu
fconno => gdfourier%fconno
flayno => gdfourier%flayno
fnumcy => gdfourier%fnumcy
foumask => gdfourier%foumask
idvar => gdfourier%idvar
fouref => gdfourier%fouref
fouwrt => gdfourier%fouwrt
nofouvar => gdfourier%nofouvar
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
fouelp => gdfourier%fouelp
founam => gdfourier%founam
fouvarnam => gdfourier%fouvarnam
fouvarnamlong => gdfourier%fouvarnamlong
fouvarunit => gdfourier%fouvarunit
foutyp => gdfourier%foutyp
!
success = .false.
ifou = 1
do i = 1, nofou
flayno(i) = 1
fconno(i) = 1
foumask(i) = 0
foutyp(i) = 'n'
fouref(i,:) = -1
fouelp(i) = 'n'
enddo
do i = 1, nofouvar
idvar(i) = 0
fouvarnam(i) = ' '
fouvarnamlong(i) = ' '
fouvarunit(i) = ' '
enddo
!
! define length of file name
!
call remove_leading_spaces(filfou, lfile)
!
errmsg = 'Times in file ' // filfou(1:lfile)
cdummy = ' '
!
linenumber = 0
fouid = 0
!
! reading file
!
! -->
20 continue
read (lunfou, '(a)',iostat=iostat) record
if (iostat/=0) then
goto 9999
endif
!
linenumber = linenumber + 1
fouid = fouid + 1
!
call str_lower(record, 132)
columns = ''
read(record,*,iostat=iostat) (columns(iveld),iveld=1,maxvld) ! TODO: cover for read errors
do nveld = maxvld,1,-1
if (len_trim(columns(nveld))>0) then
exit
endif
enddo
!
if (columns(1)(1:1)=='*' .or. nveld==0) then
fouid = fouid - 1
goto 20
endif
! <--
!
! determine array names and type (scalar or vectorial) for
! fourier analysis
!
founam(ifou) = trim(columns(1))
!
if (founam(ifou)=='wl') then
founam(ifou) = 's1 '
foutyp(ifou) = 's'
fouref(ifou,1) = fouid
elseif (founam(ifou)=='eh') then
!
! founam must be s1 to pass through s1 to fouana
! use fouelp to flag energy head
!
founam(ifou) = 's1 ' ! waterlevel
foutyp(ifou) = 's'
fouref(ifou,1) = fouid
fouelp(ifou) = 'e'
elseif (founam(ifou)=='uv') then ! cell-centre velocities, (u,v) (only vectorial here)
founam(ifou) = 'u1 '
founam(ifou+1) = 'v1 '
foutyp(ifou) = 'v'
foutyp(ifou+1) = 'v'
fouref(ifou,1) = fouid
fouref(ifou+1,1) = fouid
elseif (founam(ifou)=='qf') then ! interpolated cell-centre velocities (vector)
founam(ifou) = 'qxk ' ! ucx
foutyp(ifou) = 's'
fouref(ifou,1) = fouid
elseif (founam(ifou)=='bs') then
founam(ifou) = 'taubpu '
foutyp(ifou) = 's'
fouref(ifou,1) = fouid
elseif (founam(ifou)=='ct') then ! constituent: temperature (scalar)
if (ltem/=0) then
founam(ifou) = 'r1 '
foutyp(ifou) = 's'
fconno(ifou) = ltem
fouref(ifou,1) = fouid
else
! Exception: no temperature
goto 9999
endif
elseif (founam(ifou)=='cs') then ! constituent: salinity (scalar)
if (lsal/=0) then
founam(ifou) = 'r1 '
foutyp(ifou) = 's'
fconno(ifou) = lsal
fouref(ifou,1) = fouid
else
! Exception: no salt
goto 9999
endif
else ! constituent, anything else
read (founam(ifou)(2:2), '(i1)') fconno(ifou)
fconno(ifou) = fconno(ifou) + max(lsal, ltem)
if (fconno(ifou)>lstsc) then
write (cdummy(1:1), '(i1)') fconno(ifou) - max(lsal, ltem)
goto 9999
endif
founam(ifou) = 'r1 '
foutyp(ifou) = 's'
fouref(ifou,1) = fouid
endif
!
! read start time, stop time, number of cycles
! determine corresponding integer time step numbers and frequency
!
fmt = '(f )'
!
nopos = len_trim(columns(2))
!
if (nopos<10) then
write (fmt(3:3), '(i1)') nopos
write (fmt(4:5), '(a2)') '.0'
else
write (fmt(3:4), '(i2)') nopos
write (fmt(5:6), '(a2)') '.0'
endif
!
read (columns(2), fmt) rstart
!
ftmstr(ifou) = nint(rstart/dt)
!
! original code translates as : 'if abs(real(nint(rstart/dt)*dt-rstart) > (0.1_fp*dt)' ! RL666
if ( mod(real(rstart,fp),real(dt,fp)) > 0.1_fp*dt) then
! call prterr(lundia, 'U044', errmsg)
goto 9999
endif
!
!
if (rstart (0.1_fp*dt)' ! RL666
if ( mod(real(rstop,fp),real(dt,fp)) > 0.1_fp*dt) then
! call prterr(lundia, 'U044', errmsg)
goto 9999
endif
!
!
if (rstop>tstop) then
! call prterr(lundia, 'F006', ' ')
goto 9999
endif
!
! Fouwrt catches the end of all fourier analyses
!
fouwrt = max(fouwrt,(ftmsto(ifou)-1))
!
fmt = '(i )'
nopos = len_trim(columns(4))
!
if (nopos<10) then
write (fmt(3:3), '(i1)') nopos
else
write (fmt(3:4), '(i2)') nopos
endif
!
read (columns(4), fmt) fnumcy(ifou)
!
if (fnumcy(ifou)==0) then
foufas(ifou) = 0.
else
foufas(ifou) = 2.*pi*real(fnumcy(ifou),fp)/real(ftmsto(ifou) - ftmstr(ifou),fp)
endif
!
! read nodal amplifications and phase shifts for comparison
! with cotidal maps
!
fmt = '(f )'
nopos = len_trim(columns(5))
!
if (nopos<10) then
write (fmt(3:3), '(i1)') nopos
write (fmt(4:5), '(a2)') '.0'
else
write (fmt(3:4), '(i2)') nopos
write (fmt(4:5), '(a2)') '.0'
endif
!
read (columns(5), fmt) fknfac(ifou)
!
fmt = '(f )'
nopos = len_trim(columns(6))
!
if (nopos<10) then
write (fmt(3:3), '(i1)') nopos
write (fmt(4:5), '(a2)') '.0'
else
write (fmt(3:4), '(i2)') nopos
write (fmt(4:5), '(a2)') '.0'
endif
!
read (columns(6), fmt) fv0pu(ifou)
!
if (fv0pu(ifou)<0.) fv0pu(ifou) = fv0pu(ifou) + 360.
fv0pu(ifou) = mod(fv0pu(ifou), 360.0_fp)
!
irelp = 7
!
if (founam(ifou)(1:2)/='s1' .and. founam(ifou)(1:3)/='tau') then
fmt = '(i )'
nopos = len_trim(columns(7))
!
if (nopos<10) then
write (fmt(3:3), '(i1)') nopos
else
write (fmt(3:4), '(i2)') nopos
endif
!
read (columns(7), fmt) flayno(ifou)
if (flayno(ifou)>kmax) then
! call prterr(lundia, 'F007', ' ')
goto 9999
endif
irelp = irelp + 1
endif
!
! Elliptic parameters requested / MAX - MIN added
!
if (nveld>=irelp) then
cdummy = trim(columns(irelp))
!
! check for MAX and or MIN before Y/N
!
if (cdummy=='max' .or. fouelp(ifou)=='e') then
if (fouelp(ifou) == 'e') then
!
! Don't change fouelp. fouelp=e flags max energy head analysis
!
else
fouelp(ifou) = 'x'
endif
if (fnumcy(ifou)>0) then
fnumcy(ifou) = 0
foufas(ifou) = 0.
! call prterr(lundia, 'F008', 'max')
endif
elseif (cdummy=='min') then
if (fouelp(ifou) == 'e') then
write (message, '(3a,i0,a)') 'in file ', trim(filfou), ' line ', linenumber, &
& ': energy head in combination with "min" is not supported'
! call prterr(lundia, 'P004', trim(message))
goto 9999
else
fouelp(ifou) = 'i'
endif
if (fnumcy(ifou)>0) then
fnumcy(ifou) = 0
foufas(ifou) = 0.
! call prterr(lundia, 'F008', 'min')
endif
!
! elliptic parameters requested only for all foutyp='v'
!
elseif (foutyp(ifou)=='v') then
if (cdummy(1:1)=='n') then
fouelp(ifou) = 'n'
elseif (cdummy(1:1)=='y') then
fouelp(ifou) = 'y'
elseif (cdummy /= 'mean') then
write (message, '(3a,i0,2a)') 'in file ', trim(filfou), ' line ', linenumber, &
& ': expecting min, max, mean, yes or no, instead of ', trim(cdummy)
goto 9999
endif
else
if (cdummy /= 'mean') then
write (message, '(3a,i0,2a)') 'in file ', trim(filfou), ' line ', linenumber, &
& ': expecting min, max or mean, instead of ', trim(cdummy)
goto 9999
endif
endif
endif
!
if ((fouelp(ifou)=='x'.or.fouelp(ifou)=='e') .and. founam(ifou)=='s1') then
!
! Currently, inidryonly is only implemented for fouelp=x, founam=s1 or eh
!
if (index(record,'inidryonly') > 0) then
foumask(ifou) = 1
write (message,'(3a,i0,a)') 'in file ', trim(filfou), ' line ', linenumber, &
& ': Fourier analysis only for initially dry points'
! call prterr(lundia, 'G051', trim(message))
endif
endif
if (foutyp(ifou)=='v') then ! for the second vector component
ifou = ifou + 1
foutyp(ifou) = 'v'
fouref(ifou,1) = fouref(ifou - 1,1)
ftmstr(ifou) = ftmstr(ifou - 1)
ftmsto(ifou) = ftmsto(ifou - 1)
fnumcy(ifou) = fnumcy(ifou - 1)
flayno(ifou) = flayno(ifou - 1)
fconno(ifou) = fconno(ifou - 1)
foufas(ifou) = foufas(ifou - 1)
fknfac(ifou) = fknfac(ifou - 1)
fv0pu (ifou) = fv0pu(ifou - 1)
fouelp(ifou) = fouelp(ifou - 1)
foumask(ifou) = foumask(ifou - 1)
endif
!
ifou = ifou + 1
!
if (ifou<=nofou) goto 20
! <--
!
success = .True.
9999 continue
if (iostat>0) then
! ToDo: some exception handling....apparently something wring with the file
endif
!
! Define all variable names to be written
! Add the (start-)index ivar to fouref(..,2)
!
ivar = 0
do ifou = 1, nofou
fouref(ifou,2) = ivar + 1
if (foutyp(ifou) == 's') then
if (fouelp(ifou)=='x') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_max"
fouvarnamlong(ivar) = "maximum value"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
if (founam(ifou) == 's1') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_max_depth"
fouvarnamlong(ivar) = "maximum depth value"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
endif
if (foumask(ifou) == 1) then
write(fouvarnam (ivar ),'(2a)') trim(fouvarnam (ivar )), "_inidryonly"
write(fouvarnamlong(ivar ),'(2a)') trim(fouvarnamlong(ivar )), ", initially dry points only"
if (founam(ifou) == 's1') then
write(fouvarnam (ivar-1),'(2a)') trim(fouvarnam (ivar-1)), "_inidryonly"
write(fouvarnamlong(ivar-1),'(2a)') trim(fouvarnamlong(ivar-1)), ", initially dry points only"
endif
endif
elseif (fouelp(ifou)=='e') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_max"
fouvarnamlong(ivar) = "maximum value"
fouvarunit(ivar) = 'm'
if (foumask(ifou) == 1) then
write(fouvarnam (ivar ),'(2a)') trim(fouvarnam (ivar )), "_inidryonly"
write(fouvarnamlong(ivar ),'(2a)') trim(fouvarnamlong(ivar )), ", initially dry points only"
endif
elseif (fouelp(ifou)=='i') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_min"
fouvarnamlong(ivar) = "minimum value"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
else
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_amp"
fouvarnamlong(ivar) = "Fourier amplitude"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_phs"
fouvarnamlong(ivar) = "Fourier phase"
fouvarunit(ivar) = "degree"
endif
else
!
! foutyp=v
!
if (fouelp(ifou)=='x') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,2a)') "fourier", fouref(ifou,1), "_max_", trim(founam(ifou))
write(fouvarnamlong(ivar),'(2a)') "maximum value component ", trim(founam(ifou))
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
if (index(founam(ifou),'v') > 0) then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_max_mag"
fouvarnamlong(ivar) = "maximum value magnitude"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
endif
elseif (fouelp(ifou)=='i') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,2a)') "fourier", fouref(ifou,1), "_min_", trim(founam(ifou))
write(fouvarnamlong(ivar),'(2a)') "minimum value component ", trim(founam(ifou))
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
if (index(founam(ifou),'v') > 0) then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_min_mag"
fouvarnamlong(ivar) = "minimum value magnitude"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
endif
else
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,2a)') "fourier", fouref(ifou,1), "_amp_", trim(founam(ifou))
write(fouvarnamlong(ivar),'(2a)') "Fourier amplitude component ", trim(founam(ifou)) ! RL666 Wat is fouref
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,2a)') "fourier", fouref(ifou,1), "_phs_", trim(founam(ifou))
write(fouvarnamlong(ivar),'(2a)') "Fourier phase component ", trim(founam(ifou))
fouvarunit(ivar) = "degree"
endif
if (index(founam(ifou),'v')>0 .and. fouelp(ifou)=='y') then
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_ellip_amp"
fouvarnamlong(ivar) = "elliptic amplitude"
call setfouunit(founam(ifou), lsal, ltem, fconno(ifou), fouvarunit(ivar))
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_ellip_ecc"
fouvarnamlong(ivar) = "elliptic eccentricity"
fouvarunit(ivar) = ""
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_ellip_phs"
fouvarnamlong(ivar) = "elliptic phase"
fouvarunit(ivar) = "degree"
ivar = ivar + 1
write(fouvarnam(ivar),'(a,i3.3,a)') "fourier", fouref(ifou,1), "_ellip_inc"
fouvarnamlong(ivar) = "elliptic inclination"
fouvarunit(ivar) = ""
endif
endif
enddo
end subroutine reafou
subroutine setfouunit(founam, lsal, ltem, fconno, fouvarunit)
!
! parameters
character(*), intent(in) :: founam
integer , intent(in) :: lsal
integer , intent(in) :: ltem
integer , intent(in) :: fconno
character(*), intent(out) :: fouvarunit
!
! body
select case (founam(:2))
case ('s1')
fouvarunit = 'm'
case ('u1', 'v1')
fouvarunit = 'm/s'
case ('ta')
fouvarunit = 'N/m2'
case ('r1')
if (fconno == ltem) then
fouvarunit = 'degrees_Celsius'
elseif (fconno == lsal) then
fouvarunit = 'ppt'
else
fouvarunit = 'kg/m3'
endif
case default
fouvarunit = ''
end select
end subroutine setfouunit
subroutine fouana( ifou ,kfs ,kfst0 ,nst , rarray , &
& dps ,gdfourier ,gddimens ,umean, vmean)
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
!
! Function: - performs fourier analysis i.e. computes suma
! and sumb
! - calculates MAX or MIN value
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use m_d3ddimens
!
implicit none
!
type(gd_fourier) , pointer :: gdfourier
type(gd_dimens) , pointer :: gddimens
!
! The following list of pointer parameters is used to point inside the gdp structure
!
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
integer , dimension(:) , pointer :: foumask
real(fp) , dimension(:) , pointer :: foufas
real(fp) , dimension(:,:,:) , pointer :: fousma
real(fp) , dimension(:,:,:) , pointer :: fousmb
character(1) , dimension(:) , pointer :: fouelp
character(16) , dimension(:) , pointer :: founam
integer , pointer :: mmax
integer , pointer :: nmaxus
integer , pointer :: ndx
integer , pointer :: lnx
integer , pointer :: ndkx
integer , pointer :: lnkx
!
! Global variables
!
integer , intent(in) :: ifou !! Counter
integer , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in) :: kfs ! State of cell (0=dry,1=wet)
integer , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in) :: kfst0 ! State of cell at zero time (0=dry,1=wet)
integer , intent(in) :: nst !! Time step number
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in) :: rarray ! Array for fourier analysis
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in), optional :: umean ! Description and declaration in esm_alloc_real.f90
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in), optional :: vmean ! Description and declaration in esm_alloc_real.f90
real(prec), dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) , intent(in) :: dps
!
! Local variables
!
integer :: m ! Loop counter over MMAX
integer :: md
integer :: n ! Loop counter over NMAXUS
integer :: nd
integer :: nm ! Loop counter over NMMAX
real(fp) :: angl
real(fp) :: uuu ! umean in zeta point
real(fp) :: vvv ! vmean in zeta point
real(fp) :: utot2 ! |U|**2 = uuu**2 + vvv**2
!
!! executable statements -------------------------------------------------------
!
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
foumask => gdfourier%foumask
foufas => gdfourier%foufas
fousma => gdfourier%fousma
fousmb => gdfourier%fousmb
fouelp => gdfourier%fouelp
founam => gdfourier%founam
mmax => gddimens%mmax
nmaxus => gddimens%nmaxus
ndx => gddimens%ndx
lnx => gddimens%lnx
ndkx => gddimens%ndkx
lnkx => gddimens%lnkx
!
! Initialize for MAX = -1.0e+30 / MIN = 1.0e+30
!
if (nst==ftmstr(ifou)) then
if (fouelp(ifou)=='x' .or. fouelp(ifou)=='e') then
fousma(:, :, ifou) = -1.0e+30_fp
fousmb(:, :, ifou) = -1.0e+30_fp
elseif (fouelp(ifou)=='i') then
fousma(:, :, ifou) = 1.0e+30_fp
else
fousma(:, :, ifou) = 0.0_fp
fousmb(:, :, ifou) = 0.0_fp
endif
endif
!
! Perform fourier analysis, every timestep as long as NST value
! lies in requested time interval FTMSTR and FTMSTO
!
! The name of the variable for Fourier analysis fully specified the number of elements for the fourier loop
select case (founam(ifou))
case('s1')
nmaxus = ndx
case('u1')
nmaxus = ndkx
case('qxk')
nmaxus = lnx
case('taubpu')
nmaxus = ndx
case('r1')
nmaxus = ndkx
end select
if (nst>=ftmstr(ifou) .and. nst gdfourier%nofouvar
! initialisation
!
nofou = 0
nofouvar = 0
!
! cycle through file, while reading records
!
! -->
!
10 continue
read (lunfou, '(a)', end = 999) record
!
! reset record in smaller case characters and define contents
!
call str_lower(record, 300)
columns = ''
read(record,*,iostat=iostat) (columns(iveld),iveld=1,maxvld) ! TODO: cover for read errors
!
do nveld = maxvld,1,-1
if (len_trim(columns(nveld))>0) exit
enddo
if (columns(1)(1:1)=='*' .or. nveld==0) goto 10
!
! test for continuation record
!
! requested fourier analysis water-level
!
if (columns(1)(1:2)=='wl') then
nofou = nofou + 1
if ((index(record,'max')>0) .or. (index(record,'min')>0)) then
nofouvar = nofouvar + 1
else
!
! max: also write max depth
!
nofouvar = nofouvar + 2
endif
!
! requested fourier analysis energy head
!
elseif (columns(1)(1:2)=='eh') then
nofou = nofou + 1
nofouvar = nofouvar + 1
!
! requested fourier analysis cell-centre velocities ucx and ucy on cell circumcentres
!
elseif (columns(1)(1:2)=='uv') then
nofou = nofou + 2 !
if (index(record,'max')>0 .or. index(record,'min')>0) then ! (either min or max) for (both u and v)
nofouvar = nofouvar + 2
elseif (index(record,'y')>0) then
nofouvar = nofouvar + 8 ! both amplitude and phase) for (both u and v) + 4 elliptic parameters
else
nofouvar = nofouvar + 4 ! both amplitude and phase) for (both u and v)
endif
!
! requested fourier analysis discharge
!
elseif (columns(1)(1:2)=='qf') then
nofou = nofou + 1
if (index(record,'max')>0 .or. index(record,'min')>0) then
nofouvar = nofouvar + 3
else
nofouvar = nofouvar + 4
endif
!
! requested fourier analysis bedstress
!
elseif (columns(1)(1:2)=='bs') then ! RL666: bedstresses are scalar
nofou = nofou + 1
if (index(record,'max')>0 .or. index(record,'min')>0) then
nofouvar = nofouvar + 3
else
nofouvar = nofouvar + 4
endif
!
! requested fourier analysis temperature
!
elseif (columns(1)(1:2)=='ct') then
nofou = nofou + 1
if (index(record,'max')>0 .or. index(record,'min')>0) then
nofouvar = nofouvar + 1
else
nofouvar = nofouvar + 2
endif
!
! requested fourier analysis salinity
!
elseif (columns(1)(1:2)=='cs') then
nofou = nofou + 1
if (index(record,'max')>0 .or. index(record,'min')>0) then
nofouvar = nofouvar + 1
else
nofouvar = nofouvar + 2
endif
!
! requested fourier analysis constituent
!
elseif ((columns(1)(1:1)=='c') .and. (index('12345',columns(1)(2:2))>0)) then
nofou = nofou + 1
if (index(record,'max')>0 .or. index(record,'min')>0) then
nofouvar = nofouvar + 1
else
nofouvar = nofouvar + 2
endif
else
!
! requested fourier analysis undefined
!
goto 999
endif
!
goto 10
!
! <-- next record
!
!
success = .True.
999 continue
rewind(lunfou)
end subroutine fouini
subroutine postpr_fourier(s1,u1,ucx,ucy,constituents,taus,kfs,kfst0,dps,nst , &
& trifil ,dtsec ,versio ,namcon , &
& refdat ,hdt ,tzone ,xs, ys, xu, yu, gdfourier ,gddimens ,umean ,vmean)
use m_d3ddimens
implicit none
type(gd_fourier) , pointer :: gdfourier !< Fourier Analysis structure
type(gd_dimens) , pointer :: gddimens !< Model geometry/grid structure
double precision , intent(in), pointer :: s1(:,:) !< waterlevel
double precision , intent(in), pointer :: u1(:,:) !< face-normal velocity
double precision , intent(in), pointer :: ucx(:,:) !< cell-centre velocity x-direction
double precision , intent(in), pointer :: ucy(:,:) !< cell-centre velocity y-direction
!double precision , intent(in), pointer :: xs(:,:) !< s-point x-coordinate
!double precision , intent(in), pointer :: ys(:,:) !< s-point y-coordinate
!double precision , intent(in), pointer :: xu(:,:) !< u-point x-coordinate
!double precision , intent(in), pointer :: yu(:,:) !< u-point y-coordinate
double precision , intent(in) :: xs(:) !< s-point x-coordinate
double precision , intent(in) :: ys(:) !< s-point y-coordinate
double precision , intent(in) :: xu(:) !< u-point x-coordinate
double precision , intent(in) :: yu(:) !< u-point y-coordinate
real(fp) , intent(in) :: dtsec !< Integration time step [in seconds]
real(fp) , intent(in) :: hdt !< Half Integration time step [seconds] => gdp%gdnumeco%hdt
real(fp) , intent(in) :: tzone !< Local (FLOW) time - GMT (in hours) => gdp%gdexttim%tzone
double precision , intent(in), pointer, optional :: umean(:,:) !< depth-average cell-centre velocity x-direction
double precision , intent(in), pointer, optional :: vmean(:,:) !< depth-average cell-centre velocity y-direction
double precision , intent(in) :: dps(:,:) !< depth (inverted bottomlevel)
integer , intent(in) :: kfs(:,:) !< if cell is wet: 1, otherwise: 0
integer , intent(in) :: kfst0(:,:) !< if cell is wet at t=0: 1, otherwise: 0
integer , intent(in) :: nst !< timestep number
character(len=*) , intent(in) :: trifil !< output filename
character(len=*) , intent(in) :: versio ! Version nr. of the current package
double precision , intent(in), pointer :: constituents(:,:,:) !< array of constituents, also containing sal and tmp
character(len=*) , intent(in), dimension(:) :: namcon !< constituent names
double precision , intent(in), pointer :: taus(:,:) !< bedfriction velocity (?)
character(len=*), intent(in) :: refdat
!
! Perform analysis and write to Fourier file
!
integer :: ifou
integer , dimension(:) , pointer :: fconno
integer , dimension(:) , pointer :: flayno
integer , dimension(:) , pointer :: fnumcy
integer , pointer :: fouwrt
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
integer(pntrsize), dimension(:) , pointer :: ifoupt
integer , dimension(:) , pointer :: iofset
real(fp) , dimension(:) , pointer :: fknfac
real(fp) , dimension(:,:,:) , pointer :: foucomp
real(fp) , dimension(:) , pointer :: foufas
real(fp) , dimension(:,:,:) , pointer :: fousma
real(fp) , dimension(:,:,:) , pointer :: fousmb
real(fp) , dimension(:,:,:) , pointer :: fouvec
real(fp) , dimension(:) , pointer :: fv0pu
character(1) , dimension(:) , pointer :: fouelp
character(16) , dimension(:) , pointer :: founam
character(1) , dimension(:) , pointer :: foutyp
integer , pointer :: nmax
integer , pointer :: mmax
integer , pointer :: nlb
integer , pointer :: nub
integer , pointer :: mlb
integer , pointer :: mub
integer , pointer :: nmaxus
integer , pointer :: kmax
double precision, pointer :: fieldptr1(:,:), fieldptr2(:,:)
integer :: itdate !< Reference time in YYYYMMDD as an integer
read(refdat,*) itdate
fconno => gdfourier%fconno
flayno => gdfourier%flayno
fnumcy => gdfourier%fnumcy
fouwrt => gdfourier%fouwrt
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
ifoupt => gdfourier%ifoupt
iofset => gdfourier%iofset
fknfac => gdfourier%fknfac
foucomp => gdfourier%foucomp
foufas => gdfourier%foufas
fousma => gdfourier%fousma
fousmb => gdfourier%fousmb
fouvec => gdfourier%fouvec
fv0pu => gdfourier%fv0pu
fouelp => gdfourier%fouelp
founam => gdfourier%founam
foutyp => gdfourier%foutyp
nmax => gddimens%nmax
mmax => gddimens%mmax
nlb => gddimens%nlb
nub => gddimens%nub
mlb => gddimens%mlb
mub => gddimens%mub
nmaxus => gddimens%nmaxus
kmax => gddimens%kmax
if (nofou > 0) then
ifou = 1
do
if (ifou > nofou) exit
!
! Perform fourier analysis, every timestep as long as NST value
! lies in requested time interval FTMSTR and FTMSTO
!
if (foutyp(ifou) == 's') then ! scalar
!
! Scalar type Fourier component
!
! Get Fourier component pointer
!
select case (founam(ifou)(:2))
case ('s1') ! waterlevels
fieldptr1 => s1
case ('r1')
fieldptr1 => constituents(fconno(ifou),:,:)
case ('taubpu')
fieldptr1 => taus
case default
continue ! Unknown FourierVariable exception
end select
if (present(umean) .and. present(vmean)) then
call fouana(ifou ,kfs ,kfst0 ,nst ,fieldptr1 , & ! future use of analysis with depth average cell-centre velocities
& dps ,gdfourier ,gddimens ,umean ,vmean)
else
call fouana(ifou ,kfs ,kfst0 ,nst ,fieldptr1 , &
& dps ,gdfourier ,gddimens)
endif
ifou = ifou + 1
elseif (foutyp(ifou) == 'v') then ! vector with 2 component
!
! Vector type Fourier component: do everything for both vector components
!
select case (founam(ifou)(:2))
case ('u1') ! cell-centred velocities u and v
fieldptr1 => ucx
fieldptr2 => ucy
end select
!
if (present(umean) .and. present(vmean)) then
call fouana(ifou ,kfs ,kfst0 ,nst ,fieldptr1 , &
& dps ,gdfourier ,gddimens ,umean ,vmean)
call fouana(ifou+1 ,kfs ,kfst0 ,nst ,fieldptr2 , &
& dps ,gdfourier ,gddimens ,umean ,vmean)
else
call fouana(ifou ,kfs ,kfst0 ,nst ,fieldptr1 , &
& dps ,gdfourier ,gddimens )
call fouana(ifou+1 ,kfs ,kfst0 ,nst ,fieldptr2 , &
& dps ,gdfourier ,gddimens )
endif
!
! Determine the maximum or minimum of the magnitude of the vector
! Done for velocity (u,v), discharge (qxk,qyk) and bed shear stress (taubpu,taubpv)
!
call fouvecmax(mmax ,nmaxus ,nofou ,ifou , &
& nst ,gdfourier ,gddimens )
!
! Skip second component of vector (already done in fou_vector)
!
ifou = ifou + 2
else
!
! Incorrect Fourier type found, issue a warning
!
endif
enddo
!
! Write results of fourier analysis to data file
! only once when all fourier periods are complete
!
if (nst==fouwrt) then
call wrfou(trifil ,dtsec ,versio ,namcon , &
& itdate ,hdt ,tzone ,xs ,ys ,xu ,yu ,gdfourier ,gddimens)
endif
endif
end subroutine postpr_fourier
subroutine wrfou(trifil ,dtsec ,versio ,namcon , &
& itdate ,hdt ,tzone ,xs, ys, xu, yu, gdfourier ,gddimens)
!----- GPL ---------------------------------------------------------------------
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
!
! Function: - open fourier analysis output file
! - writes results of fourier analysis to output
! file
! - closes fourier analysis output file
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use mathconsts
use netcdf
use m_d3ddimens
use unstruc_netcdf, only: unc_write_flowgeom_filepointer
!
implicit none
!
type(gd_fourier) , pointer :: gdfourier !< Fourier Analysis structure
type(gd_dimens) , pointer :: gddimens !< Model geometry/grid structure
!
integer , pointer :: nofouvar
integer , dimension(:) , pointer :: fconno
character(1) , dimension(:) , pointer :: foutyp
character(1) , dimension(:) , pointer :: fouelp
character(16) , dimension(:) , pointer :: founam
character(50) , dimension(:) , pointer :: fouvarnam
character(50) , dimension(:) , pointer :: fouvarnamlong
character(50) , dimension(:) , pointer :: fouvarunit
integer , dimension(:,:), pointer :: fouref
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
real(fp) , dimension(:) , pointer :: foufas
integer , dimension(:) , pointer :: flayno
integer , dimension(:) , pointer :: fnumcy
integer , pointer :: iblwl
integer , pointer :: ibleh
integer , pointer :: iblcn
integer , pointer :: ibluv
integer , pointer :: iblqf
integer , pointer :: iblbs
integer , pointer :: iblep
integer , pointer :: idfile
integer , dimension(:) , pointer :: idvar
integer , pointer :: nmaxgl
integer , pointer :: mmaxgl
integer , pointer :: mmax
integer , pointer :: nmax
integer , pointer :: ndx
integer , pointer :: lnx
integer , pointer :: ndkx
integer , pointer :: lnkx
integer , pointer :: nmaxus
!
! Global variables
!
integer , intent(in) :: itdate !< Reference time in YYYYMMDD
real(fp) , intent(in) :: dtsec !< Integration time step [in seconds]
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) :: xcor !< Description and declaration in esm_alloc_real.f90
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) :: xz !< Description and declaration in esm_alloc_real.f90
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) :: ycor !< Description and declaration in esm_alloc_real.f90
real(fp) , dimension(gddimens%nlb:gddimens%nub, gddimens%mlb:gddimens%mub) :: yz !< Description and declaration in esm_alloc_real.f90
character(20) , dimension(:) :: namcon !< Description and declaration in esm_alloc_char.f90
character(5) , intent(in) :: versio !< Version nr. of the current package
character(*) , intent(in) :: trifil !< File name for FLOW NEFIS output files (tri"h/m"-"casl""labl".dat/def)
real(fp) , intent(in) :: tzone !< Local (FLOW) time - GMT (in hours) => gdp%gdexttim%tzone
real(fp) , intent(in) :: hdt !< Half Integration time step [seconds] => gdp%gdnumeco%hdt
!double precision , intent(in), pointer :: xs(:,:) !< s-point x-coordinate
!double precision , intent(in), pointer :: ys(:,:) !< s-point y-coordinate
!double precision , intent(in), pointer :: xu(:,:) !< u-point x-coordinate
!double precision , intent(in), pointer :: yu(:,:) !< u-point y-coordinate
double precision , intent(in) :: xs(:) !< s-point x-coordinate
double precision , intent(in) :: ys(:) !< s-point y-coordinate
double precision , intent(in) :: xu(:) !< u-point x-coordinate
double precision , intent(in) :: yu(:) !< u-point y-coordinate
!
! Local variables
!
integer :: filetype
integer :: iddim_nflowelem
integer :: iddim_nflowlink
integer :: iddim
integer :: ierr
!
integer, dimension(6) :: idatt_fou
!
integer :: irequest
integer :: ierror
integer :: ifou ! Local teller for fourier functions
integer :: ivar ! Local teller for fourier functions
integer :: lrid ! Length of RUNID character string
integer :: lunfou
integer , external :: newlun
integer , external :: nc_def_var
real(fp) :: freqnt ! Frequency in degrees per hour
real(fp) :: tfasto ! Stop time in minutes
real(fp) :: tfastr ! Start time in minutes
character(4) :: blnm
character(16) :: fougrp
character(20) :: namfun ! Local name for fourier function
character(30) :: namfunlong ! Local name for fourier function, including reference to the line in the fourier input file
character(256) :: filename
character(256) :: msg
integer :: idvar_xs, idvar_ys, idvar_xu, idvar_yu
integer, parameter :: REQUESTTYPE_DEFINE = 1
integer, parameter :: REQUESTTYPE_WRITE = 2
!
!! executable statements -------------------------------------------------------
!
nofouvar => gdfourier%nofouvar
fconno => gdfourier%fconno
foutyp => gdfourier%foutyp
fouelp => gdfourier%fouelp
founam => gdfourier%founam
fouvarnam => gdfourier%fouvarnam
fouvarnamlong => gdfourier%fouvarnamlong
fouvarunit => gdfourier%fouvarunit
fouref => gdfourier%fouref
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
foufas => gdfourier%foufas
flayno => gdfourier%flayno
fnumcy => gdfourier%fnumcy
iblwl => gdfourier%iblwl
ibleh => gdfourier%ibleh
iblcn => gdfourier%iblcn
ibluv => gdfourier%ibluv
iblqf => gdfourier%iblqf
iblbs => gdfourier%iblbs
iblep => gdfourier%iblep
idfile => gdfourier%idfile
idvar => gdfourier%idvar
mmax => gddimens%mmax
nmax => gddimens%nmax
lnx => gddimens%lnx
ndx => gddimens%ndx
lnkx => gddimens%lnkx
ndkx => gddimens%ndkx
nmaxus => gddimens%nmaxus
! lundia => gdp%gdinout%lundia
mmaxgl => mmax ! => gdp%gdparall%mmaxgl ! RL TODO: Hier een aparte parameter voor nodig ?
nmaxgl => nmax ! => gdp%gdparall%nmaxgl
!
fougrp = 'fou-fields'
!
if (idfile == 0) then
!filename = trifil(1:3) // 'f' // trifil(5:)
!write(filename,'(2a)') trim(filename), '.nc'
!!
!ierror = nf90_create(filename, NF90_WRITE, idfile)
ierror = nf90_create(trim(trifil), NF90_WRITE, idfile)
FouMessage = ''
if (ierror/=0) FouMessage = nf90_strerror(ierror)
! DEFINITIONS
!
! Assumption: wrimap has already created the map-NetCDF file
!
if (ierror/=0) FouMessage = nf90_strerror(ierror)
!
! Dimensions
!
ierror = nf90_def_dim(idfile, 'nFlowElemWithBnd', ndx, iddim_nflowelem) ! Flow elements (cell centres)
ierror = nf90_def_dim(idfile, 'nFlowLink', lnx, iddim_nflowlink) ! Flow links
! Define X and Y of waterlevel points :
ierr = nf90_def_var(idfile, 'FlowElem_xcc', nf90_double, (/iddim_nflowelem/), idvar_xs)
ierr = nf90_put_att(idfile, idvar_xs, 'standard_name', 'projection_x_coordinate')
! ierr = nf90_put_att(idfile, idvar_xs, 'long_name','X-coordinate s-point')
ierr = nf90_put_att(idfile, idvar_xs, 'long_name','x')
ierr = nf90_put_att(idfile, idvar_xs, 'units','m')
ierr = nf90_def_var(idfile, 'FlowElem_ycc', nf90_double, (/iddim_nflowelem/), idvar_ys)
ierr = nf90_put_att(idfile, idvar_ys, 'standard_name', 'projection_y_coordinate')
! ierr = nf90_put_att(idfile, idvar_ys, 'long_name','Y-coordinate s-point')
ierr = nf90_put_att(idfile, idvar_ys, 'long_name','y')
ierr = nf90_put_att(idfile, idvar_ys, 'units','m')
ierr = nf90_def_var(idfile, 'FlowLink_xu', nf90_double, (/ iddim_nflowlink /) , idvar_xu)
ierr = nf90_put_att(idfile, idvar_xu, 'standard_name', 'projection_x_coordinate')
! ierr = nf90_put_att(idfile, idvar_xu, 'long_name' , 'Center coordinate of net link (velocity point).')
ierr = nf90_put_att(idfile, idvar_xu, 'long_name','x')
ierr = nf90_put_att(idfile, idvar_xu, 'units','m')
ierr = nf90_def_var(idfile, 'FlowLink_yu', nf90_double, (/ iddim_nflowlink /) , idvar_yu)
ierr = nf90_put_att(idfile, idvar_yu, 'standard_name', 'projection_y_coordinate')
! ierr = nf90_put_att(idfile, idvar_yu, 'long_name' , 'Center coordinate of net link (velocity point).')
ierr = nf90_put_att(idfile, idvar_yu, 'long_name','y')
ierr = nf90_put_att(idfile, idvar_yu, 'units','m')
!
ifou = 1
iblwl = 0
ibleh = 0
iblcn = 0
ibluv = 0
iblqf = 0
iblbs = 0
iblep = 0
do ivar=1, nofouvar
if (ifou < nofou) then
if (fouref(ifou+1,2) <= ivar) then
ifou = ifou + 1
endif
endif
freqnt = foufas(ifou)*raddeg*3600.0_fp/dtsec
tfastr = real(ftmstr(ifou),fp)*dtsec/60.0_fp
tfasto = real(ftmsto(ifou),fp)*dtsec/60.0_fp
!
if (founam(ifou)(:2)=='s1') then
if (fouelp(ifou)=='e') then
ibleh = ibleh + 1
blnm = 'EH??'
write (blnm(3:4), '(i2.2)') ibleh
namfun = 'energy head'
else
iblwl = iblwl + 1
blnm = 'WL??'
write (blnm(3:4), '(i2.2)') iblwl
namfun = 'water level'
endif
endif
if (founam(ifou)(:2)=='r1') then
iblcn = iblcn + 1
blnm = 'CO??'
write (blnm(3:4), '(i2.2)') iblcn
namfun = namcon(fconno(ifou))
endif
if (founam(ifou)(:2)=='u1') then
ibluv = ibluv + 1
blnm = 'UV??'
write (blnm(3:4), '(i2.2)') ibluv
namfun = 'velocity'
endif
if (founam(ifou)(:2)=='qx') then
iblqf = iblqf + 1
blnm = 'QF??'
write (blnm(3:4), '(i2.2)') iblqf
namfun = 'unit discharge'
endif
if (founam(ifou)(:2)=='ta') then
iblbs = iblbs + 1
blnm = 'BS??'
write (blnm(3:4), '(i2.2)') iblbs
namfun = 'bed stress'
endif
write(namfunlong,'(i3.3,2a)') fouref(ifou,1), ": ", trim(namfun)
!
! Different variables on different dimensions
select case (founam(ifou))
case('s1','r1','u1')
iddim = iddim_nflowelem
case('qxk','taubpu')
iddim = iddim_nflowlink
end select
ierr = nf90_def_var(idfile,trim(fouvarnam(ivar)) , NF90_FLOAT, (/iddim/), idvar(ivar)) ! RL666: split flowlink- and flowelemen
ierr = nf90_put_att(idfile,idvar(ivar), 'long_name','Fourier analysis '//trim(namfunlong)//', '//trim(fouvarnamlong(ivar)))
ierr = nf90_put_att(idfile,idvar(ivar), 'units',fouvarunit(ivar))
ierr = nf90_put_att(idfile,idvar(ivar), 'layer_number', flayno(ifou))
ierr = nf90_put_att(idfile,idvar(ivar), 'Reference_date_in_YYYYMMDD', itdate)
ierr = nf90_put_att(idfile,idvar(ivar), 'Starttime_fourier_analysis_in_minutes_since_reference_date', tfastr)
ierr = nf90_put_att(idfile,idvar(ivar), 'Stoptime_fourier_analysis_in_minutes_since_reference_date', tfasto)
select case (founam(ifou))
case('s1','r1','u1')
ierr = nf90_put_att(idfile, idvar(ivar), 'coordinates' , 'FlowElem_xcc FlowElem_ycc')
case('qxk','taubpu')
ierr = nf90_put_att(idfile, idvar(ivar), 'coordinates' , 'FlowElem_xu FlowElem_yu')
end select
ierr = nf90_put_att(idfile,idvar(ivar), 'Number_of_cycles', fnumcy(ifou))
ierr = nf90_put_att(idfile,idvar(ivar), 'Frequency_degrees_per_hour', freqnt)
!
enddo
!
! Reset the global indices
!
iblwl = 0
ibleh = 0
iblcn = 0
ibluv = 0
iblqf = 0
iblbs = 0
iblep = 0
endif
! END DEFINITION MODE
ierr = nf90_enddef(idfile)
!ierror = nf90_put_var(idfile, idvar_xs, xs(1:ndx,1), start=(/ 1/), count = (/ndx/))
!ierror = nf90_put_var(idfile, idvar_ys, ys(1:ndx,1), start=(/ 1/), count = (/ndx/))
!ierror = nf90_put_var(idfile, idvar_xu, xu(1:lnx,1), start=(/ 1/), count = (/lnx/))
!ierror = nf90_put_var(idfile, idvar_yu, yu(1:lnx,1), start=(/ 1/), count = (/lnx/))
ierror = nf90_put_var(idfile, idvar_xs, xs(1:ndx), start=(/ 1/), count = (/ndx/))
ierror = nf90_put_var(idfile, idvar_ys, ys(1:ndx), start=(/ 1/), count = (/ndx/))
ierror = nf90_put_var(idfile, idvar_xu, xu(1:lnx), start=(/ 1/), count = (/lnx/))
ierror = nf90_put_var(idfile, idvar_yu, yu(1:lnx), start=(/ 1/), count = (/lnx/))
! START DATA MODE
!
ifou = 1
!
! Write requested fourier function output for scalar quantity
!
do
if (ifou > nofou) exit
!
if (foutyp(ifou)=='s') then
call wrfous(ifou ,dtsec ,namcon ,hdt ,tzone ,gdfourier ,gddimens )
ifou = ifou + 1
else
! call wrfous(ifou ,dtsec ,namcon ,hdt ,tzone ,gdfourier ,gddimens )
! call wrfous(ifou+1 ,dtsec ,namcon ,hdt ,tzone ,gdfourier ,gddimens )
call wrfouv(ifou ,dtsec ,hdt ,tzone ,gdfourier ,gddimens )
! Write vector as two scalars and see what happens ....
! Todo: Add the additional vector stuff, such as elliptic results etc.
!
ifou = ifou + 2
endif
enddo
!endif
!
! Close fourier output file
!
ierror = nf90_close(idfile)
end subroutine wrfou
subroutine fouvecmax(mmax ,nmaxus ,nofou , &
& ifou ,nst ,gdfourier ,gddimens )
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
!
! Function: - Determines the maximum of the different vector parameters
! 1) Velocity (u and v) [m/s]
! 2) (Unit) Discharge (qxk and qyk) [m3/m]
! 3) Bed shear stress (taubpu and taubpv) [N/m2]
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use m_d3ddimens
!
implicit none
!
type(gd_fourier),target, intent(in) :: gdfourier
type(gd_dimens), target, intent(in) :: gddimens
!
! The following list of pointer parameters is used to point inside the gdp structure
!
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
real(fp) , dimension(:,:,:) , pointer :: foucomp
real(fp) , dimension(:,:,:) , pointer :: fouvec
character(1) , dimension(:) , pointer :: fouelp
!
! Global variables
!
integer , intent(in) :: ifou !! Counter
integer , intent(in) :: mmax ! Description and declaration in esm_alloc_int.f90
integer , intent(in) :: nmaxus ! Description and declaration in esm_alloc_int.f90
integer , intent(in) :: nofou ! Description and declaration in dimens.igs
integer , intent(in) :: nst !! Time step number
!
! Local variables
!
integer :: m ! Loop counter over MMAX
integer :: n ! Loop counter over NMAXUS
real(fp) :: vecmagn ! Magnitude of vector ( =sqrt(xcomp^2 + ycomp^2) )
!
!! executable statements -------------------------------------------------------
!
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
foucomp => gdfourier%foucomp
fouvec => gdfourier%fouvec
fouelp => gdfourier%fouelp
!
vecmagn = 0.0_fp
!
! Initialize for MAX = -1.E+30 / MIN = 1.E+30
!
if (nst==ftmstr(ifou)) then
if (fouelp(ifou)=='x') then
do n = 1, nmaxus
do m = 1, mmax
fouvec(n, m, ifou) = -1.0E+30
enddo
enddo
elseif (fouelp(ifou)=='i') then
do n = 1, nmaxus
do m = 1, mmax
fouvec(n, m, ifou) = 1.0E+30
enddo
enddo
else
endif
endif
!
! For every time step between FTMSTART and FTMSTOP
!
if (nst>=ftmstr(ifou) .and. nst gdp%gdexttim%tzone
real(fp) , intent(in) :: hdt !< Half Integration time step [seconds] => gdp%gdnumeco%hdt
type(gd_fourier) , pointer :: gdfourier
type(gd_dimens) , pointer :: gddimens
!
! Local variables
!
integer :: ierror
integer :: fouvar
integer :: m ! Loop counter over MMAX
integer :: n ! Loop counter over NMAXUS
integer :: ncol ! Number of column to write to TEKAL data file
logical :: ltest ! Help variable for atan2 function test
real(fp) :: amp ! Fourier amplitude
real(fp) :: fas ! Fourier phase
real(fp) :: freqnt ! Frequency in degrees per hour
real(fp) :: shift ! Phase shift
real(fp) :: tfasto ! Stop time in minutes
real(fp) :: tfastr ! Start time in minutes
real(sp) :: defaul ! Default value
character(20) :: namfun ! Local name for fourier function
character(4) :: blnm
real(sp), dimension(:,:), allocatable, save :: glbarr2 ! global arrays in dffunctionals.f90, used in d3d, local here !!
real(sp), dimension(:,:,:), allocatable, save :: glbarr3
real(sp), dimension(:,:,:,:), allocatable, save :: glbarr4
integer, dimension(:,:), allocatable, save :: glbari2
!
!! executable statements -------------------------------------------------------
!
fconno => gdfourier%fconno
flayno => gdfourier%flayno
fnumcy => gdfourier%fnumcy
iblwl => gdfourier%iblwl
ibleh => gdfourier%ibleh
iblcn => gdfourier%iblcn
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
idvar => gdfourier%idvar
fknfac => gdfourier%fknfac
foufas => gdfourier%foufas
fousma => gdfourier%fousma
fousmb => gdfourier%fousmb
fv0pu => gdfourier%fv0pu
fouelp => gdfourier%fouelp
founam => gdfourier%founam
fouvarnam => gdfourier%fouvarnam
nofouvar => gdfourier%nofouvar
fouref => gdfourier%fouref
idfile => gdfourier%idfile
mmax => gddimens%mmax
nmax => gddimens%nmax
nmaxus => gddimens%nmaxus
ndx => gddimens%ndx
lnx => gddimens%lnx
ndkx => gddimens%ndkx
lnkx => gddimens%lnkx
mmaxgl => mmax ! => gdp%gdparall%mmaxgl ! RL TODO: Hier een aparte parameter voor nodig ?
nmaxgl => nmax ! => gdp%gdparall%nmaxgl
!
! Initialize local variables
!
defaul = -999.0_sp
!
! Frequention := 360 degree / period
! where period = [ (FTMSTO - FMTSTR) * DTSEC ] / [ FNUMCY * 3600 ]
! FOUFAS is defined in RDFOUR as
! FOUFAS = 2 * PI * FNUMCY / [(FTMSTO - FMTSTR) ]
! so FREQNT = FOUFAS * RADDEG * 3600 / DTSEC is OK
!
shift = ftmstr(ifou)*foufas(ifou)
freqnt = foufas(ifou)*raddeg*3600.0_fp/dtsec
tfastr = real(ftmstr(ifou),fp)*dtsec/60.0_fp
tfasto = real(ftmsto(ifou),fp)*dtsec/60.0_fp
!
namfun = founam(ifou)
select case (founam(ifou))
case('s1')
nmaxus = ndx
case('u1','v1')
nmaxus = ndkx
case('qxk')
nmaxus = lnx
case('taubpu')
nmaxus = lnx
case('r1')
nmaxus = ndkx
end select
if (founam(ifou)(:2)=='s1') then
if (fouelp(ifou)=='e') then
ibleh = ibleh + 1
blnm = 'EH??'
write (blnm(3:4), '(i2.2)') ibleh
namfun = 'energy head'
else
iblwl = iblwl + 1
blnm = 'WL??'
write (blnm(3:4), '(i2.2)') iblwl
namfun = 'water level'
endif
endif
if (founam(ifou)(:2)=='r1') then
iblcn = iblcn + 1
blnm = 'CO??'
write (blnm(3:4), '(i2.2)') iblcn
namfun = namcon(fconno(ifou))
endif
!
! Write data for user defined dimensions, hence NMAXUS and MMAX
! First for Maximum or Minimum
!
if (fouelp(ifou)=='x' .or. fouelp(ifou)=='i' .or. fouelp(ifou)=='e') then
if (allocated(glbarr2)) deallocate(glbarr2, stat = ierror)
allocate(glbarr2(nmaxgl,mmaxgl), stat = ierror)
glbarr2 = defaul
do n = 1, nmaxus
do m = 1, mmax
!
! Only write values unequal to initial min/max values (-/+1.0e+30)
!
if (comparereal(abs(fousma(n,m,ifou)),1.0e29_fp)==-1) then
glbarr2(n,m) = real(fousma(n,m,ifou),sp)
endif
enddo
enddo
fouvar = fouref(ifou,2)
ierror = nf90_put_var(idfile, idvar(fouvar),glbarr2(1:nmaxus,1),start=(/ 1/), count = (/nmaxus/))
if (fouelp(ifou)=='x' .and. founam(ifou)=='s1') then
!
! Write max waterdepth too
!
glbarr2 = defaul
do n = 1, nmaxus
do m = 1, mmax
!
! Only write values unequal to initial min/max values (-/+1.0e+30)
!
if (comparereal(abs(fousmb(n,m,ifou)),1.0e29_fp)==-1) then
glbarr2(n,m) = real(fousmb(n,m,ifou),sp)
endif
enddo
enddo
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar),glbarr2(1:nmaxus,1),start=(/ 1/), count = (/nmaxus/)) ! RL: split in flowlink and flownode quantities
endif
else
if (allocated(glbarr3)) deallocate(glbarr3, stat = ierror)
allocate(glbarr3(nmaxgl,mmaxgl,2), stat = ierror)
glbarr3 = defaul
!
! Write data for user defined dimensions, hence NMAXUS and MMAX
!
fouvar = fouref(ifou,2)
do n = 1, nmaxus
do m = 1, mmax
ltest = (fousma(n, m, ifou)==0.0_fp .and. fousmb(n, m, ifou)==0.0_fp)
!
! Test for active point and non-zero values
! when KCS (N,M) = 1 N > 1 and M > 1 per definition
!
if (.not.ltest) then
fousma(n, m, ifou) = fousma(n, m, ifou) &
& *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
fousmb(n, m, ifou) = fousmb(n, m, ifou) &
& *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
amp = sqrt(fousma(n, m, ifou)*fousma(n, m, ifou) &
& + fousmb(n, m, ifou)*fousmb(n, m, ifou))
fas = atan2(fousmb(n, m, ifou), fousma(n, m, ifou)) + shift
if (fnumcy(ifou)==0) then
amp = 0.5_fp*amp*cos(fas)
fas = 0.0_fp
endif
!
! Timezone correction added timezone*phase [degrees/hr].
! foufas is in [rad/timestep]
! halftimestep is in [sec/timestep]
! => timezonecorr = tzone [-] * foufas [rad/timestep] * raddeg [deg/rad] * [sec/hr] / (2 * halftimestep [sec/timestep])
!
fas = fas*raddeg + fv0pu(ifou) - tzone*foufas(ifou)*raddeg*1800.0_fp/hdt
!
! To define FAS between 0 and 360. add 720. to the MOD of
! FAS and re-use the MOD function
!
fas = mod(mod(fas, 360.0_fp) + 720.0_fp, 360.0_fp)
amp = amp/fknfac(ifou)
glbarr3(n,m,1) = real(amp,sp)
glbarr3(n,m,2) = real(fas,sp)
else
!
! Inactive point (not inside grid, can be open boundary)
! defaul instead of xz/yz needed for GPP
! '0' instead of kcs, because TEKAL does not accept '2'
!
glbarr3(n,m,1) = defaul ! amplitudes
glbarr3(n,m,2) = defaul ! phases
endif
enddo
enddo
ierror = nf90_put_var(idfile, idvar(fouvar),glbarr3(1:nmaxus,1,1),start=(/1/),count = (/nmaxus/)) ! write amplitudes
! ierror = nf90_put_var(idfile, idvar(fouvar+1),glbarr3(1:nmaxus,1,2),start=(/1/),count = (/nmaxus/)) ! write phase
ierror = nf90_put_var(idfile, idvar(fouvar+1),glbarr3(1:nmaxus,1,2),start=(/1/),count = (/nmaxus/)) ! write phase
endif
end subroutine wrfous
subroutine wrfouv(ifou ,dtsec ,hdt ,tzone ,gdfourier ,gddimens )
!-------------------------------------------------------------------------------
! $Id: fourier_analysis.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fourier_analysis.f90 $
!!--description-----------------------------------------------------------------
use precision
use mathconsts
use m_d3ddimens
use netcdf
!use datagroups
!use dfparall
!use dffunctionals
!
implicit none
!
! The following list of pointer parameters is used to point inside the gdp structure
!
integer , pointer :: nmaxgl
integer , pointer :: mmaxgl
integer , dimension(:) , pointer :: flayno
integer , dimension(:) , pointer :: fnumcy
integer , dimension(:) , pointer :: ftmsto
integer , dimension(:) , pointer :: ftmstr
integer , dimension(:) , pointer :: idvar
integer , pointer :: ibluv
integer , pointer :: iblqf
integer , pointer :: iblbs
integer , pointer :: iblep
real(fp) , dimension(:) , pointer :: fknfac
real(fp) , dimension(:) , pointer :: foufas
real(fp) , dimension(:,:,:) , pointer :: fousma
real(fp) , dimension(:,:,:) , pointer :: fousmb
real(fp) , dimension(:,:,:) , pointer :: fouvec
real(fp) , dimension(:) , pointer :: fv0pu
character(1) , dimension(:) , pointer :: fouelp
character(16) , dimension(:) , pointer :: founam
character(50) , dimension(:) , pointer :: fouvarnam
integer , pointer :: nofouvar
integer , dimension(:,:) , pointer :: fouref
integer , pointer :: idfile
integer , pointer :: mmax
integer , pointer :: nmax
integer , pointer :: nmaxus
integer , pointer :: ndx
integer , pointer :: lnx
integer , pointer :: ndkx
integer , pointer :: lnkx
!
! Global variables
!
integer , intent(in) :: ifou !< Fourier counter
real(fp) , intent(in) :: dtsec !< Integration time step [in seconds]
real(fp) , intent(in) :: tzone !< Local (FLOW) time - GMT (in hours) => gdp%gdexttim%tzone
real(fp) , intent(in) :: hdt !< Half Integration time step [seconds] => gdp%gdnumeco%hdt
type(gd_fourier) , pointer :: gdfourier
type(gd_dimens) , pointer :: gddimens
!
!
! Local variables
!
integer :: ierror
integer :: fouvar
integer :: m ! Loop counter over MMAX
integer :: md
integer :: n ! Loop counter over NMAXUS
integer :: ncol ! Number of column to write to TEKAL data file
integer :: nd
logical :: ltest ! Help variable for atan2 function test
real(fp) :: a1 ! Used in the computation of el. par.
real(fp) :: a2 ! Used in the computation of el. par.
real(fp) :: alfa ! Value of ALFAS(N,M) in radials
real(fp) :: amp ! Fourier amplitude
real(fp) :: b1 ! Used in the computation of el. par.
real(fp) :: b2 ! Used in the computation of el. par.
real(fp) :: elam ! Eleptic parameter (amplitude)
real(fp) :: elex ! Eleptic parameter (eccentricity)
real(fp) :: elfi ! Eleptic parameter (phase)
real(fp) :: elps ! Eleptic parameter (inclination)
real(fp) :: fas ! Fourier phase
real(fp) :: freqnt ! Frequency in degrees per hour
real(fp) :: fuzeta ! Help variable for fourier value in zeta point in U direction
real(fp) :: fvzeta ! Help variable for fourier value in zeta point in V direction
real(fp) :: r1 ! Used in the computation of el. par.
real(fp) :: r2 ! Used in the computation of el. par.
real(fp) :: shift ! Phase shift
real(fp) :: t1 ! Used in the computation of el. par.
real(fp) :: t2 ! Used in the computation of el. par.
real(fp) :: tfasto ! Stop time in minutes
real(fp) :: tfastr ! Start time in minutes
real(sp) :: defaul ! Default value
character(20) :: namfun ! Local name for fourier function
character(4) :: blnm
real(sp), dimension(:,:), allocatable, save :: glbarr2 ! global arrays in dffunctionals.f90, used in d3d, local here !!
real(sp), dimension(:,:,:), allocatable, save :: glbarr3
real(sp), dimension(:,:,:,:), allocatable, save :: glbarr4
integer, dimension(:,:), allocatable, save :: glbari2
!
!! executable statements -------------------------------------------------------
!
flayno => gdfourier%flayno
fnumcy => gdfourier%fnumcy
ftmsto => gdfourier%ftmsto
ftmstr => gdfourier%ftmstr
idvar => gdfourier%idvar
ibluv => gdfourier%ibluv
iblqf => gdfourier%iblqf
iblbs => gdfourier%iblbs
iblep => gdfourier%iblep
fknfac => gdfourier%fknfac
foufas => gdfourier%foufas
fousma => gdfourier%fousma
fousmb => gdfourier%fousmb
fouvec => gdfourier%fouvec
fv0pu => gdfourier%fv0pu
fouelp => gdfourier%fouelp
founam => gdfourier%founam
fouvarnam => gdfourier%fouvarnam
fouref => gdfourier%fouref
idfile => gdfourier%idfile
nofouvar => gdfourier%nofouvar
mmax => gddimens%mmax
nmax => gddimens%nmax
nmaxus => gddimens%nmaxus
ndx => gddimens%ndx
lnx => gddimens%lnx
ndkx => gddimens%ndkx
lnkx => gddimens%lnkx
mmaxgl => mmax
nmaxgl => nmax
!
! Initialize local variables
!
defaul = -999.0_sp
!
! Frequention := 360 degree / period
! where period = [ (FTMSTO - FMTSTR) * DTSEC ] / [ FNUMCY * 3600 ]
! FOUFAS is defined in RDFOUR as
! FOUFAS = 2 * PI * FNUMCY / [(FTMSTO - FMTSTR) ]
! so FREQNT = FOUFAS * RADDEG * 3600 / DTSEC is OK
!
shift = ftmstr(ifou)*foufas(ifou)
freqnt = foufas(ifou)*raddeg*3600.0_fp/dtsec
tfastr = real(ftmstr(ifou),fp)*dtsec/60.0_fp
tfasto = real(ftmsto(ifou),fp)*dtsec/60.0_fp
!
namfun = founam(ifou)
if (founam(ifou)(:2)=='u1') then
ibluv = ibluv + 1
blnm = 'UV??'
write (blnm(3:4), '(i2.2)') ibluv
namfun = 'velocity'
endif
if (founam(ifou)(:2)=='qx') then
iblqf = iblqf + 1
blnm = 'QF??'
write (blnm(3:4), '(i2.2)') iblqf
namfun = 'unit discharge'
endif
if (founam(ifou)(:2)=='ta') then
iblbs = iblbs + 1
blnm = 'BS??'
write (blnm(3:4), '(i2.2)') iblbs
namfun = 'bed stress'
endif
!
! Write data for user defined dimensions, hence NMAXUS and MMAX
!
if (fouelp(ifou)=='x' .or. fouelp(ifou)=='i') then
if (allocated(glbarr3)) deallocate(glbarr3, stat = ierror)
allocate(glbarr3(nmaxgl,mmaxgl,3), stat = ierror)
glbarr3 = defaul
do n = 1, nmaxus
do m = 1, mmax
!
! Only write values unequal to initial min/max values (-/+1.0e+30)
!
if (comparereal(abs(fousma(n,m,ifou)),1.0e29_fp)==-1) then
glbarr3(n,m,1) = real(fousma(n,m,ifou),sp)
endif
if (comparereal(abs(fousma(n,m,ifou+1)),1.0e29_fp)==-1) then
glbarr3(n,m,2) = real(fousma(n,m,ifou+1),sp)
endif
if (comparereal(abs(fouvec(n,m,ifou)),1.0e29_fp)==-1) then
glbarr3(n,m,3) = real(fouvec(n,m,ifou),sp)
endif
enddo ! over index m
enddo ! over index n
fouvar = fouref(ifou,2)
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,1), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,2), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,3), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
else
do n = 1, nmaxus
do m = 1, mmax
!
! Test FOUSMA and FOUSMB values (<> 0.) for IFOU
!
ltest = (fousma(n, m, ifou)==0.0_fp .and. fousmb(n, m, ifou)==0.0_fp)
if (.not.ltest) then
fousma(n, m, ifou) = fousma(n, m, ifou) &
& *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
fousmb(n, m, ifou) = fousmb(n, m, ifou) &
& *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
amp = sqrt(fousma(n, m, ifou)*fousma(n, m, ifou) &
& + fousmb(n, m, ifou)*fousmb(n, m, ifou))
fas = atan2(fousmb(n, m, ifou), fousma(n, m, ifou)) + shift
if (fnumcy(ifou)==0) then
amp = 0.5_fp*amp*cos(fas)
fas = 0.0_fp
endif
!
! Timezone correction added timezone*phase [degrees/hr].
! foufas is in [rad/timestep]
! halftimestep is in [sec/timestep]
! => timezonecorr = tzone [-] * foufas [rad/timestep] * raddeg [deg/rad] * [sec/hr] / (2 * halftimestep [sec/timestep])
!
fas = fas*raddeg + fv0pu(ifou) - tzone*foufas(ifou)*raddeg*1800.0_fp/hdt
!
! To define FAS between 0 and 360. add 720. to the mod
! function of FAS and re-use the mod function
!
fas = mod(mod(fas, 360.0_fp) + 720.0_fp, 360.0_fp)
amp = amp/fknfac(ifou)
fousma(n, m, ifou) = amp
fousmb(n, m, ifou) = fas
else
fousma(n, m, ifou) = 0.0_fp
fousmb(n, m, ifou) = 0.0_fp
endif
!
! Test FOUSMA and FOUSMB values (<> 0.) for IFOU + 1
!
ltest = (fousma(n, m, ifou + 1)==0.0_fp .and. &
& fousmb(n, m, ifou + 1)==0.0_fp)
if (.not.ltest) then
fousma(n, m, ifou + 1) &
& = fousma(n, m, ifou + 1) *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
fousmb(n, m, ifou + 1) &
& = fousmb(n, m, ifou + 1) *2.0_fp/(real(ftmsto(ifou) - ftmstr(ifou),fp))
amp = sqrt(fousma(n, m, ifou + 1)*fousma(n, m, ifou + 1) &
& + fousmb(n, m, ifou + 1)*fousmb(n, m, ifou + 1))
fas = atan2(fousmb(n, m, ifou + 1), fousma(n, m, ifou + 1)) &
& + shift
if (fnumcy(ifou)==0) then
amp = 0.5_fp*amp*cos(fas)
fas = 0.0_fp
endif
!
! Timezone correction added timezone*phase [degrees/hr].
! foufas is in [rad/timestep]
! halftimestep is in [sec/timestep]
! => timezonecorr = tzone [-] * foufas [rad/timestep] * raddeg [deg/rad] * [sec/hr] / (2 * halftimestep [sec/timestep])
!
fas = fas*raddeg + fv0pu(ifou) - tzone*foufas(ifou)*raddeg*1800.0_fp/hdt
!
! To define FAS between 0 and 360. add 720. to the mod
! function of FAS and re-use the mod function
!
fas = mod(mod(fas, 360.0_fp) + 720.0_fp, 360.0_fp)
amp = amp/fknfac(ifou)
fousma(n, m, ifou + 1) = amp
fousmb(n, m, ifou + 1) = fas
else
fousma(n, m, ifou + 1) = 0.0_fp
fousmb(n, m, ifou + 1) = 0.0_fp
endif
enddo
enddo
! Feed the single precision array glbarr3 to nf90_put_var
! instead of the flexible precision arrays fousma/fousmb
!
if (allocated(glbarr3)) deallocate(glbarr3, stat = ierror)
allocate(glbarr3(nmaxgl,mmaxgl,4), stat = ierror)
glbarr3 = defaul
do n = 1, nmaxus
do m = 1, mmax
!
! Only write values unequal to initial min/max values (-/+1.0e+30)
!
if (comparereal(abs(fousma(n,m,ifou)),1.0e29_fp)==-1) then
glbarr3(n,m,1) = real(fousma(n,m,ifou) , sp)
endif
if (comparereal(abs(fousmb(n,m,ifou)),1.0e29_fp)==-1) then
glbarr3(n,m,2) = real(fousmb(n,m,ifou) , sp)
endif
if (comparereal(abs(fousma(n,m,ifou+1)),1.0e29_fp)==-1) then
glbarr3(n,m,3) = real(fousma(n,m,ifou+1), sp)
endif
if (comparereal(abs(fousmb(n,m,ifou+1)),1.0e29_fp)==-1) then
glbarr3(n,m,4) = real(fousmb(n,m,ifou+1), sp)
endif
enddo
enddo
fouvar = fouref(ifou,2)
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,1), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,2), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,3), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,4), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
endif
!
! Write elliptic parameters to file (if requested)
!
if (fouelp(ifou)=='y') then
iblep = iblep + 1
blnm = 'EP??'
write (blnm(3:4), '(i2.2)') iblep
if (allocated(glbarr3)) deallocate(glbarr3, stat = ierror)
allocate(glbarr3(nmaxgl,mmaxgl,4), stat = ierror)
glbarr3 = defaul
!
do n = 1, nmaxus
do m = 1, mmax
!
! Test for active point and for FOUSMA values (<> defaul)
! for IFOU and IFOU + 1
! when KCS (N,M) = 1 N > 1 and M > 1 per definition
!
ltest = (fousma(n, m, ifou)==defaul .and. fousma(n, m, ifou + 1) &
& ==defaul)
if (.not.ltest) then
!
! Define ellips parameters
!
a1 = fousma(n, m, ifou)*cos(fousmb(n, m, ifou)/raddeg)
a2 = fousma(n, m, ifou + 1)*cos(fousmb(n, m, ifou + 1)/raddeg)
b1 = fousma(n, m, ifou)*sin(fousmb(n, m, ifou)/raddeg)
b2 = fousma(n, m, ifou + 1)*sin(fousmb(n, m, ifou + 1)/raddeg)
!
r1 = 0.5_fp*sqrt((a1 + b2)*(a1 + b2) + (a2 - b1)*(a2 - b1))
r2 = 0.5_fp*sqrt((a1 - b2)*(a1 - b2) + (a2 + b1)*(a2 + b1))
!
! Test ATAN2 input values
!
if ((a2 - b1)==0.0_fp .and. (a1 + b2)==0.0_fp) then
t1 = 0.0_fp
else
t1 = atan2((a2 - b1), (a1 + b2))
endif
!
if ((a2 + b1)==0.0_fp .and. (a1 - b2)==0.0_fp) then
t2 = 0.0_fp
else
t2 = atan2((a2 + b1), (a1 - b2))
endif
!
elam = r1 + r2
!
if ((r1 - r2)==0.0_fp .and. (r1 + r2)==0.0_fp) then
elex = 0.0_fp
elseif ((r1 + r2)==0.0_fp) then
elex = defaul
else
elex = (r1 - r2)/(r1 + r2)
endif
!
elfi = 0.5_fp*(t2 - t1)
elps = 0.5_fp*(t2 + t1)
!
! To define ELFI and ELPS between 0 and 360. add 720.
! to the mod functions and re-use the mod function
!
elfi = mod(mod(elfi*raddeg, 360.0_fp) + 720.0_fp, 360.0_fp)
elps = mod(mod(elps*raddeg, 360.0_fp) + 720.0_fp, 360.0_fp)
!
! Write to file
!
glbarr3(n,m,1) = real(elam,sp)
glbarr3(n,m,2) = real(elex,sp)
glbarr3(n,m,3) = real(elfi,sp)
glbarr3(n,m,4) = real(elps,sp)
else
glbarr3(n,m,:) = defaul
endif
enddo
enddo
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,1), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,2), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,3), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
fouvar = fouvar + 1
ierror = nf90_put_var(idfile, idvar(fouvar) , glbarr3(:,:,4), start=(/ 1, 1/), count = (/nmaxgl, mmaxgl/))
endif
end subroutine wrfouv
end module m_fourier_analysis