diff --git a/CMakeLists.txt b/CMakeLists.txt index 11d68f8e3..c373ed951 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,6 +14,7 @@ set (CMAKE_C_STANDARD 11) include (CheckCCompilerFlag) include (CheckCXXCompilerFlag) +include(GNUInstallDirs) # Check if compiler supports extra debugging flags foreach (DEBUG_FLAG "-ggdb" "-g3") @@ -727,6 +728,10 @@ if ( ${PARFLOW_HAVE_CLM} ) configure_file (pfsimulator/clm/parflow_config.F90.in ${PROJECT_BINARY_DIR}/pfsimulator/clm/parflow_config.F90) endif ( ${PARFLOW_HAVE_CLM} ) +# Install headers +install(FILES ${CMAKE_BINARY_DIR}/include/parflow_config.h + ${CMAKE_BINARY_DIR}/include/pfversion.h + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/parflow) #----------------------------------------------------------------------------- # CMAKE Subdirectories #----------------------------------------------------------------------------- diff --git a/README.md b/README.md index 30eca2474..58dbb18a5 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,14 @@ +# Fork of ParFlow pre-patched for TSMP-PDAF + +**DO NOT USE WITHOUT PDAF** + +This is a fork of ParFlow with changes added to run ParFlow in the +TSMP-PDAF framework. + +It is built via the TSMP2 repository: +https://github.com/HPSCTerrSys/TSMP2 + + # ParFlow ![ParFlow Linux CI Test](https://github.com/parflow/parflow/actions/workflows/linux.yml/badge.svg) diff --git a/cmake/modules/FindOASIS.cmake b/cmake/modules/FindOASIS.cmake index a876f5109..ca42f127a 100644 --- a/cmake/modules/FindOASIS.cmake +++ b/cmake/modules/FindOASIS.cmake @@ -22,5 +22,6 @@ if(OASIS_FOUND AND NOT TARGET OASIS3MCT::OASIS3MCT) target_link_libraries(OASIS3MCT::OASIS3MCT INTERFACE ${OASIS_Fortran_LIBRARY} ${MCT_Fortran_LIBRARY} ${MPEU_Fortran_LIBRARY} ${SCRIP_Fortran_LIBRARY}) target_link_libraries(OASIS3MCT::OASIS3MCT INTERFACE OpenMP::OpenMP_Fortran) target_link_libraries(OASIS3MCT::OASIS3MCT INTERFACE ${NETCDF_Fortran_LIBRARY}) + target_include_directories(OASIS3MCT::OASIS3MCT INTERFACE ${NETCDF_Fortran_INCLUDES}) endif() diff --git a/pfsimulator/amps/CMakeLists.txt b/pfsimulator/amps/CMakeLists.txt index 208a98097..03a8cd7a0 100644 --- a/pfsimulator/amps/CMakeLists.txt +++ b/pfsimulator/amps/CMakeLists.txt @@ -23,6 +23,10 @@ set(COMMON_SRC_FILES amps_print.c ) +list(APPEND HEADER_FILES common/amps_common.h + ${PARFLOW_AMPS_LAYER}/amps.h + ${PARFLOW_AMPS_LAYER}/amps_proto.h) + # Add RMM / Umpire wrappers if enabled if (${PARFLOW_HAVE_RMM}) list(APPEND COMMON_SRC_FILES @@ -111,7 +115,8 @@ endif() if(${PARFLOW_HAVE_OAS3}) target_link_libraries(amps PUBLIC OASIS3MCT::OASIS3MCT) -endif() + list(APPEND HEADER_FILES oas3/oas3_coupler.h oas3/oas3_external.h) +endif(${PARFLOW_HAVE_OAS3}) # ----------------------------- # Install targets and scripts @@ -123,5 +128,9 @@ set(AMPS_SCRIPTS run) string(REGEX REPLACE "([^;]+)" "${PARFLOW_AMPS_LAYER}/\\1" AMPS_SCRIPTS "${AMPS_SCRIPTS}") install(FILES ${AMPS_SCRIPTS} DESTINATION bin) +# Install headers +install(FILES ${HEADER_FILES} + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/parflow) + # Add subdirectory for tests add_subdirectory(test/src) diff --git a/pfsimulator/amps/mpi1/amps_finalize.c b/pfsimulator/amps/mpi1/amps_finalize.c index e6f02b2dc..5a16a4b82 100644 --- a/pfsimulator/amps/mpi1/amps_finalize.c +++ b/pfsimulator/amps/mpi1/amps_finalize.c @@ -68,8 +68,7 @@ int amps_Finalize() MPI_Comm_free(&s_CommNode); MPI_Comm_free(&s_CommWrite); MPI_Comm_free(&s_CommWorld); - - MPI_Finalize(); + /* MPI_Finalize(); */ } #if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS) amps_gpu_finalize(); diff --git a/pfsimulator/amps/mpi1/amps_init.c b/pfsimulator/amps/mpi1/amps_init.c index ac5e7dbee..29a8e755a 100644 --- a/pfsimulator/amps/mpi1/amps_init.c +++ b/pfsimulator/amps/mpi1/amps_init.c @@ -35,6 +35,8 @@ #include #include +MPI_Comm dacomm; /* kuw */ + /* Global flag indicating if AMPS has been initialized */ int amps_mpi_initialized = FALSE; @@ -114,10 +116,10 @@ int amps_Init(int *argc, char **argv[]) int length; #endif - MPI_Init(argc, argv); + /* MPI_Init(argc, argv); */ amps_mpi_initialized = TRUE; - MPI_Comm_dup(MPI_COMM_WORLD, &s_CommWorld); + MPI_Comm_dup(dacomm, &s_CommWorld); MPI_Comm_size(amps_CommWorld, &s_size); MPI_Comm_rank(amps_CommWorld, &s_rank); diff --git a/pfsimulator/parflow_lib/CMakeLists.txt b/pfsimulator/parflow_lib/CMakeLists.txt index eb55eafce..c332a1674 100644 --- a/pfsimulator/parflow_lib/CMakeLists.txt +++ b/pfsimulator/parflow_lib/CMakeLists.txt @@ -323,3 +323,69 @@ endif (${PARFLOW_HAVE_PDI}) if (${PARFLOW_HAVE_NETCDF}) target_link_libraries (pfsimulator NetCDF::NetCDF) endif (${PARFLOW_HAVE_NETCDF}) + +# Install headers +set (HEADER_FILES + backend_mapping.h + background.h + bc_pressure.h + char_vector.h + clustering.h + communication.h + computation.h + drv_gridmodule.h + drv_module.h + drv_tilemodule.h + f2c.h + fgetopt.h + file_versions.h + general.h + geometry.h + geostats.h + globals.h + grgeometry.h + grgeom_list.h + grgeom_octree.h + grid.h + hbt.h + hypre_dependences.h + index_space.h + info_header.h + input_database.h + kinsol_dependences.h + lb.h + llnlmath.h + llnltyps.h + logging.h + loops.h + matrix.h + metadata.h + nl_function_eval.h + n_vector.h + Parflow.hxx + parflow.h + parflow_netcdf.h + parflow_proto_f90.h + parflow_proto_f.h + parflow_proto.h + pf_cudaloops.h + pf_cudamalloc.h + pf_devices.h + pf_hypre.h + pf_kokkosloops.h + pf_kokkosmalloc.h + pf_module.h + pf_omploops.h + problem_bc.h + problem_eval.h + problem.h + region.h + solver.h + time_cycle_data.h + timing.h + vector.h + well.h + reservoir.h +) +install(FILES ${HEADER_FILES} + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/parflow) diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index 27a81f4c7..05f82d3c3 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -978,6 +978,8 @@ void PhaseRelPermFreeInstanceXtra(void); PFModule *PhaseRelPermNewPublicXtra(void); void PhaseRelPermFreePublicXtra(void); int PhaseRelPermSizeOfTempData(void); +Vector *PhaseRelPermGetAlpha(PFModule *this_module); +Vector *PhaseRelPermGetN(PFModule *this_module); typedef void (*PhaseSourceInvoke) (Vector *phase_source, int phase, Problem *problem, ProblemData *problem_data, double time); typedef PFModule *(*PhaseSourceNewPublicXtraInvoke) (int num_phases); @@ -1031,6 +1033,8 @@ typedef PFModule *(*SaturationOutputStaticInvoke) (char *file_prefix, ProblemDat /* problem_saturation.c */ void Saturation(Vector *phase_saturation, Vector *phase_pressure, Vector *phase_density, double gravity, ProblemData *problem_data, int fcn); +Vector *SaturationGetAlpha(PFModule *this_module); +Vector *SaturationGetN(PFModule *this_module); PFModule *SaturationInitInstanceXtra(Grid *grid, double *temp_data); void SaturationFreeInstanceXtra(void); PFModule *SaturationNewPublicXtra(void); @@ -1238,6 +1242,8 @@ PFModule *SolverRichardsNewPublicXtra(char *name); void SolverRichardsFreePublicXtra(void); int SolverRichardsSizeOfTempData(void); ProblemData *GetProblemDataRichards(PFModule *this_module); +PFModule *GetPhaseRelPerm(PFModule *this_module); +PFModule *GetSaturation(PFModule *this_module); Problem *GetProblemRichards(PFModule *this_module); PFModule *GetICPhasePressureRichards(PFModule *this_module); Grid *GetGrid2DRichards(PFModule *this_module); @@ -1251,6 +1257,16 @@ void AdvanceRichards(PFModule *this_module, Vector ** porosity_out, Vector ** saturation_out ); +// PDAF: this function is for initialization of OASIS only +void PseudoAdvanceRichards(PFModule *this_module, + double start_time, /* Starting time */ + double stop_time, /* Stopping time */ + PFModule *time_step_control, /* Use this module to control timestep if supplied */ + Vector * evap_trans, /* Flux from land surface model */ + Vector ** pressure_out, /* Output vars */ + Vector ** porosity_out, + Vector ** saturation_out + ); void ExportRichards(PFModule *this_module, Vector ** pressure_out, /* Output vars */ Vector ** porosity_out, diff --git a/pfsimulator/parflow_lib/problem_phase_rel_perm.c b/pfsimulator/parflow_lib/problem_phase_rel_perm.c index f0838a414..827f91ab3 100644 --- a/pfsimulator/parflow_lib/problem_phase_rel_perm.c +++ b/pfsimulator/parflow_lib/problem_phase_rel_perm.c @@ -2179,3 +2179,31 @@ int PhaseRelPermSizeOfTempData() return sz; } + +/*----------------------------------------------------------------------- + * PhaseRelPermGetAlpha + *-----------------------------------------------------------------------*/ + +Vector *PhaseRelPermGetAlpha(PFModule *this_module) +{ + PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); + + Type1 *dummy1; + dummy1 = (Type1*)(public_xtra->data); + + return (dummy1 -> alpha_values); +} + +/*----------------------------------------------------------------------- + * PhaseRelPermGetN + *-----------------------------------------------------------------------*/ + +Vector *PhaseRelPermGetN(PFModule *this_module) +{ + PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); + + Type1 *dummy1; + dummy1 = (Type1*)(public_xtra->data); + + return (dummy1 -> n_values); +} diff --git a/pfsimulator/parflow_lib/problem_saturation.c b/pfsimulator/parflow_lib/problem_saturation.c index c5dda41c0..b764e6ea4 100644 --- a/pfsimulator/parflow_lib/problem_saturation.c +++ b/pfsimulator/parflow_lib/problem_saturation.c @@ -654,6 +654,34 @@ void Saturation( } /* End switch */ } +/*-------------------------------------------------------------------------- + * SaturationGetAlpha + *--------------------------------------------------------------------------*/ + +Vector *SaturationGetAlpha(PFModule *this_module) +{ + PublicXtra *public_xtra = (PublicXtra *)PFModulePublicXtra(this_module); + + Type1 *dummy1; + dummy1 = (Type1 *)(public_xtra -> data); + + return (dummy1 -> alpha_values); +} + +/*-------------------------------------------------------------------------- + * SaturationGetN + *--------------------------------------------------------------------------*/ + +Vector *SaturationGetN(PFModule *this_module) +{ + PublicXtra *public_xtra = (PublicXtra *)PFModulePublicXtra(this_module); + + Type1 *dummy1; + dummy1 = (Type1 *)(public_xtra -> data); + + return (dummy1 -> n_values); +} + /*-------------------------------------------------------------------------- * SaturationInitInstanceXtra *--------------------------------------------------------------------------*/ diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 1b8641c2c..e1088d2b0 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -1942,10 +1942,12 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time Subgrid *subgrid; Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; double *pp, *sp, *et, *ms, *po_dat, *dz_dat; - double sw_lat = .0; - double sw_lon = .0; + /* double sw_lat = .0; */ + /* double sw_lon = .0; */ #endif +int istep = 1; + #ifdef HAVE_CLM Grid *grid = (instance_xtra->grid); Subgrid *subgrid; @@ -1967,9 +1969,9 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time int fflag, fstart, fstop; // IMF: index w/in 3D forcing array corresponding to istep int n, c; // IMF: index vars for looping over subgrid data BH: added c int ind_veg; /*BH: temporary variable to store vegetation index */ - int Stepcount = 0; /* Added for transient EvapTrans file management - NBE */ - int Loopcount = 0; /* Added for transient EvapTrans file management - NBE */ - double sw = NAN, lw = NAN, prcp = NAN, tas = NAN, u = NAN, v = NAN, patm = NAN, qatm = NAN; // IMF: 1D forcing vars (local to AdvanceRichards) + /* int Stepcount = 0; /\* Added for transient EvapTrans file management - NBE *\/ */ + /* int Loopcount = 0; /\* Added for transient EvapTrans file management - NBE *\/ */ + double sw=NAN, lw=NAN, prcp=NAN, tas=NAN, u=NAN, v=NAN, patm=NAN, qatm=NAN; // IMF: 1D forcing vars (local to AdvanceRichards) double lai[18], sai[18], z0m[18], displa[18]; /*BH: array with lai/sai/z0m/displa values for each veg class */ double *sw_data = NULL; double *lw_data = NULL; @@ -2033,37 +2035,41 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time char file_prefix[2048], file_type[2048], file_postfix[2048]; char nc_postfix[2048]; + int Stepcount = 0; /* Added for transient EvapTrans file management - NBE */ + int Loopcount = 0; /* Added for transient EvapTrans file management - NBE */ int first_tstep = 1; sprintf(file_prefix, "%s", GlobalsOutFileName); //CPS oasis definition phase #ifdef HAVE_OAS3 - int nlon = GetInt("ComputationalGrid.NX"); - int nlat = GetInt("ComputationalGrid.NY"); - double pfl_step = GetDouble("TimeStep.Value"); - double pfl_stop = GetDouble("TimingInfo.StopTime"); + /* int nlon = GetInt("ComputationalGrid.NX"); */ + /* int nlat = GetInt("ComputationalGrid.NY"); */ + /* double pfl_step = GetDouble("TimeStep.Value"); */ + /* double pfl_stop = GetDouble("TimingInfo.StopTime"); */ + /* // PDAF: getting start time */ + /* double pfl_start = GetDouble("TimingInfo.StartTime"); */ int is; - ForSubgridI(is, GridSubgrids(grid)) - { - double dx, dy; - int nx, ny, ix, iy; + /* ForSubgridI(is, GridSubgrids(grid)) */ + /* { */ + /* double dx, dy; */ + /* int nx, ny, ix, iy; */ - subgrid = GridSubgrid(grid, is); + /* subgrid = GridSubgrid(grid, is); */ - nx = SubgridNX(subgrid); - ny = SubgridNY(subgrid); + /* nx = SubgridNX(subgrid); */ + /* ny = SubgridNY(subgrid); */ - ix = SubgridIX(subgrid); - iy = SubgridIY(subgrid); + /* ix = SubgridIX(subgrid); */ + /* iy = SubgridIY(subgrid); */ - dx = SubgridDX(subgrid); - dy = SubgridDY(subgrid); + /* dx = SubgridDX(subgrid); */ + /* dy = SubgridDY(subgrid); */ - CALL_oas_pfl_define(nx, ny, dx, dy, ix, iy, sw_lon, sw_lat, nlon, nlat, - pfl_step, pfl_stop); - } + // CALL_oas_pfl_define(nx, ny, dx, dy, ix, iy, sw_lon, sw_lat, nlon, nlat, + // pfl_step, pfl_stop); + /* } */ amps_Sync(amps_CommWorld); #endif // end to HAVE_OAS3 CALL @@ -2088,6 +2094,23 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time } dt = cdt; +/* #ifdef FOR2131 */ +/* PFModuleInvokeType(SaturationInvoke, problem_saturation, */ +/* (instance_xtra -> saturation, instance_xtra -> pressure, */ +/* instance_xtra -> density, gravity, problem_data, */ +/* CALCFCN)); */ + +/* handle = InitVectorUpdate(instance_xtra -> saturation, VectorUpdateAll); */ +/* FinalizeVectorUpdate(handle); */ + +/* //sprintf(file_postfix, "pressIC.%05d",instance_xtra -> file_number); */ +/* // WritePFBinary(file_prefix, file_postfix, instance_xtra -> pressure); */ +/* //sprintf(file_postfix, "saturIC.%05d",instance_xtra -> file_number); */ +/* // WritePFBinary(file_prefix, file_postfix, instance_xtra -> saturation); */ +/* // sprintf(file_postfix, "permIC.%05d",instance_xtra -> file_number); */ +/* // WritePFBinary(file_prefix, file_postfix, ProblemDataPermeabilityX(problem_data)); */ + +/* #endif */ /* * Check to see if pressure solves are requested * start_count < 0 implies that subsurface data ONLY is requested @@ -2115,6 +2138,19 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time { if (t == ct) { + if (time_step_control) + { + PFModuleInvokeType(SelectTimeStepInvoke, time_step_control, + (&dt, &dt_info, t, problem, + problem_data)); + } + else + { + PFModuleInvokeType(SelectTimeStepInvoke, select_time_step, + (&dt, &dt_info, t, problem, + problem_data)); + } + ct += cdt; // Read in evap_trans even if not HAVE_CLM @@ -2882,8 +2918,35 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); -#endif -#ifdef HAVE_CLM + + + + + //#endif //End of call to CLM + /* NBE counter for reusing CLM input files */ + clm_next += 1; + if (clm_next > clm_skip) + { + istep = istep + 1; + clm_next = 1; + } // NBE + + //istep = istep + 1; + + EndTiming(CLMTimingIndex); + + + /* ============================================================= + * NBE: It looks like the time step isn't really scaling the CLM + * inputs, but the looping flag is working as intended as + * of 2014-04-06. + * + * It is using the different time step counter BUT then it + * isn't scaling the inputs properly. + * ============================================================= */ +#endif // End of call to CLM + + istep = istep + 1; /******************************************/ /* read transient evap trans flux file */ /******************************************/ @@ -2937,28 +3000,28 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time } - /* NBE counter for reusing CLM input files */ - clm_next += 1; - if (clm_next > clm_skip) - { - istep = istep + 1; - clm_next = 1; - } // NBE +// /* NBE counter for reusing CLM input files */ +// clm_next += 1; +// if (clm_next > clm_skip) +// { +// istep = istep + 1; +// clm_next = 1; +// } // NBE - //istep = istep + 1; +// //istep = istep + 1; - EndTiming(CLMTimingIndex); +// EndTiming(CLMTimingIndex); - /* ============================================================= - * NBE: It looks like the time step isn't really scaling the CLM - * inputs, but the looping flag is working as intended as - * of 2014-04-06. - * - * It is using the different time step counter BUT then it - * isn't scaling the inputs properly. - * ============================================================= */ -#endif +// /* ============================================================= +// * NBE: It looks like the time step isn't really scaling the CLM +// * inputs, but the looping flag is working as intended as +// * of 2014-04-06. +// * +// * It is using the different time step counter BUT then it +// * isn't scaling the inputs properly. +// * ============================================================= */ +// #endif } //Endif to check whether an entire dt is complete converged = 1; @@ -2980,18 +3043,18 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time /*******************************************************************/ if (converged) { - if (time_step_control) - { - PFModuleInvokeType(SelectTimeStepInvoke, time_step_control, - (&dt, &dt_info, t, problem, - problem_data)); - } - else - { - PFModuleInvokeType(SelectTimeStepInvoke, select_time_step, - (&dt, &dt_info, t, problem, - problem_data)); - } + // if (time_step_control) + // { + // PFModuleInvokeType(SelectTimeStepInvoke, time_step_control, + // (&dt, &dt_info, t, problem, + // problem_data)); + // } + // else + // { + // PFModuleInvokeType(SelectTimeStepInvoke, select_time_step, + // (&dt, &dt_info, t, problem, + // problem_data)); + // } PFVCopy(instance_xtra->density, instance_xtra->old_density); PFVCopy(instance_xtra->saturation, @@ -4852,6 +4915,187 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time *saturation_out = instance_xtra->saturation; } +// Compare to AdvanceRichards() in parflow v3.9.0 +void +PseudoAdvanceRichards(PFModule * this_module, double start_time, /* Starting time */ + double stop_time, /* Stopping time */ + PFModule * time_step_control, /* Use this module to control timestep if supplied */ + Vector * evap_trans, /* Flux from land surface model */ + Vector ** pressure_out, /* Output vars */ + Vector ** porosity_out, Vector ** saturation_out) +{ + //>>TSMP-PDAF internal change beginning (compare to AdvanceRichards) + printf("Pseudo richard for OAS init\n"); + //>>TSMP-PDAF internal change end + + /* PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); */ + InstanceXtra *instance_xtra = + (InstanceXtra*)PFModuleInstanceXtra(this_module); + /* Problem *problem = (public_xtra->problem); */ + + /* int max_iterations = (public_xtra->max_iterations); */ + /* int print_satur = (public_xtra->print_satur); */ + /* int print_wells = (public_xtra->print_wells); */ + + /* PFModule *problem_saturation = (instance_xtra->problem_saturation); */ + /* PFModule *phase_density = (instance_xtra->phase_density); */ + /* PFModule *select_time_step = (instance_xtra->select_time_step); */ + /* PFModule *l2_error_norm = (instance_xtra->l2_error_norm); */ + /* PFModule *nonlin_solver = (instance_xtra->nonlin_solver); */ + + /* ProblemData *problem_data = (instance_xtra->problem_data); */ + + /* int start_count = ProblemStartCount(problem); */ + /* double dump_interval = ProblemDumpInterval(problem); */ + + /* Vector *porosity = ProblemDataPorosity(problem_data); */ + /* Vector *evap_trans_sum = instance_xtra->evap_trans_sum; */ + /* Vector *overland_sum = instance_xtra->overland_sum; /\* sk: Vector of outflow at the boundary *\/ */ + + if (evap_trans == NULL) + { + evap_trans = instance_xtra->evap_trans; + } + +#ifdef HAVE_OAS3 + Grid *grid = (instance_xtra->grid); + Subgrid *subgrid; + /* Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; */ + /* double *pp, *sp, *et, *ms, *po_dat, *dz_dat; */ + double sw_lat = .0; + double sw_lon = .0; +#endif + +//>>TSMP-PDAF internal change beginning (compare to AdvanceRichards) +// #ifdef HAVE_CLM +// Grid *grid = (instance_xtra->grid); +// Subgrid *subgrid; +// Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; +// double *pp, *sp, *et, *ms, *po_dat, *dz_dat; + +// /* IMF: For CLM met forcing (local to AdvanceRichards) */ +// int istep; // IMF: counter for clm output times + +// /* NBE added for clm reuse of inputs */ +// int clm_next = 1; //NBE: Counter for reuse loop +// int clm_skip = public_xtra->clm_reuse_count; // NBE:defaults to 1 +// int clm_write_logs = public_xtra->clm_write_logs; // NBE: defaults to 1, disables log file writing if 0 +// int clm_last_rst = public_xtra->clm_last_rst; // Reuse of the RST file +// int clm_daily_rst = public_xtra->clm_daily_rst; // Daily or hourly RST files, defaults to daily + +// int fstep = INT_MIN; +// int fflag, fstart, fstop; // IMF: index w/in 3D forcing array corresponding to istep +// int n, c; // IMF: index vars for looping over subgrid data BH: added c +// int ind_veg; /*BH: temporary variable to store vegetation index */ +// double sw=NAN, lw=NAN, prcp=NAN, tas=NAN, u=NAN, v=NAN, patm=NAN, qatm=NAN; // IMF: 1D forcing vars (local to AdvanceRichards) +// double lai[18], sai[18], z0m[18], displa[18]; /*BH: array with lai/sai/z0m/displa values for each veg class */ +// double *sw_data = NULL; +// double *lw_data = NULL; +// double *prcp_data = NULL; // IMF: 2D forcing vars (SubvectorData) (local to AdvanceRichards) +// double *tas_data = NULL; +// double *u_data = NULL; +// double *v_data = NULL; +// double *patm_data = NULL; +// double *qatm_data = NULL; +// double *lai_data = NULL; +// /*BH*/ double *sai_data = NULL; +// /*BH*/ double *z0m_data = NULL; +// /*BH*/ double *displa_data = NULL; +// /*BH*/ double *veg_map_data = NULL; +// /*BH*/ /*will fail if veg_map_data is declared as int */ +// char filename[2048]; // IMF: 1D input file name *or* 2D/3D input file base name +// Subvector *sw_forc_sub, *lw_forc_sub, *prcp_forc_sub, *tas_forc_sub, *u_forc_sub, *v_forc_sub, *patm_forc_sub, *qatm_forc_sub, *lai_forc_sub, *sai_forc_sub, *z0m_forc_sub, *displa_forc_sub, *veg_map_forc_sub; /*BH: added LAI/SAI/Z0M/DISPLA/vegmap */ + +// /* Slopes */ +// Subvector *slope_x_sub, *slope_y_sub; +// double *slope_x_data, *slope_y_data; + +// /* IMF: For writing CLM output */ +// Subvector *eflx_lh_tot_sub, *eflx_lwrad_out_sub, *eflx_sh_tot_sub, +// *eflx_soil_grnd_sub, *qflx_evap_tot_sub, *qflx_evap_grnd_sub, +// *qflx_evap_soi_sub, *qflx_evap_veg_sub, *qflx_tran_veg_sub, +// *qflx_infl_sub, *swe_out_sub, *t_grnd_sub, *tsoil_sub, *irr_flag_sub, +// *qflx_qirr_sub, *qflx_qirr_inst_sub; + +// double *eflx_lh, *eflx_lwrad, *eflx_sh, *eflx_grnd, *qflx_tot, *qflx_grnd, +// *qflx_soi, *qflx_eveg, *qflx_tveg, *qflx_in, *swe, *t_g, *t_soi, *iflag, +// *qirr, *qirr_inst; +// int clm_file_dir_length; +// #endif +//>>TSMP-PDAF internal change end + + /* int any_file_dumped; */ + /* int clm_file_dumped; */ + /* int dump_files = 0; */ + + /* int retval; */ + /* int converged; */ + /* int take_more_time_steps; */ + /* int conv_failures; */ + /* int max_failures = public_xtra->max_convergence_failures; */ + + /* double t; */ + /* double dt = 0.0; */ + /* double ct = 0.0; */ + /* double cdt = 0.0; */ + /* double print_dt; */ + /* double dtmp, err_norm; */ + /* double gravity = ProblemGravity(problem); */ + + /* VectorUpdateCommHandle *handle; */ + + /* char dt_info; */ + char file_prefix[2048];//, file_type[2048], file_postfix[2048]; + /* char nc_postfix[2048]; */ + +//>>TSMP-PDAF internal change beginning (compare to AdvanceRichards) + // int Stepcount = 0; /* Added for transient EvapTrans file management - NBE */ + // int Loopcount = 0; /* Added for transient EvapTrans file management - NBE */ + // int first_tstep = 1; +//>>TSMP-PDAF internal change end + + sprintf(file_prefix, "%s", GlobalsOutFileName); + + //CPS oasis definition phase +#ifdef HAVE_OAS3 + int nlon = GetInt("ComputationalGrid.NX"); + int nlat = GetInt("ComputationalGrid.NY"); + double pfl_step = GetDouble("TimeStep.Value"); + double pfl_stop = GetDouble("TimingInfo.StopTime"); + //>>TSMP-PDAF internal change beginning (compare to AdvanceRichards) + /* double pfl_start = GetDouble("TimingInfo.StartTime"); */ + //>>TSMP-PDAF internal change end + + int is; + ForSubgridI(is, GridSubgrids(grid)) + { + double dx, dy; + int nx, ny, ix, iy; + + subgrid = GridSubgrid(grid, is); + + nx = SubgridNX(subgrid); + ny = SubgridNY(subgrid); + + ix = SubgridIX(subgrid); + iy = SubgridIY(subgrid); + + dx = SubgridDX(subgrid); + dy = SubgridDY(subgrid); + + //>>TSMP-PDAF internal change beginning (compare to AdvanceRichards) + // gw only init once + // kuw + CALL_oas_pfl_define(nx, ny, dx, dy, ix, iy, sw_lon, sw_lat, nlon, nlat, + pfl_step, pfl_stop); + // kuw end + // gw end + //>>TSMP-PDAF internal change end + } + amps_Sync(amps_CommWorld); +#endif // end to HAVE_OAS3 CALL +} + void ExportRichards(PFModule * this_module, Vector ** pressure_out, /* Output vars */ @@ -7107,6 +7351,16 @@ GetProblemDataRichards(PFModule * this_module) return(instance_xtra->problem_data); } +PFModule *GetPhaseRelPerm(PFModule *this_module){ + InstanceXtra *instance_xtra = (InstanceXtra *)PFModuleInstanceXtra(this_module); + return (instance_xtra -> phase_rel_perm); +} + +PFModule *GetSaturation(PFModule *this_module){ + InstanceXtra *instance_xtra = (InstanceXtra *)PFModuleInstanceXtra(this_module); + return (instance_xtra -> problem_saturation); +} + Problem * GetProblemRichards(PFModule * this_module) { @@ -7124,6 +7378,46 @@ GetICPhasePressureRichards(PFModule * this_module) return(instance_xtra->ic_phase_pressure); } +Vector *GetPressureRichards(PFModule *this_module) +{ + InstanceXtra *instance_xtra = + (InstanceXtra *)PFModuleInstanceXtra(this_module); + + return (instance_xtra -> pressure); +} + +Vector *GetSaturationRichards(PFModule *this_module) +{ + InstanceXtra *instance_xtra = + (InstanceXtra *)PFModuleInstanceXtra(this_module); + + return (instance_xtra -> saturation); +} + +Vector *GetDensityRichards(PFModule *this_module) +{ + InstanceXtra *instance_xtra = + (InstanceXtra *)PFModuleInstanceXtra(this_module); + + return (instance_xtra -> density); +} + +int GetEvapTransFile(PFModule *this_module) +{ + PublicXtra *public_xtra = + (PublicXtra *)PFModulePublicXtra(this_module); + + return (public_xtra -> evap_trans_file); +} + +char *GetEvapTransFilename(PFModule *this_module) +{ + PublicXtra *public_xtra = + (PublicXtra *)PFModulePublicXtra(this_module); + + return (public_xtra -> evap_trans_filename); +} + Grid * GetGrid2DRichards(PFModule * this_module) { diff --git a/pfsimulator/parflow_lib/wrf_parflow.c b/pfsimulator/parflow_lib/wrf_parflow.c index cad5b77c2..f85bde0dd 100644 --- a/pfsimulator/parflow_lib/wrf_parflow.c +++ b/pfsimulator/parflow_lib/wrf_parflow.c @@ -36,7 +36,8 @@ #include -amps_ThreadLocalDcl(PFModule *, Solver_module); +extern PFModule *Solver_module; /* defined in solver.c */ +/* amps_ThreadLocalDcl(PFModule *, Solver_module); */ amps_ThreadLocalDcl(PFModule *, solver); amps_ThreadLocalDcl(Vector *, evap_trans);