!----- 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