You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Starting with exchanges with Jason D. et al. (OWP)
Hi Joseph,
Adding the CPP flags to the SCHISM compile time sounds like an excellent idea to me. This will greatly assist us in streamlining the latest SCHISM release to our NextGen BMI development. Thank you very much for offering to take on this task.
I just quickly want to follow up on the hotstart question. Essentially, we will have no knowledge of what the source term discharges and the open water boundary levels are upon the NextGen model start time ("t0") from existing model formulations that are coupled with SCHISM. SCHISM will be able to start ingesting source and boundary data in the following hour ("t1") from other model formulations. To my understanding, the "hot_start.nc" file should provide data fields to force the model at the start time ("t0"). You've mentioned that there is another "other_hot_init.nc" that may initialize source and boundary terms at the model start time, or did I misunderstand your comment from the previous email? If there exists this "other_hot_init.nc" for our purposes, is there a certain compiling flag or param.nml variable to turn on and provide this hot_start.nc output?
As always, thank you Joseph for your continued support for the NextGen SCHISM BMI development on behalf of OWP. You and Carsten's insight for SCHISM is critical during our development phases. Have a great rest of your week.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center
On Thu, Apr 20, 2023 at 2:40 PM Y. Joseph Zhang <yjzhang@vims.edu> wrote:
Hi Jason:
I thought about these more during my CA trip. I think we can streamline the options by simply adding CPP flags (USE_NWM_BMI e.g.) without introducing additional parameters. I looked back at the nws=3 option (external atmos model), and think we can actually also replace it with USE_EXT_ATMOS.
To keep the logic straight, I think we should reuse some existing options with those CPP flags. So the users will still need to prepare source.nc if if_source=-1, with the understanding that some info won’t be required and will be overwritten by BMI. If USE_NWM_BMI is on, some reading will be bypassed, while static/geom arrays are still read in. This also helps with the hotstart issue you mentioned below; at t0 the code will read in the first record from source.nc (for elev, this is taken care of by hotstart.nc). I think assuming source/sink=0 at the restart of the run will be OK also as the init condition will be forgotten eventually.
In general, upon hotstart, SCHISM will read in previous ‘state variables’ but not the forcing time series. The latter are done in ‘other_hot_init’.
Let me know if this sounds good to me and I’ll add those CPPs if you agree. Thx.
Thank you very much for all of your comments and suggestions. This has already helped me push forward some SCHISM BMI functionality this week. I just have a few quick follow up questions regarding your answers to the previous four questions I postulated in my initial email. My questions to your responses are highlighted in yellow below:
In the case of a serial mode SCHISM simulation, why would the node variable counts "np_global" and "np" differ in this case? My understanding is that "np" is the number of resident elements within a partitioned portion of the mesh. However, when we have a single thread for your simulation, shouldn't np = np_global? Or does np_global account for resident and ghost nodes, in which ghost nodes also relate to boundary nodes?
Joseph Zhang response: I cannot reproduce your issue. With 1 MPI process for compute (scribe I/O processes are in a separate comm world), np_global=npa=np. I tested it with both OLDIO on/off and verified this in mirror.out (the beginning part of it shows the partitioning info). Anyway, np is the # of resident NODES (ne is elem) without ghost. npg is # of ghost nodes, npa=np+npg. With 1 MPI process, npg=0, so np=npa=np_global; similarly for elem/sides.
Jason Ducker response: My apologies for mixing this question up with SCHISM here Joseph. I believe this issue was actually occurring in another model (DFlowFM) and not SCHISM. I was able to confirm as well on my end that np_global=npa=np with SCHISM running on a single thread. There are no ghost nodes in this case and this simplifies the initial BMI implementation with SCHISM running in a serial mode (all variables are on the same global nodes within the domain, no local nodes based on partitioned domain assigned to varying threads).
For open water level boundaries, what internal SCHISM variable deals with the actual forcing data? My initial investigation seemed to point to the variable "ath2", which has a number of dimension dependencies (ntracers, nvrt, neta_global, 2, nthfiles2), which include vertical levels and tracers (salinity and temperature). Is this indeed the case? If so, then what dimensions physically point to water levels (given in elev2Dth.nc) for a given open boundary segment? The idea here is that we would like these boundary arrays to be flattened into a one dimensional array in the simplest manner possible without dealing with a multitude of array dimensions that are not dependent on the unstructured mesh connectivity.
Joseph Zhang response: Yes ath2() is the receiving array. In addition, the code also needs the time step used in elev2D.th.nc (‘time_step’). The global data are usually defined in Core/schism_glbl.F90. ath2() is single precision.
To reduce memory consumption, I used ath2 to receive b.c. inputs for many variables. In the case of elev, the final array is eth() found around line 1397 in schism_step.F90, where you can see the time interpolation is done from ath2 into eth.
The dimensions in ath2() are: (1) ntracers for tracers so for elev, this is 1; (2) nvrt is vertical indices, so for elev this is also 1 (no vertical structure); (3) neta_global is the lumped total number of nodes that require this input. The order is determined from info in hgrid.gr3 (b.c. part) and bctides.in (elev b.c. flag of 4 or 5); (4) ‘2’ is to account for 2 time records for time interp (i.e. store 2 records from .nc at a time); (5) nthfiles2 points to the variable type, ‘1’ is elev, ‘2’ is velocity etc.
eth(mnond_global,nope_global): 2nd index is the open bnd segment ID and 1st index is the node ID on that particular segment. This should help you flatten the array. You can get the global node ID by using iond_global(), and then use ipgl() to find local node ID (if this rank owns it); see the code near line 4934 in _step.
Mapping from ‘neta_global’ to local node IDs is more involved. You basically loop over all open bnd segments and test to see if the elev flag is 4 or 5. Maybe you can create some mapping arrays on BMI side and I can provide you the code.
Jason Ducker response: This explanation makes things much clearer with the way we would slice the "ath2" array and update the indices or entire boundary condition variable over each BMI model time step. Thank you for the source code line to see an example of how to obtain the local node ID for an index within the "ath2" array. You also raise a very interesting point regarding how water level boundaries are updated over each model iteration. If we were to bypass the "elev2D.th.nc" forcing file requirement, then the BMI must be able to update water level boundaries at the "start" and "end" time for each model timestep. Essentially, it sounds like the BMI will also need to provide "start time" and "end time" to the variable "th_time2" in order to properly update the ratio variable in line 1382 within "schism_step.F90". My assumption is that we would need to set a SCHISM compiler flag to essentially turn off the updated netcdf read statements for the variables "time_series" and "time" within the "elev2D.th.nc" file right? Since the calculation of the "eth" variable is also dependent on the time component, then it may be more ideal to just update the "ath2" and "th_time2" variables within the BMI "set_var_indices" functions. Let me know though if you think there is a simpler path here than what I proposed.
For discharge/precipitation source terms, what internal SCHISM variables deal with the actual forcing data for sources? My initial investigation pointed to the variable "ath3", which also has a number of dimension dependencies (nsources, ntracers, 2, nthfiles3). If so, then what dimensions physically point to actual discharge and precipitation source terms (in "source.nc") for a given land boundary segment?
Joseph Zhang response: Correct again but ath3(max(1,nsources,nsinks), ntracers,2,nthfiles3). Source/sink can be injected at any element, not just boundary. The IDs of injection elements are found in source.nc. The last dim ‘nthfiles3’ points to 3 types of inputs: volume source, volume sink, and source concentration. (1) nsources/nsinks: total # of source or sink elements; (2) ntracers: in 2D model you’ll only have 2 tracers (T,S) that are not really used. If nthfiles3 points to 1 or 2, this dim is 1 b/c volumes are same for all tracers; if it points to ‘3’, then this dim is 1:2 (T,S); (3) ‘2’ refers to 2 time records; (4) last dim is type explained before. To map the first index to global elem ID, use ieg_source(max(1,nsources)),ieg_sink(max(1,nsinks)). Then use iegl() to get local elem ID. The final receiving arrays are vsource(), msource().
Jason Ducker reponse: Thank you for describing the indexing for the source term. I was able to find the source code accordingly based on your description for finding the local element ID for each source term. As mentioned in my previous response to your question fixed submodule pointer to actually point to SCHISM instead of to schism-esmf #2 statement, we also appear to have temporal dependencies with updating the "ath3" source term array (at "t0" and "t1" time) and the time array "th_time3", which originate from the contents of the "source.nc" file. If we were to initiate SCHISM with a hotstart file at the beginning of the simulation, will there already be water levels and source terms available at the open and land boundary nodes respectively? The reason why I ask this is because SCHISM will have no information from the BMI about water level and source terms at "t0" (start time of simulation). Only once models begin executing to the next model time step ("t1") that SCHISM can retrieve boundary conditions fields. Therefore, water level and source term boundary arrays will have to be assumed as values of zero at the start time of the simulation unless data fields are already available at those nodes (within a hotstart file I assume). I'm curious what your take on this issue would be here for SCHISM and whether or not a hotstart file might resolve this issue.
Within the "schism_init.f90", the model attempts to initialize the bounds for a series of arrays for land/open boundary variables for "source.nc" and "elev2Dth.nc" files. In your expert opinions, could we simply just turn off the netcdf read statements for the boundary forcing data and just fill the arrays before a given model iteration? Or are there more dependencies regarding the external forcing data being initialized than my current understanding here?
Joseph Zhang response: To bypass SCHISM’s existing options, I think the best way is for me to add a new option that will essentially do nothing, and wait for external program to fill in those needed arrays. This is the approach we used for coupling with atmos model (nws=3). I don’t think there are ways to avoid any code change on SCHISM side, are there? For those 2 .nc inputs, the geometric arrays are all ‘static’ (not changing over time) and ath*() are the only dynamic part so I think you can easily fill them in with pointers.
Jason Ducker response: I completely agree with the approach for adding a SCHISM option with the param.nml file (external_bnd_forcing=1?) to tell the model to bypass external forcing files for boundary conditions and allow the user to depend on coupling software to fill the boundary array values. What's your initial guess on the level of effort this could take for implementing such a SCHISM configuration? To my knowledge, it only appears that we need to include an "#ifelse" SCHISM flag to either include or exclude netcdf inquiry commands to update arrays within "mis_subs.F90" and "schism_step.F90", but I may be missing a few pieces to the puzzle here possibly. I could try to make an attempt on source code modification possibly if you agree with my initial solution to this code approach? To my understanding, we are essentially asking the user only to specify boundary condition locations within the "elev2D.th.nc" and "source.nc" files and bypass any information regarding "time_step" "time" and "time_series"/forcing data.
We greatly appreciate the collaborations of Joseph and Carsten with the OWP development of the SCHISM BMI for the NextGen project. I think we've been making great progress so far on SCHISM and we're getting into some interesting source code topics that may be required to properly enable SCHISM BMI implementation for NextGen. If my questions require a future meeting discussion, then please let me know and I'll be happy to coordinate these efforts moving forward. Have a great weekend.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center
On Thu, Apr 13, 2023 at 12:42 PM Y. Joseph Zhang <yjzhang@vims.edu> wrote:
See below.
-Joseph
Y. Joseph Zhang
Web: schism.wiki
Office: 804 684 7466
I want to first thank you for taking a look at our SCHISM BMI development branch I sent to you previously and providing your feedback and comments to us. We greatly value your insight and expertise with the SCHISM model and its source code workflow structure. Here are the following steps we must now take to assimilate the OWP SCHISM (serial) workflow into the NextGen framework using a BMI application:
Advertise discharge/precipitation source terms (from "source.nc" file) and water level open boundaries (from "elev2Dth.nc" file) from their corresponding boundary segments using the NextGen unstructured mesh connectivity data highlighted previously in BMI functions (schism_grid_x, schism_grid_y, schism node/face/edge connectivity). This is essential to enabling the SCHISM BMI functionality for "get_var_indicies" and "set_var_indices" for the boundary features.
Bypass SCHISM requirements to ingest external boundary/source forcing data from the "source.nc" and "elev2Dth.nc" file and instead fill their respective forcing arrays with data ingested from the "set_var" or "set_var_indicies" functionality within the SCHISM BMI framework.
As mentioned below, the best approach I think is to add a new option on SCHISM side.
Test and evaluate the SCHISM BMI functions "get_var" and "set_var" within a NextGen framework test run using other hydrologic models to fill in SCHISM boundary conditions and possibly forcings within a serial model simulation for SCHISM.
In order to complete these previous tasks highlighted above, I had to complete some initial research into the SCHISM source code to understand the boundary condition setup and workflow to utilize within the BMI framework. Here are the following questions I have related to the source code for the 3 initiatives stated previously:
In the case of a serial mode SCHISM simulation, why would the node variable counts "np_global" and "np" differ in this case? My understanding is that "np" is the number of resident elements within a partitioned portion of the mesh. However, when we have a single thread for your simulation, shouldn't np = np_global? Or does np_global account for resident and ghost nodes, in which ghost nodes also relate to boundary nodes?
I cannot reproduce your issue. With 1 MPI process for compute (scribe I/O processes are in a separate comm world), np_global=npa=np. I tested it with both OLDIO on/off and verified this in mirror.out (the beginning part of it shows the partitioning info). Anyway, np is the # of resident NODES (ne is elem) without ghost. npg is # of ghost nodes, npa=np+npg. With 1 MPI process, npg=0, so np=npa=np_global; similarly for elem/sides.
For open water level boundaries, what internal SCHISM variable deals with the actual forcing data? My initial investigation seemed to point to the variable "ath2", which has a number of dimension dependencies (ntracers, nvrt, neta_global, 2, nthfiles2), which include vertical levels and tracers (salinity and temperature). Is this indeed the case? If so, then what dimensions physically point to water levels (given in elev2Dth.nc) for a given open boundary segment? The idea here is that we would like these boundary arrays to be flattened into a one dimensional array in the simplest manner possible without dealing with a multitude of array dimensions that are not dependent on the unstructured mesh connectivity.
Yes ath2() is the receiving array. In addition, the code also needs the time step used in elev2D.th.nc (‘time_step’). The global data are usually defined in Core/schism_glbl.F90. ath2() is single precision.
To reduce memory consumption, I used ath2 to receive b.c. inputs for many variables. In the case of elev, the final array is eth() found around line 1397 in schism_step.F90, where you can see the time interpolation is done from ath2 into eth.
The dimensions in ath2() are: (1) ntracers for tracers so for elev, this is 1; (2) nvrt is vertical indices, so for elev this is also 1 (no vertical structure); (3) neta_global is the lumped total number of nodes that require this input. The order is determined from info in hgrid.gr3 (b.c. part) and bctides.in (elev b.c. flag of 4 or 5); (4) ‘2’ is to account for 2 time records for time interp (i.e. store 2 records from .nc at a time); (5) nthfiles2 points to the variable type, ‘1’ is elev, ‘2’ is velocity etc.
eth(mnond_global,nope_global): 2nd index is the open bnd segment ID and 1st index is the node ID on that particular segment. This should help you flatten the array. You can get the global node ID by using iond_global(), and then use ipgl() to find local node ID (if this rank owns it); see the code near line 4934 in _step.
Mapping from ‘neta_global’ to local node IDs is more involved. You basically loop over all open bnd segments and test to see if the elev flag is 4 or 5. Maybe you can create some mapping arrays on BMI side and I can provide you the code.
For discharge/precipitation source terms, what internal SCHISM variables deal with the actual forcing data for sources? My initial investigation pointed to the variable "ath3", which also has a number of dimension dependencies (nsources, ntracers, 2, nthfiles3). If so, then what dimensions physically point to actual discharge and precipitation source terms (in "source.nc") for a given land boundary segment?
Correct again but ath3(max(1,nsources,nsinks), ntracers,2,nthfiles3). Source/sink can be injected at any element, not just boundary. The IDs of injection elements are found in source.nc. The last dim ‘nthfiles3’ points to 3 types of inputs: volume source, volume sink, and source concentration. (1) nsources/nsinks: total # of source or sink elements; (2) ntracers: in 2D model you’ll only have 2 tracers (T,S) that are not really used. If nthfiles3 points to 1 or 2, this dim is 1 b/c volumes are same for all tracers; if it points to ‘3’, then this dim is 1:2 (T,S); (3) ‘2’ refers to 2 time records; (4) last dim is type explained before. To map the first index to global elem ID, use ieg_source(max(1,nsources)),ieg_sink(max(1,nsinks)). Then use iegl() to get local elem ID. The final receiving arrays are vsource(), msource().
Within the "schism_init.f90", the model attempts to initialize the bounds for a series of arrays for land/open boundary variables for "source.nc" and "elev2Dth.nc" files. In your expert opinions, could we simply just turn off the netcdf read statements for the boundary forcing data and just fill the arrays before a given model iteration? Or are there more dependencies regarding the external forcing data being initialized than my current understanding here?
To bypass SCHISM’s existing options, I think the best way is for me to add a new option that will essentially do nothing, and wait for external program to fill in those needed arrays. This is the approach we used for coupling with atmos model (nws=3). I don’t think there are ways to avoid any code change on SCHISM side, are there? For those 2 .nc inputs, the geometric arrays are all ‘static’ (not changing over time) and ath*() are the only dynamic part so I think you can easily fill them in with pointers.
The documentation within the SCHISM manual helps us to gather the correct variables and point to the global node/element indices number to place SCHISM data in the appropriate array dimensions within the "get_var" and "set_var" functionality, so we should be able to work with that. The idea for the SCHISM BMI is that the interface needs to have the entire unstructured mesh connectivity (np_global, ne_global, ns_global) initialized at the start and then from there we pull data appropriately to and from the model arrays. However, there are likely to be hurdles for integrating the SCHISM data properly into the NextGen model engine. Therefore, I greatly appreciate the communications here for Q/A regarding SCHISM source code and its workflow logic. Thank you very much for your time and have a great rest of your week.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center
Starting with exchanges with Jason D. et al. (OWP)
Hi Joseph,
Adding the CPP flags to the SCHISM compile time sounds like an excellent idea to me. This will greatly assist us in streamlining the latest SCHISM release to our NextGen BMI development. Thank you very much for offering to take on this task.
I just quickly want to follow up on the hotstart question. Essentially, we will have no knowledge of what the source term discharges and the open water boundary levels are upon the NextGen model start time ("t0") from existing model formulations that are coupled with SCHISM. SCHISM will be able to start ingesting source and boundary data in the following hour ("t1") from other model formulations. To my understanding, the "hot_start.nc" file should provide data fields to force the model at the start time ("t0"). You've mentioned that there is another "other_hot_init.nc" that may initialize source and boundary terms at the model start time, or did I misunderstand your comment from the previous email? If there exists this "other_hot_init.nc" for our purposes, is there a certain compiling flag or param.nml variable to turn on and provide this hot_start.nc output?
As always, thank you Joseph for your continued support for the NextGen SCHISM BMI development on behalf of OWP. You and Carsten's insight for SCHISM is critical during our development phases. Have a great rest of your week.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center
On Thu, Apr 20, 2023 at 2:40 PM Y. Joseph Zhang <yjzhang@vims.edu> wrote:
Hi Jason:
I thought about these more during my CA trip. I think we can streamline the options by simply adding CPP flags (USE_NWM_BMI e.g.) without introducing additional parameters. I looked back at the nws=3 option (external atmos model), and think we can actually also replace it with USE_EXT_ATMOS.
To keep the logic straight, I think we should reuse some existing options with those CPP flags. So the users will still need to prepare source.nc if if_source=-1, with the understanding that some info won’t be required and will be overwritten by BMI. If USE_NWM_BMI is on, some reading will be bypassed, while static/geom arrays are still read in. This also helps with the hotstart issue you mentioned below; at t0 the code will read in the first record from source.nc (for elev, this is taken care of by hotstart.nc). I think assuming source/sink=0 at the restart of the run will be OK also as the init condition will be forgotten eventually.
In general, upon hotstart, SCHISM will read in previous ‘state variables’ but not the forcing time series. The latter are done in ‘other_hot_init’.
Let me know if this sounds good to me and I’ll add those CPPs if you agree. Thx.
Regards,
Joseph Zhang
(804)684 7595 (office)
SCHISM web: http://ccrm.vims.edu/schism/
From: Jason Ducker - NOAA Affiliate <jason.ducker@noaa.gov>
Sent: Friday, April 14, 2023 2:06 PM
To: Y. Joseph Zhang <yjzhang@vims.edu>
Cc: carsten.lemmen@hereon.de; Julio Zyserman - NOAA Affiliate <julio.zyserman@noaa.gov>; Camaron George - NOAA Affiliate <camaron.george@noaa.gov>
Subject: Re: SCHISM BMI development Update and Boundary Condition Questions
[EXTERNAL to VIMS received message]
Hi Joseph,
Thank you very much for all of your comments and suggestions. This has already helped me push forward some SCHISM BMI functionality this week. I just have a few quick follow up questions regarding your answers to the previous four questions I postulated in my initial email. My questions to your responses are highlighted in yellow below:
Joseph Zhang response: I cannot reproduce your issue. With 1 MPI process for compute (scribe I/O processes are in a separate comm world), np_global=npa=np. I tested it with both OLDIO on/off and verified this in mirror.out (the beginning part of it shows the partitioning info). Anyway, np is the # of resident NODES (ne is elem) without ghost. npg is # of ghost nodes, npa=np+npg. With 1 MPI process, npg=0, so np=npa=np_global; similarly for elem/sides.
Jason Ducker response: My apologies for mixing this question up with SCHISM here Joseph. I believe this issue was actually occurring in another model (DFlowFM) and not SCHISM. I was able to confirm as well on my end that np_global=npa=np with SCHISM running on a single thread. There are no ghost nodes in this case and this simplifies the initial BMI implementation with SCHISM running in a serial mode (all variables are on the same global nodes within the domain, no local nodes based on partitioned domain assigned to varying threads).
For open water level boundaries, what internal SCHISM variable deals with the actual forcing data? My initial investigation seemed to point to the variable "ath2", which has a number of dimension dependencies (ntracers, nvrt, neta_global, 2, nthfiles2), which include vertical levels and tracers (salinity and temperature). Is this indeed the case? If so, then what dimensions physically point to water levels (given in elev2Dth.nc) for a given open boundary segment? The idea here is that we would like these boundary arrays to be flattened into a one dimensional array in the simplest manner possible without dealing with a multitude of array dimensions that are not dependent on the unstructured mesh connectivity.
Joseph Zhang response: Yes ath2() is the receiving array. In addition, the code also needs the time step used in elev2D.th.nc (‘time_step’). The global data are usually defined in Core/schism_glbl.F90. ath2() is single precision.
To reduce memory consumption, I used ath2 to receive b.c. inputs for many variables. In the case of elev, the final array is eth() found around line 1397 in schism_step.F90, where you can see the time interpolation is done from ath2 into eth.
The dimensions in ath2() are: (1) ntracers for tracers so for elev, this is 1; (2) nvrt is vertical indices, so for elev this is also 1 (no vertical structure); (3) neta_global is the lumped total number of nodes that require this input. The order is determined from info in hgrid.gr3 (b.c. part) and bctides.in (elev b.c. flag of 4 or 5); (4) ‘2’ is to account for 2 time records for time interp (i.e. store 2 records from .nc at a time); (5) nthfiles2 points to the variable type, ‘1’ is elev, ‘2’ is velocity etc.
eth(mnond_global,nope_global): 2nd index is the open bnd segment ID and 1st index is the node ID on that particular segment. This should help you flatten the array. You can get the global node ID by using iond_global(), and then use ipgl() to find local node ID (if this rank owns it); see the code near line 4934 in _step.
Mapping from ‘neta_global’ to local node IDs is more involved. You basically loop over all open bnd segments and test to see if the elev flag is 4 or 5. Maybe you can create some mapping arrays on BMI side and I can provide you the code.
Jason Ducker response: This explanation makes things much clearer with the way we would slice the "ath2" array and update the indices or entire boundary condition variable over each BMI model time step. Thank you for the source code line to see an example of how to obtain the local node ID for an index within the "ath2" array. You also raise a very interesting point regarding how water level boundaries are updated over each model iteration. If we were to bypass the "elev2D.th.nc" forcing file requirement, then the BMI must be able to update water level boundaries at the "start" and "end" time for each model timestep. Essentially, it sounds like the BMI will also need to provide "start time" and "end time" to the variable "th_time2" in order to properly update the ratio variable in line 1382 within "schism_step.F90". My assumption is that we would need to set a SCHISM compiler flag to essentially turn off the updated netcdf read statements for the variables "time_series" and "time" within the "elev2D.th.nc" file right? Since the calculation of the "eth" variable is also dependent on the time component, then it may be more ideal to just update the "ath2" and "th_time2" variables within the BMI "set_var_indices" functions. Let me know though if you think there is a simpler path here than what I proposed.
For discharge/precipitation source terms, what internal SCHISM variables deal with the actual forcing data for sources? My initial investigation pointed to the variable "ath3", which also has a number of dimension dependencies (nsources, ntracers, 2, nthfiles3). If so, then what dimensions physically point to actual discharge and precipitation source terms (in "source.nc") for a given land boundary segment?
Joseph Zhang response: Correct again but ath3(max(1,nsources,nsinks), ntracers,2,nthfiles3). Source/sink can be injected at any element, not just boundary. The IDs of injection elements are found in source.nc. The last dim ‘nthfiles3’ points to 3 types of inputs: volume source, volume sink, and source concentration. (1) nsources/nsinks: total # of source or sink elements; (2) ntracers: in 2D model you’ll only have 2 tracers (T,S) that are not really used. If nthfiles3 points to 1 or 2, this dim is 1 b/c volumes are same for all tracers; if it points to ‘3’, then this dim is 1:2 (T,S); (3) ‘2’ refers to 2 time records; (4) last dim is type explained before. To map the first index to global elem ID, use ieg_source(max(1,nsources)),ieg_sink(max(1,nsinks)). Then use iegl() to get local elem ID. The final receiving arrays are vsource(), msource().
Jason Ducker reponse: Thank you for describing the indexing for the source term. I was able to find the source code accordingly based on your description for finding the local element ID for each source term. As mentioned in my previous response to your question fixed submodule pointer to actually point to
SCHISMinstead of toschism-esmf#2 statement, we also appear to have temporal dependencies with updating the "ath3" source term array (at "t0" and "t1" time) and the time array "th_time3", which originate from the contents of the "source.nc" file. If we were to initiate SCHISM with a hotstart file at the beginning of the simulation, will there already be water levels and source terms available at the open and land boundary nodes respectively? The reason why I ask this is because SCHISM will have no information from the BMI about water level and source terms at "t0" (start time of simulation). Only once models begin executing to the next model time step ("t1") that SCHISM can retrieve boundary conditions fields. Therefore, water level and source term boundary arrays will have to be assumed as values of zero at the start time of the simulation unless data fields are already available at those nodes (within a hotstart file I assume). I'm curious what your take on this issue would be here for SCHISM and whether or not a hotstart file might resolve this issue.Within the "schism_init.f90", the model attempts to initialize the bounds for a series of arrays for land/open boundary variables for "source.nc" and "elev2Dth.nc" files. In your expert opinions, could we simply just turn off the netcdf read statements for the boundary forcing data and just fill the arrays before a given model iteration? Or are there more dependencies regarding the external forcing data being initialized than my current understanding here?
Joseph Zhang response: To bypass SCHISM’s existing options, I think the best way is for me to add a new option that will essentially do nothing, and wait for external program to fill in those needed arrays. This is the approach we used for coupling with atmos model (nws=3). I don’t think there are ways to avoid any code change on SCHISM side, are there? For those 2 .nc inputs, the geometric arrays are all ‘static’ (not changing over time) and ath*() are the only dynamic part so I think you can easily fill them in with pointers.
Jason Ducker response: I completely agree with the approach for adding a SCHISM option with the param.nml file (external_bnd_forcing=1?) to tell the model to bypass external forcing files for boundary conditions and allow the user to depend on coupling software to fill the boundary array values. What's your initial guess on the level of effort this could take for implementing such a SCHISM configuration? To my knowledge, it only appears that we need to include an "#ifelse" SCHISM flag to either include or exclude netcdf inquiry commands to update arrays within "mis_subs.F90" and "schism_step.F90", but I may be missing a few pieces to the puzzle here possibly. I could try to make an attempt on source code modification possibly if you agree with my initial solution to this code approach? To my understanding, we are essentially asking the user only to specify boundary condition locations within the "elev2D.th.nc" and "source.nc" files and bypass any information regarding "time_step" "time" and "time_series"/forcing data.
We greatly appreciate the collaborations of Joseph and Carsten with the OWP development of the SCHISM BMI for the NextGen project. I think we've been making great progress so far on SCHISM and we're getting into some interesting source code topics that may be required to properly enable SCHISM BMI implementation for NextGen. If my questions require a future meeting discussion, then please let me know and I'll be happy to coordinate these efforts moving forward. Have a great weekend.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center
On Thu, Apr 13, 2023 at 12:42 PM Y. Joseph Zhang <yjzhang@vims.edu> wrote:
See below.
-Joseph
Y. Joseph Zhang
Web: schism.wiki
Office: 804 684 7466
From: Jason Ducker - NOAA Affiliate <jason.ducker@noaa.gov>
Sent: Tuesday, April 11, 2023 2:55 PM
To: carsten.lemmen@hereon.de; Y. Joseph Zhang <yjzhang@vims.edu>
Cc: Julio Zyserman - NOAA Affiliate <julio.zyserman@noaa.gov>; Camaron George - NOAA Affiliate <camaron.george@noaa.gov>
Subject: SCHISM BMI development Update and Boundary Condition Questions
[EXTERNAL to VIMS received message]
Hi Joseph and Carsten,
I want to first thank you for taking a look at our SCHISM BMI development branch I sent to you previously and providing your feedback and comments to us. We greatly value your insight and expertise with the SCHISM model and its source code workflow structure. Here are the following steps we must now take to assimilate the OWP SCHISM (serial) workflow into the NextGen framework using a BMI application:
Advertise discharge/precipitation source terms (from "source.nc" file) and water level open boundaries (from "elev2Dth.nc" file) from their corresponding boundary segments using the NextGen unstructured mesh connectivity data highlighted previously in BMI functions (schism_grid_x, schism_grid_y, schism node/face/edge connectivity). This is essential to enabling the SCHISM BMI functionality for "get_var_indicies" and "set_var_indices" for the boundary features.
Bypass SCHISM requirements to ingest external boundary/source forcing data from the "source.nc" and "elev2Dth.nc" file and instead fill their respective forcing arrays with data ingested from the "set_var" or "set_var_indicies" functionality within the SCHISM BMI framework.
In order to complete these previous tasks highlighted above, I had to complete some initial research into the SCHISM source code to understand the boundary condition setup and workflow to utilize within the BMI framework. Here are the following questions I have related to the source code for the 3 initiatives stated previously:
eth(mnond_global,nope_global): 2nd index is the open bnd segment ID and 1st index is the node ID on that particular segment. This should help you flatten the array. You can get the global node ID by using iond_global(), and then use ipgl() to find local node ID (if this rank owns it); see the code near line 4934 in _step.
Mapping from ‘neta_global’ to local node IDs is more involved. You basically loop over all open bnd segments and test to see if the elev flag is 4 or 5. Maybe you can create some mapping arrays on BMI side and I can provide you the code.
The documentation within the SCHISM manual helps us to gather the correct variables and point to the global node/element indices number to place SCHISM data in the appropriate array dimensions within the "get_var" and "set_var" functionality, so we should be able to work with that. The idea for the SCHISM BMI is that the interface needs to have the entire unstructured mesh connectivity (np_global, ne_global, ns_global) initialized at the start and then from there we pull data appropriately to and from the model arrays. However, there are likely to be hurdles for integrating the SCHISM data properly into the NextGen model engine. Therefore, I greatly appreciate the communications here for Q/A regarding SCHISM source code and its workflow logic. Thank you very much for your time and have a great rest of your week.
Cheers,
Jason Ducker
Scientific Programmer, NextGen Team
NOAA/NWS/Office of Water Prediction
National Water Center