!----- 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: meteo1.f90 43143 2015-11-19 08:26:04Z leander $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/meteo1.f90 $
module m_itdate
character(len=8) :: refdat
integer :: itdate !< should be user specified for (asc routines)
integer :: jul0, imonth0, iday0, iyear0
double precision :: Tzone ! doubling with "use m_flowtimes, only : tzone"
end module m_itdate
module timespace_parameters
! enumeration for filetypes van de providers
integer, parameter :: uniform = 1 ! kx values per tijdstap 1 dim arr uni
integer, parameter :: unimagdir = 2 ! kx values per tijdstap 1 dim arr, mag/dir transf op index 1,2 u,v
integer, parameter :: svwp = 3 ! 3 velden per tijdstap 3 dim array noint
integer, parameter :: arcinfo = 4 ! 1 veld per tijdstap 2 dim array bilin/direct
integer, parameter :: spiderweb = 5 ! 3 veld per tijdstap 3 dim array bilin/spw
integer, parameter :: curvi = 6 ! 1 veld per tijdstap 2 dim array bilin/findnm
integer, parameter :: triangulation = 7 ! 1 veld per tijdstap triang
integer, parameter :: triangulationmagdir = 8 ! 2 velden u,v per tijdstap 3 dim array triang, vectormax = 2
! op basis van windreeksen op stations mag/dir
integer, parameter :: poly_tim = 9 ! for line oriented bnd conditions, refs to uniform, fourier or harmonic
integer, parameter :: inside_polygon = 10 ! Constant value inside polygon, used for initial fields.
integer, parameter :: ncgrid = 11 ! NetCDF grid, rectangular type as arcinfo
integer, parameter :: ncflow = 12 ! NetCDF flow, with arbitrary type of input
integer, parameter :: ncwave = 14 ! NetCDF com file, with arbitrary type of input
integer, parameter :: bcascii = 17 ! .bc format as ASCII file
integer, parameter :: max_file_types = 103 ! max nr of supported types for end user in ext file.
! Enumeration for file types of sub-providers (not directly in ext file)
integer, parameter :: fourier = 101 ! period(hrs), ampl(m), phas(deg) NOTE: not directly used in ext file by users.
integer, parameter :: multiple_uni = 102 ! multiple time series, no spatial relation
integer, parameter :: qhtable = 103 ! used to link to dataprovider file
! het filetype legt vast : a) format file
! b) vectormax van grootheid / heden in file
! c) elementset waarop grootheid is gedefinieerd
! d) is daarmee bepalend voor de toepasbare interpolatiemethodes
!
integer :: mdia = 0 ! -1 ! -1 = write dia, 0 = do not write dia
! enumeration for interpolation methods of providers
integer, parameter :: justupdate = 0 ! provider just updates, another provider that
! pointers to this one does the actual interpolation
integer, parameter :: spaceandtime = 1 ! intp space and time (getval)
! keep 2 meteofields in memory
integer, parameter :: spacefirst = 2 ! first intp space (update), next intp. time (getval)
! keep 2 flowfields in memory
integer, parameter :: weightfactors = 3 ! save weightfactors, intp space and time (getval)
! keep 2 pointer- and weight sets in memory
integer, parameter :: spaceandtimeint = 100 ! spaceinterpolation and time=t1 integrated values
end module timespace_parameters
!> Module which constructs and connects the target Items for FM.
!! This is the wrapper between FM and the EC-module.
module m_meteo
use m_ec_module
use m_ec_provider
use MessageHandling
use m_itdate, only: itdate
use m_flowtimes, only : tzone
use m_wind
use m_waves
use m_ship
use m_flowexternalforcings
use unstruc_messages
use time_module
use m_observations
use string_module
implicit none
type(tEcInstance), pointer, save :: ecInstancePtr !< FM's instance of the EC-module.
character(maxMessageLen) :: message !< EC's message, to be passed to FM's log.
!
integer, dimension(:), allocatable, target :: item_tracerbnd !< dim(numtracers)
!
integer, target :: item_windx !< Unique Item id of the ext-file's 'windx' quantity's x-component.
integer, target :: item_windy !< Unique Item id of the ext-file's 'windy' quantity's y-component.
integer, target :: item_windxy_x !< Unique Item id of the ext-file's 'windxy' quantity's x-component.
integer, target :: item_windxy_y !< Unique Item id of the ext-file's 'windxy' quantity's y-component.
integer, target :: item_apwxwy_p !< Unique Item id of the ext-file's 'airpressure_windx_windy' quantity 'p'.
integer, target :: item_apwxwy_x !< Unique Item id of the ext-file's 'airpressure_windx_windy' quantity 'x'.
integer, target :: item_apwxwy_y !< Unique Item id of the ext-file's 'airpressure_windx_windy' quantity 'y'.
integer, target :: item_waterlevelbnd !< Unique Item id of the ext-file's 'waterlevelbnd' quantity's ...-component.
integer, target :: item_atmosphericpressure !< Unique Item id of the ext-file's 'atmosphericpressure' quantity
integer, target :: item_velocitybnd !< Unique Item id of the ext-file's 'velocitybnd' quantity
integer, target :: item_salinitybnd !< Unique Item id of the ext-file's 'salinitybnd' quantity
integer, target :: item_temperaturebnd !< Unique Item id of the ext-file's 'temperaturebnd' quantity
integer, target :: item_sedimentbnd !< Unique Item id of the ext-file's 'sedimentbnd' quantity
integer, target :: item_tangentialvelocitybnd !< Unique Item id of the ext-file's 'tangentialvelocitybnd' quantity
integer, target :: item_uxuyadvectionvelocitybnd !< Unique Item id of the ext-file's 'uxuyadvectionvelocitybnd'
integer, target :: item_normalvelocitybnd !< Unique Item id of the ext-file's 'normalvelocitybnd' quantity
integer, target :: item_rainfall !< Unique Item id of the ext-file's 'rainfall' quantity
integer, target :: item_qhbnd !< Unique Item id of the ext-file's 'qhbnd' quantity
integer, target :: item_shiptxy !< Unique Item id of the ext-file's 'shiptxy' quantity
integer, target :: item_movingstationtxy !< Unique Item id of the ext-file's 'movingstationtxy' quantity
integer, target :: item_pump !< Unique Item id of the ext-file's 'pump' quantity
integer, target :: item_damlevel !< Unique Item id of the ext-file's 'damlevel' quantity
integer, target :: item_gateloweredgelevel !< Unique Item id of the ext-file's 'gateloweredgelevel' quantity
integer, target :: item_generalstructure !< Unique Item id of the ext-file's 'generalstructure' quantity
integer, target :: item_hacs_humidity !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness_solarradiation' quantity 'humidity'
integer, target :: item_hacs_airtemperature !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness_solarradiation' quantity 'airtemperature'
integer, target :: item_hacs_cloudiness !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness_solarradiation' quantity 'cloudiness'
integer, target :: item_hacs_solarradiation !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness_solarradiation' quantity 'solarradiation'
integer, target :: item_hac_humidity !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness' quantity 'humidity'
integer, target :: item_hac_airtemperature !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness' quantity 'airtemperature'
integer, target :: item_hac_cloudiness !< Unique Item id of the ext-file's 'humidity_airtemperature_cloudiness' quantity 'cloudiness'
integer, target :: item_discharge_salinity_temperature_sorsin !< Unique Item id of the ext-file's 'discharge_salinity_temperature_sorsin' quantity
integer, target :: item_hrms !< Unique Item id of the ext-file's 'item_hrms' quantity
integer, target :: item_tp !< Unique Item id of the ext-file's 'item_tp' quantity
integer, target :: item_dir !< Unique Item id of the ext-file's 'item_dir' quantity
integer, target :: item_fx !< Unique Item id of the ext-file's 'item_fx' quantity
integer, target :: item_fy !< Unique Item id of the ext-file's 'item_fy' quantity
integer, target :: item_wsbu !< Unique Item id of the ext-file's 'item_wsbu' quantity
integer, target :: item_wsbv !< Unique Item id of the ext-file's 'item_wsbv' quantity
integer, target :: item_mx !< Unique Item id of the ext-file's 'item_mx' quantity
integer, target :: item_my !< Unique Item id of the ext-file's 'item_my' quantity
!
integer :: n_qhbnd !< Number of already connected qh-boundaries.
interface ec_gettimespacevalue
module procedure ec_gettimespacevalue_by_itemID
module procedure ec_gettimespacevalue_by_name
end interface ec_gettimespacevalue
contains
!> Initialize the module variables.
subroutine init_variables()
ecInstancePtr => null()
message = ' '
!
item_windx = ec_undef_int
item_windy = ec_undef_int
item_windxy_x = ec_undef_int
item_windxy_y = ec_undef_int
item_apwxwy_p = ec_undef_int
item_apwxwy_x = ec_undef_int
item_apwxwy_y = ec_undef_int
item_waterlevelbnd = ec_undef_int
item_atmosphericpressure = ec_undef_int
item_velocitybnd = ec_undef_int
item_salinitybnd = ec_undef_int
item_temperaturebnd = ec_undef_int
item_sedimentbnd = ec_undef_int
item_tangentialvelocitybnd = ec_undef_int
item_uxuyadvectionvelocitybnd = ec_undef_int
item_normalvelocitybnd = ec_undef_int
item_rainfall = ec_undef_int
item_qhbnd = ec_undef_int
item_shiptxy = ec_undef_int
item_movingstationtxy = ec_undef_int
item_pump = ec_undef_int
item_damlevel = ec_undef_int
item_gateloweredgelevel = ec_undef_int
item_generalstructure = ec_undef_int
item_hacs_humidity = ec_undef_int
item_hacs_airtemperature = ec_undef_int
item_hacs_cloudiness = ec_undef_int
item_hacs_solarradiation = ec_undef_int
item_hac_humidity = ec_undef_int
item_hac_airtemperature = ec_undef_int
item_hac_cloudiness = ec_undef_int
item_discharge_salinity_temperature_sorsin = ec_undef_int
item_hrms = ec_undef_int
item_tp = ec_undef_int
item_dir = ec_undef_int
item_fx = ec_undef_int
item_fy = ec_undef_int
item_wsbu = ec_undef_int
item_wsbv = ec_undef_int
item_mx = ec_undef_int
item_my = ec_undef_int
!
n_qhbnd = 0
!
! tracers
if ( allocated(item_tracerbnd) ) deallocate(item_tracerbnd)
allocate(item_tracerbnd(numtracers))
item_tracerbnd = ec_undef_int
end subroutine init_variables
! ==========================================================================
!> Translate FM's meteo1 'filetype' enum to EC's 'provFile' enum.
subroutine filetype_fm_to_ec(filetype, ec_filetype)
use timespace_parameters
implicit none
integer, intent(in) :: filetype
integer, intent(out) :: ec_filetype
!
select case (filetype)
case (uniform) ! 1
ec_filetype = provFile_uniform
case (unimagdir) ! 2
ec_filetype = provFile_unimagdir
case (svwp) ! 3
ec_filetype = provFile_svwp
case (arcinfo) ! 4
ec_filetype = provFile_arcinfo
case (spiderweb) ! 5
ec_filetype = provFile_spiderweb
case (curvi) ! 6
ec_filetype = provFile_curvi
case (triangulation) ! 7
ec_filetype = provFile_samples
case (triangulationmagdir) ! 8
ec_filetype = provFile_triangulationmagdir
case (poly_tim) ! 9
ec_filetype = provFile_poly_tim
case (ncgrid, ncwave) ! 11, 14
ec_filetype = provFile_netcdf
case (ncflow) ! 12
ec_filetype = provFile_undefined ! only used for timespaceinitialfield, no EC yet.
case (bcascii) ! 17
ec_filetype = provFile_bc
case (fourier) ! 101
ec_filetype = provFile_fourier
case default
ec_filetype = provFile_undefined
end select
end subroutine filetype_fm_to_ec
! ==========================================================================
!> Translate FM's meteo1 'method' enum to EC's 'interpolate' enum.
subroutine method_fm_to_ec(method, ec_method)
integer, intent(in) :: method
integer, intent(out) :: ec_method
!
select case (method)
case (0)
ec_method = interpolate_passthrough
case (1)
ec_method = interpolate_timespace
case (2)
ec_method = interpolate_spacetime
case (3)
ec_method = interpolate_spacetimeSaveWeightFactors
case (4)
ec_method = interpolate_space
case (5)
ec_method = interpolate_time
case (7)
ec_method = interpolate_time_extrapolation_ok
case default
ec_method = interpolate_unknown
end select
end subroutine method_fm_to_ec
! ==========================================================================
!> Translate FM's meteo1 'operand' enum to EC's 'operand' enum.
subroutine operand_fm_to_ec(operand, ec_operand)
character, intent(in) :: operand
integer, intent(out) :: ec_operand
!
select case (operand)
case ('O')
ec_operand = operand_replace
case ('+')
ec_operand = operand_add
case default
ec_operand = operand_undefined
end select
end subroutine operand_fm_to_ec
! ==========================================================================
!> Translate EC's 'filetype' to EC's 'convType' enum.
subroutine ec_filetype_to_conv_type(filetype, convtype)
integer, intent(in) :: filetype
integer, intent(out) :: convtype
!
select case (filetype)
case (provFile_uniform)
convtype = convType_uniform
case (provFile_fourier)
convtype = convType_fourier
case (provFile_unimagdir)
convtype = convType_unimagdir
case (provFile_svwp)
convtype = convType_undefined ! not yet implemented
case (provFile_arcinfo)
convtype = convType_arcinfo
case (provFile_spiderweb)
convtype = convType_spiderweb
case (provFile_curvi)
convtype = convType_curvi
case (provFile_samples)
convtype = convType_undefined ! not yet implemented
case (provFile_triangulationmagdir)
convtype = convType_undefined ! not yet implemented
case (provFile_poly_tim)
convtype = convType_polytim
case (provFile_netcdf)
convtype = convType_netcdf
case default
convtype = convType_undefined
end select
end subroutine ec_filetype_to_conv_type
! ==========================================================================
!> Convert quantity names as given in user input (ext file)
!! to accepted Unstruc names (as used in Fortran code)
!! Note: for old-style ext quantities, fm_name==input_name, e.g. waterlevelbnd.
!subroutine bndname_to_fm(input_name, fm_name)
! character(len=*), intent(in) :: input_name !< given by the user
! character(len=*), intent(out) :: fm_name !< known within FM
!
! character(len=256) :: tempname
!
! fm_name = input_name
! tempname = input_name
! call str_upper(tempname)
! call remove_substr(tempname,'_')
! call remove_substr(tempname,'-')
! call remove_substr(tempname,' ')
!
! select case (trim(tempname))
! case ('WATERLEVEL','VELOCITY','SALINITY','TEMPERATURE','SEDIMENT','TANGENTIALVELOCITY','NORMALVELOCITY','QH','TRACER')
! ! These are new-ext-style quantities: FM needs additional 'bnd' behind quantityid
! fm_name = trim(tempname)//'bnd'
! call str_lower(fm_name)
! end select
!end subroutine bndname_to_fm
! ==========================================================================
!> Translate EC's ext.force-file's item name to the integer EC item handle and to
!> the data pointer(s), i.e. the array that will contain the values of the target item
function fm_ext_force_name_to_ec_item(trname, qidname, &
itemPtr1, itemPtr2, itemPtr3, itemPtr4, &
dataPtr1, dataPtr2, dataPtr3, dataPtr4 ) result(success)
logical :: success
character(len=NAMTRACLEN), intent(in) :: trname
character(len=NAMTRACLEN), intent(in) :: qidname
integer, pointer :: itemPtr1, itemPtr2, itemPtr3, itemPtr4
real(hp), dimension(:), pointer :: dataPtr1, dataPtr2, dataPtr3, dataPtr4
! for tracers:
integer :: itrac
integer, external :: findname
success = .true.
itemPtr1 => null()
itemPtr2 => null()
itemPtr3 => null()
itemPtr4 => null()
dataPtr1 => null()
dataPtr2 => null()
dataPtr3 => null()
dataPtr4 => null()
select case (qidname)
case ('windx')
itemPtr1 => item_windx
dataPtr1 => wx
case ('windy')
itemPtr1 => item_windy
dataPtr1 => wy
case ('windxy')
itemPtr1 => item_windxy_x
dataPtr1 => wx
itemPtr2 => item_windxy_y
dataPtr2 => wy
case ('airpressure_windx_windy')
itemPtr1 => item_apwxwy_p
dataPtr1 => patm
itemPtr2 => item_apwxwy_x
dataPtr2 => ec_pwxwy_x
itemPtr3 => item_apwxwy_y
dataPtr3 => ec_pwxwy_y
case ('waterlevelbnd', 'neumannbnd', 'riemannbnd', 'outflowbnd')
itemPtr1 => item_waterlevelbnd
dataPtr1 => zbndz
case ('velocitybnd', 'dischargebnd', 'riemann_velocitybnd')
itemPtr1 => item_velocitybnd
dataPtr1 => zbndu
case ('salinitybnd')
itemPtr1 => item_salinitybnd
dataPtr1 => zbnds
case ('temperaturebnd')
itemPtr1 => item_temperaturebnd
dataPtr1 => zbndTM
case ('sedimentbnd')
itemPtr1 => item_sedimentbnd
dataPtr1 => zbndsd
case ('tangentialvelocitybnd')
itemPtr1 => item_tangentialvelocitybnd
dataPtr1 => zbndt
case ('uxuyadvectionvelocitybnd')
itemPtr1 => item_uxuyadvectionvelocitybnd
dataPtr1 => zbnduxy
case ('normalvelocitybnd')
itemPtr1 => item_normalvelocitybnd
dataPtr1 => zbndn
case ('atmosphericpressure')
itemPtr1 => item_atmosphericpressure
dataPtr1 => patm
case ('rainfall')
itemPtr1 => item_rainfall
dataPtr1 => rain
case ('qhbnd')
itemPtr1 => item_qhbnd
dataPtr1 => qhbndz
case ('shiptxy')
itemPtr1 => item_shiptxy
dataPtr1 => xyship
case ('movingstationtxy')
itemPtr1 => item_movingstationtxy
dataPtr1 => xyobs
case ('pump')
itemPtr1 => item_pump
!dataPtr1 => qpump
case ('damlevel')
itemPtr1 => item_damlevel
!dataPtr1 => zcdam
case ('gateloweredgelevel')
itemPtr1 => item_gateloweredgelevel
dataPtr1 => zgate
case ('generalstructure')
itemPtr1 => item_generalstructure
dataPtr1 => zcgen
case ('humidity_airtemperature_cloudiness')
itemPtr1 => item_hac_humidity
dataPtr1 => rhum
itemPtr2 => item_hac_airtemperature
dataPtr2 => tair
itemPtr3 => item_hac_cloudiness
dataPtr3 => clou
case ('humidity_airtemperature_cloudiness_solarradiation')
itemPtr1 => item_hacs_humidity
dataPtr1 => rhum
itemPtr2 => item_hacs_airtemperature
dataPtr2 => tair
itemPtr3 => item_hacs_cloudiness
dataPtr3 => clou
itemPtr4 => item_hacs_solarradiation
dataPtr4 => qrad
case ('discharge_salinity_temperature_sorsin')
itemPtr1 => item_discharge_salinity_temperature_sorsin
dataPtr1 => qstss
case ('hrms')
itemPtr1 => item_hrms
dataPtr1 => hwav
case ('tp', 'tps', 'rtp')
itemPtr1 => item_tp
dataPtr1 => twav
case ('dir')
itemPtr1 => item_dir
dataPtr1 => phiwav
case ('fx')
itemPtr1 => item_fx
dataPtr1 => sxwav
case ('fy')
itemPtr1 => item_fy
dataPtr1 => sywav
case ('wsbu')
itemPtr1 => item_wsbu
dataPtr1 => sbxwav
case ('wsbv')
itemPtr1 => item_wsbv
dataPtr1 => sbywav
case ('mx')
itemPtr1 => item_mx
dataPtr1 => mxwav
case ('my')
itemPtr1 => item_my
dataPtr1 => mywav
case ('tracerbnd')
! get tracer (boundary) number
itrac = findname(numtracers, trnames, trname)
itemPtr1 => item_tracerbnd(itrac)
dataPtr1 => bndtr(itrac)%z
case default
call mess(LEVEL_FATAL, 'm_meteo::fm_ext_force_name_to_ec_item: Unsupported quantity specified in ext-file (construct target field): '//qidname)
success = .false.
end select
end function fm_ext_force_name_to_ec_item
! ==========================================================================
!> Construct and initialize a new Instance of the EC-module.
subroutine initialize_ec_module()
use m_sferic
implicit none
! FM re-initialize call: First destroy the EC-module instance.
if (associated(ecInstancePtr)) then
if (.not. ecFreeInstance(ecInstancePtr)) then
message = dumpECMessageStack(LEVEL_WARN,callback_msg)
end if
end if
! FM initialize call or second phase of re-initialize call.
if (.not. associated(ecInstancePtr)) then
call init_variables()
if (.not. ecCreateInstance(ecInstancePtr)) then
message = dumpECMessageStack(LEVEL_WARN,callback_msg)
end if
end if
if (jsferic == 1) then
ecInstancePtr%coordsystem = EC_COORDS_SFERIC
else
ecInstancePtr%coordsystem = EC_COORDS_CARTHESIAN
endif
end subroutine initialize_ec_module
! ==========================================================================
!> Helper function for creating and initializing a target Item.
function createItem(instancePtr, itemId, quantityId, elementSetId, fieldId) result(success)
logical :: success !< function status
type(tEcInstance), pointer :: instancePtr !<
integer, intent(inout) :: itemId !< Unique Item id.
integer, intent(inout) :: quantityId !< Unique Quantity id.
integer, intent(inout) :: elementSetId !< Unique ElementSet id.
integer, intent(inout) :: fieldId !< Unique Field id.
!
success = .true.
if (itemId == ec_undef_int) then
itemId = ecCreateItem(ecInstancePtr)
success = ecSetItemRole(instancePtr, itemId, itemType_target)
if (success) success = ecSetItemQuantity(instancePtr, itemId, quantityId)
if (success) success = ecSetItemElementSet(instancePtr, itemId, elementSetId)
if (success) success = ecSetItemTargetField(instancePtr, itemId, fieldId)
end if
end function createItem
! ==========================================================================
!> Helper function for initializing a Converter.
function initializeConverter(instancePtr, converterId, convtype, operand, method, srcmask) result(success)
logical :: success !< function status
type(tEcInstance), pointer :: instancePtr !<
integer :: converterId !<
integer :: convtype !<
integer :: operand !<
integer :: method !<
type (tEcMask), optional :: srcmask !<
!
success = ecSetConverterType(instancePtr, converterId, convtype)
if (success) success = ecSetConverterOperand(instancePtr, converterId, operand)
if (success) success = ecSetConverterInterpolation(instancePtr, converterId, method)
if (present(srcmask)) then
if (success) success = ecSetConverterMask(instancePtr, converterId, srcmask)
end if
end function initializeConverter
! ==========================================================================
!> Helper function for initializing a Connection.
function initializeConnection(instancePtr, connectionId, sourceItemId, targetItemId) result(success)
logical :: success !< function status
type(tEcInstance), pointer :: instancePtr !<
integer, intent(inout) :: connectionId !<
integer, intent(inout) :: sourceItemId !<
integer, intent(inout) :: targetItemId !<
!
success = ecAddConnectionSourceItem(instancePtr, connectionId, sourceItemId)
if (success) success = ecAddConnectionTargetItem(instancePtr, connectionId, targetItemId)
if (success) success = ecAddItemConnection(instancePtr, targetItemId, connectionId)
end function initializeConnection
! ==========================================================================
!> Helper function for Connection initialization.
function checkFileType(actualfiletype, requiredfiletype, name) result(success)
logical :: success !< function status
integer, intent(in) :: actualfiletype !< EC-module's filetype enumeration.
integer, intent(in) :: requiredfiletype !< EC-module's filetype enumeration.
character(*), intent(in) :: name !< Name for the target Quantity.
!
success = .true.
if (.not. actualfiletype == requiredfiletype) then
message = 'm_meteo::checkFileType: Unsupported filetype for quantity '//name//'.'
success = .false.
end if
end function checkFileType
! ==========================================================================
!> Replacement function for FM's meteo1 'addtimespacerelation' function.
logical function ec_addtimespacerelation(name, x, y, mask, vectormax, filename, filetype, method, operand, &
xyen, z, targetIndex, forcingfile, srcmaskfile, dtnodal, quiet)
use m_ec_support, only: ecSupportFindFileReader ! TODO: Refactor this private data access (UNST-703).
use m_ec_filereader_read, only: ecParseARCinfoMask
use m_flow, only: kmx
character(len=*), intent(inout) :: name !< Name for the target Quantity, possibly compounded with a tracer name.
real(hp), dimension(:), intent(in) :: x !< Array of x-coordinates for the target ElementSet.
real(hp), dimension(:), intent(in) :: y !< Array of y-coordinates for the target ElementSet.
integer, intent(in) :: vectormax !< Vector max (length of data values at each element location).
integer, dimension(:), intent(in) :: mask !< Array of masking values for the target ElementSet.
character(len=*), intent(in) :: filename !< File name of meteo data file.
integer, intent(in) :: filetype !< FM's filetype enumeration.
integer, intent(in) :: method !< FM's method enumeration.
character(len=1), intent(in) :: operand !< FM's operand enumeration.
real(hp), optional, intent(in) :: xyen(:,:) !< FM's distance tolerance / cellsize of ElementSet.
real(hp), dimension(:), optional, target, intent(in) :: z !< FM's array of z/sigma coordinates
integer, optional, intent(in) :: targetIndex !< target position or rank of (complete!) vector in target array
character(len=*), optional, intent(in) :: forcingfile !< file containing the forcing data for pli-file 'filename'
character(len=*), optional, intent(in) :: srcmaskfile !< file containing mask applicable to the arcinfo source data
real(hp), optional, intent(in) :: dtnodal !< update interval for nodal factors
logical, optional, intent(in) :: quiet !< When .true., in case of errors, do not write the errors to screen/dia at the end of the routine.
!
integer :: ec_filetype !< EC-module's provFile_ enumeration.
integer :: ec_convtype !< EC-module's convType_ enumeration.
integer :: ec_method !< EC-module's interpolate_ enumeration.
integer :: ec_operand !< EC-module's operand_ enumeration.
!
integer :: fileReaderId !< Unique FileReader id.
integer :: quantityId !< Unique Quantity id.
integer :: elementSetId !< Unique ElementSet id.
integer :: fieldId !< Unique Field id.
integer :: fieldId_2 !< Unique Field id.
integer :: fieldId_3 !< Unique Field id.
integer :: fieldId_4 !< Unique Field id.
integer :: converterId !< Unique Converter id.
integer :: connectionId !< Unique Connection id.
integer :: sourceItemId !< Unique source item id.
integer :: sourceItemId_2 !< Unique additional source item id.
integer :: sourceItemId_3 !< Unique additional third source item id.
integer :: sourceItemId_4 !< Unique additional fourth source item id.
integer :: i !< loop counter
integer :: iostat
!
character(len=maxnamelen) :: sourceItemName !< name of source item (as created by provider)
character(len=maxnamelen) :: target_name !< Unstruc target name derived from user-specified name
integer, pointer :: targetItemPtr1 => null() !< pointer to the target item id
integer, pointer :: targetItemPtr2 => null() !< pointer to optional second target item id (e.g. in case of windxy)
integer, pointer :: targetItemPtr3 => null() !< pointer to optional third target item id (e.g. in case of spiderweb)
integer, pointer :: targetItemPtr4 => null() !< pointer to optional fourth target item id (e.g. in case of hacs)
real(hp), dimension(:), pointer :: dataPtr1 => null() !< Pointer to FM's 1D data arrays.
real(hp), dimension(:), pointer :: dataPtr2 => null() !< Pointer to FM's optional extra 1D data array (e.g. in case of windxy)
real(hp), dimension(:), pointer :: dataPtr3 => null() !< Pointer to FM's optional third 1D data array (e.g. in case of spiderweb)
real(hp), dimension(:), pointer :: dataPtr4 => null() !< Pointer to FM's optional fourth 1D data array (e.g. in case of hacs)
type(tEcFileReader) , pointer :: fileReaderPtr => null() !<
logical :: success
logical :: quiet_
integer :: itrac !< tracer index
integer :: n_rows, n_cols
character(len=NAMTRACLEN) :: trname, qidname
integer, external :: findname
type (tEcMask) :: srcmask
character(len=maxMessageLen) :: messages(15)
ec_addtimespacerelation = .false.
if (present(quiet)) then
quiet_ = quiet
else
quiet_ = .false. ! Default: print errors at the end of routine, if no success
end if
! ========================================================
! Translate FM's enumerations to EC-module's enumerations.
! ========================================================
call filetype_fm_to_ec(filetype, ec_filetype)
if (ec_filetype == provFile_undefined) then
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype.')
return
end if
call method_fm_to_ec(method, ec_method)
if (ec_method == interpolate_unknown) then
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported method.')
return
end if
call operand_fm_to_ec(operand, ec_operand)
if (ec_operand == operand_undefined) then
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported operand.')
return
end if
! =================================================
! Convert ext file names to accepted Unstruc names.
! =================================================
! Name conversion: (targetname=qidname==name for all names, except name=tracerbndfoo --> qidname=tracerbnd)
call get_tracername(name, trname, qidname)
target_name = qidname
call clearECMessage()
! ============================================================
! Construct the FileReader, which constructs the source Items.
! ============================================================
fileReaderId = ecCreateFileReader(ecInstancePtr)
fileReaderPtr => ecSupportFindFileReader(ecInstancePtr, fileReaderId) ! TODO: Refactor this private data access (UNST-703).
fileReaderPtr%vectormax = vectormax
if (present(forcingfile)) then
if (present(dtnodal)) then
success = ecSetFileReaderProperties(ecInstancePtr, fileReaderId, ec_filetype, filename, itdate, tzone, ec_second, name, forcingfile=forcingfile, dtnodal=dtnodal)
else
success = ecSetFileReaderProperties(ecInstancePtr, fileReaderId, ec_filetype, filename, itdate, tzone, ec_second, name, forcingfile=forcingfile)
end if
!message = dumpECMessageStack(LEVEL_WARN,callback_msg)
if (.not. success) then
goto 1234
end if
if (ecAtLeastOnePointIsCorrection) then ! TODO: Refactor this shortcut (UNST-180).
ecAtLeastOnePointIsCorrection = .false. ! TODO: Refactor this shortcut (UNST-180).
ec_addtimespacerelation = .true.
return
end if
else
if (present(dtnodal)) then
success = ecSetFileReaderProperties(ecInstancePtr, fileReaderId, ec_filetype, filename, itdate, tzone, ec_second, name, dtnodal=dtnodal)
else
success = ecSetFileReaderProperties(ecInstancePtr, fileReaderId, ec_filetype, filename, itdate, tzone, ec_second, name)
end if
if (.not. success) then
! message = ecGetMessage()
! message = dumpECMessageStack(LEVEL_WARN,callback_msg)
! NOTE: do all error dumping (if any) at the end of this routine at label 1234
! NOTE: in relation to WAVE: all calling WAVE-related routines now pass quiet=.true. to this addtimespace routine.
! When running online with WAVE and the first WAVE calculation is after the first DFlowFM calculation,
! this message will be generated. This must be a warning: notify the user that DFlowFM is going to do
! a calculation with zero wave values. This message should be written every time step, until proper
! wave data is available. The user has to check whether this behaviour is as expected.
goto 1234
end if
end if
! ==============================
! Construct the target Quantity.
! ==============================
quantityId = ecCreateQuantity(ecInstancePtr)
if (.not. ecSetQuantity(ecInstancePtr, quantityId, target_name, ' ', vectormax)) then ! TODO: VectorMax feels like a target Field property. (UNST-720)
goto 1234
end if
! ================================
! Construct the target ElementSet.
! ================================
elementSetId = ecCreateElementSet(ecInstancePtr)
if (ec_filetype /= provFile_spiderweb) then
success = ecSetElementSetType(ecInstancePtr, elementSetId, elmSetType_cartesian)
if (success) success = ecSetElementSetXArray(ecInstancePtr, elementSetId, x)
if (success) success = ecSetElementSetYArray(ecInstancePtr, elementSetId, y)
if (success) success = ecSetElementSetMaskArray(ecInstancePtr, elementSetId, mask)
if (success) success = ecSetElementSetNumberOfCoordinates(ecInstancePtr, elementSetId, size(x))
else
success = ecSetElementSetType(ecInstancePtr, elementSetId, elmSetType_spheric)
if (success) success = ecSetElementSetLatitudeArray(ecInstancePtr, elementSetId, x)
if (success) success = ecSetElementSetLongitudeArray(ecInstancePtr, elementSetId, y)
if (success) success = ecSetElementSetMaskArray(ecInstancePtr, elementSetId, mask)
if (success) success = ecSetElementSetNumberOfCoordinates(ecInstancePtr, elementSetId, size(x))
end if
if (present(xyen)) then
if (success) success = ecSetElementSetXyen(ecInstancePtr, elementSetId, xyen)
end if
if (present(z)) then ! 3D
if (success) success = ecSetElementSetZArray(ecInstancePtr, elementSetId, z, .true.)
if (success) success = ecSetElementSetItype3D(ecInstancePtr, elementSetID, 0) ! sigma layers
! add 3D settings if needed
if (ec_filetype == provFile_poly_tim .and. (target_name == 'salinitybnd' .or. target_name == 'temperaturebnd' .or. target_name == 'tracerbnd')) then
if (success) success = ecSetElementSetMaskArray(ecInstancePtr, elementSetId, mask)
if (success) success = ecSetElementSetNumberOfCoordinates(ecInstancePtr, elementSetId, size(x))
end if
end if
if (.not. success) then
goto 1234
end if
! ==============================================
! Construct the target field and the target item
! ==============================================
! determine which target item (id) will be created, and which FM data array has to be used
if (.not. fm_ext_force_name_to_ec_item(trname, qidname, &
targetItemPtr1, targetItemPtr2, targetItemPtr3, targetItemPtr4, &
dataPtr1 , dataPtr2 , dataPtr3 , dataPtr4 )) then
return
end if
! Create the field and the target item, and if needed additional ones.
fieldId = ecCreateField(ecInstancePtr)
success = ecSetField1dArray(ecInstancePtr, fieldId, dataPtr1)
if (success) success = createItem(ecInstancePtr, targetItemPtr1, quantityId, elementSetId, fieldId)
if (associated(targetItemPtr2)) then
! second field (e.g. for 'windxy')
fieldId_2 = ecCreateField(ecInstancePtr)
if (success) success = ecSetField1dArray(ecInstancePtr, fieldId_2, dataPtr2)
if (success) success = createItem(ecInstancePtr, targetItemPtr2, quantityId, elementSetId, fieldId_2)
end if
if (associated(targetItemPtr3)) then
! third field (e.g. for 'airpressure_windx_windy'
fieldId_3 = ecCreateField(ecInstancePtr)
if (success) success = ecSetField1dArray(ecInstancePtr, fieldId_3, dataPtr3)
if (success) success = createItem(ecInstancePtr, targetItemPtr3, quantityId, elementSetId, fieldId_3)
end if
if (associated(targetItemPtr4)) then
! third field (e.g. for 'humidity_airtemperatur_cloudiness_solarradiation'
fieldId_4 = ecCreateField(ecInstancePtr)
if (success) success = ecSetField1dArray(ecInstancePtr, fieldId_4, dataPtr4)
if (success) success = createItem(ecInstancePtr, targetItemPtr4, quantityId, elementSetId, fieldId_4)
end if
if (.not. success) then
goto 1234
end if
! ==========================
! Construct a new Converter.
! ==========================
call ec_filetype_to_conv_type(ec_filetype, ec_convtype)
if (ec_convtype == convType_undefined) then
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported converter.')
return
end if
converterId = ecCreateConverter(ecInstancePtr)
select case(target_name)
case ('shiptxy', 'movingstationtxy', 'discharge_salinity_temperature_sorsin', 'pump', 'damlevel', 'gateloweredgelevel', 'generalstructure')
! for the FM 'target' arrays, the index is provided by the caller
if (.not. present(targetIndex)) then
message = 'Internal program error: missing targetIndex for quantity '''//trim(target_name)
call mess(LEVEL_ERROR, message)
return
end if
success = initializeConverter(ecInstancePtr, converterId, ec_convtype, operand_replace_element, ec_method)
if (success) success = ecSetConverterElement(ecInstancePtr, converterId, targetIndex)
case ('qhbnd')
! count qh boundaries
n_qhbnd = n_qhbnd + 1
success = initializeConverter(ecInstancePtr, converterId, ec_convtype, operand_replace_element, interpolate_passthrough)
! Each qhbnd polytim file replaces exactly one element in the target data array.
! Converter will put qh value in target_array(n_qhbnd)
if (success) success = ecSetConverterElement(ecInstancePtr, converterId, n_qhbnd)
case ('windx', 'windy', 'windxy', 'atmosphericpressure', 'airpressure_windx_windy')
if (present(srcmaskfile)) then
if (ec_filetype == provFile_arcinfo .or. ec_filetype == provFile_curvi) then
if (.not.ecParseARCinfoMask(srcmaskfile, srcmask, fileReaderPtr)) then
write (msgbuf, '(a,a,a,a,a)') 'Error while reading mask file ''', trim(srcmaskfile),'''.'
call err_flush()
return
endif
if (.not.initializeConverter(ecInstancePtr, converterId, ec_convtype, ec_operand, ec_method, srcmask=srcmask)) then
write (msgbuf, '(a,a,a,a,a)') 'Error while setting mask to converter (file=''', trim(srcmaskfile), ''', associated with meteo file ''', trim(filename), '''.'
call err_flush()
return
end if
end if
else
success = initializeConverter(ecInstancePtr, converterId, ec_convtype, ec_operand, ec_method)
end if
case default
success = initializeConverter(ecInstancePtr, converterId, ec_convtype, ec_operand, ec_method)
end select
if (.not. success) then
goto 1234
end if
! ================================================================
! Construct a new Connection, and connect source and target Items.
! ================================================================
connectionId = ecCreateConnection(ecInstancePtr)
if (.not. ecSetConnectionConverter(ecInstancePtr, connectionId, converterId)) then
goto 1234
end if
! determine the source item's name
! note 1: this can be determined (and be improved) when creating the file reader
! note 2: the source item's name is set in the select case switch below. In some cases
! of this switch ('special cases') the source-target connections is established
! immediatly, and sourceItemName is NOT set.
! So the generic 'connect source and target' statements after the switch are
! only executed if sourceItemName IS set.
!
sourceItemName = ' '
sourceItemId = 0
sourceItemId_2 = 0
sourceItemId_3 = 0
sourceItemId_4 = 0
select case (target_name)
case ('shiptxy' , 'movingstationtxy', 'discharge_salinity_temperature_sorsin')
if (.not. checkFileType(ec_filetype, provFile_uniform, target_name)) then
return
end if
! the file reader will have created an item called 'uniform_item'
sourceItemName = 'uniform_item'
case ('pump','generalstructure','damlevel','gateloweredgelevel')
if (checkFileType(ec_filetype, provFile_uniform, target_name)) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'uniform_item')
if (sourceItemId==ec_undef_int) then
! Add something to the EC message stack about missing source item
return
endif
if (.not.ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)) return
else if (checkFileType(ec_filetype, provFile_fourier, target_name)) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'period')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'magnitude')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'phase')
if ((sourceItemId==ec_undef_int) .or. (sourceItemId_2==ec_undef_int) .or. (sourceItemId_3==ec_undef_int)) then
! Add something to the EC message stack about missing source item
return
else
if (.not.ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)) return
if (.not.ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_2)) return
if (.not.ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_3)) return
endif
else
! Add something to the EC message stack about mismatching filetype bla bla
return
endif
if (.not.ecAddConnectionTargetItem(ecInstancePtr, connectionId, targetItemPtr1)) return
if (.not.ecAddItemConnection(ecInstancePtr, targetItemPtr1, connectionId)) return
case ('velocitybnd', 'dischargebnd', 'waterlevelbnd', 'salinitybnd', 'tracerbnd', &
'neumannbnd', 'riemannbnd', 'riemann_velocitybnd', 'outflowbnd', &
'temperaturebnd', 'sedimentbnd', 'tangentialvelocitybnd', 'uxuyadvectionvelocitybnd', &
'normalvelocitybnd', 'qhbnd')
if (.not. checkFileType(ec_filetype, provFile_poly_tim, target_name)) then
return
end if
! the file reader will have created an item called 'polytim_item'
sourceItemName = 'polytim_item'
case ('rainfall')
! the name of the source item depends on the file reader
if (ec_filetype == provFile_uniform) then
sourceItemName = 'uniform_item'
else if (ec_filetype == provFile_netcdf) then
sourceItemName = 'precipitation'
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity rainfall.')
return
end if
case ('hrms', 'tp', 'tps', 'rtp', 'dir', 'fx', 'fy', 'wsbu', 'wsbv', 'mx', 'my')
! the name of the source item created by the file reader will be the same as the ext.force. quant name
sourceItemName = target_name
case ('atmosphericpressure')
if (ec_filetype == provFile_arcinfo) then
! the arc-info file contains 'wind_p'
sourceItemName = 'wind_p'
else if (ec_filetype == provFile_netcdf) then
! the arc-info file contains 'air_pressure', which is also the standard_name
sourceItemName = 'air_pressure'
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity wind_p.')
return
end if
case ('windx')
! the name of the source item depends on the file reader
if (ec_filetype == provFile_arcinfo) then
sourceItemName = 'wind_u'
else if (ec_filetype == provFile_uniform) then
sourceItemName = 'uniform_item'
else if (ec_filetype == provFile_netcdf) then
sourceItemName = 'eastward_wind'
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity windx.')
return
end if
case ('windy')
! the name of the source item depends on the file reader
if (ec_filetype == provFile_arcinfo) then
sourceItemName = 'wind_v'
else if (ec_filetype == provFile_uniform) then
sourceItemName = 'uniform_item'
else if (ec_filetype == provFile_netcdf) then
sourceItemName = 'northward_wind'
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity windy.')
return
end if
case ('windxy')
! special case: m:n converter, (for now) handle here in case switch
if (ec_filetype == provFile_unimagdir) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'uniform_item')
success = (sourceItemId /= ec_undef_int)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId )
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_x)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_y)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_x, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_y, connectionId)
if (.not. success) then
goto 1234
end if
else if (ec_filetype == provFile_uniform) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'uniform_item')
success = (sourceItemId /= ec_undef_int)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId )
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_x)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_y)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_x, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_y, connectionId)
if (.not. success) then
goto 1234
end if
else if (ec_filetype == provFile_netcdf) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'eastward_wind')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'northward_wind')
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_2)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_x)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_windxy_y)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_x, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_windxy_y, connectionId)
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity windxy.')
return
end if
case ('airpressure_windx_windy')
! special case: m:n converter, (for now) handle seperately
if (ec_filetype == provFile_curvi) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_1')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_2')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_3')
else if (ec_filetype == provFile_spiderweb) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'windspeed')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'winddirection')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'p_drop')
else if (ec_filetype == provFile_netcdf) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'air_pressure')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'eastward_wind')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'northward_wind')
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity airpressure_windx_windy.')
return
end if
if (sourceItemId == ec_undef_int .or. sourceItemId_2 == ec_undef_int .or. sourceItemId_3 == ec_undef_int) then
goto 1234
end if
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_2)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_3)
if (ec_filetype == provFile_curvi .or. ec_filetype == provFile_netcdf) then
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_p)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_x)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_y)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_p, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_x, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_y, connectionId)
else if (ec_filetype == provFile_spiderweb) then
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_x)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_y)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_apwxwy_p)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_x, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_y, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_apwxwy_p, connectionId)
end if
if (.not. success) then
goto 1234
end if
case ('humidity_airtemperature_cloudiness')
! special case: m:n converter, (for now) handle seperately
if (ec_filetype == provFile_curvi .or. ec_filetype == provFile_uniform) then
if (ec_filetype == provFile_curvi) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_1')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_2')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_3')
if (sourceItemId == ec_undef_int .or. sourceItemId_2 == ec_undef_int .or. sourceItemId_3 == ec_undef_int) then
goto 1234
end if
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_2)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_3)
else if (ec_filetype == provFile_uniform) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'uniform_item')
if (sourceItemId == ec_undef_int) then
goto 1234
end if
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
end if
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hac_humidity)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hac_airtemperature)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hac_cloudiness)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hac_humidity, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hac_airtemperature, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hac_cloudiness, connectionId)
if (.not. success) then
goto 1234
end if
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity humidity_airtemperature_cloudiness.')
return
end if
case ('humidity_airtemperature_cloudiness_solarradiation')
! special case: m:n converter, (for now) handle seperately
if (ec_filetype == provFile_curvi .or. ec_filetype == provFile_uniform) then
if (ec_filetype == provFile_curvi) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_1')
sourceItemId_2 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_2')
sourceItemId_3 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_3')
sourceItemId_4 = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'curvi_source_item_4')
if (sourceItemId == ec_undef_int .or. sourceItemId_2 == ec_undef_int .or. &
sourceItemId_3 == ec_undef_int .or. sourceItemId_4 == ec_undef_int ) then
goto 1234
end if
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_2)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_3)
if (success) success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId_4)
else if (ec_filetype == provFile_uniform) then
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, 'uniform_item')
if (sourceItemId == ec_undef_int) then
goto 1234
end if
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
end if
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hacs_humidity)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hacs_airtemperature)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hacs_cloudiness)
if (success) success = ecAddConnectionTargetItem(ecInstancePtr, connectionId, item_hacs_solarradiation)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hacs_humidity, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hacs_airtemperature, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hacs_cloudiness, connectionId)
if (success) success = ecAddItemConnection(ecInstancePtr, item_hacs_solarradiation, connectionId)
if (.not. success) then
goto 1234
end if
else
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported filetype for quantity humidity_airtemperature_cloudiness_solarradiation.')
return
end if
case default
call mess(LEVEL_FATAL, 'm_meteo::ec_addtimespacerelation: Unsupported quantity specified in ext-file (connect source and target): '//target_name)
return
end select
if (sourceItemName /= ' ') then
! not a special case, connect source and target
sourceItemId = ecFindItemInFileReader(ecInstancePtr, fileReaderId, sourceItemName)
if (sourceItemId == ec_undef_int) then
goto 1234
end if
if (.not. initializeConnection(ecInstancePtr, connectionId, sourceItemId , targetItemPtr1)) then
goto 1234
end if
end if
ec_addtimespacerelation = .true.
return
! Error handling.
1234 continue
ec_addtimespacerelation = .false.
! message = ecGetMessage()
if (.not. quiet_) then
! TODO: AvD: I'd rather have a full message stack that will combine EC + meteo + dflowfm, and any caller may print any pending messages.
! For now: Print the EC message stack here, and leave the rest to the caller.
message = dumpECMessageStack(LEVEL_WARN, callback_msg)
! Leave this concluding message for the caller to print or not. (via getmeteoerror())
end if
message = 'm_meteo::ec_addtimespacerelation: Error while initializing '''//trim(name)//''' from file: '''//trim(filename)//''''
if (present(forcingfile)) then
message = trim(message)//' ('''//trim(forcingfile)//''')'
endif
end function ec_addtimespacerelation
! ==========================================================================
! ==========================================================================
!> Replacement function for FM's meteo1 'gettimespacevalue' function.
function ec_gettimespacevalue_by_itemID(instancePtr, itemId, timesteps, target_array) result(success)
logical :: success !< function status
type(tEcInstance), pointer :: instancePtr !< intent(in)
integer, intent(inout) :: itemID !< unique Item id
real(hp), intent(in) :: timesteps !< get data corresponding to this number of timesteps since FM's refdate
real(hp), dimension(:), target, optional, intent(inout) :: target_array !< kernel's data array for the requested values
character(len=maxMessageLen) :: messages(15)
integer :: iostat, i
!
if (itemId == ec_undef_int) then
success = .true.
return
end if
!
call clearECMessage()
success = ecGetValues(instancePtr, itemId, timesteps, target_array)
if (.not. success) then
message = dumpECMessageStack(LEVEL_WARN,callback_msg)
end if
end function ec_gettimespacevalue_by_itemID
! ==========================================================================
!> Convenience wrapper around ec_gettimespacevalue_by_itemID.
function ec_gettimespacevalue_by_name(instancePtr, group_name, timesteps) result(success)
logical :: success !< function status
type(tEcInstance), pointer :: instancePtr !< intent(in)
character(len=*), intent(in) :: group_name !< unique group name
real(hp), intent(in) :: timesteps !< get data corresponding to this number of timesteps since FM's refdate
!
success = .true.
!
if (trim(group_name) == 'humidity_airtemperature_cloudiness') then
success = ec_gettimespacevalue_by_itemID(instancePtr, item_hac_humidity , timesteps)
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_hac_airtemperature, timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_hac_cloudiness , timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
end if
if (trim(group_name) == 'humidity_airtemperature_cloudiness_solarradiation') then
success = ec_gettimespacevalue_by_itemID(instancePtr, item_hacs_humidity , timesteps)
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_hacs_airtemperature, timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_hacs_cloudiness , timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_hacs_solarradiation, timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
end if
if (trim(group_name) == 'airpressure_windx_windy') then
success = ec_gettimespacevalue_by_itemID(instancePtr, item_apwxwy_p, timesteps)
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_apwxwy_x, timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
!if (success) success = ec_gettimespacevalue_by_itemID(instancePtr, item_apwxwy_y, timesteps) ! Due to removal of up-to-date check in ec_item.f90::ecItemGetValues.
end if
end function ec_gettimespacevalue_by_name
end module m_meteo
! ==========================================================================
!>
module timespace_read
!!--copyright-------------------------------------------------------------------
! Copyright (c) 2005, WL | Delft Hydraulics. All rights reserved.
!!--disclaimer------------------------------------------------------------------
! This code is part of the Delft3D software system. WL|Delft Hydraulics has
! developed c.q. manufactured this code to its best ability and according to the
! state of the art. Nevertheless, there is no express or implied warranty as to
! this software whether tangible or intangible. In particular, there is no
! express or implied warranty as to the fitness for a particular purpose of this
! software, whether tangible or intangible. The intellectual property rights
! related to this software code remain with WL|Delft Hydraulics at all times.
! For details on the licensing agreement, we refer to the Delft3D software
! license and any modifications to this license, if applicable. These documents
! are available upon request.
!!--version information---------------------------------------------------------
! $Author$
! $Date$
! $Revision$
!!--description-----------------------------------------------------------------
!
!!--pseudo code and references--------------------------------------------------
!
! Stef.Hummel@WlDelft.nl
! Herman.Kernkamp@WlDelft.nl
! Adri.Mourits@WlDelft.nl
!
!!--declarations----------------------------------------------------------------
use precision
implicit none
integer, parameter :: maxnamelen = 256
double precision, parameter :: dmiss_default = -999.0_fp ! Default missing value in meteo arrays
double precision, parameter :: xymiss = -999.0_fp ! Default missing value in elementset
character(300), target :: errormessage = ' ' ! When an error occurs, a message is set in message.
! function getmeteoerror returns the message
double precision :: pi ! pi
double precision :: d2r ! degrees to radials
double precision :: r2d ! degrees to radials
double precision, private, parameter :: earthrad = 6378137.0_fp ! Mathworld, IUGG
contains
!
!
! ==========================================================================
!>
!> Parses an UDUnit-conventions datetime unit string.
!! TODO: replace this by calling C-API from UDUnits(-2).
function parse_ud_timeunit(timeunitstr, iunit, iyear, imonth, iday, ihour, imin, isec) result(ierr)
character(len=*), intent(in) :: timeunitstr !< Time unit by UDUnits conventions, e.g. 'seconds since 2012-01-01 00:00:00.0 +0000'.
integer, intent(out) :: iunit !< Unit in seconds, i.e. 'hours since..' has iunit=3600.
integer, intent(out) :: iyear !< Year in reference datetime.
integer, intent(out) :: imonth !< Month in reference datetime.
integer, intent(out) :: iday !< Day in reference datetime.
integer, intent(out) :: ihour !< Hour in reference datetime.
integer, intent(out) :: imin !< Minute in reference datetime.
integer, intent(out) :: isec !< Seconds in reference datetime.
integer :: ierr !< Error status, only 0 when successful.
integer :: i, n, ifound, iostat
character(len=7) :: unitstr
ierr = 0
unitstr = ' '
n = len_trim(timeunitstr)
ifound = 0
do i = 1,n
if (timeunitstr(i:i) == ' ') then ! First space found
if (timeunitstr(i+1:min(n, i+5)) == 'since') then
unitstr = timeunitstr(1:i-1)
ifound = 1
else
ierr = 1
end if
exit ! Found or error, look no further.
end if
end do
if (ifound == 1) then
select case(trim(unitstr))
case('seconds')
iunit = 1
case('minutes')
iunit = 60
case('hours')
iunit = 3600
case('days')
iunit = 86400
case('weeks')
iunit = 604800
case default
iunit = -1
end select
read (timeunitstr(i+7:n), '(I4,1H,I2,1H,I2,1H,I2,1H,I2,1H,I2)', iostat = iostat) iyear, imonth, iday, ihour, imin, isec
end if
end function parse_ud_timeunit
end module timespace_read
!
!
!
! ==========================================================================
! ==========================================================================
! ==========================================================================
!>
!> Deze module doet ruimte/tijdinterpolatie
!! Voor een gegeven quantity met ruimtedefinitie in een elementset,
!! worden de bijdragen van alle dataproviders aan die quantity gesommeerd.
!! Hierbij heeft iedere dataprovider een eigen tijd/ruimtedefinitie.
!! Zitten meerdere quantities of dezelfde tijd/ruimtedefinitie dan hoeft de tijd/ruimteinterpolatie
!! maar 1 keer uitgevoerd te worden.
!! De gevraagde grootheid moet dan niet als scalair maar als vector aangeboden worden.
module timespace_data
!!--copyright-------------------------------------------------------------------
! Copyright (c) 2007, WL | Delft Hydraulics. All rights reserved.
!!--disclaimer------------------------------------------------------------------
! This code is part of the Delft3D software system. WL|Delft Hydraulics has
! developed c.q. manufactured this code to its best ability and according to the
! state of the art. Nevertheless, there is no express or implied warranty as to
! this software whether tangible or intangible. In particular, there is no
! express or implied warranty as to the fitness for a particular purpose of this
! software, whether tangible or intangible. The intellectual property rights
! related to this software code remain with WL|Delft Hydraulics at all times.
! For details on the licensing agreement, we refer to the Delft3D software
! license and any modifications to this license, if applicable. These documents
! are available upon request.
!!--version information---------------------------------------------------------
! $Author$
! $Date$
! $Revision$
!!--description-----------------------------------------------------------------
!
!!--pseudo code and references--------------------------------------------------
!
! Stef.Hummel@WlDelft.nl
! Herman.Kernkamp@WlDelft.nl
! Adri.Mourits@WlDelft.nl
!
!!--declarations----------------------------------------------------------------
use precision
use timespace_read
use timespace_parameters
implicit none
double precision :: timelast = -1d10 ! time of most recent value requested
! if time =< timelast, no updates
double precision :: t01ini = -1d10 ! initial time for dataproviders t0 and t1 fields
! AvD: NOTE
! De pointers in alle onderstaande types worden puur gebruikt om dynamisch
! te kunnen alloceren. In Fortran 95 mag je namelijk geen allocatables in
! user-defined types opnemen. In Fortran 2003 mag dit wel, dus waarom
! binnenkort niet overstappen?
! Naar allocatables mag je ook pointeren (xyen => provider%xyen), en verder
! gebruiken we uberhaupt geen pointer(omleg-)functionaliteit. Performance
! schijnt ook slechter te zijn van pointers.
! allocables hoef je ook niet te nullifyen om de allocated check te laten
! slagen. Dit geldt wel voor de associated check van pointers.
contains
!
!
! ==========================================================================
!>
!> Read the next quantity block that is found in a file.
!! The (external forcing) file is opened elsewhere and read block-by-block
!! by consecutive calls to this routine.
subroutine readprovider(minp,qid,filename,filetype,method,operand,transformcoef,ja,smask)
! globals
integer, intent(in) :: minp !< File handle to already opened input file.
integer, intent(out) :: filetype !< File type of current quantity.
integer, intent(out) :: method !< Time-interpolation method for current quantity.
character (len=*), intent(out) :: filename !< Name of data file for current quantity.
character (len=*), intent(out) :: qid !< Identifier of current quantity (i.e., 'waterlevelbnd')
character (len=1), intent(out) :: operand !< Operand w.r.t. previous data ('O'verride or '+'Append)
double precision, intent(out) :: transformcoef(:) !< Transformation coefficients
integer, intent(out) :: ja !< Whether a block was successfully read or not.
character (len=*), intent(out), optional :: smask !< Name of mask-file applied to source arcinfo meteo-data
! locals
character (len=maxnamelen) :: rec, keywrd
integer :: l1, l2, i, jaopt, k
character (len=maxnamelen) :: generalkeywrd(25)
if (minp == 0) then
ja = 0
return
end if
keywrd = 'QUANTITY'
call zoekja(minp,rec,trim(keywrd), ja)
if (ja .eq. 1) then
l1 = index(rec,'=') + 1
call checkForSpacesInProvider(rec, l1, l2) ! l2 = l1 + #spaces after the equal-sign
read(rec(l2:),'(a)',err=990) qid
else
return
end if
keywrd = 'FILENAME'
call zoekja(minp,rec,trim(keywrd), ja)
if (ja .eq. 1) then
l1 = index(rec,'=') + 1
call checkForSpacesInProvider(rec, l1, l2) ! l2 = l1 + #spaces after the equal-sign
read(rec(l2:),'(a)',err=990) filename
else
return
end if
if (present(smask)) then ! todo: shouldn't this argument be compulsory ? .....
keywrd = 'SOURCEMASK'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) smask
else
smask = ''
end if
endif
keywrd = 'FILETYPE'
call zoekja(minp,rec,trim(keywrd), ja)
if (ja .eq. 1) then
l1 = index(rec,'=') + 1
call checkForSpacesInProvider(rec, l1, l2) ! l2 = l1 + #spaces after the equal-sign
read(rec(l2:),*, err=990) filetype
else
return
end if
keywrd = 'METHOD'
method = spaceandtime ! default : spaceandtime
call zoekja(minp,rec,trim(keywrd), ja)
if (ja .eq. 1) then
l1 = index(rec,'=') + 1
call checkForSpacesInProvider(rec, l1, l2) ! l2 = l1 + #spaces after the equal-sign
read(rec(l2:),*, err=990) method
else
return
end if
keywrd = 'OPERAND'
OPERAND = 'O' ! hk : default =O
call zoekja(minp,rec,trim(keywrd), ja)
if (ja .eq. 1) then
l1 = index(rec,'=') + 1
call checkForSpacesInProvider(rec, l1, l2) ! l2 = l1 + #spaces after the equal-sign
read(rec(l2:l2),'(a1)', err=990) operand
else
return
end if
transformcoef = -999d0
keywrd = 'VALUE'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(1)
end if
keywrd = 'FACTOR'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(2)
end if
keywrd = 'IFRCTYP'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(3)
end if
keywrd = 'AVERAGINGTYPE'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(4)
end if
keywrd = 'RELATIVESEARCHCELLSIZE'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(5)
end if
keywrd = 'EXTRAPOLTOL'
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(6)
end if
! constant keywrd = 'DISCHARGE'/'SALINITY'/'TEMPERATURE' removed, now always via time series, in future also via new ext [discharge]
keywrd = 'AREA' ! Area for source-sink pipe
call zoekopt(minp, rec, trim(keywrd), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(4)
end if
generalkeywrd( 1) = 'widthleftW1' !< ! generalstructure this and following: see Sobek manual
generalkeywrd( 2) = 'levelleftZb1'
generalkeywrd( 3) = 'widthleftWsdl'
generalkeywrd( 4) = 'levelleftZbsl'
generalkeywrd( 5) = 'widthcenter'
generalkeywrd( 6) = 'levelcenter'
generalkeywrd( 7) = 'widthrightWsdr'
generalkeywrd( 8) = 'levelrightZbsr'
generalkeywrd( 9) = 'widthrightW2'
generalkeywrd(10) = 'levelrightZb2'
generalkeywrd(11) = 'gateheight'
generalkeywrd(12) = 'gateheightintervalcntrl'
generalkeywrd(13) = 'pos_freegateflowcoeff'
generalkeywrd(14) = 'pos_drowngateflowcoeff'
generalkeywrd(15) = 'pos_freeweirflowcoeff'
generalkeywrd(16) = 'pos_drownweirflowcoeff'
generalkeywrd(17) = 'pos_contrcoeffreegate'
generalkeywrd(18) = 'neg_freegateflowcoeff'
generalkeywrd(19) = 'neg_drowngateflowcoeff'
generalkeywrd(20) = 'neg_freeweirflowcoeff'
generalkeywrd(21) = 'neg_drownweirflowcoeff'
generalkeywrd(22) = 'neg_contrcoeffreegate'
generalkeywrd(23) = 'extraresistance'
generalkeywrd(24) = 'dynstructext'
generalkeywrd(25) = 'gatedoorheight'
do k = 1,25 ! generalstructure
call zoekopt(minp, rec, trim(generalkeywrd(k)), jaopt)
if (jaopt == 1) then
read (rec,*) transformcoef(k)
end if
enddo
return
990 call readerror('reading '//trim(keywrd)//' but getting ', rec, minp)
end subroutine readprovider
!
!
! ==========================================================================
!>
subroutine checkForSpacesInProvider(rec, eqsign, eqsignsp)
! I/O
character (len=256), intent(in) :: rec !< Name of record that includes the keyword and record
integer, intent(in) :: eqsign !< Location of the equal-sign in the entire record string
integer, intent(out) :: eqsignsp !< Location of the equal-sign plus first spaces after equal-sign
! Locals
integer :: i ! Counter
eqsignsp = eqsign
do i = 0, 256-eqsign
if (rec(eqsign+i:eqsign+i) .eq. ' ') then
eqsignsp = eqsignsp + 1
else
exit
end if
enddo
end subroutine checkForSpacesInProvider
!
!
! ==========================================================================
!>
subroutine read1polylin(minp,xs,ys,ns)
integer :: minp
double precision :: xs(:)
double precision :: ys(:)
integer :: ns
character (len=maxnamelen) :: rec
integer :: k
ns = 0
10 read(minp,'(a)',end = 999) rec
if (rec(1:1) == '*' ) goto 10
read(minp,'(a)',end = 999) rec
read(rec ,* ,err = 888) ns
do k = 1,ns
read(minp,'(a)',end = 999) rec
read(rec ,* ,err = 777) xs(k), ys(k)
enddo
call doclose(minp)
return
999 call eoferror(minp)
return
888 call readerror('reading nrows but getting ', rec, minp)
return
777 call readerror('reading x, y but getting ', rec, minp)
return
end subroutine read1polylin
!
!
! ==========================================================================
!>
subroutine settimespacerefdat(refda, jul00, tz, timjan)
use m_itdate
character (len=8) :: refda
integer :: jul00
double precision :: tz, timjan
integer, external :: julday
integer :: juljan
refdat = refda
read (refdat,*) itdate
read(refdat(1:4),*) iyear0
read(refdat(5:6),*) imonth0
read(refdat(7:8),*) iday0
jul0 = julday(imonth0,iday0,iyear0)
jul00 = jul0
Tzone = tz
juljan = julday(1,1,iyear0)
timjan = (jul0 - juljan)*24.d0
end subroutine settimespacerefdat
!
!
! ==========================================================================
!>
function getmeteoerror( ) result(retval)
implicit none
character(300), pointer :: retval
retval => errormessage
end function getmeteoerror
!
!
! ==========================================================================
!>
subroutine meteo_tidepotential(jul0, TIME , xz , yz , TIDEP, ndx, dstart, dstop , eps) ! call schrama's routines on reduced set
use m_sferic
integer :: jul0, ndx ! interpolate results in ndx
double precision :: time, dstart, dstop , eps
double precision :: xz(ndx), yz(ndx), tidep(ndx), xx(4), yy(4), DAREA, DLENGTH, DLENMX
double precision, allocatable, save :: xz2(:,:), yz2(:,:), td2(:,:), self(:,:), avhs(:,:), area(:,:)
double precision :: xmn, xmx, ymn, ymx, di, dj, f11,f21,f12,f22
integer, save :: ndx2, INI = 0, jaself = 1
integer :: i,j,n,ierr, m1,m2,n1,n2
integer, save :: i1,i2,j1,j2
if (INI == 0 ) then
INI = 1
XMN = 1D30 ; YMN = 1D30 ; XMX = -1D30 ; YMX = -1D30
DO I = 1,ndx
xmn = min(xz(i),xmn)
xmx = max(xz(i),xmx)
ymn = min(yz(i),ymn)
ymx = max(yz(i),ymx)
ENDDO
i1 = floor(xmn); i2 = floor(xmx) + 1
j1 = floor(ymn); j2 = floor(ymx) + 1
IF (ALLOCATED (XZ2) ) DEALLOCATE (XZ2,YZ2,TD2)
allocate ( xz2(i1:i2,j1:j2), stat=ierr) ! tot aerr
allocate ( yz2(i1:i2,j1:j2), stat=ierr)
allocate ( td2(i1:i2,j1:j2), stat=ierr)
if (jaself == 1) then
if (allocated(self) ) deallocate ( self, avhs, area )
allocate ( self(i1:i2,j1:j2), stat=ierr)
allocate (avhs(i1:i2,j1:j2), stat=ierr)
allocate ( area(i1:i2,j1:j2), stat=ierr)
do i = i1,i2
do j = j1,j2
xx(1) = dble(i)-0.5d0 ; yy(1) = dble(j)-0.5d0
xx(2) = dble(i)+0.5d0 ; yy(2) = dble(j)-0.5d0
xx(3) = dble(i)+0.5d0 ; yy(3) = dble(j)+0.5d0
xx(4) = dble(i)-0.5d0 ; yy(4) = dble(j)+0.5d0
call dAREAN( XX, YY, 4, DAREA, DLENGTH, DLENMX )
area(i,j) = darea
enddo
enddo
endif
do i = i1,i2
do j = j1,j2
xz2(i,j) = i*dg2rd
yz2(i,j) = j*dg2rd
enddo
enddo
ndx2 = (i2-i1+1)*( j2-j1+1)
end if
call tforce( jul0, TIME , xz2 , yz2 , Td2, ndx2, dstart, dstop , eps)
if (jaself == 1) then
call aggregatewaterlevels(avhs, area, i1,i2,j1,j2 )
call selfattraction(avhs, self, i1,i2,j1,j2 )
endif
do n = 1,ndx
m1 = floor(xz(n)) ; m2 = m1+1
n1 = floor(yz(n)) ; n2 = n1+1
di = xz(n) - m1
dj = yz(n) - n1
f11 = (1d0-di)*(1d0-dj)
f21 = ( di)*(1d0-dj)
f22 = ( di)*( dj)
f12 = (1d0-di)*( dj)
if (jaself == 1) then
tidep(n) = ( td2(m1,n1) + self(m1,n1) )*f11 + &
( td2(m2,n1) + self(m1,n1) )*f21 + &
( td2(m2,n2) + self(m1,n1) )*f22 + &
( td2(m1,n2) + self(m1,n1) )*f12
else
tidep(n) = td2(m1,n1) *f11 + &
td2(m2,n1) *f21 + &
td2(m2,n2) *f22 + &
td2(m1,n2) *f12
endif
enddo
end subroutine meteo_tidepotential
subroutine selfattraction(avhs, self, i1, i2, j1, j2)
implicit none
integer :: i1,i2,j1,j2
double precision :: avhs(i1:i2,j1:j2), self(i1:i2,j1:j2)
integer :: i,j
self = 0d0
do j = j1,j2
do i = i1,i2
enddo
enddo
end subroutine selfattraction
subroutine aggregatewaterlevels(avhs, area, i1,i2,j1,j2 )
use m_flow
use m_flowgeom
implicit none
integer :: i1,i2,j1,j2, k,i,j
double precision :: avhs(i1:i2,j1:j2), area(i1:i2,j1:j2)
avhs = 0d0
do k = 1,ndx
i = nint(xz(k))
j = nint(yz(k))
avhs(i,j) = avhs(i,j) + vol1(k)
enddo
do j = j1,j2
do i = i1,i2
avhs(i,j) = avhs(i,j) / area(i,j)
enddo
enddo
end subroutine aggregatewaterlevels
!
!
! ==========================================================================
!>
subroutine tforce( jul0, TIME , xzeta , yzeta , TIDEP, IDIM1, dstart, dstop , eps)
!
! ====================================================================
!
! Programmer E.J.O. Schrama
!
! Original URL: https://repos.deltares.nl/repos/simona/bo_omgeving/simona/src/waqua/waqpro/routines/wastfr.f
! $Revision: 1850 $, $Date: 2008-04-18 09:19:37 +0200 (Fri, 18 Apr 2008) $
!
! Version 1.1 Date 22-05-2008 c81402: extended for evaluation of
! tidal forces on grids (AVe,
! VORtech)
! Version 1.0 Date 24-01-2008 initial version
!
! Copyright (c) "E.J.O. Schrama".
! Permission to copy or distribute this software or documentation
! in hard copy or soft copy granted only by written license
! obtained from "Rijkswaterstaat".
! All rights reserved. No part of this publication may be
! reproduced, stored in a retrieval system (e.g., in memory, disk,
! or core) or be transmitted by any means, electronic, mechanical,
! photocopy, recording, or otherwise, without written permission
! from the publisher.
!
! ********************************************************************
!
! DESCRIPTION
!
! Computes the tidal potential for each active grid point
!
! ********************************************************************
!
! COMMON BLOCKS
!
implicit none
!
! **********************************************************************
!
! INPUT / OUTPUT PARAMETERS
!
integer idim1, luhar, jul0
double precision :: rmjdat, dstart, dstop, eps, TIME
double precision xzeta(idim1), yzeta(idim1), tidep(idim1)
DOUBLE PRECISION, allocatable, SAVE :: tideuc(:,:,:), tideus(:,:,:) ! (idim1, 0:3,2:3),
INTEGER, SAVE :: IRC = 0
character*256 harfil
CHARACTER(len=40), dimension(484) :: RECS ! ZAT IN FILE 'HARMONICS'
!
! dstart i starting Doodson number
! dstop i stopping Doodson number
! eps i tolerance level for tidal force formula
! harfil i file with tidal harmonics
! idim1 i first dimension of fullbox array (nmax for SIMONA)
! i second dimension of fullbox array (mmax+6 for SIMONA)
! irc i input parameter for tforce
! irc = 0: initialisation phase
! irc = 1: simulation phase
! luhar i logical unit number to read file with tidal harmonics
! name i character string containing the name of the
! the calling subroutine.
! Only used for error messages.
! rmjdat i modified julian day (24-jan-2008 0:00 UTC : 54489.00000)
! tidep o tidal potential
! tideuc i o cosine component of tidal potential (array)
! tideus i o sine component of tidal potential (array)
! xzeta i latitude (in radians) of grid-points in physiscal plane
! yzeta i longitude (in radians) of grid-points in physiscal plane
!
! ********************************************************************
!
! LOCAL PARAMETERS
!
integer maxdat, maxfld, idebug, i1, i1dbg, i2dbg, N
double precision :: pi, g, rmu, re, d2r, reps
parameter (idebug=0, i1dbg=0, i2dbg=0)
parameter (maxdat=500) ! maximal # records in table
parameter (maxfld=7) ! maximal # fields in table
parameter (pi = 3.14159265358979, re = 6378137d0, &
d2r = pi/180d0, rmu = 3.9860044d14, &
g = rmu / re / re, reps = 1d-5 )
integer ntable, nskip
integer itable(maxdat,maxfld)
double precision :: amps(maxdat), plsmin(6), rklove(3), rhlove(3), &
factor(2:3), pol1(0:3,2:3), cm1(0:3), sm1(0:3)
integer i,j,nr,nq,mq, IERR
integer kk(10 )
double precision :: fnm, pnm, har, argum, argfct, dtab1, dtab2, &
dtab, rlslat, rlslon, rlat, rlong, potent
double precision :: elmnts(6), can(maxdat), san(maxdat)
double precision :: cansum(0:3,2:3), sansum(0:3,2:3)
character*80 record
logical permnt, first
DOUBLE PRECISION, SAVE :: FACTORIAL(0:6)
!
! amps table with scaled amplitudes for selected tidal components
! argfct multiplication factor needed to compute argument
! argum argument for time-dependent harmonic components ca,sa
! can table with scaled harmonic components
! dcos(argument) * amp(i)
! cansum selected sum of elements of can for fixed mq,nq
! cm1 cosine-component of potential
! d2r conversion factor pi/180
! dtab Doodson number: dtab1 + dtab2 / 1000d0
! dtab1 first 3 digits of Doodson number
! dtab2 second 3 digits of Doodson number
! elmnts array needed for calculation of can, san
! factor table with multiplication factors needed to compute
! tidal potential
! fnm global normalization factor (fnm) found in cartwright 1993
! needed for pol1
! g gravity acceleration
! har amplitude of tidal component
! i loop variable
! i1 loop variable
! idebug flag wheter debug output is printed
! itable CTE table read from harmonics.table with the selected CTE
! lines
! j loop variable
! kk array of fields for CTE table
! maxdat maximum number of CTE lines
! maxfld maximum number of fields for each CTE line
! mq order of cm1, sm1 and pol1 (table(i,1))
! nq degree of Legendre polynomial (table(i,7))
! nr number of CTE line
! nskip number of skipped CTE lines in harmonics.table
! ntable number of selected CTE lines
! permnt flag whether tidal line should be skipped
! pi 3.14159265358979
! plsmin array with +/- signs for all fields
! pnm output value of legpol1
! pol1 Legendre polynomial array
! potent tidal potential
! re radius of the earth
! record character array
! reps real value 1d-5 used for inequalities
! rhlove Love numbers describing the indirect potential effects
! rklove Love numbers describing the geometric-radial effects
! rlat northern latitude in radians
! rlong eastern longitude in radians
! rlslat previous value of rlat
! rlslon previous value of rlong
! rmu gravitational constant (3.9860044d14)
! san table with scaled harmonic components
! dsin(argument) * amp(i)
! sansum selected sum of elements of san for fixed mq,nq
! sm1 sine-component of potential
!
! ********************************************************************
!
! I / O
!
! harfil = WAQPRO**/harmonics.table I CTE table (unit luhar)
!
! ********************************************************************
!
! SUBROUTINES CALLED
!
!
! ASTROL This copied from richard's subroutine astrol, in goes the
! modified Julian date, out comes an array of six double
! precision variables used for Doodson number computations
! LEGPOL1 compute unnormalized associated legendre polynomials up to
! degree 3 and order 3
!
! ********************************************************************
!
! ERROR MESSAGES
!
! ********************************************************************
!
! PSEUDO CODE
!
! save static variables
!
! if (first call) then
! initialise
! read CTE-table (table, amps, ntable)
! if (Dstart <= d(i) <= Dstop and H(i) >= eps) then
! read tidal line
! amps(i) = har * g * factor
! end if
!
! computation and storage of tideuc,tideus at each grid-point:
! do (all grid-points)
! do nq = 2, 3
! do mq = 0, nq
! if (changed latitude) then
! update pol1(m,n) (call legpol1)
! end if
! if (changed longitude) then
! update cm1(m) = cos(), sm1(m) = sin()
! end if
! compute tideuc(nm,mq,nq), tideus(nm,mq,nq)
! store tideuc,tideus
! enddo
! enddo
! enddo
! end
!
! update elements-table (call astrol)
! compute arrays can, san:
! do (all tidal components)
! compute argum
! can(i) = dcos(argum) * amps(i)
! san(i) = dsin(argum) * amps(i)
! enddo
! compute arrays cansum, sansum:
! do nq = 2, 3
! do mq = 0, nq
! do (all tidal components)
! if ( itable(i,7).eq.nq .and. itable(i,1).eq.mq) then
! cansum(mq, nq) = cansum(mq, nq) + can(i)
! sansum(mq, nq) = sansum(mq, nq) + san(i)
! end if
! enddo
! enddo
! enddo
!
! computation of the tidal potential at each grid-point:
! do (all grid-points)
! do nq = 2, 3
! do mq = 0, nq
! potent = potent + tideuc(nm, mq, nq) * cansum(mq,nq) +
! tideus(nm, mq, nq) * sansum(mq,nq)
! enddo
! enddo
! tidep(n,m) = potent
! enddo
!
! ********************************************************************
!
! DATA STATEMENTS
!
data plsmin / +1d0, +1d0, +1d0, +1d0, -1d0, +1d0 /
!
save itable, ntable, pol1, cm1, sm1, amps
!
!c ====================================================================
!
! --- initialisation phase
!
! --- compute modified Julian date (rmjdat)
! rmjdat = itdjul - 2400001.0 + time / (24 * 60)
rmjdat = jul0 - 2400001.0 + time / (24 * 60)
if (irc.eq.0 ) then
IRC = 1
FACTORIAL(0) = 1D0
FACTORIAL(1) = 1D0
FACTORIAL(2) = 2D0
FACTORIAL(3) = 6D0
FACTORIAL(4) = 24D0
FACTORIAL(5) = 120D0
FACTORIAL(6) = 720D0
if (allocated (tideuc) ) deallocate(tideuc, tideus)
ALLOCATE ( tideuc(0:3,2:3,IDIM1), STAT = IERR); tideuc = 0d0
ALLOCATE ( tideus(0:3,2:3,IDIM1), STAT = IERR); tideus = 0d0
call iniharmonics(recs)
if (idebug.ge.10) then
write(6,*) '*** Start reading the harmonics ' // &
'table ***'
write(6,*)
end if
!
! --- k and h love numbers for degree 2 and 3
!
rklove(1) = 0d0
rklove(2) = 0.303d0
rklove(3) = 0.0937d0
rhlove(1) = 0d0
rhlove(2) = 0.612d0
rhlove(3) = 0.293d0
!
do nq = 2, 3
factor(nq) = (1d0 + rklove(nq) - rhlove(nq))
enddo
!
ntable = 0
nskip = 0
N = 0
10 continue
N = N + 1
IF (N > 484) GOTO 20
RECORD = RECS(N)
! read(luhar,'(a)',end=20) record
if (idebug.ge.10) write(6,*) record
if (record(1:1).eq.'%') go to 10
read(record,*) (kk(i),i = 1, 7), har
!
! in the CTE tables there is a null line for theoretic
! --- reasons, never use this line to compute the tidal
! potential.
!
permnt = kk(1).eq.0 .and. kk(2).eq.0 .and. &
kk(3).eq.0 .and. kk(4).eq.0 .and. &
kk(5).eq.0 .and. kk(6).eq.0
!
! --- dtab is the doodson number for a table entry,
! select the lines where dstart <= dtab <= dstop
!
dtab1 = kk(1)*100d0 + (kk(2)+5d0)*10d0 + (kk(3)+5d0)
dtab2 = (kk(4)+5d0)*100d0 + (kk(5)+5d0)*10d0 + (kk(6)+5d0)
dtab = dtab1 + dtab2 / 1000d0
if (.not.permnt .and. abs(har).ge.eps .and. dstart.le.dtab .and. dtab.le.dstop ) then
ntable = ntable + 1
if (ntable.gt.maxdat) then
irc = -2
return
end if
do i = 1, 7
itable(ntable,i) = kk(i)
end do
nq = kk(7)
amps(ntable) = har * g * factor(nq)
else
nskip = nskip + 1
end if
go to 10
20 continue
! rewind(luhar)
rlslat = -9999d0
rlslon = -9999d0
if (idebug.ge.10) then
write(6,*) 'ntable = ',ntable,' nskip = ',nskip
end if
!
! --- storage of uc, us
!
do 90 i1 = 1, idim1
rlat = yzeta(i1)
rlong = xzeta(i1)
! compute legendre polynomials times the global
! --- normalization factor
! (fnm) found in cartwright 1993, also compute the
! astronomical elements
!
if (abs(rlat-rlslat).gt.reps) then
do nq = 2, 3
do mq = 0, nq
fnm = 2d0 / dble(2*nq+1) * factorial(nq+mq) / factorial(nq-mq)
fnm = sqrt(1d0 / (2d0 * pi * fnm)) * ((-1d0)**mq)
call legpol1( rlat, nq, mq, pnm )
pol1(mq, nq) = fnm * pnm
enddo
enddo
end if
if (abs(rlong-rlslon).gt.reps ) then
do mq = 0, 3
cm1(mq) = +cos(dble(mq) * rlong)
sm1(mq) = +sin(dble(mq) * rlong)
enddo
end if
!
! --- compute arrays with tideuc, tideus
!
do nq = 2, 3
do mq = 0, nq
if (mod(nq+mq, 2).eq.0) then
tideuc(mq, nq, I1) = +cm1(mq) * pol1(mq, nq)
tideus(mq, nq, I1) = -sm1(mq) * pol1(mq, nq)
else
tideuc(mq, nq, I1) = +sm1(mq) * pol1(mq, nq)
tideus(mq, nq, I1) = +cm1(mq) * pol1(mq, nq)
end if
if (idebug.ge.2 .and. i1.eq.i1dbg) then
write(6,123) 'mq,nq=',mq,nq, &
': rlat=',rlat,', rlong=',rlong, &
', pol1=',pol1(mq,nq),', cm1=',cm1(mq),',', &
'sm1=',sm1(mq), &
', tideuc=',tideuc(mq,nq,I1), &
', tideus=',tideus(mq,nq,I1)
123 format(a,2i2,4(a,f8.4),a,/,12x,3(a,f8.4))
end if
enddo
enddo
!
! --- and finally save the last known latitude and
! longitude
!
rlslat = rlat
rlslon = rlong
90 continue
end if
!
! --- end of initialisation phase
!
! ------------------------------------------------------------------------
!
! --- update elements-table
!
call astrol( rmjdat, elmnts )
!
! --- compute tabels can, san
!
do i = 1, ntable
argum = 0d0
do j = 1, 6
argfct = dble(itable(i, j))
argum = argum + argfct * elmnts(j) * plsmin(j)
end do
! argum = mod(argum, 360d0)
! if (argum.lt.0d0) argum = argum + 360d0
argum = argum * d2r
can(i) = cos(argum) * amps(i)
san(i) = sin(argum) * amps(i)
enddo
!
! --- compute tables cansum, sansum
!
do nq = 2, 3
do mq = 0, nq
cansum(mq, nq) = 0d0
sansum(mq, nq) = 0d0
do i = 1, ntable
if ( itable(i,7).eq.nq .and. itable(i,1).eq.mq) then
cansum(mq, nq) = cansum(mq, nq) + can(i)
sansum(mq, nq) = sansum(mq, nq) + san(i)
end if
enddo
enddo
enddo
!
! --- computation of the tidal potential at each grid-point:
!
do 190 i1 = 1, idim1
potent = 0d0
do nq = 2, 3
do mq = 0, nq
potent = potent + tideuc(mq, nq, I1) * cansum(mq,nq)
potent = potent + tideus(mq, nq, I1) * sansum(mq,nq)
! if (idebug.ge.2 .and. i1.eq.i1dbg) then
! write(6,'(a,2i2,2(2(a,f10.6)))') 'mq,nq=',mq,nq, &
! ': term uc=',tideuc(i1,mq,nq),'*', &
! cansum(mq,nq), &
! ', us=',tideus(i1,mq,nq),'*', &
! sansum(mq,nq)
!end if
enddo
enddo
tidep(i1) = potent
190 continue
! if (idebug.ge.1 .and. i1dbg.ge.1) write(6,*) 'tidep=', tidep(i1dbg)
end subroutine tforce
!
!
! ==========================================================================
!>
subroutine astrol( mjdate,six )
! ====================================================================
!
! Programmer R. D. Ray
!
! Version 1.0 Date dec. 1990 initial version
!
! ********************************************************************
!
! DESCRIPTION
!
! This copied from richard's subroutine astrol, in goes the
! modified Julian date, out comes an array of six double precision
! variables used for Doodson number computations
!
! Computes the basic astronomical mean longitudes s, h, p, N.
! Note N is not N', i.e. N is decreasing with time.
! These formulae are for the period 1990 - 2010, and were derived
! by David Cartwright (personal comm., Nov. 1990).
! TIME is UTC in decimal MJD.
! All longitudes returned in degrees.
! R. D. Ray Dec. 1990
!
! Non-vectorized version.
!
! ********************************************************************
!
! COMMON BLOCKS
!
implicit none
!
! ********************************************************************
!
! INPUT / OUTPUT PARAMETERS
!
double precision :: six(6),mjdate
!
! mjdate i modified julian day (24-jan-2008 0:00 UTC : 54489.00000)
! six o array of six double precision variables used for Doodson
! number computations
! see also Cartwright 1993, summer school lecture notes,
! page 108
! six(1) (tau) mean time angle in lunar days
! six(2) (q) mean longitude of moon
! six(3) (q') mean longitude of sun
! six(4) (p) mean longitude of lunar perigee
! six(5) (N) mean longitude of ascending lunar node
! six(6) (p') mean longitude of the Sun at perihelion
!
! ********************************************************************
!
! LOCAL PARAMETERS
!
! --- constant values:
!
double precision :: circle
parameter (CIRCLE=360.0D0)
!
! circle number of degrees in a circle
!
! --- variables:
!
double precision :: T,TIME,UT
integer i
!
! T translated time: TIME - 51544.4993D0
! TIME input time (mjdate)
! UT fractional part of mjdate: (mjdate - int(mjdate))
!
!=======================================================================
!
! --- start of code
!
TIME = mjdate
T = TIME - 51544.4993D0 ! reference to 2000/1/1 1200 o'clock
!
! --- perform translations using translation table of symbols:
!
! nr Cartwright Doodson Brown
! 1 tau tau 360*t-D+180
! 2 q s L
! 3 q' h L'
! 4 p p \overline(\omega)
! 5 N -N' \Omega
! 6 p' p_1 \overline{\omega}'
!
six(2) = 218.3164D0 + 13.17639648D0 * T
six(3) = 280.4661D0 + 0.98564736D0 * T
six(4) = 83.3535D0 + 0.11140353D0 * T
six(5) = 125.0445D0 - 0.05295377D0 * T
six(6) = 282.9384D0 + 0.0000471d0 * T
!
! --- get them in the right quadrant
!
do i = 2, 6
six(i) = mod(six(i), circle)
if (six(i).lt.0d0) six(i) = six(i) + circle
end do
!
! argument 1 in the doodson number denotes the mean lunar time.
! According to equation 13 and the inline remark after equation 14
! it is computed by
! alpha_G = 360 * "Universal time in fractions of a day" + q'(T) - 180
! tau = alpha_G - q
!
UT = (mjdate - int(mjdate))
six(1) = 360d0 * UT + six(3) - 180d0 - six(2)
end subroutine astrol
!
!
! ==========================================================================
!>
subroutine legpol1( theta,n,m,pnm )
! ====================================================================
!
! Programmer E. Schrama
!
! ********************************************************************
!
! DESCRIPTION
!
! A routine to compute unnormalized associated legendre polynomials
! up to degree 3 and order 3.
!
! ********************************************************************
!
! COMMON BLOCKS
!
implicit none
!
! ********************************************************************
!
! INPUT / OUTPUT PARAMETERS
!
integer n,m
double precision :: theta,pnm
!
! m i degree of Legendre polynomial
! n i order of Legendre polynomial
! pnm o value of Legendre polynomial
! theta i phase
!
! ********************************************************************
!
! LOCAL PARAMETERS
!
double precision :: cp,sp
!
! cp cos(theta)
! sp sin(theta)
!
!=======================================================================
!
pnm = 1d38
cp = cos( theta )
sp = sin( theta )
!
! --- I think this comes from (Lambeck,1988), what again are the rules for
! obtaining associated Legendre functions?
!
if (n.eq.0 ) then
if (m.eq.0 ) pnm = 1d0
ELSE if (n.eq.1 ) then
if (m.eq.0 ) pnm = sp
if (m.eq.1 ) pnm = cp
ELSE if (n.eq.2 ) then
if (m.eq.0 ) pnm = 1.5d0*sp*sp - 0.5d0
if (m.eq.1 ) pnm = 3.0d0*sp*cp
if (m.eq.2 ) pnm = 3.0d0*cp*cp
ELSE IF (n.eq.3 ) then
if (m.eq.0 ) pnm = 2.5d0*sp*sp*sp - 1.5d0*sp
if (m.eq.1 ) pnm = cp * (7.5d0*sp*sp - 1.5d0)
if (m.eq.2 ) pnm = 15d0*cp*cp*sp
if (m.eq.3 ) pnm = 15d0*cp*cp*cp
end if
end subroutine legpol1
!
!
! ==========================================================================
!>
SUBROUTINE INIHARMONICS(RECS)
CHARACTER(len=40), dimension(484) :: RECS ! ZAT IN FILE 'HARMONICS'
!%refsys 2000
!%mjd0 47893.00000000
!%mjd1 55196.00000000
!%dmjd .11410938
!%ndata 64000
!%gmearth 3.9860044d14
!%reearth 6378137
RECS( 1) = ' 0 0 0 0 0 0 2 -.31459'
RECS( 2) = ' 0 0 0 0 1 0 2 .02793'
RECS( 3) = ' 0 0 0 0 2 0 2 -.00028'
RECS( 4) = ' 0 0 0 2 1 0 2 .00004'
RECS( 5) = ' 0 0 1 0 -1 -1 2 -.00004'
RECS( 6) = ' 0 0 1 0 0 -1 2 -.00492'
RECS( 7) = ' 0 0 1 0 0 1 2 .00026'
RECS( 8) = ' 0 0 1 0 1 -1 2 .00004'
RECS( 9) = ' 0 0 2 -2 -1 0 2 .00002'
RECS( 10) = ' 0 0 2 -2 0 0 2 -.00031'
RECS( 11) = ' 0 0 2 0 0 0 2 -.03099'
RECS( 12) = ' 0 0 2 0 0 -2 2 -.00012'
RECS( 13) = ' 0 0 2 0 1 0 2 .00076'
RECS( 14) = ' 0 0 2 0 2 0 2 .00017'
RECS( 15) = ' 0 0 3 0 0 -1 2 -.00181'
RECS( 16) = ' 0 0 3 0 1 -1 2 .00003'
RECS( 17) = ' 0 0 4 0 0 -2 2 -.00007'
RECS( 18) = ' 0 1 -3 1 -1 1 2 .00002'
RECS( 19) = ' 0 1 -3 1 0 1 2 -.00029'
RECS( 20) = ' 0 1 -2 -1 -2 0 2 .00002'
RECS( 21) = ' 0 1 -2 -1 -1 0 2 .00006'
RECS( 22) = ' 0 1 -2 1 -1 0 2 .00048'
RECS( 23) = ' 0 1 -2 1 0 0 2 -.00673'
RECS( 24) = ' 0 1 -2 1 1 0 2 .00044'
RECS( 25) = ' 0 1 -1 -1 0 1 2 -.00021'
RECS( 26) = ' 0 1 -1 0 0 0 2 .00019'
RECS( 27) = ' 0 1 -1 1 0 -1 2 .00005'
RECS( 28) = ' 0 1 0 -1 -2 0 2 -.00003'
RECS( 29) = ' 0 1 0 -1 -1 0 2 .00231'
RECS( 30) = ' 0 1 0 -1 0 0 2 -.03518'
RECS( 31) = ' 0 1 0 -1 1 0 2 .00228'
RECS( 32) = ' 0 1 0 1 0 0 2 .00188'
RECS( 33) = ' 0 1 0 1 1 0 2 .00077'
RECS( 34) = ' 0 1 0 1 2 0 2 .00021'
RECS( 35) = ' 0 1 1 -1 0 -1 2 .00018'
RECS( 36) = ' 0 1 2 -1 0 0 2 .00049'
RECS( 37) = ' 0 1 2 -1 1 0 2 .00024'
RECS( 38) = ' 0 1 2 -1 2 0 2 .00004'
RECS( 39) = ' 0 1 3 -1 0 -1 2 .00002'
RECS( 40) = ' 0 2 -4 2 0 0 2 -.00011'
RECS( 41) = ' 0 2 -3 0 0 1 2 -.00039'
RECS( 42) = ' 0 2 -3 0 1 1 2 .00002'
RECS( 43) = ' 0 2 -2 0 -1 0 2 -.00042'
RECS( 44) = ' 0 2 -2 0 0 0 2 -.00584'
RECS( 45) = ' 0 2 -2 0 1 0 2 .00037'
RECS( 46) = ' 0 2 -2 2 0 0 2 .00004'
RECS( 47) = ' 0 2 -1 -2 0 1 2 -.00004'
RECS( 48) = ' 0 2 -1 -1 0 0 2 .00003'
RECS( 49) = ' 0 2 -1 0 0 -1 2 .00007'
RECS( 50) = ' 0 2 -1 0 0 1 2 -.00020'
RECS( 51) = ' 0 2 -1 0 1 1 2 -.00004'
RECS( 52) = ' 0 2 0 -2 -1 0 2 .00015'
RECS( 53) = ' 0 2 0 -2 0 0 2 -.00288'
RECS( 54) = ' 0 2 0 -2 1 0 2 .00019'
RECS( 55) = ' 0 2 0 0 0 0 2 -.06660'
RECS( 56) = ' 0 2 0 0 1 0 2 -.02761'
RECS( 57) = ' 0 2 0 0 2 0 2 -.00258'
RECS( 58) = ' 0 2 0 0 3 0 2 .00006'
RECS( 59) = ' 0 2 1 -2 0 -1 2 .00003'
RECS( 60) = ' 0 2 1 0 0 -1 2 .00023'
RECS( 61) = ' 0 2 1 0 1 -1 2 .00006'
RECS( 62) = ' 0 2 2 -2 0 0 2 .00020'
RECS( 63) = ' 0 2 2 -2 1 0 2 .00008'
RECS( 64) = ' 0 2 2 0 2 0 2 .00003'
RECS( 65) = ' 0 3 -5 1 0 1 2 -.00002'
RECS( 66) = ' 0 3 -4 1 0 0 2 -.00018'
RECS( 67) = ' 0 3 -3 -1 0 1 2 -.00007'
RECS( 68) = ' 0 3 -3 1 0 1 2 -.00011'
RECS( 69) = ' 0 3 -3 1 1 1 2 -.00005'
RECS( 70) = ' 0 3 -2 -1 -1 0 2 -.00009'
RECS( 71) = ' 0 3 -2 -1 0 0 2 -.00092'
RECS( 72) = ' 0 3 -2 -1 1 0 2 .00006'
RECS( 73) = ' 0 3 -2 1 0 0 2 -.00242'
RECS( 74) = ' 0 3 -2 1 1 0 2 -.00100'
RECS( 75) = ' 0 3 -2 1 2 0 2 -.00009'
RECS( 76) = ' 0 3 -1 -1 0 1 2 -.00013'
RECS( 77) = ' 0 3 -1 -1 1 1 2 -.00004'
RECS( 78) = ' 0 3 -1 0 0 0 2 .00007'
RECS( 79) = ' 0 3 -1 0 1 0 2 .00003'
RECS( 80) = ' 0 3 -1 1 0 -1 2 .00003'
RECS( 81) = ' 0 3 0 -3 0 0 2 -.00023'
RECS( 82) = ' 0 3 0 -3 1 -1 2 .00003'
RECS( 83) = ' 0 3 0 -3 1 1 2 .00003'
RECS( 84) = ' 0 3 0 -1 0 0 2 -.01275'
RECS( 85) = ' 0 3 0 -1 1 0 2 -.00529'
RECS( 86) = ' 0 3 0 -1 2 0 2 -.00050'
RECS( 87) = ' 0 3 0 1 2 0 2 .00005'
RECS( 88) = ' 0 3 0 1 3 0 2 .00002'
RECS( 89) = ' 0 3 1 -1 0 -1 2 .00011'
RECS( 90) = ' 0 3 1 -1 1 -1 2 .00004'
RECS( 91) = ' 0 4 -4 0 0 0 2 -.00008'
RECS( 92) = ' 0 4 -4 2 0 0 2 -.00006'
RECS( 93) = ' 0 4 -4 2 1 0 2 -.00003'
RECS( 94) = ' 0 4 -3 0 0 1 2 -.00014'
RECS( 95) = ' 0 4 -3 0 1 1 2 -.00006'
RECS( 96) = ' 0 4 -2 -2 0 0 2 -.00011'
RECS( 97) = ' 0 4 -2 0 0 0 2 -.00204'
RECS( 98) = ' 0 4 -2 0 1 0 2 -.00084'
RECS( 99) = ' 0 4 -2 0 2 0 2 -.00008'
RECS(100) = ' 0 4 -1 -2 0 1 2 -.00003'
RECS(101) = ' 0 4 -1 0 0 -1 2 .00003'
RECS(102) = ' 0 4 0 -2 0 0 2 -.00169'
RECS(103) = ' 0 4 0 -2 1 0 2 -.00070'
RECS(104) = ' 0 4 0 -2 2 0 2 -.00007'
RECS(105) = ' 1 -4 0 3 -1 0 2 .00014'
RECS(106) = ' 1 -4 0 3 0 0 2 .00075'
RECS(107) = ' 1 -4 1 1 0 1 2 -.00003'
RECS(108) = ' 1 -4 2 1 -1 0 2 .00036'
RECS(109) = ' 1 -4 2 1 0 0 2 .00194'
RECS(110) = ' 1 -4 3 1 0 -1 2 .00015'
RECS(111) = ' 1 -4 4 -1 -1 0 2 .00007'
RECS(112) = ' 1 -4 4 -1 0 0 2 .00037'
RECS(113) = ' 1 -4 5 -1 0 -1 2 .00004'
RECS(114) = ' 1 -3 -1 2 0 1 2 -.00009'
RECS(115) = ' 1 -3 0 0 -2 0 2 -.00004'
RECS(116) = ' 1 -3 0 2 -2 0 2 -.00003'
RECS(117) = ' 1 -3 0 2 -1 0 2 .00125'
RECS(118) = ' 1 -3 0 2 0 0 2 .00664'
RECS(119) = ' 1 -3 1 0 0 1 2 -.00011'
RECS(120) = ' 1 -3 1 1 0 0 2 -.00007'
RECS(121) = ' 1 -3 1 2 0 -1 2 .00010'
RECS(122) = ' 1 -3 2 0 -2 0 2 -.00004'
RECS(123) = ' 1 -3 2 0 -1 0 2 .00151'
RECS(124) = ' 1 -3 2 0 0 0 2 .00801'
RECS(125) = ' 1 -3 2 2 0 0 2 -.00007'
RECS(126) = ' 1 -3 3 0 -1 -1 2 .00010'
RECS(127) = ' 1 -3 3 0 0 -1 2 .00054'
RECS(128) = ' 1 -3 4 -2 -1 0 2 .00005'
RECS(129) = ' 1 -3 4 -2 0 0 2 .00024'
RECS(130) = ' 1 -3 4 0 0 0 2 -.00008'
RECS(131) = ' 1 -3 4 0 1 0 2 .00003'
RECS(132) = ' 1 -2 -2 1 -2 0 2 -.00004'
RECS(133) = ' 1 -2 -2 3 0 0 2 -.00016'
RECS(134) = ' 1 -2 -1 1 -1 1 2 -.00007'
RECS(135) = ' 1 -2 -1 1 0 1 2 -.00042'
RECS(136) = ' 1 -2 0 -1 -3 0 2 -.00004'
RECS(137) = ' 1 -2 0 -1 -2 0 2 -.00019'
RECS(138) = ' 1 -2 0 1 -2 0 2 -.00029'
RECS(139) = ' 1 -2 0 0 0 1 2 .00004'
RECS(140) = ' 1 -2 0 1 -1 0 2 .00947'
RECS(141) = ' 1 -2 0 1 0 0 2 .05019'
RECS(142) = ' 1 -2 0 3 0 0 2 -.00014'
RECS(143) = ' 1 -2 1 -1 0 1 2 -.00009'
RECS(144) = ' 1 -2 1 0 -1 0 2 -.00005'
RECS(145) = ' 1 -2 1 0 0 0 2 -.00027'
RECS(146) = ' 1 -2 1 1 -1 -1 2 .00007'
RECS(147) = ' 1 -2 1 1 0 -1 2 .00046'
RECS(148) = ' 1 -2 2 -1 -2 0 2 -.00005'
RECS(149) = ' 1 -2 2 -1 -1 0 2 .00180'
RECS(150) = ' 1 -2 2 -1 0 0 2 .00953'
RECS(151) = ' 1 -2 2 1 0 0 2 -.00055'
RECS(152) = ' 1 -2 2 1 1 0 2 .00017'
RECS(153) = ' 1 -2 3 -1 -1 -1 2 .00008'
RECS(154) = ' 1 -2 3 -1 0 -1 2 .00044'
RECS(155) = ' 1 -2 3 1 0 -1 2 -.00004'
RECS(156) = ' 1 -2 4 -1 0 0 2 -.00012'
RECS(157) = ' 1 -1 -2 0 -2 0 2 -.00012'
RECS(158) = ' 1 -1 -2 2 -1 0 2 -.00014'
RECS(159) = ' 1 -1 -2 2 0 0 2 -.00079'
RECS(160) = ' 1 -1 -1 0 -1 1 2 -.00011'
RECS(161) = ' 1 -1 -1 0 0 1 2 -.00090'
RECS(162) = ' 1 -1 -1 1 0 0 2 .00004'
RECS(163) = ' 1 -1 0 0 -2 0 2 -.00152'
RECS(164) = ' 1 -1 0 0 -1 0 2 .04946'
RECS(165) = ' 1 -1 0 0 0 0 2 .26216'
RECS(166) = ' 1 -1 0 2 -1 0 2 .00005'
RECS(167) = ' 1 -1 0 2 0 0 2 -.00169'
RECS(168) = ' 1 -1 0 2 1 0 2 -.00028'
RECS(169) = ' 1 -1 1 0 -1 -1 2 .00008'
RECS(170) = ' 1 -1 1 0 0 -1 2 .00076'
RECS(171) = ' 1 -1 2 -2 0 0 2 -.00015'
RECS(172) = ' 1 -1 2 0 -1 0 2 .00010'
RECS(173) = ' 1 -1 2 0 0 0 2 -.00343'
RECS(174) = ' 1 -1 2 0 1 0 2 .00075'
RECS(175) = ' 1 -1 2 0 2 0 2 .00005'
RECS(176) = ' 1 -1 3 0 0 -1 2 -.00022'
RECS(177) = ' 1 -1 4 -2 0 0 2 -.00007'
RECS(178) = ' 1 0 -3 1 0 1 2 -.00009'
RECS(179) = ' 1 0 -2 1 -1 0 2 -.00044'
RECS(180) = ' 1 0 -2 1 0 0 2 -.00193'
RECS(181) = ' 1 0 -1 0 0 0 2 .00005'
RECS(182) = ' 1 0 -1 1 0 1 2 .00010'
RECS(183) = ' 1 0 0 -1 -2 0 2 .00012'
RECS(184) = ' 1 0 0 -1 -1 0 2 -.00137'
RECS(185) = ' 1 0 0 -1 0 0 2 -.00741'
RECS(186) = ' 1 0 0 1 -1 0 2 .00059'
RECS(187) = ' 1 0 0 1 0 0 2 -.02062'
RECS(188) = ' 1 0 0 1 1 0 2 -.00414'
RECS(189) = ' 1 0 0 1 2 0 2 .00012'
RECS(190) = ' 1 0 1 0 0 0 2 .00012'
RECS(191) = ' 1 0 1 1 0 -1 2 -.00013'
RECS(192) = ' 1 0 2 -1 -1 0 2 .00011'
RECS(193) = ' 1 0 2 -1 0 0 2 -.00394'
RECS(194) = ' 1 0 2 -1 1 0 2 -.00087'
RECS(195) = ' 1 0 3 -1 0 -1 2 -.00017'
RECS(196) = ' 1 0 3 -1 1 -1 2 -.00004'
RECS(197) = ' 1 1 -4 0 0 2 2 .00029'
RECS(198) = ' 1 1 -3 0 -1 1 2 -.00006'
RECS(199) = ' 1 1 -3 0 0 1 2 .00713'
RECS(200) = ' 1 1 -2 0 -2 0 2 .00010'
RECS(201) = ' 1 1 -2 0 -1 0 2 -.00137'
RECS(202) = ' 1 1 -2 0 0 0 2 .12199'
RECS(203) = ' 1 1 -2 0 0 2 2 -.00007'
RECS(204) = ' 1 1 -2 2 0 0 2 -.00018'
RECS(205) = ' 1 1 -2 2 1 0 2 -.00004'
RECS(206) = ' 1 1 -1 0 0 -1 2 -.00102'
RECS(207) = ' 1 1 -1 0 0 1 2 -.00288'
RECS(208) = ' 1 1 -1 0 1 1 2 .00008'
RECS(209) = ' 1 1 0 -2 -1 0 2 -.00007'
RECS(210) = ' 1 1 0 0 -2 0 2 -.00005'
RECS(211) = ' 1 1 0 0 -1 0 2 .00730'
RECS(212) = ' 1 1 0 0 0 0 2 -.36872'
RECS(213) = ' 1 1 0 0 1 0 2 -.05002'
RECS(214) = ' 1 1 0 0 2 0 2 .00108'
RECS(215) = ' 1 1 1 0 0 -1 2 -.00292'
RECS(216) = ' 1 1 1 0 1 -1 2 -.00005'
RECS(217) = ' 1 1 2 -2 0 0 2 -.00018'
RECS(218) = ' 1 1 2 -2 1 0 2 -.00005'
RECS(219) = ' 1 1 2 0 0 -2 2 -.00007'
RECS(220) = ' 1 1 2 0 0 0 2 -.00525'
RECS(221) = ' 1 1 2 0 1 0 2 .00020'
RECS(222) = ' 1 1 2 0 2 0 2 .00010'
RECS(223) = ' 1 1 3 0 0 -1 2 -.00030'
RECS(224) = ' 1 2 -3 1 0 1 2 -.00017'
RECS(225) = ' 1 2 -2 -1 -1 0 2 -.00012'
RECS(226) = ' 1 2 -2 1 -1 0 2 .00012'
RECS(227) = ' 1 2 -2 1 0 0 2 -.00394'
RECS(228) = ' 1 2 -2 1 1 0 2 -.00078'
RECS(229) = ' 1 2 -1 -1 0 1 2 -.00013'
RECS(230) = ' 1 2 -1 0 0 0 2 .00011'
RECS(231) = ' 1 2 0 -1 -1 0 2 .00060'
RECS(232) = ' 1 2 0 -1 0 0 2 -.02062'
RECS(233) = ' 1 2 0 -1 1 0 2 -.00409'
RECS(234) = ' 1 2 0 -1 2 0 2 .00008'
RECS(235) = ' 1 2 0 1 0 0 2 .00032'
RECS(236) = ' 1 2 0 1 1 0 2 .00020'
RECS(237) = ' 1 2 0 1 2 0 2 .00012'
RECS(238) = ' 1 2 1 -1 0 -1 2 .00011'
RECS(239) = ' 1 2 2 -1 0 0 2 .00008'
RECS(240) = ' 1 2 2 -1 1 0 2 .00006'
RECS(241) = ' 1 3 -4 2 0 0 2 -.00006'
RECS(242) = ' 1 3 -3 0 0 1 2 -.00023'
RECS(243) = ' 1 3 -3 0 1 1 2 -.00004'
RECS(244) = ' 1 3 -2 0 -1 0 2 -.00011'
RECS(245) = ' 1 3 -2 0 0 0 2 -.00342'
RECS(246) = ' 1 3 -2 0 1 0 2 -.00067'
RECS(247) = ' 1 3 -1 0 0 -1 2 .00007'
RECS(248) = ' 1 3 0 -2 -1 0 2 .00004'
RECS(249) = ' 1 3 0 -2 0 0 2 -.00169'
RECS(250) = ' 1 3 0 -2 1 0 2 -.00034'
RECS(251) = ' 1 3 0 0 0 0 2 -.01128'
RECS(252) = ' 1 3 0 0 1 0 2 -.00723'
RECS(253) = ' 1 3 0 0 2 0 2 -.00151'
RECS(254) = ' 1 3 0 0 3 0 2 -.00010'
RECS(255) = ' 1 3 1 0 0 -1 2 .00004'
RECS(256) = ' 1 4 -4 1 0 0 2 -.00011'
RECS(257) = ' 1 4 -3 -1 0 1 2 -.00004'
RECS(258) = ' 1 4 -2 -1 0 0 2 -.00054'
RECS(259) = ' 1 4 -2 -1 1 0 2 -.00010'
RECS(260) = ' 1 4 -2 1 0 0 2 -.00041'
RECS(261) = ' 1 4 -2 1 1 0 2 -.00026'
RECS(262) = ' 1 4 -2 1 2 0 2 -.00005'
RECS(263) = ' 1 4 0 -3 0 0 2 -.00014'
RECS(264) = ' 1 4 0 -1 0 0 2 -.00216'
RECS(265) = ' 1 4 0 -1 1 0 2 -.00138'
RECS(266) = ' 1 4 0 -1 2 0 2 -.00029'
RECS(267) = ' 2 -4 0 4 0 0 2 .00018'
RECS(268) = ' 2 -4 2 2 0 0 2 .00077'
RECS(269) = ' 2 -4 3 2 0 -1 2 .00006'
RECS(270) = ' 2 -4 4 0 0 0 2 .00048'
RECS(271) = ' 2 -4 5 0 0 -1 2 .00006'
RECS(272) = ' 2 -3 0 3 -1 0 2 -.00007'
RECS(273) = ' 2 -3 0 3 0 0 2 .00180'
RECS(274) = ' 2 -3 1 1 0 1 2 -.00009'
RECS(275) = ' 2 -3 1 3 0 -1 2 .00004'
RECS(276) = ' 2 -3 2 1 -1 0 2 -.00017'
RECS(277) = ' 2 -3 2 1 0 0 2 .00467'
RECS(278) = ' 2 -3 3 1 0 -1 2 .00035'
RECS(279) = ' 2 -3 4 -1 -1 0 2 -.00003'
RECS(280) = ' 2 -3 4 -1 0 0 2 .00090'
RECS(281) = ' 2 -3 5 -1 0 -1 2 .00010'
RECS(282) = ' 2 -2 -2 4 0 0 2 -.00006'
RECS(283) = ' 2 -2 -1 2 0 1 2 -.00022'
RECS(284) = ' 2 -2 0 0 -2 0 2 -.00010'
RECS(285) = ' 2 -2 0 2 -1 0 2 -.00060'
RECS(286) = ' 2 -2 0 2 0 0 2 .01601'
RECS(287) = ' 2 -2 1 0 0 1 2 -.00027'
RECS(288) = ' 2 -2 1 1 0 0 2 -.00017'
RECS(289) = ' 2 -2 1 2 0 -1 2 .00025'
RECS(290) = ' 2 -2 2 0 -1 0 2 -.00072'
RECS(291) = ' 2 -2 2 0 0 0 2 .01932'
RECS(292) = ' 2 -2 3 -1 0 0 2 -.00004'
RECS(293) = ' 2 -2 3 0 -1 -1 2 -.00005'
RECS(294) = ' 2 -2 3 0 0 -1 2 .00130'
RECS(295) = ' 2 -2 4 -2 0 0 2 .00059'
RECS(296) = ' 2 -2 4 0 0 -2 2 .00005'
RECS(297) = ' 2 -2 5 -2 0 -1 2 .00005'
RECS(298) = ' 2 -1 -2 1 -2 0 2 -.00010'
RECS(299) = ' 2 -1 -2 3 0 0 2 -.00039'
RECS(300) = ' 2 -1 -1 1 -1 1 2 .00003'
RECS(301) = ' 2 -1 -1 1 0 1 2 -.00102'
RECS(302) = ' 2 -1 0 -1 -2 0 2 -.00046'
RECS(303) = ' 2 -1 0 1 -2 0 2 .00007'
RECS(304) = ' 2 -1 0 0 0 1 2 .00009'
RECS(305) = ' 2 -1 0 1 -1 0 2 -.00451'
RECS(306) = ' 2 -1 0 1 0 0 2 .12099'
RECS(307) = ' 2 -1 1 -1 0 1 2 -.00023'
RECS(308) = ' 2 -1 1 0 0 0 2 -.00065'
RECS(309) = ' 2 -1 1 1 -1 -1 2 -.00004'
RECS(310) = ' 2 -1 1 1 0 -1 2 .00113'
RECS(311) = ' 2 -1 2 -1 -1 0 2 -.00086'
RECS(312) = ' 2 -1 2 -1 0 0 2 .02298'
RECS(313) = ' 2 -1 2 1 0 0 2 .00010'
RECS(314) = ' 2 -1 2 1 1 0 2 -.00008'
RECS(315) = ' 2 -1 3 -1 -1 -1 2 -.00004'
RECS(316) = ' 2 -1 3 -1 0 -1 2 .00106'
RECS(317) = ' 2 0 -3 2 0 1 2 -.00008'
RECS(318) = ' 2 0 -2 0 -2 0 2 -.00028'
RECS(319) = ' 2 0 -2 2 -1 0 2 .00007'
RECS(320) = ' 2 0 -2 2 0 0 2 -.00190'
RECS(321) = ' 2 0 -1 0 -1 1 2 .00005'
RECS(322) = ' 2 0 -1 0 0 1 2 -.00217'
RECS(323) = ' 2 0 -1 1 0 0 2 .00010'
RECS(324) = ' 2 0 0 0 -2 0 2 .00033'
RECS(325) = ' 2 0 0 0 -1 0 2 -.02358'
RECS(326) = ' 2 0 0 0 0 0 2 .63194'
RECS(327) = ' 2 0 0 2 0 0 2 .00036'
RECS(328) = ' 2 0 0 2 1 0 2 .00013'
RECS(329) = ' 2 0 1 0 -1 -1 2 -.00004'
RECS(330) = ' 2 0 1 0 0 -1 2 .00192'
RECS(331) = ' 2 0 2 -2 0 0 2 -.00036'
RECS(332) = ' 2 0 2 0 0 0 2 .00072'
RECS(333) = ' 2 0 2 0 1 0 2 -.00036'
RECS(334) = ' 2 0 2 0 2 0 2 .00012'
RECS(335) = ' 2 0 3 0 0 -1 2 .00005'
RECS(336) = ' 2 1 -3 1 0 1 2 -.00022'
RECS(337) = ' 2 1 -2 1 -1 0 2 .00021'
RECS(338) = ' 2 1 -2 1 0 0 2 -.00466'
RECS(339) = ' 2 1 -1 -1 0 1 2 -.00007'
RECS(340) = ' 2 1 -1 0 0 0 2 .00010'
RECS(341) = ' 2 1 0 -1 -1 0 2 .00065'
RECS(342) = ' 2 1 0 -1 0 0 2 -.01786'
RECS(343) = ' 2 1 0 1 -1 0 2 -.00008'
RECS(344) = ' 2 1 0 1 0 0 2 .00447'
RECS(345) = ' 2 1 0 1 1 0 2 .00197'
RECS(346) = ' 2 1 0 1 2 0 2 .00028'
RECS(347) = ' 2 1 2 -1 0 0 2 .00085'
RECS(348) = ' 2 1 2 -1 1 0 2 .00041'
RECS(349) = ' 2 1 2 -1 2 0 2 .00005'
RECS(350) = ' 2 2 -4 0 0 2 2 .00070'
RECS(351) = ' 2 2 -3 0 0 1 2 .01719'
RECS(352) = ' 2 2 -2 0 -1 0 2 .00066'
RECS(353) = ' 2 2 -2 0 0 0 2 .29401'
RECS(354) = ' 2 2 -2 2 0 0 2 .00004'
RECS(355) = ' 2 2 -1 0 0 -1 2 -.00246'
RECS(356) = ' 2 2 -1 0 0 1 2 .00062'
RECS(357) = ' 2 2 -1 0 1 1 2 -.00004'
RECS(358) = ' 2 2 0 0 -1 0 2 -.00103'
RECS(359) = ' 2 2 0 0 0 0 2 .07992'
RECS(360) = ' 2 2 0 0 1 0 2 .02382'
RECS(361) = ' 2 2 0 0 2 0 2 .00259'
RECS(362) = ' 2 2 1 0 0 -1 2 .00063'
RECS(363) = ' 2 2 2 -2 0 0 2 .00004'
RECS(364) = ' 2 2 2 0 0 0 2 .00053'
RECS(365) = ' 2 3 -3 1 0 1 2 .00003'
RECS(366) = ' 2 3 -2 -1 -1 0 2 .00006'
RECS(367) = ' 2 3 -2 -1 0 0 2 .00004'
RECS(368) = ' 2 3 -2 1 0 0 2 .00085'
RECS(369) = ' 2 3 -2 1 1 0 2 .00037'
RECS(370) = ' 2 3 -2 1 2 0 2 .00004'
RECS(371) = ' 2 3 0 -1 -1 0 2 -.00009'
RECS(372) = ' 2 3 0 -1 0 0 2 .00447'
RECS(373) = ' 2 3 0 -1 1 0 2 .00195'
RECS(374) = ' 2 3 0 -1 2 0 2 .00022'
RECS(375) = ' 2 3 0 1 0 0 2 -.00003'
RECS(376) = ' 2 4 -3 0 0 1 2 .00005'
RECS(377) = ' 2 4 -2 0 0 0 2 .00074'
RECS(378) = ' 2 4 -2 0 1 0 2 .00032'
RECS(379) = ' 2 4 -2 0 2 0 2 .00003'
RECS(380) = ' 2 4 0 -2 0 0 2 .00037'
RECS(381) = ' 2 4 0 -2 1 0 2 .00016'
RECS(382) = ' 2 4 0 0 0 0 2 .00117'
RECS(383) = ' 2 4 0 0 1 0 2 .00101'
RECS(384) = ' 2 4 0 0 2 0 2 .00033'
RECS(385) = ' 2 4 0 0 3 0 2 .00005'
RECS(386) = ' 0 0 0 1 0 0 3 -.00021'
RECS(387) = ' 0 0 2 -1 0 0 3 -.00004'
RECS(388) = ' 0 1 -2 0 0 0 3 .00004'
RECS(389) = ' 0 1 0 0 -1 0 3 .00019'
RECS(390) = ' 0 1 0 0 0 0 3 -.00375'
RECS(391) = ' 0 1 0 0 1 0 3 -.00059'
RECS(392) = ' 0 1 0 0 2 0 3 .00005'
RECS(393) = ' 0 2 -2 1 0 0 3 -.00012'
RECS(394) = ' 0 2 0 -1 0 0 3 -.00061'
RECS(395) = ' 0 2 0 -1 1 0 3 -.00010'
RECS(396) = ' 0 3 -2 0 0 0 3 -.00010'
RECS(397) = ' 0 3 0 -2 0 0 3 -.00007'
RECS(398) = ' 0 3 0 0 0 0 3 -.00031'
RECS(399) = ' 0 3 0 0 1 0 3 -.00019'
RECS(400) = ' 0 3 0 0 2 0 3 -.00004'
RECS(401) = ' 0 4 0 -1 0 0 3 -.00008'
RECS(402) = ' 0 4 0 -1 1 0 3 -.00005'
RECS(403) = ' 1 -4 0 2 0 0 3 .00006'
RECS(404) = ' 1 -4 2 0 0 0 3 .00006'
RECS(405) = ' 1 -3 0 1 -1 0 3 .00014'
RECS(406) = ' 1 -3 0 1 0 0 3 .00035'
RECS(407) = ' 1 -3 2 -1 0 0 3 .00006'
RECS(408) = ' 1 -2 0 0 -2 0 3 .00004'
RECS(409) = ' 1 -2 0 0 -1 0 3 .00051'
RECS(410) = ' 1 -2 0 0 0 0 3 .00128'
RECS(411) = ' 1 -2 0 2 0 0 3 .00008'
RECS(412) = ' 1 -2 2 0 0 0 3 .00011'
RECS(413) = ' 1 -1 0 -1 0 0 3 -.00007'
RECS(414) = ' 1 -1 0 1 -1 0 3 -.00009'
RECS(415) = ' 1 -1 0 1 0 0 3 .00065'
RECS(416) = ' 1 -1 0 1 1 0 3 -.00009'
RECS(417) = ' 1 -1 2 -1 0 0 3 .00013'
RECS(418) = ' 1 0 0 0 -1 0 3 -.00059'
RECS(419) = ' 1 0 0 0 0 0 3 .00399'
RECS(420) = ' 1 0 0 0 1 0 3 -.00052'
RECS(421) = ' 1 1 -2 1 0 0 3 .00004'
RECS(422) = ' 1 1 0 -1 -1 0 3 -.00003'
RECS(423) = ' 1 1 0 -1 0 0 3 .00022'
RECS(424) = ' 1 1 0 -1 1 0 3 -.00003'
RECS(425) = ' 1 1 0 1 0 0 3 .00008'
RECS(426) = ' 1 1 0 1 1 0 3 .00003'
RECS(427) = ' 1 2 -2 0 0 0 3 .00005'
RECS(428) = ' 1 2 0 0 -1 0 3 -.00005'
RECS(429) = ' 1 2 0 0 0 0 3 .00146'
RECS(430) = ' 1 2 0 0 1 0 3 .00059'
RECS(431) = ' 1 2 0 0 2 0 3 .00005'
RECS(432) = ' 1 3 -2 1 0 0 3 .00005'
RECS(433) = ' 1 3 0 -1 0 0 3 .00024'
RECS(434) = ' 1 3 0 -1 1 0 3 .00010'
RECS(435) = ' 1 4 -2 0 0 0 3 .00004'
RECS(436) = ' 1 4 0 0 0 0 3 .00005'
RECS(437) = ' 1 4 0 0 1 0 3 .00005'
RECS(438) = ' 2 -4 2 1 0 0 3 -.00006'
RECS(439) = ' 2 -3 0 2 0 0 3 -.00019'
RECS(440) = ' 2 -3 2 0 -1 0 3 -.00003'
RECS(441) = ' 2 -3 2 0 0 0 3 -.00019'
RECS(442) = ' 2 -2 0 1 -1 0 3 -.00018'
RECS(443) = ' 2 -2 0 1 0 0 3 -.00107'
RECS(444) = ' 2 -2 2 -1 -1 0 3 -.00003'
RECS(445) = ' 2 -2 2 -1 0 0 3 -.00020'
RECS(446) = ' 2 -1 0 0 -2 0 3 .00004'
RECS(447) = ' 2 -1 0 0 -1 0 3 -.00066'
RECS(448) = ' 2 -1 0 0 0 0 3 -.00389'
RECS(449) = ' 2 -1 0 2 0 0 3 .00007'
RECS(450) = ' 2 -1 2 0 0 0 3 .00010'
RECS(451) = ' 2 0 -2 1 0 0 3 .00005'
RECS(452) = ' 2 0 0 -1 -1 0 3 .00004'
RECS(453) = ' 2 0 0 -1 0 0 3 .00022'
RECS(454) = ' 2 0 0 1 -1 0 3 -.00003'
RECS(455) = ' 2 0 0 1 0 0 3 .00059'
RECS(456) = ' 2 0 0 1 1 0 3 .00011'
RECS(457) = ' 2 0 2 -1 0 0 3 .00011'
RECS(458) = ' 2 1 0 0 -1 0 3 -.00021'
RECS(459) = ' 2 1 0 0 0 0 3 .00359'
RECS(460) = ' 2 1 0 0 1 0 3 .00067'
RECS(461) = ' 2 2 -2 1 0 0 3 .00004'
RECS(462) = ' 2 2 0 -1 0 0 3 .00019'
RECS(463) = ' 2 2 0 -1 1 0 3 .00004'
RECS(464) = ' 2 3 -2 0 0 0 3 .00004'
RECS(465) = ' 2 3 0 0 0 0 3 .00033'
RECS(466) = ' 2 3 0 0 1 0 3 .00021'
RECS(467) = ' 2 3 0 0 2 0 3 .00004'
RECS(468) = ' 2 4 0 -1 0 0 3 .00005'
RECS(469) = ' 3 -2 0 2 0 0 3 -.00036'
RECS(470) = ' 3 -2 2 0 0 0 3 -.00036'
RECS(471) = ' 3 -1 0 1 -1 0 3 .00012'
RECS(472) = ' 3 -1 0 1 0 0 3 -.00210'
RECS(473) = ' 3 -1 2 -1 0 0 3 -.00039'
RECS(474) = ' 3 0 -2 2 0 0 3 .00005'
RECS(475) = ' 3 0 0 0 -1 0 3 .00043'
RECS(476) = ' 3 0 0 0 0 0 3 -.00765'
RECS(477) = ' 3 1 -2 1 0 0 3 .00011'
RECS(478) = ' 3 1 0 -1 0 0 3 .00043'
RECS(479) = ' 3 1 0 1 0 0 3 -.00016'
RECS(480) = ' 3 1 0 1 1 0 3 -.00007'
RECS(481) = ' 3 2 0 0 -1 0 3 .00004'
RECS(482) = ' 3 2 0 0 0 0 3 -.00100'
RECS(483) = ' 3 2 0 0 1 0 3 -.00044'
RECS(484) = ' 3 2 0 0 2 0 3 -.00005'
END SUBROUTINE INIHARMONICS
end module timespace_data
!
!
!
! ==========================================================================
! ==========================================================================
! ==========================================================================
!>
Module M_arcuv ! plotbuitenbeentje
implicit none
double precision, allocatable :: arcuv(:,:,:)
End module M_arcuv
!
!
!
! ==========================================================================
! ==========================================================================
! ==========================================================================
!>
module m_spiderweb ! plot spiderweb
implicit none
double precision, allocatable :: spw(:,:,:)
end module m_spiderweb
!
!
!
! ==========================================================================
! ==========================================================================
! ==========================================================================
!>
!!--copyright-------------------------------------------------------------------
! Copyright (c) 2003, WL | Delft Hydraulics. All rights reserved.
!!--disclaimer------------------------------------------------------------------
! This code is part of the Delft3D software system. WL|Delft Hydraulics has
! developed c.q. manufactured this code to its best ability and according to the
! state of the art. Nevertheless, there is no express or implied warranty as to
! this software whether tangible or intangible. In particular, there is no
! express or implied warranty as to the fitness for a particular purpose of this
! software, whether tangible or intangible. The intellectual property rights
! related to this software code remain with WL|Delft Hydraulics at all times.
! For details on the licensing agreement, we refer to the Delft3D software
! license and any modifications to this license, if applicable. These documents
! are available upon request.
!!--version information---------------------------------------------------------
! $Author$
! $Date$
! $Revision$
!!--description-----------------------------------------------------------------
! NONE
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
module timespace_triangle
use precision
use timespace_data
use m_alloc
implicit none
integer :: nsold ! nr of samples in previous triangulation
integer :: numtri
integer , allocatable, dimension(:, :) :: indx
double precision, allocatable, dimension(:) :: xcent
double precision, allocatable, dimension(:) :: ycent
interface triint
module procedure triint_z1D
module procedure triint_z2D
module procedure triint_z3D
end interface triint
interface get_extend
module procedure get_extend1D
module procedure get_extend2D
end interface get_extend
interface find_nearest
module procedure find_nearest1D
module procedure find_nearest2D
module procedure find_nearest1D_missing_value
module procedure find_nearest2D_missing_value
end interface find_nearest
contains
!
!
! ==========================================================================
!>
subroutine pinpok(xl, yl, n, x, y, inside)
! Author: H. Kernkamp
implicit none
double precision , intent(in) :: xl, yl ! point under consideration
integer , intent(in) :: n
double precision, dimension(n), intent(in) :: x, y ! polygon(n)
integer , intent(out) :: inside
integer :: i, i1, i2, np, rechts
double precision :: rl, rm, x1, x2, y1, y2
if (n .le. 2) then
inside = 1
else
np = 0
5 continue
np = np + 1
if (np .le. n) then
if ( x(np) .ne. dmiss_default) goto 5
end if
np = np - 1
inside = 0
rechts = 0
i = 0
10 continue
i1 = mod(i,np) + 1
i2 = mod(i1,np) + 1
x1 = x(i1)
x2 = x(i2)
y1 = y(i1)
y2 = y(i2)
if (xl .ge. min(x1,x2) .and. xl .le. max(x1,x2) ) then
if (xl .eq. x1 .and. yl .eq. y1 .or. & ! tussen of op lijnstuk
(x1 .eq. x2 .and. & ! op punt 1
yl .ge. min(y1,y2) .and. yl .le. max(y1,y2) ) .or. & ! op verticale lijn
(yl .eq. y1 .and. y1 .eq. y2) ) then ! op horizontale lijn
inside = 1
return
else if (x1 .ne. x2) then ! scheve lijn
rl = ( xl - x1 ) / ( x2 - x1 )
rm = ( y1 - yl ) + rl * ( y2 - y1 )
if (rm .eq. 0) then ! op scheve lijn
inside = 1
return
else if (rm .gt. 0d0) then ! onder scheve lijn
if (xl .eq. x1 .or. xl .eq. x2) then
if (x1 .gt. xl .or. x2 .gt. xl) then
rechts = rechts + 1
end if
end if
inside = 1 - inside
end if
end if
end if
i = i + 1
if (i .lt. np) goto 10
if (mod(rechts,2) .ne. 0) inside = 1 - inside
end if
end subroutine pinpok
!
!
! ==========================================================================
!>
! This subroutine interpolates one unstructured dataset xss, yss, zss, kcss, nss to another x, y, z, kcs, nx
! It is the only one in this module that is of practical interest to the meteo module.
! The rest of the subroutines in this module are assisting this one.
! JDLA = 1 (re)triangulates
subroutine triint_z2D( xss, yss, zss, kcsss, nss, &
x , y , z , kcs , kx , mnx, jdla , indxn, wfn )
implicit none
! Global variables
integer, intent(in) :: nss ! Dimension of samples
double precision, dimension(:), intent(in) :: xss ! samples
double precision, dimension(:), intent(in) :: yss
double precision, dimension(:), intent(in) :: zss ! dimension: nss*kx
integer , dimension(:), intent(in) :: kcsss ! samples mask
integer, intent(in) :: mnx ! Dimension of grid
integer, intent(in) :: kx ! vectormax
double precision, dimension(:), intent(in) :: x ! grid
double precision, dimension(:), intent(in) :: y
double precision, dimension(:,:), intent(out) :: z ! dimension: nx*kx
integer , dimension(:), intent(in) :: kcs ! grid mask
integer, intent(in) :: jdla ! refresh delauney yes /no
integer , optional :: indxn(:,:) ! if present get weightfactors and indices
double precision, optional :: wfn (:,:)
call triint_z1D( xss, yss, zss, kcsss, nss, &
x , y , z , kcs , kx , mnx, jdla, indxn, wfn )
end subroutine triint_z2D
!
!
! ==========================================================================
!>
subroutine triint_z3D( xss, yss, zss, kcsss, nss, &
x , y , z , kcs , kx , mnx, jdla, indxn, wfn )
implicit none
! Global variables
integer, intent(in) :: nss ! Dimension of samples
double precision, dimension(:), intent(in) :: xss ! samples
double precision, dimension(:), intent(in) :: yss
double precision, dimension(:), intent(in) :: zss ! dimension: nss*kx
integer , dimension(:), intent(in) :: kcsss ! samples mask
integer, intent(in) :: mnx ! Dimension of grid
integer, intent(in) :: kx ! vectormax
double precision, dimension(:), intent(in) :: x ! grid
double precision, dimension(:), intent(in) :: y
double precision, dimension(:,:,:), intent(out) :: z ! dimension: nx*kx
integer , dimension(:), intent(in) :: kcs ! grid mask
integer, intent(in) :: jdla ! refresh delauney yes /no
integer , optional :: indxn(:,:) ! if present get weightfactors and indices
double precision, optional :: wfn (:,:)
call triint_z1D( xss, yss, zss, kcsss, nss, &
x , y , z , kcs , kx , mnx, jdla, indxn, wfn )
end subroutine triint_z3D
!
!
! ==========================================================================
!>
subroutine triint_z1D( xss, yss, zss, kcsss, nss, &
x , y , z , kcs , kx , mnx, jdla, indxn, wfn )
implicit none
! Global variables
integer, intent(in) :: nss ! Dimension of samples
double precision, dimension(:), intent(in) :: xss ! samples
double precision, dimension(:), intent(in) :: yss
double precision, dimension(:), intent(in) :: zss ! dimension: nss*kx
integer , dimension(:), intent(in) :: kcsss ! samples mask
integer, intent(in) :: mnx ! Dimension of grid
integer, intent(in) :: kx ! vectormax
double precision, dimension(:), intent(in) :: x ! grid
double precision, dimension(:), intent(in) :: y
double precision, dimension(kx*mnx), intent(out):: z ! dimension: mnx*kx
integer , dimension(:), intent(in) :: kcs ! grid mask
integer, intent(in) :: jdla ! refresh delauney yes /no
integer , optional :: indxn(:,:) ! if present get weightfactors and indices
double precision, optional :: wfn (:,:)
! Local variables
double precision, dimension(8) :: x_set
double precision, dimension(8) :: y_set
integer , dimension(8) :: kcs_set = 1
double precision, dimension(4) :: x_extr
double precision, dimension(4) :: y_extr
double precision, dimension(4) :: z_extr
double precision, dimension(3) :: zp
integer , dimension(3) :: indxp
double precision, dimension(:), allocatable :: xs
double precision, dimension(:), allocatable :: ys
double precision, dimension(:), allocatable :: zs
integer , dimension(:), allocatable :: kcss
integer :: ns
integer :: k, n, jgetw, ierr ! , MOUT
logical :: extra = .false. ! nu even niet
!! executable statements -------------------------------------------------------
!
! JDLA=1, DO DE LAUNEY
! JSLO=1, ALSO SLOPES RD4
if (nss<1) then
return
end if
call realloc(xs,nss+8,1)
call realloc(ys,nss+8,1)
call realloc(zs,nss+8,1)
call realloc(kcss,nss+8,1)
ns = 0
do k = 1,nss
if (kcsss(k) == 1) then
ns = ns + 1
xs(ns) = xss(k)
ys(ns) = yss(k)
do n = 1,kx
zs(kx*(ns-1)+n) = zss(kx*(k-1)+n)
enddo
kcss(ns) = 1
end if
enddo
if (extra) then
call get_extend(mnx, x, y, kcs, x_set(1:4), y_set(1:4))
call get_extend(ns, xs, ys, kcss, x_set(5:8), y_set(5:8))
call get_extend(8, x_set, y_set, kcs_set, x_extr, y_extr)
call extrapolate(ns, xs, ys, zs, kcss, 4, x_extr, y_extr, z_extr)
xs(ns + 1:ns + 4) = x_extr
ys(ns + 1:ns + 4) = y_extr
zs(ns + 1:ns + 4) = z_extr
ns = ns + 4
end if
if (jdla==1) then
! call dlauny(xs, ys, ns)
call DLAUN(XS,YS,NS,1,ierr)
end if
jgetw = 0 ! niets met gewichten, doe interpolatie
if ( present(indxn) .and. jdla .eq. 1) jgetw = 1 ! haal gewichten doe interpolatie , gebruik gewichten
if ( present(indxn) .and. jdla .eq. 0) jgetw = 2 ! doe interpolatie , gebruik gewichten
do n = 1,mnx
if (kcs(n) .eq. 1) then
if (jgetw .le. 1) then
call findtri_indices_weights (x(n),y( n), xs, ys, ns, zp, indxp) ! zoeken bij 0 en 1
end if
if (jgetw .eq. 1) then ! zetten bij 1
do k = 1,3
indxn(k,n) = indxp(k)
wfn(k,n) = zp(k)
enddo
else if (jgetw .eq. 2) then ! halen bij 2, je hoeft niet te zoeken
do k = 1,3
indxp(k) = indxn(k,n)
zp(k) = wfn(k,n)
enddo
end if
! en altijd interpoleren
do k = 1,kx ! over vectormax loop
if (indxp(1)==0 .or. indxp(2)==0 .or. indxp(3)==0 ) then
! z(mnx*(k-1) + n) = -999
else
z(mnx*(k-1) + n) = zp(1)*zs(kx*(indxp(1)-1)+k) + zp(2)*zs(kx*(indxp(2)-1)+k) + zp(3)*zs(kx*(indxp(3)-1)+k)
end if
enddo
end if
enddo
deallocate(xs)
deallocate(ys)
deallocate(zs)
deallocate(kcss)
end subroutine triint_z1D
!
!
! ==========================================================================
!>
subroutine findtri_indices_weights(xp, yp, xs, ys, ns, zp, indxp)
implicit none
! Global variables
double precision, intent(in) :: xp ! for this point
double precision, intent(in) :: yp
integer , intent(in) :: ns
double precision, dimension(ns), intent(in) :: xs ! on this set
double precision, dimension(ns), intent(in) :: ys
integer , dimension(3) , intent(out) :: indxp ! find indices to set
double precision, dimension(3) , intent(out) :: zp ! and corresponding weightfactors
! Local variables
integer :: k
integer :: k1
integer :: k2, n3
integer :: intri
integer :: nroldfind, nrfind
double precision :: xtmax
double precision :: xtmin
double precision :: ytmax
double precision :: ytmin
double precision, dimension(3) :: xt
double precision, dimension(3) :: yt
!
!
data nroldfind/0/
!
!! executable statements -------------------------------------------------------
!
!
indxp = 0
n3 = 3
5 continue
if (nroldfind/=0) then
k1 = max(1, nroldfind - 200)
k2 = min(numtri, nroldfind + 200)
else
k1 = 1
k2 = numtri
end if
!
do k = k1, k2
xt(1) = xs(indx(1, k))
xt(2) = xs(indx(2, k))
xt(3) = xs(indx(3, k))
yt(1) = ys(indx(1, k))
yt(2) = ys(indx(2, k))
yt(3) = ys(indx(3, k))
xtmax = max(xt(1), max(xt(2), xt(3)))
ytmax = max(yt(1), max(yt(2), yt(3)))
xtmin = min(xt(1), min(xt(2), xt(3)))
ytmin = min(yt(1), min(yt(2), yt(3)))
if (xp>=xtmin .and. xp<=xtmax .and. yp>=ytmin .and. yp<=ytmax) then
call pinpok(xp ,yp, n3, xt, yt, intri)
if (intri==1) then
nrfind = k
nroldfind = nrfind
indxp(1) = indx(1, k)
indxp(2) = indx(2, k)
indxp(3) = indx(3, k)
call linweight(xt ,yt ,xp ,yp, zp)
! write(*,*) xp, yp, k, indxp(1), indxp(2), indxp(3)
return
end if
end if
enddo
if (nroldfind/=0) then
nroldfind = 0
goto 5
end if
end subroutine findtri_indices_weights
!
!
! ==========================================================================
!>
subroutine linweight(xt ,yt ,xp ,yp, zp)
double precision, intent(in) :: xp ! for this point
double precision, intent(in) :: yp
double precision, dimension(3) :: xt ! in this triangle
double precision, dimension(3) :: yt
double precision, dimension(3) , intent(out) :: zp ! the weightfactors are...
double precision :: a11, a12, a21, a22, b1, b2, det
zp = 0
a11 = xt(2) - xt(1)
a21 = yt(2) - yt(1)
a12 = xt(3) - xt(1)
a22 = yt(3) - yt(1)
b1 = xp - xt(1)
b2 = yp - yt(1)
det = a11*a22 - a12*a21
if (abs(det)<1E-9) then
return
end if
!
zp(2) = ( a22*b1 - a12*b2)/det
zp(3) = ( -a21*b1 + a11*b2)/det
zp(1) = 1d0 - zp(2) - zp(3)
end subroutine linweight
!
!
! ==========================================================================
!>
subroutine linear(x ,y ,z ,xp ,yp , &
& zp ,jslo ,slo )
use precision
implicit none
!
!
! COMMON variables
!
double precision :: dmiss
data dmiss /-999d0/
!
! Global variables
!
integer, intent(in) :: jslo
double precision, intent(out) :: slo
double precision :: xp
double precision :: yp
double precision :: zp
double precision, dimension(3) :: x
double precision, dimension(3) :: y
double precision, dimension(3), intent(in) :: z
!
!
! Local variables
!
double precision :: a11
double precision :: a12
double precision :: a21
double precision :: a22
double precision :: a31
double precision :: a32
double precision :: b1
double precision :: b2
double precision :: det
double precision :: r3
double precision :: rlam
double precision :: rmhu
double precision :: x3
double precision :: xn
double precision :: xy
double precision :: y3
double precision :: yn
double precision :: z3
double precision :: zn
!
!
!! executable statements -------------------------------------------------------
!
!
!
!
!
zp = dmiss
a11 = x(2) - x(1)
a21 = y(2) - y(1)
a12 = x(3) - x(1)
a22 = y(3) - y(1)
b1 = xp - x(1)
b2 = yp - y(1)
!
det = a11*a22 - a12*a21
if (abs(det)<1E-12) then ! Jan Mooiman 07-01-2015
return
end if
!
rlam = (a22*b1 - a12*b2)/det
rmhu = ( - a21*b1 + a11*b2)/det
!
zp = z(1) + rlam*(z(2) - z(1)) + rmhu*(z(3) - z(1))
if (jslo==1) then
a31 = z(2) - z(1)
a32 = z(3) - z(1)
x3 = (a21*a32 - a22*a31)
y3 = -(a11*a32 - a12*a31)
z3 = (a11*a22 - a12*a21)
r3 = sqrt(x3*x3 + y3*y3 + z3*z3)
if (r3/=0) then
xn = x3/r3
yn = y3/r3
zn = z3/r3
xy = sqrt(xn*xn + yn*yn)
if (zn/=0) then
slo = abs(xy/zn)
else
slo = dmiss
end if
else
slo = dmiss
end if
end if
end subroutine linear
!
!
! ==========================================================================
!>
subroutine minmax_h(x, n, xmin, xmax ) ! BEPAAL MINIMUM EN MAXIMUM VAN EEN EENDIMENSIONALE ARRAY
use precision
implicit none
! Global variables
integer, intent(in) :: n
double precision, dimension(n), intent(in) :: x
double precision :: xmax
double precision :: xmin
integer :: i
xmin = 1E30
xmax = -1E30
do i = 1, n
xmin = min(xmin, x(i))
xmax = max(xmax, x(i))
enddo
end subroutine minmax_h
!
!
! ==========================================================================
!>
subroutine get_extend2D(n, m, x, y, kcs, x_ext, y_ext)
double precision, dimension(:,:) :: x
double precision, dimension(:,:) :: y
integer , dimension(:,:) :: kcs
integer :: n
integer :: m
double precision, dimension(:) :: x_ext
double precision, dimension(:) :: y_ext
call get_extend1D(n*m, x, y, kcs, x_ext, y_ext)
end subroutine get_extend2D
!
!
! ==========================================================================
!>
subroutine get_extend1D(n, x, y, kcs, x_ext, y_ext)
integer :: n
double precision, dimension(n) :: x
double precision, dimension(n) :: y
integer , dimension(n) :: kcs
double precision, dimension(4) :: x_ext
double precision, dimension(4) :: y_ext
double precision :: x_min
double precision :: x_max
double precision :: x_dist
double precision :: y_min
double precision :: y_max
double precision :: y_dist
integer :: index
x_min = 1E30
x_max = -1E30
y_min = 1E30
y_max = -1E30
do index = 1, n
if(kcs(index) == 1) then
if(x_min > x(index)) then
x_min = x(index)
end if
if(x_max < x(index)) then
x_max = x(index)
end if
if(y_min > y(index)) then
y_min = y(index)
end if
if(y_max < y(index)) then
y_max = y(index)
end if
end if
enddo
x_dist = x_max - x_min
y_dist = y_max - y_min
x_min = x_min - 0.01d0*x_dist
x_max = x_max + 0.01d0*x_dist
y_min = y_min - 0.01d0*y_dist
y_max = y_max + 0.01d0*y_dist
x_ext(1) = x_min
y_ext(1) = y_min
x_ext(2) = x_min
y_ext(2) = y_max
x_ext(3) = x_max
y_ext(3) = y_max
x_ext(4) = x_max
y_ext(4) = y_min
end subroutine get_extend1D
!
!
! ==========================================================================
!>
subroutine extrapolate(n, x, y, z, kcs, n_extr, x_extr, y_extr, z_extr)
integer :: n
double precision, dimension(n) :: x
double precision, dimension(n) :: y
double precision, dimension(n) :: z
integer , dimension(n) :: kcs
integer :: n_extr
double precision, dimension(n_extr), target :: x_extr
double precision, dimension(n_extr), target :: y_extr
double precision, dimension(n_extr), target :: z_extr
integer :: i_extr
integer :: i_min
double precision, pointer :: x_a
double precision, pointer :: y_a
double precision, pointer :: z_a
double precision :: dist_min
dist_min = 1E30
i_min = 0
do i_extr = 1, n_extr
x_a => x_extr(i_extr)
y_a => y_extr(i_extr)
z_a => z_extr(i_extr)
call find_nearest(n, x, y, z, kcs, x_a, y_a, i_min, dist_min)
z_a = z(i_min)
enddo
end subroutine extrapolate
!
!
! ==========================================================================
!>
subroutine find_nearest2D(n, m, x, y, kcs, x_a, y_a, n_min, m_min, dist_min)
use precision
integer :: n
integer :: m
double precision, dimension(:,:) :: x
double precision, dimension(:,:) :: y
integer , dimension(:,:) :: kcs
integer :: n_min
integer :: m_min
integer :: i_min
double precision :: x_a
double precision :: y_a
double precision :: dist_min
call find_nearest1D(n*m, x, y, kcs, x_a, y_a, i_min, dist_min)
m_min = i_min/n
n_min = i_min - (m_min * n)
m_min = m_min + 1
end subroutine find_nearest2D
!
!
! ==========================================================================
!>
subroutine find_nearest2D_missing_value(n, m, x, y, z, kcs, x_a, y_a, n_min, m_min, dist_min)
use precision
integer :: n
integer :: m
double precision, dimension(:,:) :: x
double precision, dimension(:,:) :: y
double precision, dimension(:,:) :: z
integer , dimension(:,:) :: kcs
integer :: n_min
integer :: m_min
integer :: i_min
double precision :: x_a
double precision :: y_a
double precision :: dist_min
call find_nearest1D_missing_value(n*m, x, y, z, kcs, x_a, y_a, i_min, dist_min)
m_min = i_min/n
n_min = i_min - (m_min * n)
m_min = m_min + 1
end subroutine find_nearest2D_missing_value
!
!
! ==========================================================================
!>
subroutine find_nearest1D(n, x, y, kcs, x_a, y_a, i_min, dist_min)
use precision
integer :: n
double precision, dimension(n) :: x
double precision, dimension(n) :: y
integer , dimension(n) :: kcs
integer :: i
integer :: i_min
double precision :: x_a
double precision :: y_a
double precision :: dist
double precision :: dist_min
dist_min = 1E30
i_min = 0
do i = 1, n
if(kcs(i) == 1) then
dist = (x(i)-x_a)**2 + (y(i)-y_a)**2
if(dist < dist_min) then
dist_min = dist
i_min = i
end if
end if
enddo
dist_min = sqrt(dist_min)
end subroutine find_nearest1D
!
!
! ==========================================================================
!>
subroutine find_nearest1D_missing_value(n, x, y, z, kcs, x_a, y_a, i_min, dist_min)
use precision
integer :: n
double precision, dimension(n) :: x
double precision, dimension(n) :: y
double precision, dimension(n) :: z
integer , dimension(n) :: kcs
integer :: i
integer :: i_min
double precision :: x_a
double precision :: y_a
double precision :: dist
double precision :: dist_min
dist_min = 1E30
i_min = 0
do i = 1, n
if(kcs(i) == 1) then
dist = (x(i)-x_a)**2 + (y(i)-y_a)**2
if((dist < dist_min).and.(z(i)/=-999d0)) then
dist_min = dist
i_min = i
end if
end if
enddo
dist_min = sqrt(dist_min)
end subroutine find_nearest1D_missing_value
!
!
! ==========================================================================
!>
subroutine xxpolyint( xs, ys, zs ,kcs, ns, & ! interpolate in a polyline like way
x , y ,z ,kc , kx , mnx, jintp, xyen, indxn, wfn)
implicit none
! Global variables
integer, intent(in) :: ns !< Dimension of polygon OR LINE BOUNDARY
double precision, dimension(:), intent(in) :: xs !< polyline point coordinates
double precision, dimension(:), intent(in) :: ys
double precision, dimension(:), intent(in) :: zs !< Values at all points. Dimension: ns*kx
integer, dimension(:), intent(in) :: kcs !< polyline mask
integer, intent(in) :: mnx !< Dimension of target points
integer, intent(in) :: kx !< #values at each point (vectormax)
double precision, dimension(:), intent(in) :: x !< Grid points (where to interpolate to)
double precision, dimension(:), intent(in) :: y
double precision, dimension(kx*mnx), intent(out) :: z !< Output array for interpolated values. Dimension: mnx*kx
integer , dimension(:), intent(in) :: kc !< Target (grid) points mask
integer, intent(in) :: jintp !< (Re-)interpolate if 1 (otherwise use index weights)
double precision, dimension(:,:), intent(in) :: xyen !< cellsize / tol
integer, dimension(:,:), intent(inout), optional :: indxn !< pli segment is identified by its first node nr.
double precision, dimension(:,:), intent(inout), optional :: wfn !< If present, get weight index and factor
! locals
double precision:: wL, wR
integer :: m, k, kL, kR, jgetw
jgetw = 0 ! niets met gewichten, doe interpolatie
if ( present(indxn) .and. jintp .eq. 1) jgetw = 1 ! haal gewichten doe interpolatie , gebruik gewichten
if ( present(indxn) .and. jintp .eq. 0) jgetw = 2 ! doe interpolatie , gebruik gewichten
do m = 1, mnx
if (jgetw .le. 1) then
!call polyindexweight( x(m), y(m), xs, ys, kcs, ns, xyen(:,m), k1, rl) ! interpolate in a polyline like way
call polyindexweight( x(m), y(m), xyen(1,m), xyen(2,m), xs, ys, kcs, ns, kL, wL, kR, wR) ! interpolate in a polyline like way
!call findtri_indices_weights (x(n),y( n), xs, ys, ns, zp, indxp) ! zoeken bij 0 en 1
if (jgetw .eq. 1) then ! zetten bij 1
indxn(1,m) = kL
wfn(1,m) = wL
indxn(2,m) = kR
wfn(2,m) = wR
end if
elseif (jgetw .eq. 2) then ! halen bij 2, je hoeft niet te zoeken
kL = indxn(1,m)
wL = wfn(1,m)
kR = indxn(2,m)
wR = wfn(2,m)
end if
! Now do the actual interpolation of data zs -> z
if (kL > 0) then
if (kR > 0) then
do k = 1,kx
z(kx*(m-1)+k) = wL*zs(kx*(kL-1)+k) + wR*zs(kx*(kR-1)+k)
end do
else ! Just left point
do k = 1,kx
z(kx*(m-1)+k) = wL*zs(kx*(kL-1)+k)
end do
end if
else if (kR > 0) then
do k = 1,kx
z(kx*(m-1)+k) = wR*zs(kx*(kR-1)+k)
end do
end if
enddo
end subroutine xxpolyint
!
! ==========================================================================
!>
!subroutine polyindexweight( xe, ye, xs, ys, kcs, ns, xyen, k1, rl) ! interpolate in a polyline like way
!
! ! Global variables
! integer , intent(in) :: ns ! Dimension of polygon OR LINE BOUNDARY
! double precision, dimension(:), intent(in) :: xs ! polygon
! double precision, dimension(:), intent(in) :: ys
! integer, dimension(:), intent(in) :: kcs ! polygon mask
! double precision :: xyen(:)
! double precision :: xe, ye, rl
!
!
! integer :: ja1, ja2, k, km, k1, k2
! double precision:: x1,x2,y1,y2,dis,xn,yn,dx,dy
! double precision:: dism, dis1, dis2, rl1, rl2, dbdistance
!
!
! dism = 1e30
! do k = 1, ns
! dis = DbdISTANCE( Xe,Ye,XS(K),YS(K) )
! if (dis < dism) then
! dism = dis
! km = k
! end if
! enddo
!
! k1 = 0
!
! if (km == 1) then
! x1 = xs(km ); y1 = ys(km )
! x2 = xs(km+1); y2 = ys(km+1)
! call LINEDISQ(Xe,Ye,X1,Y1,X2,Y2,JA1,DIS1,XN,YN,RL)
! if (ja1 == 1) then
! if (dis1 < rdis) k1 = km
! end if
! else if (km == ns) then
! x1 = xs(km-1); y1 = ys(km-1)
! x2 = xs(km ); y2 = ys(km )
! call LINEDISQ(Xe,Ye,X1,Y1,X2,Y2,JA1,DIS1,XN,YN,RL)
! if (ja1 == 1) then
! if (dis1 < rdis) k1 = km-1
! end if
! else
! x1 = xs(km-1); y1 = ys(km-1)
! x2 = xs(km) ; y2 = ys(km)
! call LINEDISQ(Xe,Ye,X1,Y1,X2,Y2,JA1,DIS1,XN,YN,RL1)
! x1 = xs(km) ; y1 = ys(km)
! x2 = xs(km+1); y2 = ys(km+1)
! call LINEDISQ(Xe,Ye,X1,Y1,X2,Y2,JA2,DIS2,XN,YN,RL2)
! if (ja1 == 1) then ! if on line 1
! if (dis1 < rdis) then
! k1 = km-1 ; rl = rl1
! end if
! else if (ja2 == 1) then
! if (dis2 < rdis) then
! k1 = km ; rl = rl2
! end if
! else ! niet op een van beiden, maar wel in de buurt, uitwerken. Nu dus alleen convexe randen
! end if
! end if
!
!end subroutine polyindexweight
!
!
! ==========================================================================
!>
!> Selects the index of the polyline segment that intersects with line e--en
!! with the intersection closest to point e.
!! The search range is thus from e to en, and not a distance rdis as before.
!! The normal direction is now
!! defined by e--en and not normal to the polyline. Also, *all* polyline
!! segments are checked, not the closest based on dbdistance of pli points.
subroutine polyindexweight( xe, ye, xen, yen, xs, ys, kcs, ns, kL, wL, kR, wR)
USE M_sferic
! Global variables
integer , intent(in) :: ns !< Dimension of polygon OR LINE BOUNDARY
double precision, intent(in) :: xs(:) !< polygon
double precision, intent(in) :: ys(:)
integer , intent(in) :: kcs(:) !< polygon mask
double precision, intent(in) :: xe, ye !
double precision, intent(in) :: xen, yen !< in input uitstekers, on output SL and CRP
integer , intent(out) :: kL !< Index of left nearest polyline point (with kcs==1!)
double precision, intent(out) :: wL !< Relative weight of left nearest polyline point.
integer , intent(out) :: kR !< Index of right nearest polyline point (with kcs==1!)
double precision, intent(out) :: wR !< Relative weight of right nearest polyline point.
double precision, external :: dbdistance
integer :: k, km, klastvalid, JACROS
double precision :: dis, disM, disL, disR !, rl1, rl2,
double precision :: SL,SM,SMM,SLM,XCR,YCR,CRP,CRPM,DEPS
DISM = huge(DISM)
kL = 0 ! Default: No valid point found
kR = 0 ! idem
wL = 0d0
wR = 0d0
km = 0
crpm = 0
disL = 0d0
disR = 0d0
DEPS = 1d-3
do k = 1, ns-1
call CROSS(xe, ye, xen, yen, xs(k), ys(k), xs(k+1), ys(k+1), JACROS,SL,SM,XCR,YCR,CRP)
if (SL >= 0d0 .AND. SL <= 1D0 .AND. SM > -DEPS .AND. SM < 1.0D0+DEPS) then ! instead of jacros==1, solves firmijn's problem
DIS = DBDISTANCE(XE,YE, XCR, YCR)
if (DIS < DISM) then ! Found a better intersection point
DISM = DIS
km = k
SMM = SM
SLM = SL
CRPM = CRP
end if
end if
enddo
if (km > 0) then
dis = dbdistance(xs(km), ys(km), xs(km+1), ys(km+1)) ! Length of this segment.
! Find nearest valid polyline point left of the intersection (i.e.: kcs(kL) == 1)
disL = SMM*dis
do k = km,1,-1
if (kcs(k) == 1) then
kL = k
exit ! Valid point on the left (distance was already included in disL)
else if (k > 1) then
disL = disL + dbdistance(xs(k-1), ys(k-1), xs(k), ys(k)) ! Add entire length of this segment.
end if
end do
! Find nearest valid polyline point right of the intersection (i.e.: kcs(kR) == 1)
disR = (1d0-SMM)*dis
do k = km+1,ns
if (kcs(k) == 1) then
kR = k
exit ! Valid point on the left (distance was already included in disL)
else if (k < ns) then
disR = disR + dbdistance(xs(k), ys(k), xs(k+1), ys(k+1)) ! Add entire length of this segment.
end if
end do
end if
if (kL /= 0 .and. kR /= 0) then
wL = disR/(disL+disR)
wR = 1d0-wL
else if (kL /= 0) then
wL = 1d0
else if (kR /= 0) then
wR = 1d0
end if
end subroutine polyindexweight
!
!
! ==========================================================================
!>
SUBROUTINE LINEDISq(X3,Y3,X1,Y1,X2,Y2,JA,DIS,XN,YN,rl) ! = dlinesdis2
integer :: ja
DOUBLE PRECISION :: X1,Y1,X2,Y2,X3,Y3,DIS,XN,YN
DOUBLE PRECISION :: R2,RL,X21,Y21,X31,Y31,getdx,getdy,dbdistance
! korste afstand tot lijnelement tussen eindpunten
JA = 0
!X21 = getdx(x1,y1,x2,y2)
!Y21 = getdy(x1,y1,x2,y2)
call getdxdy(x1,y1,x2,y2,x21,y21)
!X31 = getdx(x1,y1,x3,y3)
!Y31 = getdy(x1,y1,x3,y3)
call getdxdy(x1,y1,x3,y3,x31,y31)
R2 = dbdistance(x2,y2,x1,y1)
R2 = R2*R2
IF (R2 .NE. 0) THEN
RL = (X31*X21 + Y31*Y21) / R2
IF (0d0 .LE. RL .AND. RL .LE. 1d0) then
JA = 1
end if
XN = X1 + RL*(x2-x1)
YN = Y1 + RL*(y2-y1)
DIS = dbdistance(x3,y3,xn,yn)
end if
RETURN
END subroutine LINEDISq
end module timespace_triangle ! met leading dimensions 3 of 4
!
!
!
! ==========================================================================
! ==========================================================================
! ==========================================================================
!>
module timespace
!!--copyright-------------------------------------------------------------------
! Copyright (c) 2006, WL | Delft Hydraulics. All rights reserved.
!!--disclaimer------------------------------------------------------------------
! This code is part of the Delft3D software system. WL|Delft Hydraulics has
! developed c.q. manufactured this code to its best ability and according to the
! state of the art. Nevertheless, there is no express or implied warranty as to
! this software whether tangible or intangible. In particular, there is no
! express or implied warranty as to the fitness for a particular purpose of this
! software, whether tangible or intangible. The intellectual property rights
! related to this software code remain with WL|Delft Hydraulics at all times.
! For details on the licensing agreement, we refer to the Delft3D software
! license and any modifications to this license, if applicable. These documents
! are available upon request.
!!--version information---------------------------------------------------------
! $Author$
! $Date$
! $Revision$
!!--description-----------------------------------------------------------------
!
! Read time series in five possible formats:
! uniform : Delft3D-FLOW format: time, uniform windspeed, direction and pressure
! space varying : Delft3D-FLOW format: time and fields of patm, windx, windy
! on Delft3D-FLOW m,n grid
! arcinfo : time and fields on own equidistant grid
! spiderweb : time and fields of patm, windspeed, direction op spiderweb grid
! curvi : time and fields on own curvilinear grid
!
! Main calls from Delft3D-FLOW:
! readmd -> rdmeteo:
! initmeteo : allocate meteo structure for this domain
! adddataprovider : allocate and initialized an input quantity
! with specified format
! checkmeteo : check whether input is available for the complete
! time interval
! trisol -> incmeteo:
! meteoupdate : prepare meteo data for the current time
! getmeteoval : return meteo data for the current time and position
! use optional m and n parameters to speed up in case of curvi
! getspiderval : same as getmeteoval for spiderweb data
!
! gdp_dealloc:
! deallocmeteo
!
! Additional calls:
! getmeteoerror : returns a string containing an error message
! to be used in case success = false for a main call
! meteogetpaver : returns the average atmospheric pressure read
! meteogetpcorr : returns whether pressure correction is switched on on
! the boundaries
!
!!--pseudo code and references--------------------------------------------------
!
! Stef.Hummel@WlDelft.nl
! Herman.Kernkamp@WlDelft.nl
! Adri.Mourits@WlDelft.nl
!
!!--declarations----------------------------------------------------------------
use precision
use timespace_data
use timespace_triangle
implicit none
contains
!
!
! ==========================================================================
!>
!> this function selects points (kc = 1) that can receive data from the provider in file =filename
!! All points have an allowable 'search range', defined by a line from x,y
!! to xyen(1,) to xyen(2,). Generally, the points in xyen are endpoints of
!! rrtol times a perpendicular vector to edge links.
subroutine selectelset( filename, filetype, x, y, xyen, kc, mnx, ki, num, usemask, rrtolrel)
use MessageHandling
implicit none
! arguments
integer , intent(in) :: mnx !< dimension of quantity
double precision, intent(in) :: x(:) !< x of elset of all possible points in model
double precision, intent(in) :: y(:) !< y of elset
double precision, intent(in) :: xyen(:,:) !< Points on opposite edges of elementset
integer , intent(inout) :: kc(:) !< kcs of elset, allowable kandidates have 1, eg. points with less links than edges
integer , intent(out) :: ki(:) !< Returned indices of allowable points (in x/y) that fall near provided data
integer :: num !< nr of points served bij this provider
character(*), intent(in) :: filename ! file name for meteo data file
integer , intent(in) :: filetype ! spw, arcinfo, uniuvp etc
logical, intent(in) :: usemask !< Whether to use the mask array kc, or not (allows you to keep kc, but disable it for certain quantities, for example salinitybnd).
double precision, intent(in), optional :: rrtolrel !< Optional, a more strict rrtolerance value than the global rrtol. selectelset will succeed if cross SL value <= rrtolrel
! locals
double precision, allocatable :: xs (:) ! temporary array to hold polygon
double precision, allocatable :: ys (:) !
integer , allocatable :: kcs(:) !
double precision :: wL, wR
integer :: kL, kR, minp, ns, m
integer :: JACROS
double precision :: SL,SM,XCR,YCR,CRP
num = 0
ki = 0
if (filetype == poly_tim) then
if (.not. allocated(xs)) then
call realloc(xs,10000)
call realloc(ys,10000)
call realloc(kcs,10000)
end if
kcs = 1 ! todo make this safe
call oldfil(minp, filename)
call read1polylin(minp,xs,ys,ns)
do m = 1,mnx
if (iabs(kc(m)) == 1) then ! point is a possible candidate for a line boundary
call polyindexweight( x(m), y(m), xyen(1,m), xyen(2,m), xs, ys, kcs, ns, kL, wL, kR, wR)
if (kL > 0 .or. kR > 0) then
if (present(rrtolrel)) then
! x,y -> xyen =approx D + 2*rrtol * D
! This bnd requests a more strict tolerance than the global rrtol, namely: D + 2*rrtolb * D, so:
call CROSS(x(m), y(m), xyen(1,m), xyen(2,m), xs(kL), ys(kL), xs(kR), ys(kR), JACROS,SL,SM,XCR,YCR,CRP)
if (SL > rrtolrel) then
! More strict rrtolrel check failed, so do not accept this node.
cycle
end if
end if
if (usemask .and. kc(m) .eq. -1 ) then
errormessage = 'Boundary location already claimed; Overlap with other bnds?'
return
else
num = num + 1
ki(num) = m
if (usemask) then ! If we don't use the mask, also don't administer this opened bnd location (e.g. for salinitybnd)
kc(m) = -1 ! this tells you this point is already claimed by some bnd
end if
end if
end if
end if
enddo
write(msgbuf,'(a,a,a,i0,a)') 'boundary: ''', trim(filename), ''' opened ', num, ' cells.'
call msg_flush()
deallocate(xs, ys, kcs)
else ! en andere behandelmethodes voor andere mogelijkheden.
end if
end subroutine selectelset
!
!
! ==========================================================================
!>
subroutine selectelset_internal_links( filename, filetype, xz, yz, ln, lnx, keg, numg ) ! find links cut by polyline filetype 9
implicit none
character(len=*), intent(in) :: filename ! file name for meteo data file
integer , intent(in) :: filetype ! spw, arcinfo, uniuvp etc
double precision, intent(in) :: xz (:)
double precision, intent(in) :: yz (:)
integer , intent(in) :: ln (:,:)
integer , intent(in) :: lnx
integer , intent(out) :: keg(:)
integer :: numg
integer :: isec
integer :: minp, np, L, k1, k2, ja
double precision, allocatable :: xp(:) , yp(:)
double precision :: xa, ya, xb, yb,xm, ym, CRPM
numg = 0
if (filetype == poly_tim) then
call realloc(xp,100000)
call realloc(yp,100000)
call oldfil(minp, filename)
call read1polylin(minp,xp,yp,np)
do L = 1,lnx
k1 = ln(1,L) ; k2 = ln(2,L)
xa = xz(k1) ; ya = yz(k1)
xb = xz(k2) ; yb = yz(k2)
call CROSSPOLY(xa,ya,xb,yb,xp,yp,np,XM,YM,CRPM,JA,isec)
if (ja == 1) then
numg = numg + 1
if (crpm > 0) then
keg(numg) = -L
else
keg(numg) = L
end if
end if
enddo
deallocate(xp,yp)
end if
end subroutine selectelset_internal_links
!
!
! ==========================================================================
!> Combine a newly computed (external forcings-)value with an existing one, based on the operand type.
subroutine operate(a,b,operand)
use precision
implicit none
double precision, intent(inout) :: a !< Current value, will be updated based on b and operand.
double precision, intent(in) :: b !< New value, to be combined with existing value a.
character(len=1), intent(in) :: operand !< Operand type, valid values: 'O', 'A', '+', '*', 'X', 'N'.
! b = factor*b + offset ! todo doorplussen
if (operand == 'O') then ! Override, regardless of what was specified before
a = b
else if (operand == 'A') then ! Add, means: only if nothing was specified before
if (a == dmiss_default ) then
a = b
end if
else if (a .ne. dmiss_default) then ! algebra only if not missing
if (operand == '+') then
a = a + b
else if (operand == '*' ) then
a = a * b
else if (operand == 'X' ) then
a = max(a,b)
else if (operand == 'N' ) then
a = min(a,b)
end if
end if
end subroutine operate
!
!
! ==========================================================================
!>
function timespaceinitialfield(xu, yu, zu, nx, filename, filetype, method, operand, transformcoef, iprimpos) result(success) !
use m_kdtree2
use m_samples
use m_netw
use m_flowgeom, only : xz, yz, ln2lne, Ln, Lnx
use m_partitioninfo
use unstruc_netcdf
use m_flowexternalforcings, only: qid
use M_INTERPOLATIONSETTINGS
use m_missing
implicit none
logical :: success
integer, intent(in) :: nx
double precision, intent(in) :: xu(nx)
double precision, intent(in) :: yu(nx)
double precision, intent(out) :: zu(nx)
character(*), intent(in) :: filename ! file name for meteo data file
integer , intent(in) :: filetype ! spw, arcinfo, uniuvp etc
integer , intent(in) :: method ! time/space interpolation method
! 4 : inside polygon
! 5 : triangulation
! 6 : averaging
! 7 : index triangulation
! 8 : smoothing
! 9 : internal diffusion
character(1), intent(in) :: operand ! override, add
double precision, intent(in) :: transformcoef(:) !< Transformation coefficients
integer , intent(in) :: iprimpos ! only needed for averaging, position of primitive variables in network
! 1 = u point, cellfacemid, 2 = zeta point, cell centre, 3 = netnode
double precision, allocatable :: zh(:)
integer :: ierr
integer :: minp0, inside, npli,k, jdla, mout
double precision, allocatable :: xx(:,:), yy(:,:)
integer , allocatable :: nnn (:)
double precision, allocatable :: xxx(:), yyy(:)
integer :: n6 , L, Lk, ja, n, nn, kk, n1, n2, i
integer :: mpol
integer :: ierror
integer :: jakdtree=1
character(len=1), external :: get_dirsep
character(len=5) :: sd
logical :: file_exists
success = .false.
minp0 = 0
allocate( zh(nx) , stat=ierr)
zh = dmiss_default
call aerr( 'zh(nx)', ierr, nx)
if (filetype .ne. ncflow) then
call oldfil(minp0, filename)
end if
if (method == 4) then ! polyfil
call savepol()
call reapol(minp0, 0)
inside = -1
do k=1,nx
call dbpinpol(xu(k), yu(k), inside)
if (inside == 1) then
call operate(zu(k), transformcoef(1), operand)
zh(k) = zu(k)
end if
enddo
call restorepol()
else if (method == 5 .or. method == 6) then ! triangulation todo
if (filetype == ncflow) then
call read_flowsamples_from_netcdf(filename, qid, ierr)
else
call reasam(minp0, 0)
end if
if (method == 5) then
jdla = 1
call triinterp2(xu,yu,zh,nx,jdla)
else if (method == 6) then ! and this only applies to flow-link data
if ( transformcoef(4).ne.DMISS ) then
iav = int(transformcoef(4))
end if
if ( transformcoef(5).ne.DMISS ) then
rcel = transformcoef(5)
end if
if (iprimpos == 1) then ! primitime position = velocitypoint, cellfacemid
n6 = 4
allocate( xx(n6,lnx), yy(n6,lnx), nnn(lnx) )
do L = 1,lnx
xx(1,L) = xzw(ln(1,L))
yy(1,L) = yzw(ln(1,L))
xx(3,L) = xzw(ln(2,L))
yy(3,L) = yzw(ln(2,L))
Lk = ln2lne(L)
xx(2,L) = xk(kn(1,Lk))
yy(2,L) = yk(kn(1,Lk))
xx(4,L) = xk(kn(2,Lk))
yy(4,L) = yk(kn(2,Lk))
enddo
nnn = 4 ! array nnn
else if (iprimpos == 2) then ! primitime position = waterlevelpoint, cell centre
n6 = maxval(netcell%n)
allocate( xx(n6,nx), yy(n6,nx), nnn(nx) )
do n = 1,nx
nnn(n) = netcell(n)%n
do nn = 1, nnn(n)
xx(nn,n) = xk(netcell(n)%nod(nn))
yy(nn,n) = yk(netcell(n)%nod(nn))
enddo
enddo
else if (iprimpos == 3) then ! primitime position = netnode, cell corner
n6 = 2*maxval(nmk) ! 2: safe upper bound
allocate( xx(n6,numk), yy(n6,numk), nnn(numk), xxx(n6), yyy(n6) )
do k = 1,numk
! get the celllist
call make_dual_cell(k, n6, rcel, xxx, yyy, nnn(k))
do i=1,nnn(k)
xx(i,k) = xxx(i)
yy(i,k) = yyy(i)
enddo
enddo
DEALLOCATE(xxx,yyy)
end if
if ( jakdtree.eq.1 ) then
! initialize kdtree
call build_kdtree(treeglob,Ns,xs,ys,ierror)
if ( ierror.ne.0 ) then
! disable kdtree
call delete_kdtree2(treeglob)
jakdtree = 0
end if
end if
call averaging2(1,ns,xs,ys,zs,ipsam,xu,yu,zh,nx,xx,yy,n6,nnn,jakdtree)
deallocate(xx,yy,nnn)
iav = 1
rcel = RCEL_DEFAULT
if ( jakdtree.eq.1 ) then
call delete_kdtree2(treeglob)
end if
end if
do k=1,nx
if ( zh(k).ne.dmiss_default ) then
call operate(zu(k), zh(k), operand)
zh(k) = zu(k)
end if
end do
call delsam(0)
end if
success = .true.
call doclose(minp0)
N1 = index (trim(filename) , get_dirsep() , .true.)
! fix for Linux-prepared input on Windows
if ( N1.eq.0 ) then
N1 = index(trim(filename), char(47), .true.)
end if
sd = ''
if (jampi == 1) then
sd = '_'//trim(sdmn)
end if
N2 = INDEX (trim(filename) , '.' , .true.)
if (n2 == 0) then
n2 = len_trim(filename)
else
n2 = n2 -1
end if
call newfil(mout, 'DFM_interpreted_values_'//trim(filename(n1+1:n2))//trim(sd)//'.xyz')
do k = 1,nx
if (zh(k) .ne. dmiss_default) then
write(mout,*) xu(k), yu(k), zu(k)
end if
enddo
call doclose(mout)
if (allocated (zh) ) deallocate(zh)
end function timespaceinitialfield
!
!
! ==========================================================================
!>
function timespaceinitialfield_int(xz, yz, zz, nx, filename, filetype, method, operand, transformcoef) result(success) ! deze subroutine moet veralgemeniseerd en naar meteo module
implicit none
logical :: success
integer, intent(in) :: nx
double precision, intent(in) :: xz(nx)
double precision, intent(in) :: yz(nx)
integer , intent(out) :: zz(nx)
character(*), intent(in) :: filename ! file name for meteo data file
integer , intent(in) :: filetype ! spw, arcinfo, uniuvp etc
integer , intent(in) :: method ! time/space interpolation method
character(1), intent(in) :: operand ! file name for meteo data file
double precision, intent(in) :: transformcoef(:) !< Transformation coefficients
double precision, allocatable :: xpli(:), ypli(:)
integer :: maxpli = 10000
integer :: minp0, inside, npli,k, jdla = 1
success = .false.
call oldfil(minp0, filename)
if (filetype == inside_polygon) then ! polyfil
allocate(xpli(maxpli), ypli(maxpli))
call read1polylin(minp0,xpli,ypli,npli)
do k=1,nx
call pinpok(xz(k), yz(k), npli, xpli, ypli, inside)
if (inside == 1) then
if (operand == '+') then
zz(k) = zz(k) + transformcoef(1)
else
zz(k) = transformcoef(1)
end if
end if
enddo
deallocate(xpli, ypli)
else if (filetype == arcinfo) then ! arcinfo bilinear todo
else if (filetype == triangulation) then ! triangulation todo
end if
success = .true.
end function timespaceinitialfield_int
end module timespace