!----- 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: unstruc_bmi.F90 43233 2015-11-24 21:17:21Z baart_f $ ! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/unstruc_bmi.F90 $ #ifdef HAVE_CONFIG_H #include "config.h" #endif module bmi use iso_c_binding use unstruc_api use unstruc_display, only: jaGUI ! this should be removed when jaGUI = 0 by default use m_partitioninfo use m_flow use m_waves use m_ship use m_wind use m_flowgeom use network_data use m_flowexternalforcings use m_partitioninfo use m_sobekdfm use m_transport implicit none ! Define some global constants ! We can't return an array of pointers with an undertimened length so ! we fix it to something, in this case 100. character(*), parameter :: PREFIX = "FM" !DEC$ ATTRIBUTES DLLEXPORT :: PREFIX integer(c_int) :: MAXNAMES = 100 integer(c_int), BIND(C, name="MAXSTRLEN") :: MAXSTRLEN = 1024 !DEC$ ATTRIBUTES DLLEXPORT :: MAXSTRLEN integer(c_int), BIND(C, name="MAXDIMS") :: MAXDIMS = 6 !DEC$ ATTRIBUTES DLLEXPORT :: MAXDIMS ! TODO, get rid of these custom variables.... real(c_double), target, allocatable, save :: uabs(:) real(c_double), target, allocatable, save :: froude(:) real(c_double), target, allocatable, save :: sed1(:) real(c_double), target, allocatable, save :: const_t(:,:) !< Placeholder array for returning constituents (transposed: dim(numconst,ndkx)) integer, private :: iconst integer, external :: findname integer(c_int), parameter :: var_count_compound = 6 ! pumps, gates, weirs, sourcesinks, observations, crosssections ! TODO: AvD: temp, as long as this is not templated contains ! Control subroutine main() bind(C, name="main") !DEC$ ATTRIBUTES DLLEXPORT :: main ! Somehow intel fortran compiler expects a main routine in the dll, it is required since interactor is used (and win calls) end subroutine main !> The initialize() function accepts a string argument that !! gives the name (and path) of its "main input file", called !! a configuration file. This function should perform all tasks !! that are to take place before entering the model's time loop. integer(c_int) function initialize(c_config_file) result(c_iresult) bind(C, name="initialize") !DEC$ ATTRIBUTES DLLEXPORT :: initialize use iso_c_binding, only: c_char use unstruc_model use unstruc_files use m_partitioninfo #ifdef HAVE_MPI use mpi #endif character(kind=c_char),intent(in) :: c_config_file(MAXSTRLEN) character(len=strlen(c_config_file)) :: config_file ! Extra local variables integer :: inerr ! number of the initialisation error logical :: mpi_initd integer(c_int), target, allocatable, save :: x(:,:) type(c_ptr) :: xptr integer :: i,j,k c_iresult = 0 ! TODO: is this return value BMI-compliant? #ifdef HAVE_MPI call mpi_initialized(mpi_initd, inerr) if (.not. mpi_initd) then ja_mpi_init_by_fm = 1 call mpi_init(inerr) else ja_mpi_init_by_fm = 0 end if call mpi_comm_rank(DFM_COMM_DFMWORLD,my_rank,inerr) call mpi_comm_size(DFM_COMM_DFMWORLD,numranks,inerr) if ( numranks.le.1 ) then jampi = 0 end if ! make domain number string as soon as possible write(sdmn, '(I4.4)') my_rank open (unit = 2, file = "mpi_"//sdmn//".out") write(2, *) 'my_rank =', my_rank close(2) #else numranks=1 #endif ! do this until default has changed jaGUI = 0 ! Store the name config_file = char_array_to_string(c_config_file, strlen(c_config_file)) ! Now we can initialize with the config_file ! It would be great if there would be 1 init function that takes an optional configuration file name ! TODO: check why these are needed to avoid a segfault KNX = 8 MXB = 10 MAXLAN = 500 MAXPOL = MAXLAN ! call start() ! call resetFullFlowModel() !call loadmodel(config_file) call inidia(config_file(1:(len_trim(config_file)-len_trim('.mdu')))) CALL INIDAT() call api_loadmodel(config_file) c_iresult = flowinit() time_user = tstart_user ! Just terminate if we get an error.... ! if (inerr > 0) stop ! initialize = 0 end function initialize subroutine get_version_string(c_version_string) bind(C, name="get_version_string") !DEC$ ATTRIBUTES DLLEXPORT :: get_version_string use iso_c_binding, only: c_char use unstruc_version_module, only: unstruc_version character(kind=c_char), intent(out) :: c_version_string(MAXSTRLEN) character(len=MAXSTRLEN) :: name name = unstruc_version c_version_string = string_to_char_array(trim(name), len(trim(name))) end subroutine get_version_string !! USER TIMESTEP !! !! NOTE: for greater flexibility the time step routines below are split into three !! init_, run_ and finalize_ routines. Make sure to call all of them in proper order. !> Initializes a 'user' timestep (outer time level, for meteo updates and I/O). !! !! The actual time-stepping will have to be done by dfm_init/run/finalize_computational_timestep. !! @see dfm_init_computational_timestep function dfm_init_user_timestep(timetarget) result(iresult) bind(C, name="dfm_init_user_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_init_user_timestep use m_flowtimes, only: time_user, dt_user use dfm_error use MessageHandling real(c_double), intent(in) :: timetarget !< Target time. For now, this is assumed to be equal to upcoming next user time. If not, errorstatus returned. integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. character*(MAXSTRLEN) :: msg iresult = DFM_NOERR if (timetarget - time_user /= dt_user) then ! We don't support yet a changing dt_user (may affect input/output, meteo, etc.) iresult = DFM_INVALIDTARGETTIME write(msg,'(a,f8.3,a,f8.3,a)') 'Mismatch of requested step (',timetarget - time_user,') and the specified user timestep (',dt_user,').' call mess(LEVEL_INFO,msg) goto 888 end if call flow_init_usertimestep(iresult) return 888 continue end function dfm_init_user_timestep !> Runs the actual 'user' timestep (outer time level, preceded by meteo updates and to be followed by I/O). function dfm_run_user_timestep() result(iresult) bind(C, name="dfm_run_user_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_run_user_timestep integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. integer :: dummykey call flow_run_usertimestep(dummykey, iresult) end function dfm_run_user_timestep !> Finalizes a 'user' timestep (outer time level, typically timers and I/O). function dfm_finalize_user_timestep() result(iresult) bind(C, name="dfm_finalize_user_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_finalize_user_timestep integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. call flow_finalize_usertimestep(iresult) end function dfm_finalize_user_timestep !! COMPUTATIONAL TIMESTEP !! !! NOTE: for greater flexibility the time step routines below are split into three !! init_, run_ and finalize_ routines. Make sure to call all of them in proper order. !> Initializes a single computational timestep. function dfm_init_computational_timestep(timetarget, dtpredict) result(iresult) bind(C, name="dfm_init_computational_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_init_computational_timestep use m_flowtimes, only: dts use dfm_error use MessageHandling real(c_double), intent(in) :: timetarget !< Target time, resulting timestep may (will generally) be smaller. For now, this is assumed to be equal to upcoming next user time. If not, errorstatus returned. real(c_double), intent(out) :: dtpredict !< The predicted computational timestep, based on stability criteria. Pass this value (or smaller) on to run_computational_timestep. integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. character*(MAXSTRLEN) :: msg iresult = DFM_NOERR if (timetarget /= time_user) then ! We don't yet support another targettime than upcoming time_user iresult = DFM_INVALIDTARGETTIME write(msg,'(a,f8.3,a,f8.3,a)') 'Mismatch of requested step (',timetarget - time_user,') and the specified user timestep (',dt_user,').' call mess(LEVEL_INFO,msg) goto 888 end if if (time0+dts > tstop_user) then dts = tstop_user - time0 endif call flow_init_single_timestep(iresult) dtpredict = dts return 888 continue end function dfm_init_computational_timestep !> Runs the actual single computational timestep. !! !! Should be preceded by an init_computational_timestep, and followed by finalize_computational_timestep. function dfm_run_computational_timestep(dtactual) result(iresult) bind(C, name="dfm_run_computational_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_run_computational_timestep use m_flowtimes, only: dts real(c_double), intent(inout) :: dtactual !< The actual timestep to be used (may have changed upon return, if time setbacks have occurred) integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. integer :: dummykey dts = dtactual call flow_run_single_timestep(dummykey, iresult) ! TODO: AvD/ JN: do we need a nonzero iresult on SETBACK/other errors? YES dtactual = dts end function dfm_run_computational_timestep !> Finalizes a single computational timestep. !! !! Should be called immediately after a run_computational_timestep. function dfm_finalize_computational_timestep() result(iresult) bind(C, name="dfm_finalize_computational_timestep") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_finalize_computational_timestep integer(c_int) :: iresult !< Result status, DFM_NOERR(=0) if successful. call flow_finalize_single_timestep(iresult) end function dfm_finalize_computational_timestep integer function update(dt) bind(C, name="update") !DEC$ ATTRIBUTES DLLEXPORT :: update use messagehandling use iso_c_binding, only: c_double real(c_double), value, intent(in) :: dt integer :: key, ierr ! The time loop seems to be located in unstruc->flow ! It is important that we can simulate up to a time set from the outside ! We might have to set time_user or dt_user call flow_run_sometimesteps(dt, ierr) update = ierr end function update subroutine update_until(t) bind(C, name="update_until") use iso_c_binding, only: c_double real(c_double),intent(in) :: t ! Calls update(t-tnow) end subroutine update_until integer function finalize() bind(C, name="finalize") !DEC$ ATTRIBUTES DLLEXPORT :: finalize use m_partitioninfo if ( jampi.eq.1 ) then ! finalize before exit call partition_finalize() end if ! clean up, close files finalize = 0 end function finalize ! Attributes subroutine get_start_time(t) bind(C, name="get_start_time") !DEC$ ATTRIBUTES DLLEXPORT :: get_start_time use m_flowtimes use iso_c_binding, only: c_double real(c_double), intent(out) :: t t = tstart_user end subroutine get_start_time subroutine get_end_time(t) bind(C, name="get_end_time") !DEC$ ATTRIBUTES DLLEXPORT :: get_end_time use m_flowtimes use iso_c_binding, only: c_double real(c_double), intent(out) :: t t = tstop_user end subroutine get_end_time subroutine get_time_step(dt) bind(C, name="get_time_step") !DEC$ ATTRIBUTES DLLEXPORT :: get_time_step use iso_c_binding, only: c_double real(c_double), intent(out) :: dt dt = dt_user end subroutine get_time_step subroutine get_current_time(t) bind(C, name="get_current_time") !DEC$ ATTRIBUTES DLLEXPORT :: get_current_time use iso_c_binding, only: c_double real(c_double), intent(out) :: t t = time1 end subroutine get_current_time subroutine dfm_compute_1d2d_coefficients() bind(C, name="dfm_compute_1d2d_coefficients") !DEC$ ATTRIBUTES DLLEXPORT :: dfm_compute_1d2d_coefficients call compute_1d2d_coefficients() end subroutine dfm_compute_1d2d_coefficients subroutine get_time_units(unit) bind(C, name="get_time_units") ! returns unit string for model time, e.g. ‘days since 1970-01-01' character(kind=c_char), intent(in) :: unit(*) end subroutine get_time_units subroutine get_n_attributes(n) bind(C, name="get_n_attributes") use iso_c_binding, only: c_char, c_ptr, c_loc, c_int, c_null_char integer(c_int), intent(out) :: n n = 0 ! expose attributes, properties, parameters that are settable by users or can be used for calibration ! This is not part of bmi, but it allows us to easily loop over attributes: ! get_n_attributes(n) ! do i=1,n ! get_attribute_name(i, name) ! ! dosomething with name ! end end subroutine get_n_attributes subroutine get_attribute_name(i, c_att_name) bind(C, name="get_attribute_name") use iso_c_binding, only: c_char, c_ptr, c_loc, c_int, c_null_char integer(c_int), intent(in) :: i character(kind=c_char), intent(out) :: c_att_name(MAXSTRLEN) character(len=MAXSTRLEN) :: name name = 'some attribute name' c_att_name = string_to_char_array(trim(name), len(trim(name))) ! get name of attribute i end subroutine get_attribute_name subroutine get_attribute_type(c_att_name, c_type) bind(C, name="get_attribute_type") use iso_c_binding, only: c_char, c_ptr, c_null_char character(kind=c_char), intent(in) :: c_att_name(*) character(kind=c_char), intent(out) :: c_type(MAXSTRLEN) ! Use one of the following types ! BMI datatype C datatype NumPy datatype ! BMI_STRING char* S< ! BMI_INT int int16 ! BMI_DOUBLE double float64 ! The fortran name of the attribute name character(len=strlen(c_att_name)) :: att_name character(len=MAXSTRLEN) :: type_name ! Store the name att_name = char_array_to_string(c_att_name, strlen(c_att_name)) ! Lookup the type of the variable type_name = 'BMI_INT' c_type = string_to_char_array(trim(type_name), len(trim(type_name))) end subroutine get_attribute_type subroutine get_double_attribute(c_att_name, value) bind(C, name="get_double_attribute") use iso_c_binding, only: c_double, c_char character(kind=c_char), intent(in) :: c_att_name(*) real(c_double), intent(out) :: value ! The fortran name of the attribute name character(len=strlen(c_att_name)) :: att_name ! Store the name att_name = char_array_to_string(c_att_name, strlen(c_att_name)) ! Look up the value of att_name value = -1.0d0 end subroutine get_double_attribute subroutine get_int_attribute(c_att_name, value) bind(C, name="get_int_attribute") use iso_c_binding, only: c_double, c_char character(kind=c_char), intent(in) :: c_att_name(*) integer(c_int), intent(out) :: value ! The fortran name of the attribute name character(len=strlen(c_att_name)) :: att_name ! Store the name att_name = char_array_to_string(c_att_name, strlen(c_att_name)) ! Look up the value of att_name value = -1 end subroutine get_int_attribute subroutine get_string_attribute(c_att_name, c_value) bind(C, name="get_string_attribute") !DEC$ ATTRIBUTES DLLEXPORT :: get_string_attribute use iso_c_binding, only: c_char use m_flowtimes character(kind=c_char), intent(in) :: c_att_name(*) character(kind=c_char), intent(out) :: c_value(MAXSTRLEN) ! The fortran name of the attribute name character(len=strlen(c_att_name)) :: att_name character(len=MAXSTRLEN) :: value ! Store the name att_name = char_array_to_string(c_att_name, strlen(c_att_name)) ! Lookup the type of the variable select case(att_name) case ("refdat") value = refdat case default value = '' end select c_value = string_to_char_array(trim(value), len(trim(value))) end subroutine get_string_attribute ! Variables subroutine get_var_count(c_var_count) bind(C, name="get_var_count") ! non-BMI !DEC$ ATTRIBUTES DLLEXPORT :: get_var_count use iso_c_binding, only: c_ptr, c_int integer(c_int), intent(out) :: c_var_count include "bmi_get_var_count.inc" ! plus extra local vars c_var_count = var_count + var_count_compound + numconst + 5 end subroutine get_var_count subroutine get_var_name(var_index, c_var_name) bind(C, name="get_var_name") !DEC$ ATTRIBUTES DLLEXPORT :: get_var_name use iso_c_binding, only: c_char, c_int integer(kind=c_int), value, intent(in) :: var_index character(kind=c_char), intent(out) :: c_var_name(MAXSTRLEN) integer :: iloc character(MAXSTRLEN) :: var_name include "bmi_get_var_count.inc" include "bmi_get_var_name.inc" ! Var ordering: ! 1:var_count = autogenerated ! var_count+1:var_count+var_count_compound = compound types (pumps, etc.) ! var_count+var_count_compound+1:var_count+var_count_compound+numconst = constituents (tracers, salt, etc.) ! var_count+var_count_compound+numconst+1:end = remaining hand-crafted variables select case(var_index) case(var_count + 1) var_name = "pumps" case(var_count + 2) var_name = "weirs" case(var_count + 3) var_name = "gates" case(var_count + 4) var_name = "sourcesinks" case(var_count + 5) var_name = "observations" case(var_count + 6) var_name = "crosssections" end select if (var_index > var_count + var_count_compound .and. var_index <= var_count + var_count_compound + numconst) then var_name = const_names(var_index - var_count - var_count_compound) end if iloc = var_index - (var_count + var_count_compound + numconst) select case (iloc) case (1) var_name = "uabs" case (2) var_name = "sed1" case (3) var_name = "netelemnode" case(4) var_name = "flowelemnode" case(5) var_name = "kn" end select c_var_name = string_to_char_array(trim(var_name), len(trim(var_name))) end subroutine get_var_name subroutine get_var_names(names) bind(C, name="get_var_names") use iso_c_binding, only: c_char, c_ptr character(kind=c_char), dimension(MAXNAMES), intent(out) :: names(*) ! I can't get this to work..... ! http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-bridges-using-iso-c-binding ! The way we do it is to use a C_PTR array to point to strings. For example: ! CHARACTER(LEN=100), DIMENSION(numStrings), TARGET :: stringArray ! TYPE(C_PTR), DIMENSION(numStrings) :: stringPtrs ! then we set our strings in stringArray, remembering to null-terminate them such as: ! DO ns = 1, numStrings ! stringArray(ns) = "My String"//C_NULL_CHAR ! stringPtrs(ns) = C_LOC(stringArray(ns)) ! END DO ! and pass stringPtrs to the C function. ! The C function has the interface: ! void stringFunc(int *numStrings, char **stringArray) { ! int i; ! for(i=0;i<*numStrings;++i) { ! printf("%s\n",stringArray[i]); ! } ! } end subroutine get_var_names subroutine get_input_var_names(names) bind(C, name="get_input_var_names") use iso_c_binding, only: c_char, c_ptr character(kind=c_char), dimension(MAXNAMES), intent(out) :: names(*) !type(c_ptr), dimension(:) :: names end subroutine get_input_var_names subroutine get_output_var_names(names) bind(C, name="get_output_var_names") use iso_c_binding, only: c_char character(kind=c_char), dimension(MAXNAMES), intent(out) :: names(*) end subroutine get_output_var_names subroutine get_var_type(c_var_name, c_type) bind(C, name="get_var_type") !DEC$ ATTRIBUTES DLLEXPORT :: get_var_type character(kind=c_char), intent(in) :: c_var_name(*) character(kind=c_char), intent(out) :: c_type(MAXSTRLEN) character(len=MAXSTRLEN) :: type_name, var_name ! Use one of the following types ! BMI datatype C datatype NumPy datatype ! BMI_STRING char* S< ! BMI_INT int int16 ! BMI_LONG long int int32 ! BMI_FLOAT float float32 ! BMI_DOUBLE double float64 var_name = char_array_to_string(c_var_name, strlen(c_var_name)) include "bmi_get_var_type.inc" select case(var_name) case("netelemnode") type_name = "int" case("flowelemnode") type_name = "int" end select iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then type_name = "double" end if c_type = string_to_char_array(trim(type_name), len(trim(type_name))) end subroutine get_var_type subroutine get_var_role(c_var_name, role) bind(C, name="get_var_role") character(kind=c_char), intent(in) :: c_var_name(*) character(kind=c_char), intent(out) :: role(*) ! Roles: ! BMI_INPUT ! BMI_OUTPUT ! BMI_INPUTOUTPUT end subroutine get_var_role subroutine get_var_units( c_var_name, unit ) bind(C, name="get_var_units") character(kind=c_char), intent(in) :: c_var_name(*) character(kind=c_char), intent(out) :: unit(*) end subroutine get_var_units !> Returns the rank of a variable, i.e., its dimensionality. !! !! rank = count(shape) !! @see get_var_shape subroutine get_var_rank(c_var_name, rank) bind(C, name="get_var_rank") !DEC$ ATTRIBUTES DLLEXPORT :: get_var_rank use iso_c_binding, only: c_int, c_char character(kind=c_char), intent(in) :: c_var_name(*) integer(c_int), intent(out) :: rank ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) include "bmi_get_var_rank.inc" select case(var_name) case("netelemnode") rank = 2 case("flowelemnode") rank = 2 case("pumps", "weirs", "gates", "sourcesinks", "observations", "crosssections") ! Compound vars: shape = [numobj, numfields_per_obj] rank = 2 end select iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then rank = 1 return end if end subroutine get_var_rank !> Returns the shape of a variable, i.e., an array with length equal to this variables's rank. !! !! count(shape) = rank !! @see get_var_rank subroutine get_var_shape(c_var_name, shape) bind(C, name="get_var_shape") !DEC$ ATTRIBUTES DLLEXPORT :: get_var_shape use iso_c_binding, only: c_int, c_char, c_loc use m_flowgeom use network_data use m_observations, only: numobs, nummovobs, MAXNUMVALOBS2D, MAXNUMVALOBS3D, MAXNUMVALOBS3Dw use m_crosssections, only: ncrs, maxnval character(kind=c_char), intent(in) :: c_var_name(*) integer(c_int), intent(inout) :: shape(MAXDIMS) character(len=strlen(c_var_name)) :: var_name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) shape = (/0, 0, 0, 0, 0, 0/) select case(var_name) case("netelemnode") shape(1) = nump shape(2) = get_net_elem_max_nodes() return case("flowelemnode") shape(1) = ndx2d shape(2) = get_flow_elem_max_nodes() return ! Compounds: case("pumps") shape(1) = npumpsg shape(2) = 1 return case("weirs") shape(1) = nweirgen shape(2) = 2 return case("gates") shape(1) = ngategen shape(2) = 5 return case("sourcesinks") shape(1) = numsrc shape(2) = 3 return case("observations") shape(1) = numobs+nummovobs shape(2) = MAXNUMVALOBS2D+MAXNUMVALOBS3D+MAXNUMVALOBS3Dw return case("crosssections") shape(1) = ncrs shape(2) = maxnval return end select include "bmi_get_var_shape.inc" iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then shape(1) = ndkx return end if end subroutine get_var_shape subroutine get_2d_int(c_var_name, xptr) bind(C, name="get_2d_int") !DEC$ ATTRIBUTES DLLEXPORT :: get_2d_int use iso_c_binding use m_flow use m_flowgeom use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(out) :: xptr integer(c_int), target, allocatable, save :: x(:,:) integer :: i,j,k ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case("netelemnode") if (allocated(x)) deallocate(x) ! Deallocate if already allocated. allocate(x(nump, get_net_elem_max_nodes())) ! initialize to 0 x = -1 do k=1,nump do i=1,netcell(k)%n x(k,i) = netcell(k)%nod(i) end do end do xptr = c_loc(x) case("flowelemnode") if (allocated(x)) deallocate(x) ! Deallocate if already allocated. allocate(x(ndx, get_flow_elem_max_nodes())) ! initialize to 0 x = -1 do k=1,ndx do i=1,size(nd(k)%nod,1) x(k,i) = nd(k)%nod(i) end do end do xptr = c_loc(x) case("kn") xptr = c_loc(kn) end select end subroutine get_2d_int !> Gets a pointer to a model variable, given its name. !! subroutine get_var(c_var_name, x) bind(C, name="get_var") !DEC$ ATTRIBUTES DLLEXPORT :: get_var ! Return a pointer to the variable use iso_c_binding, only: c_double, c_char, c_loc use iso_c_utils use m_flow use m_flowgeom use m_sediment use m_flowparameters use network_data use m_sobekdfm use m_alloc use string_module character(kind=c_char), intent(in) :: c_var_name(*) !< Variable name. May be slash separated string "name/item/field": then get_compound_field is called. type(c_ptr), intent(inout) :: x integer(c_int), target, allocatable, save :: xi(:,:) integer :: i, j, k ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name character(len=strlen(c_var_name)) :: tmp_var_name character(len=strlen(c_var_name)) :: varset_name, item_name, field_name !< For parsing compound variable names. ! Store the name var_name = char_array_to_string(c_var_name) ! Please be conservative in adding variables here. Most variables ! can be computed outside. ! You can generate extra variables here. select case(var_name) case("uabs") if (.not.allocated(uabs)) allocate(uabs(size(ucx))) do i=1, size(ucx) uabs(i) = dble(sqrt( ucx(i)*ucx(i) + ucy(i)*ucy(i) )) end do x = c_loc(uabs) return case("sed1") ! sediment concentration, fraction1 if (.not.allocated(sed1)) allocate(sed1(size(ucx))) if(allocated(sed)) then do i=1, size(ucx) sed1(i) = sed(1, i) end do else sed1 = 0 end if x = c_loc(sed1) return case("netelemnode") if (allocated(xi)) deallocate(xi) ! Deallocate if already allocated. allocate(xi(nump, get_net_elem_max_nodes())) ! TODO transpose ! 1.... ndx2d, ndx2d+1...ndxi, ndxi+1...ndx ! 2d flow , 1d flow , boundary ! initialize to 0 xi = -1 do k=1,nump do i=1,netcell(k)%n xi(k,i) = netcell(k)%nod(i) end do end do x = c_loc(xi) return case("flowelemnode") if (allocated(xi)) deallocate(xi) ! Deallocate if already allocated. allocate(xi(get_flow_elem_max_nodes(), ndx2D)) ! initialize to 0 xi = -1 do k=1,ndx2D do i=1,size(nd(k)%nod,1) if (allocated(nd(k)%nod)) then xi(i, k) = nd(k)%nod(i) end if end do end do x = c_loc(xi) return case("kn") x = c_loc(kn) return end select ! Try to parse variable name as slash-separated id (e.g., 'weirs/Lith/crest_level') tmp_var_name = var_name call str_token(tmp_var_name, varset_name, DELIMS='/') ! Check for valid group/set name (e.g. 'observations') select case(varset_name) case ("pumps", "weirs", "gates", "sourcesinks", "observations", "crosssections") ! A valid group name, now parse the location id first... call str_token(tmp_var_name, item_name, DELIMS='/') if (len_trim(item_name) > 0) then ! A valid item name, now parse the field name... call str_token(tmp_var_name, field_name, DELIMS='/') if (len_trim(field_name) > 0) then ! Finally, a field_name was found, call the compound getter and return directly. call get_compound_field(string_to_char_array(varset_name), & string_to_char_array(item_name), & string_to_char_array(field_name), & x) return end if end if case ("controllabledam") call str_token(tmp_var_name, field_name, DELIMS='/') select case(field_name) case ("damlevel") x = c_loc(zcdam(1)) end select ! If we return here, the var_name was no valid triplet 'weirs/Lith/crest_level', ! so continue below with all default variables. continue end select ! ! check automatically generated interface include "bmi_get_var.inc" ! TODO: AvD: add returns to all auto generated cases to avoid unnecessary fall-through iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then call realloc(const_t, (/ ndkx, numconst /), keepExisting = .true.) do k=1,ndkx const_t(k, iconst) = constituents(iconst, k) end do x = c_loc(const_t(:,iconst)) return end if end subroutine get_var subroutine set_var(c_var_name, xptr) bind(C, name="set_var") !DEC$ ATTRIBUTES DLLEXPORT :: set_var ! Return a pointer to the variable use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), value, intent(in) :: xptr real(c_double), pointer :: x_0d_double_ptr real(c_double), pointer :: x_1d_double_ptr(:) real(c_double), pointer :: x_2d_double_ptr(:,:) real(c_double), pointer :: x_3d_double_ptr(:,:,:) integer(c_int), pointer :: x_0d_int_ptr integer(c_int), pointer :: x_1d_int_ptr(:) integer(c_int), pointer :: x_2d_int_ptr(:,:) integer(c_int), pointer :: x_3d_int_ptr(:,:,:) real(c_float), pointer :: x_0d_float_ptr real(c_float), pointer :: x_1d_float_ptr(:) real(c_float), pointer :: x_2d_float_ptr(:,:) real(c_float), pointer :: x_3d_float_ptr(:,:,:) ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name integer :: i ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) include "bmi_set_var.inc" iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then call c_f_pointer(xptr, x_1d_double_ptr, (/ ndkx /)) do i=1,ndkx constituents(iconst, i) = x_1d_double_ptr(i) end do return end if end subroutine set_var subroutine set_var_slice(c_var_name, c_start, c_count, xptr) bind(C, name="set_var_slice") !DEC$ ATTRIBUTES DLLEXPORT :: set_var_slice ! Return a pointer to the variable use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer integer(c_int), intent(in) :: c_start(*) integer(c_int), intent(in) :: c_count(*) character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), value, intent(in) :: xptr integer :: i real(c_double), pointer :: x_0d_double_ptr real(c_double), pointer :: x_1d_double_ptr(:) real(c_double), pointer :: x_2d_double_ptr(:,:) real(c_double), pointer :: x_3d_double_ptr(:,:,:) integer(c_int), pointer :: x_0d_int_ptr integer(c_int), pointer :: x_1d_int_ptr(:) integer(c_int), pointer :: x_2d_int_ptr(:,:) integer(c_int), pointer :: x_3d_int_ptr(:,:,:) real(c_float), pointer :: x_0d_float_ptr real(c_float), pointer :: x_1d_float_ptr(:) real(c_float), pointer :: x_2d_float_ptr(:,:) real(c_float), pointer :: x_3d_float_ptr(:,:,:) ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) call c_f_pointer(xptr, x_1d_double_ptr, (/ c_count(1) /)) ! custom overrides select case(var_name) case("ucx") !cell = netcell(index) ! !do edgeIndex = 1, cell%n ! ! calculate edge angle ! linkX1 = XK(KN(1,edgeIndex)) ! linkX2 = XK(KN(2,edgeIndex)) ! linkY1 = YK(KN(1,edgeIndex)) ! linkY2 = YK(KN(2,edgeIndex)) ! ! angle = atan((linkY2 - linkY1) / (linkX2 - linkX1)) ! ! u0(cell%lin(edgeIndex)) = value * cos(angle) ! u1(cell%lin(edgeIndex)) = value * cos(angle) !enddo ! convert it to the velocity increment on cell interfaces ! update u0 - velocity at cell edges ! a) calculate weights based on distances to the cell edges ! b) project current velocity component to the cell edge tangential velocity ! c) add resulting velocity increment to the every cell edge velocity ! for every link: ! 1. A - angle of the link, relative to the X axis (or calculate it from the link coordinates = atan((y2-y1)/(x2-x1)) ! 2. for Y component: V_link += Uy * sin(A) ! 3. for X component: V_link += Ux * cos(A) !ucx(index + 1) = value !case("ucy") ! convert it to the velocity increment on cell interfaces !ucy(index + 1) = value case("unorm") ! u1(index + 1) = value u1(c_start(1) + 1 : c_start(1) + c_count(1)) = x_1d_double_ptr(1:c_count(1)) return case("zk") do i = 1, c_count(1) call update_land(c_start(1), x_1d_double_ptr(i)) enddo !zkdropstep = value - zk(index + 1) !call dropland(xz(index + 1), yz(index + 1), 1) return end select include 'bmi_set_var_slice.inc' iconst = findname(numconst, const_names, var_name) if (iconst /= 0) then call c_f_pointer(xptr, x_1d_double_ptr, (/c_count(1)/)) do i=1,c_count(1) constituents(iconst, c_start(1)+i) = x_1d_double_ptr(i) end do return end if end subroutine set_var_slice !> Gets data for a specific field for a specific item from a set-variable of compound values. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! a single pump selected by its name, and the actual data by its field name. !! !! To convert the returned pointer to an actual data type in your calling application, !! use get_compound_field_type to determine the correct data type. subroutine get_compound_field(c_var_name, c_item_name, c_field_name, x) bind(C, name="get_compound_field") !DEC$ ATTRIBUTES DLLEXPORT :: get_compound_field ! Return a pointer to the compound's field variable use iso_c_binding, only: c_double, c_char, c_loc use iso_c_utils use m_flowexternalforcings use m_observations use m_crosssections character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' character(kind=c_char), intent(in) :: c_item_name(*) !< Name of a single item's index/location, e.g., 'Pump01' character(kind=c_char), intent(in) :: c_field_name(*) !< Name of the field to get, e.g., 'capacity' type(c_ptr), intent(inout) :: x !< Pointer (by reference) to requested value data, NULL if not available. integer :: item_index ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: item_name character(len=MAXSTRLEN) :: field_name ! Store the name var_name = char_array_to_string(c_var_name) item_name = char_array_to_string(c_item_name) field_name = char_array_to_string(c_field_name) select case(var_name) ! PUMPS case("pumps") call getStructureIndex('pumps', item_name, item_index) if (item_index == 0) then return endif select case(field_name) case("capacity") x = c_loc(qpump(item_index)) return end select ! WEIRS case("weirs") call getStructureIndex('weirs', item_name, item_index) if (item_index == 0) then return endif select case(field_name) case("crest_level") x = c_loc(zcgen((item_index-1)*3+1)) return case("lat_contr_coeff") ! TODO: RTC: AvD: get this from weir params return end select ! GATES case("gates") call getStructureIndex('gates', item_name, item_index) if (item_index == 0) then return endif select case(field_name) case("sill_level") x = c_loc(zcgen((item_index-1)*3+1)) return case("door_height") ! TODO: RTC: AvD: get this from gate/genstru params. return case("lower_edge_level") x = c_loc(zcgen((item_index-1)*3+2)) return case("opening_width") x = c_loc(zcgen((item_index-1)*3+3)) return case("horizontal_opening_direction") ! TODO: RTC: AvD: get this from gate/genstru params return end select ! SOURCE-SINKS case("sourcesinks") call getStructureIndex('sourcesinks', item_name, item_index) if (item_index == 0) then return endif select case(field_name) case("discharge") x = c_loc(qstss((item_index-1)*3+1)) return case("change_in_salinity") x = c_loc(qstss((item_index-1)*3+2)) return case("change_in_temperature") x = c_loc(qstss((item_index-1)*3+3)) return end select ! OBSERVATION STATIONS case("observations") call getObservationIndex(item_name, item_index) if (item_index == 0) then return end if select case(field_name) case("water_level") x = c_loc(valobs(IPNT_S1, item_index)) return case("water_depth") x = c_loc(valobs(IPNT_HS, item_index)) return case("salinity") x = c_loc(valobs(IPNT_SA1, item_index)) return case default ! TODO: AvD: error to warn for unimplemented feature? return end select ! MONITORING CROSSSECTIONS case("crosssections") call getCrosssectionIndex(item_name, item_index) if (item_index == 0) then return end if select case(field_name) case("discharge") x = c_loc(crs(item_index)%sumvalcur(IPNT_Q1C)) return case("velocity") x = c_loc(crs(item_index)%sumvalcur(IPNT_U1A)) return case("water_level") x = c_loc(crs(item_index)%sumvalcur(IPNT_S1A)) return case("water_depth") x = c_loc(crs(item_index)%sumvalcur(IPNT_HUA)) return case default ! TODO: AvD: error to warn for unimplemented feature? return end select end select ! var_name ! check automatically generated interface ! TODO: AvD: include "bmi_get_compound_field.inc" end subroutine get_compound_field !> Sets the value for a specific field for a specific item in a set-variable of compound values. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! a single pump selected by its name, and the actual data by its field name. !! !! The input value enters as a generic pointer, and will be converted to the required data type, e.g., double. !! If necessary, use get_compound_field_type and get_compound_field_shape to determine which input is expected. subroutine set_compound_field(c_var_name, c_item_name, c_field_name, xptr) bind(C, name="set_compound_field") !DEC$ ATTRIBUTES DLLEXPORT :: set_compound_field use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use iso_c_utils use unstruc_messages character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' character(kind=c_char), intent(in) :: c_item_name(*) !< Name of a single item's index/location, e.g., 'Pump01' character(kind=c_char), intent(in) :: c_field_name(*) !< Name of the field to get, e.g., 'capacity' type(c_ptr), value, intent(in) :: xptr !< Pointer (by value) to the C-compatible value data to be set. real(c_double), pointer :: x_0d_double_ptr integer :: item_index integer :: iostat ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: item_name character(len=MAXSTRLEN) :: field_name ! Store the name var_name = char_array_to_string(c_var_name) item_name = char_array_to_string(c_item_name) field_name = char_array_to_string(c_field_name) ! Debugging printing only: guess that it's a scalar double value, for now. call c_f_pointer(xptr, x_0d_double_ptr) write(msgbuf, '(6a,f20.6,a)', iostat=iostat) 'set_compound_field for ', trim(var_name), '(', trim(item_name), ')::', trim(field_name), ', will be set to value = ', x_0d_double_ptr, '.' call dbg_flush() ! TODO: AvD: include "bmi_set_compound_field.inc" select case(var_name) ! PUMPS case("pumps") call getStructureIndex('pumps', item_name, item_index) select case(field_name) case("capacity") call c_f_pointer(xptr, x_0d_double_ptr) qpump(item_index) = x_0d_double_ptr return end select ! WEIRS case("weirs") call getStructureIndex('weirs', item_name, item_index) select case(field_name) case("crest_level") call c_f_pointer(xptr, x_0d_double_ptr) zcgen((item_index-1)*3+1) = x_0d_double_ptr return case("lat_contr_coeff") ! TODO: RTC: AvD: set this in weir params return end select ! GATES case("gates") call getStructureIndex('gates', item_name, item_index) select case(field_name) case("sill_level") call c_f_pointer(xptr, x_0d_double_ptr) zcgen((item_index-1)*3+1) = x_0d_double_ptr return case("door_height") ! TODO: RTC: AvD: set this in gate/genstru params return case("lower_edge_level") call c_f_pointer(xptr, x_0d_double_ptr) zcgen((item_index-1)*3+2) = x_0d_double_ptr return case("opening_width") call c_f_pointer(xptr, x_0d_double_ptr) zcgen((item_index-1)*3+3) = x_0d_double_ptr return case("horizontal_opening_direction") ! TODO: RTC: AvD: set this once it's used return end select ! SOURCE-SINKS case("sourcesinks") call getStructureIndex('sourcesinks', item_name, item_index) if (item_index == 0) then return endif select case(field_name) case("discharge") call c_f_pointer(xptr, x_0d_double_ptr) qstss((item_index-1)*3+1) = x_0d_double_ptr return case("change_in_salinity") call c_f_pointer(xptr, x_0d_double_ptr) qstss((item_index-1)*3+2) = x_0d_double_ptr return case("change_in_temperature") call c_f_pointer(xptr, x_0d_double_ptr) qstss((item_index-1)*3+3) = x_0d_double_ptr return end select ! NOTE: observations and crosssections are read-only! end select end subroutine set_compound_field !> Gets the name for a specific field in a compound variable. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! each pump has the same list of fields, each with a name, for example 'capacity'. !! Use the returned field name to call other get_compound_field_* routines. subroutine get_compound_field_name(c_var_name, c_field_index, c_field_name) bind(C, name="get_compound_field_name") !DEC$ ATTRIBUTES DLLEXPORT :: get_compound_field_name use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use iso_c_utils character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' integer(kind=c_int), intent(in) :: c_field_index !< Index of the field in the compound variable's field list character(kind=c_char), intent(out) :: c_field_name(MAXSTRLEN) !< Returned name of the field to , e.g., 'capacity' integer :: field_index ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: field_name ! Store the name var_name = char_array_to_string(c_var_name) field_index = c_field_index select case(var_name) ! PUMPS case("pumps") select case(field_index) case(1) field_name = "capacity" end select ! WEIRS case("weirs") select case(field_index) case(1) field_name = "crest_level" case(2) field_name = "lat_contr_coeff" end select ! GATES case("gates") select case(field_index) case(1) field_name = "sill_level" case(2) field_name = "door_height" case(3) field_name = "lower_edge_level" case(4) field_name = "opening_width" case(5) field_name = "horizontal_opening_direction" end select ! SOURCE-SINKS case("sourcesinks") select case(field_index) case(1) field_name = "discharge" case(2) field_name = "change_in_salinity" case(3) field_name = "change_in_temperature" end select ! OBSERVATION STATIONS case("observations") select case(field_index) case(1) field_name = "water_level" case(2) field_name = "water_depth" case(3) field_name = "salinity" case default return end select ! MONITORING CROSSSECTIONS case("crosssections") select case(field_index) case(1) field_name = "discharge" case(2) field_name = "velocity" case(3) field_name = "water_level" case(4) field_name = "water_depth" case default ! TODO: AvD: error to warn for unimplemented feature? return end select end select ! var_name c_field_name = string_to_char_array(trim(field_name)) end subroutine get_compound_field_name !> Gets the type for a specific field in a compound variable. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! each pump has the same list of fields, each with a type, for example 'double'. !! Use the returned field type to prepare for calls to the set/get_compound_field routines. subroutine get_compound_field_type(c_var_name, c_field_name, c_type_name) bind(C, name="get_compound_field_type") !DEC$ ATTRIBUTES DLLEXPORT :: get_compound_field_type use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use iso_c_utils character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' character(kind=c_char), intent(in) :: c_field_name(*) !< Name of the field, e.g., 'capacity' character(kind=c_char), intent(out) :: c_type_name(MAXSTRLEN) !< Returned type of the field, e.g., 'double' ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: field_name character(len=MAXSTRLEN) :: type_name ! Store the name var_name = char_array_to_string(c_var_name) field_name = char_array_to_string(c_field_name) type_name = "double" c_type_name = string_to_char_array(trim(type_name)) end subroutine get_compound_field_type !> Gets the rank for a specific field in a compound variable. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! each pump has the same list of fields, each with a rank, typically 0 (scalar). !! Use the returned field rank to prepare for calls to the set/get_compound_field routines. subroutine get_compound_field_rank(c_var_name, c_field_name, rank) bind(C, name="get_compound_field_rank") !DEC$ ATTRIBUTES DLLEXPORT :: get_compound_field_rank use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use iso_c_utils character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' character(kind=c_char), intent(in) :: c_field_name(*) !< Name of the field, e.g., 'capacity' integer(kind=c_int), intent(out) :: rank !< Returned rank of the field, e.g., 1 ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: field_name ! Store the name var_name = char_array_to_string(c_var_name) field_name = char_array_to_string(c_field_name) rank = 0 ! for all scalar vars/fields now (unless pumps start to have multiple stages) end subroutine get_compound_field_rank !> Gets the shape for a specific field in a compound variable. !! !! For example, all pumping stations in a model may be collected in a single compound variable set, named 'pumps', !! each pump has the same list of fields, each with a shape, which is now only for completeness: shape(1:rank) = 0. !! Use the returned field shape to prepare for calls to the set/get_compound_field routines. subroutine get_compound_field_shape(c_var_name, c_field_name, shape) bind(C, name="get_compound_field_shape") !DEC$ ATTRIBUTES DLLEXPORT :: get_compound_field_shape use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use iso_c_utils character(kind=c_char), intent(in) :: c_var_name(*) !< Name of the set variable, e.g., 'pumps' character(kind=c_char), intent(in) :: c_field_name(*) !< Name of the field, e.g., 'capacity' integer(kind=c_int), intent(inout) :: shape(MAXDIMS) !< Returned shape of the field, e.g., [1] ! The fortran name of the attribute name character(len=MAXSTRLEN) :: var_name character(len=MAXSTRLEN) :: field_name ! Store the name var_name = char_array_to_string(c_var_name) field_name = char_array_to_string(c_field_name) shape = 0 ! All fields now still scalar: rank=0, shape irrelevant. end subroutine get_compound_field_shape subroutine get_1d_double(c_var_name, x) bind(C, name="get_1d_double") !DEC$ ATTRIBUTES DLLEXPORT :: get_1d_double use iso_c_binding, only: c_double, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(inout) :: x ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call get_var(c_var_name, x) end select end subroutine get_1d_double subroutine get_1d_int(c_var_name, x) bind(C, name="get_1d_int") !DEC$ ATTRIBUTES DLLEXPORT :: get_1d_int use iso_c_binding, only: c_int, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(inout) :: x character(len=strlen(c_var_name)) :: var_name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call get_var(c_var_name, x) end select end subroutine get_1d_int subroutine get_2d_double(c_var_name, xptr) bind(C, name="get_2d_double") !DEC$ ATTRIBUTES DLLEXPORT :: get_2d_double use iso_c_binding, only: c_double, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(inout) :: xptr real(c_double), target, allocatable, save :: x(:,:) integer :: i ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) if (allocated(x)) deallocate(x) select case(var_name) case("flowelemcontour_x") allocate(x(ndx, get_flow_elem_max_contour())) x=0 ! Set this to nan? do i=1,ndx x(i,1:size(nd(i)%x)) = nd(i)%x enddo xptr = c_loc(x) case("flowelemcontour_y") allocate(x(ndx, get_flow_elem_max_contour())) x=0 ! I would like to set nans here. but how? do i=1,ndx x(i,1:size(nd(i)%x)) = nd(i)%y enddo xptr = c_loc(x) end select end subroutine get_2d_double subroutine set_1d_double_at_index(c_var_name, index, value) bind(C, name="set_1d_double_at_index") !DEC$ ATTRIBUTES DLLEXPORT :: set_1d_double_at_index use m_flow use network_data use m_flowparameters use iso_c_binding, only: c_double, c_char, c_int character(kind=c_char), intent(in) :: c_var_name(*) !< Name of a BMI-available variable (should be associated with a 1d double array) integer(c_int), value, intent(in) :: index !< Position in array where to set the value. 0-based, pass by value. real(c_double), value, intent(in) :: value !< The new value to set in the array. character(len=strlen(c_var_name)) :: var_name type (tface) :: cell integer :: edgeIndex double precision :: linkX1, linkX2, linkY1, linkY2 double precision :: angle real(c_double), target :: valuet type(c_ptr) :: xptr integer(c_int) :: size1(1) var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case("s1") s1(index + 1) = value case("ucx") cell = netcell(index) do edgeIndex = 1, cell%n ! calculate edge angle linkX1 = XK(KN(1,edgeIndex)) linkX2 = XK(KN(2,edgeIndex)) linkY1 = YK(KN(1,edgeIndex)) linkY2 = YK(KN(2,edgeIndex)) angle = atan((linkY2 - linkY1) / (linkX2 - linkX1)) u0(cell%lin(edgeIndex)) = value * cos(angle) u1(cell%lin(edgeIndex)) = value * cos(angle) enddo ! convert it to the velocity increment on cell interfaces ! update u0 - velocity at cell edges ! a) calculate weights based on distances to the cell edges ! b) project current velocity component to the cell edge tangential velocity ! c) add resulting velocity increment to the every cell edge velocity ! for every link: ! 1. A - angle of the link, relative to the X axis (or calculate it from the link coordinates = atan((y2-y1)/(x2-x1)) ! 2. for Y component: V_link += Uy * sin(A) ! 3. for X component: V_link += Ux * cos(A) !ucx(index + 1) = value case("ucy") ! convert it to the velocity increment on cell interfaces !ucy(index + 1) = value case("unorm") u1(index + 1) = value case("zk") call update_land(index + 1, value) !zkdropstep = value - zk(index + 1) !call dropland(xz(index + 1), yz(index + 1), 1) case("fstuw") fstuw(index + 1) = value case("froer") froer(index + 1) = value case default ! Fall back to any of the auto-generated BMI vars in set_var_slice (with slice count just 1) size1(1) = 1 valuet = value xptr = c_loc(valuet) ! To pass on this subroutine argument to set_var_slice as a C pointer call set_var_slice(c_var_name, (/ index /), size1, xptr) end select end subroutine set_1d_double_at_index subroutine set_1d_double(c_var_name, x) bind(C, name="set_1d_double") !DEC$ ATTRIBUTES DLLEXPORT :: set_1d_double use iso_c_binding, only: c_double, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(in) :: x ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call set_var(c_var_name, x) end select end subroutine set_1d_double subroutine set_2d_double(c_var_name, x) bind(C, name="set_2d_double") !DEC$ ATTRIBUTES DLLEXPORT :: set_2d_double use iso_c_binding, only: c_double, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(in) :: x ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call set_var(c_var_name, x) end select end subroutine set_2d_double subroutine set_1d_int(c_var_name, x) bind(C, name="set_1d_int") !DEC$ ATTRIBUTES DLLEXPORT :: set_1d_int use iso_c_binding, only: c_int, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(in) :: x ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call set_var(c_var_name, x) end select end subroutine set_1d_int subroutine set_2d_int(c_var_name, x) bind(C, name="set_2d_int") !DEC$ ATTRIBUTES DLLEXPORT :: set_2d_int use iso_c_binding, only: c_int, c_char, c_loc use m_flow use network_data character(kind=c_char), intent(in) :: c_var_name(*) type(c_ptr), intent(in) :: x ! The fortran name of the attribute name character(len=strlen(c_var_name)) :: var_name ! Store the name var_name = char_array_to_string(c_var_name, strlen(c_var_name)) select case(var_name) case default call set_var(c_var_name, x) end select end subroutine set_2d_int ! TODO: ! grid routines ! Utility functions, move these to interop module ! Make functions pure so they can be used as input arguments. integer(c_int) pure function strlen(char_array) character(c_char), intent(in) :: char_array(MAXSTRLEN) integer :: inull, i strlen = 0 do i = 1, size(char_array) if (char_array(i) .eq. C_NULL_CHAR) then strlen = i-1 exit end if end do end function strlen pure function char_array_to_string(char_array, length) integer(c_int), intent(in) :: length character(c_char),intent(in) :: char_array(length) character(len=length) :: char_array_to_string integer :: i do i = 1, length char_array_to_string(i:i) = char_array(i) enddo end function char_array_to_string pure function string_to_char_array(string, length) character(len=length), intent(in) :: string integer(c_int),intent(in) :: length character(kind=c_char,len=1) :: string_to_char_array(length+1) integer :: i do i = 1, length string_to_char_array(i) = string(i:i) enddo string_to_char_array(length+1) = C_NULL_CHAR end function string_to_char_array ! TODO move to network_data... integer function get_net_elem_max_nodes() use network_data integer :: k ! Determine max nr. of vertices in NetElems (netcells) get_net_elem_max_nodes = 0 do k=1,nump get_net_elem_max_nodes = max(get_net_elem_max_nodes, netcell(k)%n) end do end function get_net_elem_max_nodes integer function get_flow_elem_max_nodes() use network_data integer :: k ! Determine max nr. of flow nodes get_flow_elem_max_nodes = 0 do k=1,ndx get_flow_elem_max_nodes = max(get_flow_elem_max_nodes, size(nd(k)%nod)) end do end function get_flow_elem_max_nodes integer function get_flow_elem_max_contour() use network_data integer :: k ! Determine max nr. of flow contour points (TODO: nodes+1?) get_flow_elem_max_contour = 0 do k=1,ndx get_flow_elem_max_contour = max(get_flow_elem_max_contour, size(nd(k)%x)) end do end function get_flow_elem_max_contour ! Geometry dll functions subroutine triang(cptr_sx, cptr_sy, cptr_sv, c_numS, cptr_dx, cptr_dy, c_numD, cptr_res) bind(C, name="triang") !DEC$ ATTRIBUTES DLLEXPORT :: triang use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use unstruc_model use m_samples use m_missing implicit none ! parameters type(c_ptr), intent(in) :: cptr_sx ! samples x, y, values type(c_ptr), intent(in) :: cptr_sy type(c_ptr), intent(in) :: cptr_sv integer(c_int), intent(in) :: c_numS ! num samples type(c_ptr), intent(in) :: cptr_dx ! destinations x, y type(c_ptr), intent(in) :: cptr_dy integer(c_int), intent(in) :: c_numD ! num destination points type(c_ptr), intent(inout) :: cptr_res ! return values (ptr to double array) ! local variables integer :: numS integer :: numD integer :: jdla = 1 real(c_double), pointer :: ptr(:) real(c_double), pointer :: dx(:) real(c_double), pointer :: dy(:) real(c_double), pointer :: dRes(:) numS = c_numS numD = c_numD ! (re)allocate sample arrays if (allocated(XS)) then deallocate(XS,YS,ZS) end if allocate(XS(numS), YS(numS), ZS(numS)) ! copy ptr's to fortran arrays call c_f_pointer(cptr_sx, ptr, (/numS/)) XS(:) = ptr call c_f_pointer(cptr_sy, ptr, (/numS/)) YS(:) = ptr call c_f_pointer(cptr_sv, ptr, (/numS/)) ZS(:) = ptr call c_f_pointer(cptr_dx, dx, (/numD/)) call c_f_pointer(cptr_dy, dy, (/numD/)) call c_f_pointer(cptr_res, dRes, (/numD/)) ! set max stuff NS = numS NSMAX = numS ! assign 'missing value' to all elements of dRes dRes = DMISS ! call triangulate call triinterp2(dx, dy, dRes, numD, jdla) end subroutine triang subroutine averaging(cptr_sx, cptr_sy, cptr_sv, c_nums, cptr_cx, cptr_cy, cptr_cxx, cptr_cyy, cptr_cnp, c_numc, c_n6, cptr_res, cptr_meth, cptr_nmin, cptr_csize) bind(C, name="averaging") !DEC$ ATTRIBUTES DLLEXPORT :: averaging use iso_c_binding, only: c_double, c_char, c_loc, c_f_pointer use unstruc_model use m_missing use M_INTERPOLATIONSETTINGS use m_kdtree2 implicit none ! parameters type(c_ptr), intent(in) :: cptr_sx ! samples x, y, values type(c_ptr), intent(in) :: cptr_sy type(c_ptr), intent(in) :: cptr_sv integer(c_int), intent(in) :: c_nums ! number of samples type(c_ptr), intent(in) :: cptr_cx ! destination cell center x, y type(c_ptr), intent(in) :: cptr_cy type(c_ptr), intent(in) :: cptr_cxx ! destination cell corner x, y type(c_ptr), intent(in) :: cptr_cyy type(c_ptr), intent(in) :: cptr_cnp ! destination cell corner array lengths integer(c_int), intent(in) :: c_numc ! number of destination cells integer(c_int), intent(in) :: c_n6 ! max. cell corner array length type(c_ptr), intent(inout) :: cptr_res ! return values (ptr to double array) integer(c_int), intent(in) :: cptr_meth ! averaging method integer(c_int), intent(in) :: cptr_nmin ! minimum nr of samples for avaraging real(c_double), intent(in) :: cptr_csize ! relative search cell size ! local variables real(c_double), pointer :: sx(:) real(c_double), pointer :: sy(:) real(c_double), pointer :: svtmp(:) integer :: nums real(c_double), pointer :: cx(:) real(c_double), pointer :: cy(:) real(c_double), pointer :: cxtmp(:) real(c_double), pointer :: cytmp(:) integer, pointer :: cnp(:) integer :: numc integer :: n6 real(c_double), pointer :: res(:) double precision, allocatable :: sv(:,:) integer, allocatable :: ipsam(:) double precision, allocatable :: cz(:,:) double precision, allocatable :: cxx(:,:) double precision, allocatable :: cyy(:,:) integer :: meth integer :: nmin double precision :: csize integer :: i, j, k, IAVtmp, NUMMINtmp, INTTYPEtmp, ierr double precision :: RCELtmp ! cache interpolation settings IAVtmp = IAV NUMMINtmp = NUMMIN INTTYPEtmp = INTERPOLATIONTYPE RCELtmp = RCEL ! assign ranges and settings nums = c_nums numc = c_numc n6 = c_n6 meth = cptr_meth nmin = cptr_nmin csize = cptr_csize !assign pointers call c_f_pointer(cptr_sx, sx, (/nums/)) call c_f_pointer(cptr_sy, sy, (/nums/)) call c_f_pointer(cptr_sv, svtmp, (/nums/)) call c_f_pointer(cptr_cx, cx, (/numc/)) call c_f_pointer(cptr_cy, cy, (/numc/)) call c_f_pointer(cptr_cxx, cxtmp, (/n6*numc/)) call c_f_pointer(cptr_cyy, cytmp, (/n6*numc/)) call c_f_pointer(cptr_cnp, cnp, (/numc/)) call c_f_pointer(cptr_res, res, (/numc/)) !allocate & copy to 2d arrays allocate(sv(1, nums), ipsam(nums), cz(1,numc), cxx(n6, numc), cyy(n6, numc)) sv(1,:) = svtmp(:) ipsam(:) = 1 k = 1 do i = 1, numc cz(1,i) = DMISS do j=1,n6 cxx(j, i) = cxtmp(k) cyy(j, i) = cytmp(k) k = k + 1 enddo enddo if(meth > 0 .and. meth < 8) then IAV = meth else goto 1234 endif if(nmin > 0) then NUMMIN = nmin else goto 1234 endif if(csize > 0 .and. csize < 10) then RCEL = csize else goto 1234 endif INTERPOLATIONTYPE = 2 call build_kdtree(treeglob, nums, sx, sy, ierr) call averaging2(1, nums, sx, sy, sv, ipsam, cx, cy, cz, numc, cxx, cyy, n6, cnp, 1) call delete_kdtree2(treeglob) !copy values back res(:) = cz(1,:) 1234 continue !unroll & cleanup IAV = IAVtmp NUMMIN = NUMMINtmp INTERPOLATIONTYPE = INTTYPEtmp RCEL = RCELtmp deallocate(sv, ipsam, cz, cxx, cyy) end subroutine averaging ! Further custom api functions subroutine find_cells(c_net_file, c_numCells, c_maxPerCell, cptr_netElemNode) bind(C, name="find_cells") !DEC$ ATTRIBUTES DLLEXPORT :: find_cells use iso_c_binding, only: c_int, c_char, c_ptr use network_data use unstruc_files use m_netw character(kind=c_char),intent(in) :: c_net_file(MAXSTRLEN) integer(c_int), intent(out) :: c_numCells integer(c_int), intent(out) :: c_maxPerCell type(c_ptr), intent(inout) :: cptr_netElemNode ! return values (ptr to 2d int array) integer, pointer :: netElemNode(:) character(len=strlen(c_net_file)) :: net_file integer :: numk_read, numl_read, istat integer :: maxNodes, ci, ni, i call resetFullFlowModel() call INIDAT() ! read grid filename from c-ptr net_file = char_array_to_string(c_net_file, strlen(c_net_file)) ! read grid from netcdf call loadNetwork(net_file, istat, 0) nump = 0 ! find cells call FINDCELLS(0) maxNodes = 0 ! find max num of vertices per cell do ci=1, nump ! nump = num cells maxNodes = max(maxNodes, netcell(ci)%n) end do ! allocate 1d array (for 2d data) allocate(netElemNode(nump*maxNodes)) ! fill NetElemNode array i=1 do ci=1, nump do ni=1, maxNodes if (ni <= netcell(ci)%n) then netElemNode(i) = netcell(ci)%nod(ni) else netElemNode(i) = -1 endif i=i+1 end do end do ! copy to c variables cptr_netElemNode = c_loc(netElemNode(1)) ! TODO: Fedor/Tiemen: please review whether this still works (was fix for GNU Fortran "Error: Argument 'netelemnode' to 'c_loc' at (1) must be an associated scalar POINTER") c_numCells = nump c_maxPerCell = maxNodes end subroutine find_cells !> snap polylines to mesh subroutine get_snapped_feature(c_feature_type, c_Nin, cptr_xin, cptr_yin, c_Nout, cptr_xout, cptr_yout, cptr_feature_ids, c_ierror) bind(C, name="get_snapped_feature") !DEC$ ATTRIBUTES DLLEXPORT :: get_snapped_feature use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_f_pointer use m_snappol use m_missing implicit none character(kind=c_char), intent(in) :: c_feature_type(MAXSTRLEN) !< feature type ('thindam') integer(c_int), intent(in) :: c_Nin !< input feature array length type(c_ptr), intent(in) :: cptr_xin, cptr_yin !< input feature coordinates integer(c_int), intent(out) :: c_Nout !< output array length type(c_ptr), intent(inout) :: cptr_xout, cptr_yout !< output feature coordinates type(c_ptr), intent(inout) :: cptr_feature_ids !< output feature ids integer(c_int), intent(out) :: c_ierror !< ierror (1) or not (0) character(len=strlen(c_feature_type)) :: feature_type real(c_double), pointer :: ptr(:) ! temporary pointer double precision, dimension(:), target, allocatable, save :: xout, yout !< memory leak integer, dimension(:), target, allocatable, save :: feature_ids !< memory leak double precision, dimension(:), allocatable :: xin, yin c_ierror = 1 ! read feature type from c-ptr feature_type = char_array_to_string(c_feature_type, strlen(c_feature_type)) ! allocate allocate(xin(c_Nin), yin(c_Nin)) ! copy pointers to fortran array call c_f_pointer(cptr_xin, ptr, (/c_Nin/)) xin(:) = ptr call c_f_pointer(cptr_yin, ptr, (/c_Nin/)) yin(:) = ptr select case( feature_type ) case("thindam", "fixedweir", "crosssection", "gate", "weir", "pump") call snappol(c_Nin, xin, yin, DMISS, c_Nout, xout, yout, feature_ids, c_ierror) case("obspoint") call snappnt(c_Nin, xin, yin, DMISS, c_Nout, xout, yout, feature_ids, c_ierror) case DEFAULT call snapbnd(feature_type, c_Nin, xin, yin, DMISS, c_Nout, xout, yout, feature_ids, c_ierror) end select if ( c_ierror.ne.0 ) goto 1234 cptr_xout = c_loc(xout) cptr_yout = c_loc(yout) cptr_feature_ids = c_loc(feature_ids) c_ierror = 0 1234 continue if ( allocated(xin) ) deallocate(xin) if ( allocated(yin) ) deallocate(yin) return end subroutine get_snapped_feature subroutine write_netgeom(c_net_file) bind(C, name="write_netgeom") !DEC$ ATTRIBUTES DLLEXPORT :: write_netgeom use iso_c_binding, only: c_int, c_char, c_ptr use network_data use unstruc_netcdf character(kind=c_char),intent(in) :: c_net_file(MAXSTRLEN) character(len=strlen(c_net_file)) :: netgeom_file ! read grid filename from c-ptr netgeom_file = char_array_to_string(c_net_file, strlen(c_net_file)) if ( netstat.ne.NETSTAT_OK ) then call findcells(0) end if call unc_write_net(netgeom_file, 1, 0) end subroutine write_netgeom subroutine write_partition_metis(c_net_file, c_npart, c_jacontiguous) bind(C, name="write_partition_metis") !DEC$ ATTRIBUTES DLLEXPORT :: write_partition_metis use iso_c_binding, only: c_int, c_char, c_ptr use network_data use unstruc_netcdf use m_partitioninfo character(kind=c_char), intent(in) :: c_net_file(MAXSTRLEN) integer(c_int), intent(in) :: c_npart integer(c_int), intent(in) :: c_jacontiguous character(len=strlen(c_net_file)) :: netfilename integer :: npart integer :: jacontiguous if ( netstat.ne.NETSTAT_OK ) then call findcells(0) end if if(c_npart < 2) then return else npart=c_npart endif if(c_jacontiguous < 1) then jacontiguous = 0 else jacontiguous = 1 endif call partition_METIS_to_idomain(npart, jacontiguous) ndomains = npart call generate_partition_pol_from_idomain() netfilename = char_array_to_string(c_net_file, strlen(c_net_file)) if(ndomains > 1) then call partition_write_domains(netfilename,0) endif end subroutine write_partition_metis subroutine write_partition_pol(c_net_file, c_pol_file) bind(C, name="write_partition_pol") !DEC$ ATTRIBUTES DLLEXPORT :: write_partition_pol use iso_c_binding, only: c_int, c_char, c_ptr use network_data use unstruc_netcdf use m_partitioninfo use m_polygon use unstruc_files character(kind=c_char), intent(in) :: c_net_file(MAXSTRLEN) character(kind=c_char), intent(in) :: c_pol_file(MAXSTRLEN) character(len=strlen(c_net_file)) :: netfilename character(len=strlen(c_pol_file)) :: polfilename integer :: minp if ( netstat.ne.NETSTAT_OK ) then call findcells(0) end if polfilename = char_array_to_string(c_pol_file, strlen(c_pol_file)) call newfil(minp, polfilename) call reapol(minp, 0) if (npl > 1) then ! use polygons call generate_partitioning_from_pol() end if netfilename = char_array_to_string(c_net_file, strlen(c_net_file)) if(ndomains > 1) then call partition_write_domains(netfilename,0) endif end subroutine write_partition_pol end module bmi