From 0f99f27bd65b85400db686e8dee1ca7422015212 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Thu, 6 Apr 2023 11:07:42 +0200 Subject: [PATCH 01/19] Added soil ice fraction parameter to Richards solver --- pfsimulator/amps/oas3/oas3_coupler.h | 6 ++--- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 24 ++++++++++++------- pfsimulator/amps/oas3/receive_fld2_clm.F90 | 5 ++-- .../parflow_lib/kinsol_nonlin_solver.c | 3 ++- pfsimulator/parflow_lib/nl_function_eval.c | 17 +++++++------ pfsimulator/parflow_lib/nl_function_eval.h | 2 ++ pfsimulator/parflow_lib/parflow_proto.h | 12 +++++----- .../parflow_lib/problem_phase_rel_perm.c | 3 ++- .../parflow_lib/richards_jacobian_eval.c | 4 ++-- pfsimulator/parflow_lib/solver_richards.c | 14 +++++++---- 10 files changed, 55 insertions(+), 35 deletions(-) diff --git a/pfsimulator/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index 94113032a..aaf7e47c9 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -28,7 +28,7 @@ void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep, double *porosity, double *dz); #define receive_fld2_clm receive_fld2_clm_ -#define CALL_receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ - receive_fld2_clm(evap_trans, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) -void receive_fld2_clm(double *evap_trans, double *topo, int *ix, int*iy, +#define CALL_receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ + receive_fld2_clm(evap_trans, ice_fraction, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) +void receive_fld2_clm(double *evap_trans, double *ice_fraction, double *topo, int *ix, int*iy, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep); diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 6ee3b78bc..938cac7fb 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -176,16 +176,18 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n deallocate(pressure_3d) end subroutine send_fld2_clm - subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') + subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') !----------------------------------------------- ! Receives evapotranspiration fluxes from eCLM. !----------------------------------------------- - integer, intent(in) :: ix, iy, & ! starting coordinate of the subgrid (UNUSED) - nx, ny, nz, & ! subgrid dimensions - nx_f, ny_f ! dimensions of ET subvector - real(kind=8), intent(in) :: pstep, & ! current model time in hours - topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) - real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! evapotranspiration [1/hrs] + integer, intent(in) :: ix, iy, & ! starting coordinate of the subgrid (UNUSED) + nx, ny, nz, & ! subgrid dimensions + nx_f, ny_f ! dimensions of ET subvector + real(kind=8), intent(in) :: pstep, & ! current model time in hours + topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) + real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! evapotranspiration [1/hrs] + real(kind=8), intent(inout) :: ice_fraction((nx+2)*(ny+2)*(nz+2)) ! soil ice fraction + ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the ! extra "+2" terms account for the ghost nodes/halo points @@ -196,9 +198,12 @@ subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, ps integer :: z ! subsurface level (z=nz topmost layer, z=1 deepest layer) integer :: top_z_level(nx,ny) ! topmost z level of active ParFlow cells real(kind=8), allocatable :: evap_trans_3d(:,:,:) ! Root ET fluxes received from eCLM [1/hrs] + real(kind=8), allocatable :: ice_frac_3d(:,:,:) ! Soil ice fraction received from eCLM [unitless; values from 0 to 1] ! Receive ET fluxes from eCLM allocate(evap_trans_3d(nx,ny,nlevsoi)) + allocate(ice_frac_3d(nx,ny,nlevsoi)) + seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) @@ -211,13 +216,14 @@ subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, ps do k = 1, nlevsoi ! eCLM: 1=topmost layer, nlevsoi=deepest layer z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) - evap_trans(l) = evap_trans_3d(i,j,k) + evap_trans(l) = evap_trans_3d(i,j,k) + ice_fraction(l) = ice_frac_3d(i,j,k) end do end if end do end do - deallocate(evap_trans_3d) + deallocate(evap_trans_3d, ice_frac_3d) end subroutine receive_fld2_clm function get_top_z_level(nx, ny, nz, topo) diff --git a/pfsimulator/amps/oas3/receive_fld2_clm.F90 b/pfsimulator/amps/oas3/receive_fld2_clm.F90 index 4143cfb83..4ef3acadd 100644 --- a/pfsimulator/amps/oas3/receive_fld2_clm.F90 +++ b/pfsimulator/amps/oas3/receive_fld2_clm.F90 @@ -1,4 +1,4 @@ -SUBROUTINE receive_fld2_clm(evap_trans,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) +SUBROUTINE receive_fld2_clm(evap_trans,ice_fraction,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) !---------------------------------------------------------------------------- ! @@ -46,7 +46,7 @@ SUBROUTINE receive_fld2_clm(evap_trans,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) REAL(KIND=8), INTENT(IN) :: pstep ! Parflow model time-step in hours REAL(KIND=8), INTENT(IN) :: topo((nx+2)*(ny+2)*(nz+2)) ! mask (0 for inactive, 1 for active) REAL(KIND=8), INTENT(INOUT) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! source/sink (1/T) - +REAL(KIND=8), INTENT(INOUT) :: ice_fraction((nx+2)*(ny+2)*(nz+2)) ! soil ice fraction ! All vecotrs from parflow on grid w/ ghost nodes for current proc !Local Variables @@ -102,6 +102,7 @@ SUBROUTINE receive_fld2_clm(evap_trans,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) ENDDO ! evap_trans = 0. !CPS initialize for masking + ice_fraction = 1. DO i = 1, nx DO j = 1, ny DO k = 1, nlevsoil diff --git a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c index 83b8033e8..b9cb61065 100644 --- a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c +++ b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c @@ -197,7 +197,7 @@ void PrintFinalStats( * KinsolNonlinSolver *--------------------------------------------------------------------------*/ -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity) +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction) { PFModule *this_module = ThisPFModule; PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); @@ -246,6 +246,7 @@ int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, V StatePrecond(current_state) = precond; StateEvapTrans(current_state) = evap_trans; /*sk*/ StateOvrlBcFlx(current_state) = ovrl_bc_flx; /*sk*/ + StateIceFraction(current_state) = ice_fraction; StateXvel(current_state) = x_velocity; //jjb StateYvel(current_state) = y_velocity; //jjb StateZvel(current_state) = z_velocity; //jjb diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index bf86c8c8b..e866bd481 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -84,8 +84,9 @@ void KINSolFunctionEval( Vector *old_density = StateOldDensity(((State*)current_state)); double dt = StateDt(((State*)current_state)); double time = StateTime(((State*)current_state)); - Vector *evap_trans = StateEvapTrans(((State*)current_state)); - Vector *ovrl_bc_flx = StateOvrlBcFlx(((State*)current_state)); + Vector *evap_trans = StateEvapTrans(((State*)current_state)); + Vector *ovrl_bc_flx = StateOvrlBcFlx(((State*)current_state)); + Vector *ice_fraction = StateIceFraction(((State*)current_state)); /* velocity vectors jjb */ Vector *x_velocity = StateXvel(((State*)current_state)); @@ -97,7 +98,7 @@ void KINSolFunctionEval( PFModuleInvokeType(NlFunctionEvalInvoke, nl_function_eval, (pressure, fval, problem_data, saturation, old_saturation, density, old_density, dt, time, old_pressure, evap_trans, - ovrl_bc_flx, x_velocity, y_velocity, z_velocity)); + ovrl_bc_flx, ice_fraction, x_velocity, y_velocity, z_velocity)); return; } @@ -117,9 +118,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double dt, /* Time step size */ double time, /* New time value */ Vector * old_pressure, - Vector * evap_trans, /*sk sink term from land surface model*/ - Vector * ovrl_bc_flx, /*sk overland flow boundary fluxes*/ - Vector * x_velocity, /* velocity vectors jjb */ + Vector * evap_trans, /*sk sink term from land surface model*/ + Vector * ovrl_bc_flx, /*sk overland flow boundary fluxes*/ + Vector * ice_fraction, /* soil ice fraction */ + Vector * x_velocity, /* velocity vectors jjb */ Vector * y_velocity, Vector * z_velocity) { @@ -538,7 +540,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ * phase source values */ PFModuleInvokeType(PhaseRelPermInvoke, rel_perm_module, - (rel_perm, pressure, density, gravity, problem_data, + (rel_perm, pressure, density, ice_fraction, gravity, problem_data, CALCFCN)); @@ -562,6 +564,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ permx_sub = VectorSubvector(permeability_x, is); permy_sub = VectorSubvector(permeability_y, is); permz_sub = VectorSubvector(permeability_z, is); + /* @RMM added to provide access to x/y slopes */ x_ssl_sub = VectorSubvector(x_ssl, is); y_ssl_sub = VectorSubvector(y_ssl, is); diff --git a/pfsimulator/parflow_lib/nl_function_eval.h b/pfsimulator/parflow_lib/nl_function_eval.h index 88210329a..b72f9f2e6 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.h +++ b/pfsimulator/parflow_lib/nl_function_eval.h @@ -51,6 +51,7 @@ typedef struct { double *outflow; /*sk*/ Vector *evap_trans; /*sk*/ + Vector *ice_fraction; Vector *ovrl_bc_flx; /*sk*/ Vector *x_velocity; //jjb @@ -83,6 +84,7 @@ typedef struct { #define StateOutflow(state) ((state)->outflow) /*sk*/ #define StateEvapTrans(state) ((state)->evap_trans) /*sk*/ #define StateOvrlBcFlx(state) ((state)->ovrl_bc_flx) /*sk*/ +#define StateIceFraction(state) ((state)->ice_fraction) #define StateXvel(state) ((state)->x_velocity) //jjb #define StateYvel(state) ((state)->y_velocity) //jjb #define StateZvel(state) ((state)->z_velocity) //jjb diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index cddf5b5b3..3c598a7e7 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -293,14 +293,14 @@ PFModule *InputRFNewPublicXtra(char *geom_name); void InputRFFreePublicXtra(void); int InputRFSizeOfTempData(void); -typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction); typedef PFModule *(*NonlinSolverInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); /* kinsol_nonlin_solver.c */ int KINSolInitPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtemp1, N_Vector vtemp2, void *nl_function, double uround, long int *nfePtr, void *current_state); int KINSolCallPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtem, N_Vector ftem, void *nl_function, double uround, long int *nfePtr, void *current_state); void PrintFinalStats(FILE *out_file, long int *integer_outputs_now, long int *integer_outputs_total); -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction); PFModule *KinsolNonlinSolverInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); void KinsolNonlinSolverFreeInstanceXtra(void); PFModule *KinsolNonlinSolverNewPublicXtra(void); @@ -433,12 +433,12 @@ void FreeTempVector(Vector *vector); /* new_endpts.c */ void NewEndpts(double *alpha, double *beta, double *pp, int *size_ptr, int n, double *a_ptr, double *b_ptr, double *cond_ptr, double ereps); -typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NlFunctionEvalInitInstanceXtraInvoke) (Problem *problem, Grid *grid, double *temp_data); /* nl_function_eval.c */ void KINSolFunctionEval(int size, N_Vector pressure, N_Vector fval, void *current_state); -void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *NlFunctionEvalInitInstanceXtra(Problem *problem, Grid *grid, double *temp_data); void NlFunctionEvalFreeInstanceXtra(void); PFModule *NlFunctionEvalNewPublicXtra(char *name); @@ -909,11 +909,11 @@ PFModule *PhaseMobilityNewPublicXtra(int num_phases); void PhaseMobilityFreePublicXtra(void); int PhaseMobilitySizeOfTempData(void); -typedef void (*PhaseRelPermInvoke) (Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, double gravity, ProblemData *problem_data, int fcn); +typedef void (*PhaseRelPermInvoke) (Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_fraction, double gravity, ProblemData *problem_data, int fcn); typedef PFModule *(*PhaseRelPermInitInstanceXtraInvoke) (Grid *grid, double *temp_data); /* problem_phase_rel_perm.c */ -void PhaseRelPerm(Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, double gravity, ProblemData *problem_data, int fcn); +void PhaseRelPerm(Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_fraction, double gravity, ProblemData *problem_data, int fcn); PFModule *PhaseRelPermInitInstanceXtra(Grid *grid, double *temp_data); void PhaseRelPermFreeInstanceXtra(void); PFModule *PhaseRelPermNewPublicXtra(void); diff --git a/pfsimulator/parflow_lib/problem_phase_rel_perm.c b/pfsimulator/parflow_lib/problem_phase_rel_perm.c index bfbd696a5..f492642a2 100644 --- a/pfsimulator/parflow_lib/problem_phase_rel_perm.c +++ b/pfsimulator/parflow_lib/problem_phase_rel_perm.c @@ -406,7 +406,8 @@ void PhaseRelPerm( Vector * phase_rel_perm, /* Vector of return rel. perms. at each block */ Vector * phase_pressure, /* Vector of pressures at each block */ - Vector * phase_density, /* Vector of densities at each block */ + Vector * phase_density, /* Vector of densities at each block */ + Vector * ice_fraction, /* Vector of soil ice fraction at each block */ double gravity, /* Magnitude of gravity in neg. z direction */ ProblemData *problem_data, /* Contains geometry info for the problem */ int fcn) /* Flag determining what to calculate diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index 9db749277..f709ffe12 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -507,11 +507,11 @@ void RichardsJacobianEval( /* Calculate rel_perm and rel_perm_der */ PFModuleInvokeType(PhaseRelPermInvoke, rel_perm_module, - (rel_perm, pressure, density, gravity, problem_data, + (rel_perm, pressure, density, NULL, gravity, problem_data, CALCFCN)); PFModuleInvokeType(PhaseRelPermInvoke, rel_perm_module, - (rel_perm_der, pressure, density, gravity, problem_data, + (rel_perm_der, pressure, density, NULL, gravity, problem_data, CALCDER)); /* Calculate contributions from second order derivatives and gravity */ diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 3c9c86f9b..d2cf3e3e3 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -248,6 +248,7 @@ typedef struct { Vector *x_velocity; /* vector containing x-velocity face values */ Vector *y_velocity; /* vector containing y-velocity face values */ Vector *z_velocity; /* vector containing z-velocity face values */ + Vector *ice_fraction; /* soil ice fraction from eCLM */ #ifdef HAVE_CLM /* RM: vars for pf printing of clm output */ Vector *eflx_lh_tot; /* total LH flux from canopy height to atmosphere [W/m^2] */ @@ -1571,8 +1572,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time #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; + Subvector *p_sub, *s_sub, *et_sub, *ice_frac_sub, *m_sub, *po_sub, *dz_sub; + double *pp, *sp, *et, *ice_frac, *ms, *po_dat, *dz_dat; double sw_lat = .0; double sw_lon = .0; #endif @@ -1759,6 +1760,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time p_sub = VectorSubvector(instance_xtra->pressure, is); s_sub = VectorSubvector(instance_xtra->saturation, is); et_sub = VectorSubvector(evap_trans, is); + ice_frac_sub = VectorSubvector(instance_xtra->ice_fraction, is); m_sub = VectorSubvector(instance_xtra->mask, is); po_sub = VectorSubvector(porosity, is); dz_sub = VectorSubvector(instance_xtra->dz_mult, is); @@ -1775,6 +1777,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time sp = SubvectorData(s_sub); pp = SubvectorData(p_sub); et = SubvectorData(et_sub); + ice_frac = SubvectorData(et_sub); ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); @@ -1782,11 +1785,13 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t,po_dat,dz_dat); amps_Sync(amps_CommWorld); - CALL_receive_fld2_clm(et, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); + CALL_receive_fld2_clm(et, ice_frac, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); } amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); + handle = InitVectorUpdate(instance_xtra->ice_fraction, VectorUpdateAll); + FinalizeVectorUpdate(handle); #endif // end to HAVE_OAS3 CALL // IMF: Added to include CLM dumps in file_number updating. @@ -2880,7 +2885,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time instance_xtra->ovrl_bc_flx, instance_xtra->x_velocity, instance_xtra->y_velocity, - instance_xtra->z_velocity)); + instance_xtra->z_velocity, + instance_xtra->ice_fraction)); if (retval != 0) { From 5e9ddec249be3b4e7ddebb4f27e567d35068791e Mon Sep 17 00:00:00 2001 From: kvrigor Date: Thu, 6 Apr 2023 13:32:47 +0200 Subject: [PATCH 02/19] Initialized ice_fraction vector to 1.0 --- pfsimulator/parflow_lib/kinsol_nonlin_solver.c | 2 +- pfsimulator/parflow_lib/parflow_proto.h | 4 ++-- pfsimulator/parflow_lib/solver_richards.c | 11 +++++++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c index b9cb61065..6cb850749 100644 --- a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c +++ b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c @@ -197,7 +197,7 @@ void PrintFinalStats( * KinsolNonlinSolver *--------------------------------------------------------------------------*/ -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction) +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ice_fraction, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity) { PFModule *this_module = ThisPFModule; PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index 3c598a7e7..48c9eda1f 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -293,14 +293,14 @@ PFModule *InputRFNewPublicXtra(char *geom_name); void InputRFFreePublicXtra(void); int InputRFSizeOfTempData(void); -typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction); +typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NonlinSolverInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); /* kinsol_nonlin_solver.c */ int KINSolInitPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtemp1, N_Vector vtemp2, void *nl_function, double uround, long int *nfePtr, void *current_state); int KINSolCallPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtem, N_Vector ftem, void *nl_function, double uround, long int *nfePtr, void *current_state); void PrintFinalStats(FILE *out_file, long int *integer_outputs_now, long int *integer_outputs_total); -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *ice_fraction); +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *KinsolNonlinSolverInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); void KinsolNonlinSolverFreeInstanceXtra(void); PFModule *KinsolNonlinSolverNewPublicXtra(void); diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index d2cf3e3e3..d5c204445 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -248,7 +248,7 @@ typedef struct { Vector *x_velocity; /* vector containing x-velocity face values */ Vector *y_velocity; /* vector containing y-velocity face values */ Vector *z_velocity; /* vector containing z-velocity face values */ - Vector *ice_fraction; /* soil ice fraction from eCLM */ + Vector *ice_fraction; /* soil ice fraction from eCLM */ #ifdef HAVE_CLM /* RM: vars for pf printing of clm output */ Vector *eflx_lh_tot; /* total LH flux from canopy height to atmosphere [W/m^2] */ @@ -850,6 +850,9 @@ SetupRichards(PFModule * this_module) instance_xtra->evap_trans = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->evap_trans, 0.0); + instance_xtra->ice_fraction = NewVectorType(grid, 1, 1, vector_cell_centered); + InitVectorAll(instance_xtra->ice_fraction, 1.0); + if (public_xtra->evap_trans_file) { //sprintf(filename, "%s", public_xtra->evap_trans_filename); @@ -1777,7 +1780,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time sp = SubvectorData(s_sub); pp = SubvectorData(p_sub); et = SubvectorData(et_sub); - ice_frac = SubvectorData(et_sub); + ice_frac = SubvectorData(ice_frac_sub); ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); @@ -2883,10 +2886,10 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time instance_xtra->old_pressure, evap_trans, instance_xtra->ovrl_bc_flx, + instance_xtra->ice_fraction, instance_xtra->x_velocity, instance_xtra->y_velocity, - instance_xtra->z_velocity, - instance_xtra->ice_fraction)); + instance_xtra->z_velocity)); if (retval != 0) { From 7f9570e05afab5a8bb9fe6b2b15a24c541d6817a Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sun, 14 May 2023 17:04:16 +0200 Subject: [PATCH 03/19] Added cpl field PFL_ICE_FRAC and added option to log ice fraction to NetCDF output --- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 6 +++++- pfsimulator/parflow_lib/solver_richards.c | 22 +++++++++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 938cac7fb..fd69fdb53 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -30,7 +30,7 @@ module oas_pfl_mod ! Local variables used by OASIS integer :: comp_id, rank, ierror - integer :: soilliq_id, psi_id, et_id ! coupling field IDs + integer :: soilliq_id, psi_id, et_id, ice_frac_id ! coupling field IDs ! TODO: Get these values from eCLM instead of hardcoding it here integer :: nlevsoi = 20 ! Number of hydrologically active soil layers in eCLM @@ -111,6 +111,9 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, call oasis_def_var (et_id, "PFL_ET", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') + call oasis_def_var (ice_frac_id, "PFL_ICE_FRAC", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) + if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') + var_nodims(2) = nlevgrnd call oasis_def_var (psi_id, "PFL_PSI", part_id, var_nodims, OASIS_Out, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_PSI') @@ -206,6 +209,7 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) + call oasis_get(ice_frac_id, seconds_elapsed, ice_frac_3d, ierror) ! Save ET fluxes to ParFlow evap_trans vector evap_trans = 0. diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index d5c204445..80c752884 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -186,6 +186,7 @@ typedef struct { int write_netcdf_satur; /* write saturations? */ int write_netcdf_evaptrans; /* write evaptrans? */ int write_netcdf_evaptrans_sum; /* write evaptrans_sum? */ + int write_netcdf_ice_frac; /* write soil ice fraction? */ int write_netcdf_overland_sum; /* write overland_sum? */ int write_netcdf_overland_bc_flux; /* write overland_bc_flux? */ int write_netcdf_mask; /* write mask? */ @@ -851,7 +852,7 @@ SetupRichards(PFModule * this_module) InitVectorAll(instance_xtra->evap_trans, 0.0); instance_xtra->ice_fraction = NewVectorType(grid, 1, 1, vector_cell_centered); - InitVectorAll(instance_xtra->ice_fraction, 1.0); + InitVectorAll(instance_xtra->ice_fraction, 0.0); if (public_xtra->evap_trans_file) { @@ -3097,6 +3098,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time || public_xtra->write_netcdf_satur || public_xtra->write_netcdf_evaptrans || public_xtra->write_netcdf_evaptrans_sum + || public_xtra->write_netcdf_ice_frac || public_xtra->write_netcdf_overland_sum || public_xtra->write_netcdf_overland_bc_flux) { @@ -3271,6 +3273,14 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time any_file_dumped = 1; } + if (public_xtra->write_netcdf_ice_frac) + { + sprintf(nc_postfix, "%05d", instance_xtra->file_number); + WritePFNC(file_prefix, nc_postfix, t, instance_xtra->ice_fraction, + public_xtra->numVarTimeVariant, "ice_fraction", 3, + false, public_xtra->numVarIni); + any_file_dumped = 1; + } if (public_xtra->print_evaptrans_sum || public_xtra->write_silo_evaptrans_sum @@ -5371,6 +5381,15 @@ SolverRichardsNewPublicXtra(char *name) } public_xtra->write_netcdf_evaptrans = switch_value; + sprintf(key, "NetCDF.WriteIceFrac"); + switch_name = GetStringDefault(key, "False"); + switch_value = NA_NameToIndexExitOnError(switch_na, switch_name, key); + if (switch_value == 1) + { + public_xtra->numVarTimeVariant++; + } + public_xtra->write_netcdf_ice_frac = switch_value; + sprintf(key, "NetCDF.WriteEvapTransSum"); switch_name = GetStringDefault(key, "False"); switch_value = NA_NameToIndexExitOnError(switch_na, switch_name, key); @@ -5460,6 +5479,7 @@ SolverRichardsNewPublicXtra(char *name) if (public_xtra->write_netcdf_press || public_xtra->write_netcdf_satur || public_xtra->write_netcdf_evaptrans || public_xtra->write_netcdf_evaptrans_sum + || public_xtra->write_netcdf_ice_frac || public_xtra->write_netcdf_overland_sum || public_xtra->write_netcdf_overland_bc_flux) From ab036f33089b46b157535c3c30e5f70a0ffe4512 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 13 Jun 2023 17:32:27 +0200 Subject: [PATCH 04/19] Fixed incorrect dimensions of ice_fraction --- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index fd69fdb53..cf4cff560 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -111,14 +111,13 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, call oasis_def_var (et_id, "PFL_ET", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') - call oasis_def_var (ice_frac_id, "PFL_ICE_FRAC", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) - if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') - var_nodims(2) = nlevgrnd call oasis_def_var (psi_id, "PFL_PSI", part_id, var_nodims, OASIS_Out, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_PSI') call oasis_def_var (soilliq_id, "PFL_SOILLIQ", part_id, var_nodims, OASIS_Out, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_SOILLIQ') + call oasis_def_var (ice_frac_id, "PFL_ICE_FRAC", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) + if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') call oasis_enddef ( ierror ) if (ierror /= 0) call oasis_abort (comp_id, 'oas_pfl_define', 'oasis_enddef failed') @@ -190,9 +189,8 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! evapotranspiration [1/hrs] real(kind=8), intent(inout) :: ice_fraction((nx+2)*(ny+2)*(nz+2)) ! soil ice fraction - - ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the - ! extra "+2" terms account for the ghost nodes/halo points + ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the + ! extra "+2" terms account for the ghost nodes/halo points ! Local variables integer :: seconds_elapsed ! current model time in seconds @@ -205,7 +203,7 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, ! Receive ET fluxes from eCLM allocate(evap_trans_3d(nx,ny,nlevsoi)) - allocate(ice_frac_3d(nx,ny,nlevsoi)) + allocate(ice_frac_3d(nx,ny,nlevgrnd)) seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) @@ -217,10 +215,10 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, do i = 1, nx do j = 1, ny if (top_z_level(i,j) > 0) then ! ***** Subsurface level indexing convention ***** - do k = 1, nlevsoi ! eCLM: 1=topmost layer, nlevsoi=deepest layer + do k = 1, nlevgrnd ! eCLM: 1=topmost layer, nlevsoi=deepest layer z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) - evap_trans(l) = evap_trans_3d(i,j,k) + if (k <= nlevsoi) evap_trans(l) = evap_trans_3d(i,j,k) ice_fraction(l) = ice_frac_3d(i,j,k) end do end if From 5a1f5e968cb997d5a807fcb74364e53952dae288 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sun, 18 Jun 2023 21:20:59 +0200 Subject: [PATCH 05/19] Added NetCDF key for ice_fraction --- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 3 +- pfsimulator/parflow_lib/solver_richards.c | 1 + .../parflow_lib/write_parflow_netcdf.c | 42 +++++++++++++++++++ 3 files changed, 45 insertions(+), 1 deletion(-) diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index cf4cff560..9e3170eb0 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -211,6 +211,7 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, ! Save ET fluxes to ParFlow evap_trans vector evap_trans = 0. + ice_fraction = 0. top_z_level = get_top_z_level(nx, ny, nz, topo) do i = 1, nx do j = 1, ny @@ -218,7 +219,7 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, do k = 1, nlevgrnd ! eCLM: 1=topmost layer, nlevsoi=deepest layer z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) - if (k <= nlevsoi) evap_trans(l) = evap_trans_3d(i,j,k) + if (k <= nlevsoi) evap_trans(l) = evap_trans_3d(i,j,k) ice_fraction(l) = ice_frac_3d(i,j,k) end do end if diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 80c752884..43dcffd44 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -4095,6 +4095,7 @@ TeardownRichards(PFModule * this_module) FreeVector(instance_xtra->y_velocity); FreeVector(instance_xtra->z_velocity); FreeVector(instance_xtra->evap_trans); + FreeVector(instance_xtra->ice_fraction); if (instance_xtra->evap_trans_sum) { diff --git a/pfsimulator/parflow_lib/write_parflow_netcdf.c b/pfsimulator/parflow_lib/write_parflow_netcdf.c index b9707a193..ade43f07a 100644 --- a/pfsimulator/parflow_lib/write_parflow_netcdf.c +++ b/pfsimulator/parflow_lib/write_parflow_netcdf.c @@ -948,6 +948,48 @@ int LookUpInventory(char * varName, varNCData **myVarNCData, int *netCDFIDs) return evaptransVarID; } + if (strcmp(varName, "ice_fraction") == 0) + { + *myVarNCData = malloc(sizeof(varNCData)); + (*myVarNCData)->varName = varName; + (*myVarNCData)->ncType = NC_DOUBLE; + (*myVarNCData)->dimSize = 4; + (*myVarNCData)->dimIDs = malloc((*myVarNCData)->dimSize * sizeof(int)); + (*myVarNCData)->dimIDs[0] = netCDFIDs[1]; + (*myVarNCData)->dimIDs[1] = netCDFIDs[2]; + (*myVarNCData)->dimIDs[2] = netCDFIDs[3]; + (*myVarNCData)->dimIDs[3] = netCDFIDs[4]; + int iceFractionVarID; + int res = nc_def_var(netCDFIDs[0], varName, (*myVarNCData)->ncType, (*myVarNCData)->dimSize, + (*myVarNCData)->dimIDs, &iceFractionVarID); + if (res != NC_ENAMEINUSE) + { + char *switch_name; + char key[IDB_MAX_KEY_LEN]; + char *default_val = "None"; + sprintf(key, "NetCDF.Chunking"); + switch_name = GetStringDefault(key, "None"); + if (strcmp(switch_name, default_val) != 0) + { + size_t chunksize[(*myVarNCData)->dimSize]; + chunksize[0] = 1; + chunksize[1] = GetInt("NetCDF.ChunkZ"); + chunksize[2] = GetInt("NetCDF.ChunkY"); + chunksize[3] = GetInt("NetCDF.ChunkX"); + nc_def_var_chunking(netCDFIDs[0], iceFractionVarID, NC_CHUNKED, chunksize); + } + if (enable_netcdf_compression) { + nc_def_var_deflate(netCDFIDs[0],iceFractionVarID,0,1,compression_level); + } + + } + if (res == NC_ENAMEINUSE) + { + res = nc_inq_varid(netCDFIDs[0], varName, &iceFractionVarID); + } + return iceFractionVarID; + } + if (strcmp(varName, "evaptrans_sum") == 0) { *myVarNCData = malloc(sizeof(varNCData)); From f4539d0ed3994eebbd03f4c354127591ffb9e8aa Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 29 Aug 2023 10:39:40 +0200 Subject: [PATCH 06/19] Changed received field to ice impedance --- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 20 ++--- pfsimulator/parflow_lib/solver_richards.c | 85 ++++++++++++++----- .../parflow_lib/write_parflow_netcdf.c | 2 +- 3 files changed, 74 insertions(+), 33 deletions(-) diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 9e3170eb0..7cb9b34fb 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -30,7 +30,7 @@ module oas_pfl_mod ! Local variables used by OASIS integer :: comp_id, rank, ierror - integer :: soilliq_id, psi_id, et_id, ice_frac_id ! coupling field IDs + integer :: soilliq_id, psi_id, et_id, ice_impedance_id ! coupling field IDs ! TODO: Get these values from eCLM instead of hardcoding it here integer :: nlevsoi = 20 ! Number of hydrologically active soil layers in eCLM @@ -116,7 +116,7 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_PSI') call oasis_def_var (soilliq_id, "PFL_SOILLIQ", part_id, var_nodims, OASIS_Out, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_SOILLIQ') - call oasis_def_var (ice_frac_id, "PFL_ICE_FRAC", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) + call oasis_def_var (ice_impedance_id, "PFL_ICE_IMPEDANCE", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ET') call oasis_enddef ( ierror ) @@ -178,7 +178,7 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n deallocate(pressure_3d) end subroutine send_fld2_clm - subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') + subroutine receive_fld2_clm(evap_trans, ice_impedance, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') !----------------------------------------------- ! Receives evapotranspiration fluxes from eCLM. !----------------------------------------------- @@ -188,7 +188,7 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, real(kind=8), intent(in) :: pstep, & ! current model time in hours topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! evapotranspiration [1/hrs] - real(kind=8), intent(inout) :: ice_fraction((nx+2)*(ny+2)*(nz+2)) ! soil ice fraction + real(kind=8), intent(inout) :: ice_impedance((nx+2)*(ny+2)*(nz+2))! soil ice impedance ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the ! extra "+2" terms account for the ghost nodes/halo points @@ -199,19 +199,19 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, integer :: z ! subsurface level (z=nz topmost layer, z=1 deepest layer) integer :: top_z_level(nx,ny) ! topmost z level of active ParFlow cells real(kind=8), allocatable :: evap_trans_3d(:,:,:) ! Root ET fluxes received from eCLM [1/hrs] - real(kind=8), allocatable :: ice_frac_3d(:,:,:) ! Soil ice fraction received from eCLM [unitless; values from 0 to 1] + real(kind=8), allocatable :: ice_impedance_3d(:,:,:) ! Soil ice impedance received from eCLM [unitless; values from 0 to 1] ! Receive ET fluxes from eCLM allocate(evap_trans_3d(nx,ny,nlevsoi)) - allocate(ice_frac_3d(nx,ny,nlevgrnd)) + allocate(ice_impedance_3d(nx,ny,nlevgrnd)) seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) - call oasis_get(ice_frac_id, seconds_elapsed, ice_frac_3d, ierror) + call oasis_get(ice_impedance_id, seconds_elapsed, ice_impedance_3d, ierror) ! Save ET fluxes to ParFlow evap_trans vector evap_trans = 0. - ice_fraction = 0. + ice_impedance = 0. top_z_level = get_top_z_level(nx, ny, nz, topo) do i = 1, nx do j = 1, ny @@ -220,13 +220,13 @@ subroutine receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) if (k <= nlevsoi) evap_trans(l) = evap_trans_3d(i,j,k) - ice_fraction(l) = ice_frac_3d(i,j,k) + ice_impedance(l) = ice_impedance_3d(i,j,k) end do end if end do end do - deallocate(evap_trans_3d, ice_frac_3d) + deallocate(evap_trans_3d, ice_impedance_3d) end subroutine receive_fld2_clm function get_top_z_level(nx, ny, nz, topo) diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 43dcffd44..bdc0643c8 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -185,10 +185,10 @@ typedef struct { int write_netcdf_press; /* write pressures? */ int write_netcdf_satur; /* write saturations? */ int write_netcdf_evaptrans; /* write evaptrans? */ - int write_netcdf_evaptrans_sum; /* write evaptrans_sum? */ - int write_netcdf_ice_frac; /* write soil ice fraction? */ - int write_netcdf_overland_sum; /* write overland_sum? */ - int write_netcdf_overland_bc_flux; /* write overland_bc_flux? */ + int write_netcdf_evaptrans_sum; /* write evaptrans_sum? */ + int write_netcdf_ice_impedance; /* write soil ice impedance? */ + int write_netcdf_overland_sum; /* write overland_sum? */ + int write_netcdf_overland_bc_flux; /* write overland_bc_flux? */ int write_netcdf_mask; /* write mask? */ int write_netcdf_mannings; /* write mask? */ int write_netcdf_subsurface; /* write subsurface? */ @@ -249,7 +249,7 @@ typedef struct { Vector *x_velocity; /* vector containing x-velocity face values */ Vector *y_velocity; /* vector containing y-velocity face values */ Vector *z_velocity; /* vector containing z-velocity face values */ - Vector *ice_fraction; /* soil ice fraction from eCLM */ + Vector *ice_impedance; /* soil ice impedance from eCLM */ #ifdef HAVE_CLM /* RM: vars for pf printing of clm output */ Vector *eflx_lh_tot; /* total LH flux from canopy height to atmosphere [W/m^2] */ @@ -851,8 +851,8 @@ SetupRichards(PFModule * this_module) instance_xtra->evap_trans = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->evap_trans, 0.0); - instance_xtra->ice_fraction = NewVectorType(grid, 1, 1, vector_cell_centered); - InitVectorAll(instance_xtra->ice_fraction, 0.0); + instance_xtra->ice_impedance = NewVectorType(grid, 1, 1, vector_cell_centered); + InitVectorAll(instance_xtra->ice_impedance, 0.0); if (public_xtra->evap_trans_file) { @@ -1377,6 +1377,27 @@ SetupRichards(PFModule * this_module) any_file_dumped = 1; } + /*----------------------------------------------------------------- + * Print out Subsurface data in NetCDF? + *-----------------------------------------------------------------*/ + // if (public_xtra->write_netcdf_ice_impedance) + // { + // sprintf(nc_postfix, "%05d", instance_xtra->file_number); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBx(problem_data), + // public_xtra->numVarTimeVariant, "FBx", 3, true, + // public_xtra->numVarIni); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBy(problem_data), + // public_xtra->numVarTimeVariant, "FBy", 3, true, + // public_xtra->numVarIni); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBz(problem_data), + // public_xtra->numVarTimeVariant, "FBz", 3, true, + // public_xtra->numVarIni); + // any_file_dumped = 1; + // } + /*----------------------------------------------------------------- * Print out Slopes in NetCDF? *-----------------------------------------------------------------*/ @@ -1576,8 +1597,11 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time #ifdef HAVE_OAS3 Grid *grid = (instance_xtra->grid); Subgrid *subgrid; - Subvector *p_sub, *s_sub, *et_sub, *ice_frac_sub, *m_sub, *po_sub, *dz_sub; - double *pp, *sp, *et, *ice_frac, *ms, *po_dat, *dz_dat; + Subvector *p_sub, *s_sub, *et_sub, *ice_impedance_sub, *m_sub, *po_sub, *dz_sub; + Vector *FBx = ProblemDataFBx(problem_data); + Vector *FBy = ProblemDataFBy(problem_data); + Vector *FBz = ProblemDataFBz(problem_data); + double *pp, *sp, *et, *ice_impedance, *ms, *po_dat, *dz_dat; double sw_lat = .0; double sw_lon = .0; #endif @@ -1764,7 +1788,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time p_sub = VectorSubvector(instance_xtra->pressure, is); s_sub = VectorSubvector(instance_xtra->saturation, is); et_sub = VectorSubvector(evap_trans, is); - ice_frac_sub = VectorSubvector(instance_xtra->ice_fraction, is); + ice_impedance_sub = VectorSubvector(instance_xtra->ice_impedance, is); + //ice_impedance_sub = VectorSubvector(FBx, is); m_sub = VectorSubvector(instance_xtra->mask, is); po_sub = VectorSubvector(porosity, is); dz_sub = VectorSubvector(instance_xtra->dz_mult, is); @@ -1781,7 +1806,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time sp = SubvectorData(s_sub); pp = SubvectorData(p_sub); et = SubvectorData(et_sub); - ice_frac = SubvectorData(ice_frac_sub); + ice_impedance = SubvectorData(ice_impedance_sub); ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); @@ -1789,13 +1814,17 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t,po_dat,dz_dat); amps_Sync(amps_CommWorld); - CALL_receive_fld2_clm(et, ice_frac, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); + CALL_receive_fld2_clm(et, ice_impedance, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); } amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); - handle = InitVectorUpdate(instance_xtra->ice_fraction, VectorUpdateAll); + handle = InitVectorUpdate(instance_xtra->ice_impedance, VectorUpdateAll); + //handle = InitVectorUpdate(FBx, VectorUpdateAll); FinalizeVectorUpdate(handle); + //PFVCopy(FBx, FBy); + //PFVCopy(FBx, FBz); + #endif // end to HAVE_OAS3 CALL // IMF: Added to include CLM dumps in file_number updating. @@ -2887,7 +2916,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time instance_xtra->old_pressure, evap_trans, instance_xtra->ovrl_bc_flx, - instance_xtra->ice_fraction, + instance_xtra->ice_impedance, instance_xtra->x_velocity, instance_xtra->y_velocity, instance_xtra->z_velocity)); @@ -3098,7 +3127,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time || public_xtra->write_netcdf_satur || public_xtra->write_netcdf_evaptrans || public_xtra->write_netcdf_evaptrans_sum - || public_xtra->write_netcdf_ice_frac + || public_xtra->write_netcdf_ice_impedance || public_xtra->write_netcdf_overland_sum || public_xtra->write_netcdf_overland_bc_flux) { @@ -3273,12 +3302,24 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time any_file_dumped = 1; } - if (public_xtra->write_netcdf_ice_frac) + if (public_xtra->write_netcdf_ice_impedance) { sprintf(nc_postfix, "%05d", instance_xtra->file_number); - WritePFNC(file_prefix, nc_postfix, t, instance_xtra->ice_fraction, - public_xtra->numVarTimeVariant, "ice_fraction", 3, + WritePFNC(file_prefix, nc_postfix, t, instance_xtra->ice_impedance, + public_xtra->numVarTimeVariant, "ice_impedance", 3, false, public_xtra->numVarIni); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBx(problem_data), + // public_xtra->numVarTimeVariant, "FBx", 3, true, + // public_xtra->numVarIni); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBy(problem_data), + // public_xtra->numVarTimeVariant, "FBy", 3, true, + // public_xtra->numVarIni); + // WritePFNC(file_prefix, nc_postfix, t, + // ProblemDataFBz(problem_data), + // public_xtra->numVarTimeVariant, "FBz", 3, true, + // public_xtra->numVarIni); any_file_dumped = 1; } @@ -4095,7 +4136,7 @@ TeardownRichards(PFModule * this_module) FreeVector(instance_xtra->y_velocity); FreeVector(instance_xtra->z_velocity); FreeVector(instance_xtra->evap_trans); - FreeVector(instance_xtra->ice_fraction); + FreeVector(instance_xtra->ice_impedance); if (instance_xtra->evap_trans_sum) { @@ -5382,14 +5423,14 @@ SolverRichardsNewPublicXtra(char *name) } public_xtra->write_netcdf_evaptrans = switch_value; - sprintf(key, "NetCDF.WriteIceFrac"); + sprintf(key, "NetCDF.WriteIceImpedance"); switch_name = GetStringDefault(key, "False"); switch_value = NA_NameToIndexExitOnError(switch_na, switch_name, key); if (switch_value == 1) { public_xtra->numVarTimeVariant++; } - public_xtra->write_netcdf_ice_frac = switch_value; + public_xtra->write_netcdf_ice_impedance = switch_value; sprintf(key, "NetCDF.WriteEvapTransSum"); switch_name = GetStringDefault(key, "False"); @@ -5480,7 +5521,7 @@ SolverRichardsNewPublicXtra(char *name) if (public_xtra->write_netcdf_press || public_xtra->write_netcdf_satur || public_xtra->write_netcdf_evaptrans || public_xtra->write_netcdf_evaptrans_sum - || public_xtra->write_netcdf_ice_frac + || public_xtra->write_netcdf_ice_impedance || public_xtra->write_netcdf_overland_sum || public_xtra->write_netcdf_overland_bc_flux) diff --git a/pfsimulator/parflow_lib/write_parflow_netcdf.c b/pfsimulator/parflow_lib/write_parflow_netcdf.c index ade43f07a..59e9e0429 100644 --- a/pfsimulator/parflow_lib/write_parflow_netcdf.c +++ b/pfsimulator/parflow_lib/write_parflow_netcdf.c @@ -948,7 +948,7 @@ int LookUpInventory(char * varName, varNCData **myVarNCData, int *netCDFIDs) return evaptransVarID; } - if (strcmp(varName, "ice_fraction") == 0) + if (strcmp(varName, "ice_impedance") == 0) { *myVarNCData = malloc(sizeof(varNCData)); (*myVarNCData)->varName = varName; From 3c3cc36a132819d41a95cd0477fd27285e490faa Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 29 Aug 2023 12:12:05 +0200 Subject: [PATCH 07/19] Renamed ice fraction variables to ice impedance --- pfsimulator/amps/oas3/oas3_coupler.h | 6 +-- pfsimulator/amps/oas3/receive_fld2_clm.F90 | 6 +-- pfsimulator/parflow_lib/parflow_proto.h | 12 +++--- .../parflow_lib/problem_phase_rel_perm.c | 2 +- pfsimulator/parflow_lib/solver_richards.c | 40 ------------------- .../parflow_lib/write_parflow_netcdf.c | 12 +++--- 6 files changed, 19 insertions(+), 59 deletions(-) diff --git a/pfsimulator/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index aaf7e47c9..dc1bbd3ae 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -28,7 +28,7 @@ void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep, double *porosity, double *dz); #define receive_fld2_clm receive_fld2_clm_ -#define CALL_receive_fld2_clm(evap_trans, ice_fraction, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ - receive_fld2_clm(evap_trans, ice_fraction, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) -void receive_fld2_clm(double *evap_trans, double *ice_fraction, double *topo, int *ix, int*iy, +#define CALL_receive_fld2_clm(evap_trans, ice_impedance, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ + receive_fld2_clm(evap_trans, ice_impedance, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) +void receive_fld2_clm(double *evap_trans, double *ice_impedance, double *topo, int *ix, int*iy, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep); diff --git a/pfsimulator/amps/oas3/receive_fld2_clm.F90 b/pfsimulator/amps/oas3/receive_fld2_clm.F90 index 4ef3acadd..96781ccb5 100644 --- a/pfsimulator/amps/oas3/receive_fld2_clm.F90 +++ b/pfsimulator/amps/oas3/receive_fld2_clm.F90 @@ -1,4 +1,4 @@ -SUBROUTINE receive_fld2_clm(evap_trans,ice_fraction,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) +SUBROUTINE receive_fld2_clm(evap_trans,ice_impedance,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep) !---------------------------------------------------------------------------- ! @@ -46,7 +46,7 @@ SUBROUTINE receive_fld2_clm(evap_trans,ice_fraction,topo,ix,iy,nx,ny,nz,nx_f,ny_ REAL(KIND=8), INTENT(IN) :: pstep ! Parflow model time-step in hours REAL(KIND=8), INTENT(IN) :: topo((nx+2)*(ny+2)*(nz+2)) ! mask (0 for inactive, 1 for active) REAL(KIND=8), INTENT(INOUT) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! source/sink (1/T) -REAL(KIND=8), INTENT(INOUT) :: ice_fraction((nx+2)*(ny+2)*(nz+2)) ! soil ice fraction +REAL(KIND=8), INTENT(INOUT) :: ice_impedance((nx+2)*(ny+2)*(nz+2))! soil ice fraction ! All vecotrs from parflow on grid w/ ghost nodes for current proc !Local Variables @@ -102,7 +102,7 @@ SUBROUTINE receive_fld2_clm(evap_trans,ice_fraction,topo,ix,iy,nx,ny,nz,nx_f,ny_ ENDDO ! evap_trans = 0. !CPS initialize for masking - ice_fraction = 1. + ice_impedance = 1. DO i = 1, nx DO j = 1, ny DO k = 1, nlevsoil diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index 48c9eda1f..64e894676 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -293,14 +293,14 @@ PFModule *InputRFNewPublicXtra(char *geom_name); void InputRFFreePublicXtra(void); int InputRFSizeOfTempData(void); -typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NonlinSolverInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); /* kinsol_nonlin_solver.c */ int KINSolInitPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtemp1, N_Vector vtemp2, void *nl_function, double uround, long int *nfePtr, void *current_state); int KINSolCallPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtem, N_Vector ftem, void *nl_function, double uround, long int *nfePtr, void *current_state); void PrintFinalStats(FILE *out_file, long int *integer_outputs_now, long int *integer_outputs_total); -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *KinsolNonlinSolverInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); void KinsolNonlinSolverFreeInstanceXtra(void); PFModule *KinsolNonlinSolverNewPublicXtra(void); @@ -433,12 +433,12 @@ void FreeTempVector(Vector *vector); /* new_endpts.c */ void NewEndpts(double *alpha, double *beta, double *pp, int *size_ptr, int n, double *a_ptr, double *b_ptr, double *cond_ptr, double ereps); -typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NlFunctionEvalInitInstanceXtraInvoke) (Problem *problem, Grid *grid, double *temp_data); /* nl_function_eval.c */ void KINSolFunctionEval(int size, N_Vector pressure, N_Vector fval, void *current_state); -void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_fraction, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *NlFunctionEvalInitInstanceXtra(Problem *problem, Grid *grid, double *temp_data); void NlFunctionEvalFreeInstanceXtra(void); PFModule *NlFunctionEvalNewPublicXtra(char *name); @@ -909,11 +909,11 @@ PFModule *PhaseMobilityNewPublicXtra(int num_phases); void PhaseMobilityFreePublicXtra(void); int PhaseMobilitySizeOfTempData(void); -typedef void (*PhaseRelPermInvoke) (Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_fraction, double gravity, ProblemData *problem_data, int fcn); +typedef void (*PhaseRelPermInvoke) (Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_impedance, double gravity, ProblemData *problem_data, int fcn); typedef PFModule *(*PhaseRelPermInitInstanceXtraInvoke) (Grid *grid, double *temp_data); /* problem_phase_rel_perm.c */ -void PhaseRelPerm(Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_fraction, double gravity, ProblemData *problem_data, int fcn); +void PhaseRelPerm(Vector *phase_rel_perm, Vector *phase_pressure, Vector *phase_density, Vector *ice_impedance, double gravity, ProblemData *problem_data, int fcn); PFModule *PhaseRelPermInitInstanceXtra(Grid *grid, double *temp_data); void PhaseRelPermFreeInstanceXtra(void); PFModule *PhaseRelPermNewPublicXtra(void); diff --git a/pfsimulator/parflow_lib/problem_phase_rel_perm.c b/pfsimulator/parflow_lib/problem_phase_rel_perm.c index f492642a2..5716d8475 100644 --- a/pfsimulator/parflow_lib/problem_phase_rel_perm.c +++ b/pfsimulator/parflow_lib/problem_phase_rel_perm.c @@ -407,7 +407,7 @@ void PhaseRelPerm( Vector * phase_rel_perm, /* Vector of return rel. perms. at each block */ Vector * phase_pressure, /* Vector of pressures at each block */ Vector * phase_density, /* Vector of densities at each block */ - Vector * ice_fraction, /* Vector of soil ice fraction at each block */ + Vector * ice_impedance, /* Vector of soil ice impedance at each block */ double gravity, /* Magnitude of gravity in neg. z direction */ ProblemData *problem_data, /* Contains geometry info for the problem */ int fcn) /* Flag determining what to calculate diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index bdc0643c8..afe38aa8d 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -1377,27 +1377,6 @@ SetupRichards(PFModule * this_module) any_file_dumped = 1; } - /*----------------------------------------------------------------- - * Print out Subsurface data in NetCDF? - *-----------------------------------------------------------------*/ - // if (public_xtra->write_netcdf_ice_impedance) - // { - // sprintf(nc_postfix, "%05d", instance_xtra->file_number); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBx(problem_data), - // public_xtra->numVarTimeVariant, "FBx", 3, true, - // public_xtra->numVarIni); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBy(problem_data), - // public_xtra->numVarTimeVariant, "FBy", 3, true, - // public_xtra->numVarIni); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBz(problem_data), - // public_xtra->numVarTimeVariant, "FBz", 3, true, - // public_xtra->numVarIni); - // any_file_dumped = 1; - // } - /*----------------------------------------------------------------- * Print out Slopes in NetCDF? *-----------------------------------------------------------------*/ @@ -1598,9 +1577,6 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time Grid *grid = (instance_xtra->grid); Subgrid *subgrid; Subvector *p_sub, *s_sub, *et_sub, *ice_impedance_sub, *m_sub, *po_sub, *dz_sub; - Vector *FBx = ProblemDataFBx(problem_data); - Vector *FBy = ProblemDataFBy(problem_data); - Vector *FBz = ProblemDataFBz(problem_data); double *pp, *sp, *et, *ice_impedance, *ms, *po_dat, *dz_dat; double sw_lat = .0; double sw_lon = .0; @@ -1789,7 +1765,6 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time s_sub = VectorSubvector(instance_xtra->saturation, is); et_sub = VectorSubvector(evap_trans, is); ice_impedance_sub = VectorSubvector(instance_xtra->ice_impedance, is); - //ice_impedance_sub = VectorSubvector(FBx, is); m_sub = VectorSubvector(instance_xtra->mask, is); po_sub = VectorSubvector(porosity, is); dz_sub = VectorSubvector(instance_xtra->dz_mult, is); @@ -1820,10 +1795,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); handle = InitVectorUpdate(instance_xtra->ice_impedance, VectorUpdateAll); - //handle = InitVectorUpdate(FBx, VectorUpdateAll); FinalizeVectorUpdate(handle); - //PFVCopy(FBx, FBy); - //PFVCopy(FBx, FBz); #endif // end to HAVE_OAS3 CALL @@ -3308,18 +3280,6 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time WritePFNC(file_prefix, nc_postfix, t, instance_xtra->ice_impedance, public_xtra->numVarTimeVariant, "ice_impedance", 3, false, public_xtra->numVarIni); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBx(problem_data), - // public_xtra->numVarTimeVariant, "FBx", 3, true, - // public_xtra->numVarIni); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBy(problem_data), - // public_xtra->numVarTimeVariant, "FBy", 3, true, - // public_xtra->numVarIni); - // WritePFNC(file_prefix, nc_postfix, t, - // ProblemDataFBz(problem_data), - // public_xtra->numVarTimeVariant, "FBz", 3, true, - // public_xtra->numVarIni); any_file_dumped = 1; } diff --git a/pfsimulator/parflow_lib/write_parflow_netcdf.c b/pfsimulator/parflow_lib/write_parflow_netcdf.c index 59e9e0429..e90c4dc69 100644 --- a/pfsimulator/parflow_lib/write_parflow_netcdf.c +++ b/pfsimulator/parflow_lib/write_parflow_netcdf.c @@ -959,9 +959,9 @@ int LookUpInventory(char * varName, varNCData **myVarNCData, int *netCDFIDs) (*myVarNCData)->dimIDs[1] = netCDFIDs[2]; (*myVarNCData)->dimIDs[2] = netCDFIDs[3]; (*myVarNCData)->dimIDs[3] = netCDFIDs[4]; - int iceFractionVarID; + int iceImpedanceVarID; int res = nc_def_var(netCDFIDs[0], varName, (*myVarNCData)->ncType, (*myVarNCData)->dimSize, - (*myVarNCData)->dimIDs, &iceFractionVarID); + (*myVarNCData)->dimIDs, &iceImpedanceVarID); if (res != NC_ENAMEINUSE) { char *switch_name; @@ -976,18 +976,18 @@ int LookUpInventory(char * varName, varNCData **myVarNCData, int *netCDFIDs) chunksize[1] = GetInt("NetCDF.ChunkZ"); chunksize[2] = GetInt("NetCDF.ChunkY"); chunksize[3] = GetInt("NetCDF.ChunkX"); - nc_def_var_chunking(netCDFIDs[0], iceFractionVarID, NC_CHUNKED, chunksize); + nc_def_var_chunking(netCDFIDs[0], iceImpedanceVarID, NC_CHUNKED, chunksize); } if (enable_netcdf_compression) { - nc_def_var_deflate(netCDFIDs[0],iceFractionVarID,0,1,compression_level); + nc_def_var_deflate(netCDFIDs[0],iceImpedanceVarID,0,1,compression_level); } } if (res == NC_ENAMEINUSE) { - res = nc_inq_varid(netCDFIDs[0], varName, &iceFractionVarID); + res = nc_inq_varid(netCDFIDs[0], varName, &iceImpedanceVarID); } - return iceFractionVarID; + return iceImpedanceVarID; } if (strcmp(varName, "evaptrans_sum") == 0) From d4e0cab8e0e95be6fcfe77fb2829a7091c5ec121 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Tue, 29 Aug 2023 14:42:49 +0200 Subject: [PATCH 08/19] Inject ice impedance values to flow barrier vectors --- pfsimulator/parflow_lib/kinsol_nonlin_solver.c | 4 ++-- pfsimulator/parflow_lib/nl_function_eval.c | 17 ++++++++++++----- pfsimulator/parflow_lib/nl_function_eval.h | 12 ++++++------ 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c index 6cb850749..e3a7ed4f0 100644 --- a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c +++ b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c @@ -197,7 +197,7 @@ void PrintFinalStats( * KinsolNonlinSolver *--------------------------------------------------------------------------*/ -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ice_fraction, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity) +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ice_impedance, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity) { PFModule *this_module = ThisPFModule; PublicXtra *public_xtra = (PublicXtra*)PFModulePublicXtra(this_module); @@ -246,7 +246,7 @@ int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, V StatePrecond(current_state) = precond; StateEvapTrans(current_state) = evap_trans; /*sk*/ StateOvrlBcFlx(current_state) = ovrl_bc_flx; /*sk*/ - StateIceFraction(current_state) = ice_fraction; + StateIceImpedance(current_state) = ice_impedance; StateXvel(current_state) = x_velocity; //jjb StateYvel(current_state) = y_velocity; //jjb StateZvel(current_state) = z_velocity; //jjb diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index e866bd481..a0450da6b 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -86,7 +86,7 @@ void KINSolFunctionEval( double time = StateTime(((State*)current_state)); Vector *evap_trans = StateEvapTrans(((State*)current_state)); Vector *ovrl_bc_flx = StateOvrlBcFlx(((State*)current_state)); - Vector *ice_fraction = StateIceFraction(((State*)current_state)); + Vector *ice_impedance = StateIceImpedance(((State*)current_state)); /* velocity vectors jjb */ Vector *x_velocity = StateXvel(((State*)current_state)); @@ -98,7 +98,7 @@ void KINSolFunctionEval( PFModuleInvokeType(NlFunctionEvalInvoke, nl_function_eval, (pressure, fval, problem_data, saturation, old_saturation, density, old_density, dt, time, old_pressure, evap_trans, - ovrl_bc_flx, ice_fraction, x_velocity, y_velocity, z_velocity)); + ovrl_bc_flx, ice_impedance, x_velocity, y_velocity, z_velocity)); return; } @@ -120,7 +120,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Vector * old_pressure, Vector * evap_trans, /*sk sink term from land surface model*/ Vector * ovrl_bc_flx, /*sk overland flow boundary fluxes*/ - Vector * ice_fraction, /* soil ice fraction */ + Vector * ice_impedance, /* soil ice impedance */ Vector * x_velocity, /* velocity vectors jjb */ Vector * y_velocity, Vector * z_velocity) @@ -177,10 +177,17 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Subvector *z_mult_sub; //@RMM double *z_mult_dat; //@RMM -/* @RMM Flow Barrier / Boundary values */ +#ifdef HAVE_ECLM + amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); + Vector *FBx = ice_impedance; + Vector *FBy = ice_impedance; + Vector *FBz = ice_impedance; +#else + /* @RMM Flow Barrier / Boundary values */ Vector *FBx = ProblemDataFBx(problem_data); Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); +#endif Subvector *FBx_sub, *FBy_sub, *FBz_sub; //@RMM double *FBx_dat=NULL, *FBy_dat=NULL, *FBz_dat=NULL; //@RMM @@ -540,7 +547,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ * phase source values */ PFModuleInvokeType(PhaseRelPermInvoke, rel_perm_module, - (rel_perm, pressure, density, ice_fraction, gravity, problem_data, + (rel_perm, pressure, density, ice_impedance, gravity, problem_data, CALCFCN)); diff --git a/pfsimulator/parflow_lib/nl_function_eval.h b/pfsimulator/parflow_lib/nl_function_eval.h index b72f9f2e6..b792ab0e4 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.h +++ b/pfsimulator/parflow_lib/nl_function_eval.h @@ -51,7 +51,7 @@ typedef struct { double *outflow; /*sk*/ Vector *evap_trans; /*sk*/ - Vector *ice_fraction; + Vector *ice_impedance; Vector *ovrl_bc_flx; /*sk*/ Vector *x_velocity; //jjb @@ -66,7 +66,7 @@ typedef struct { *--------------------------------------------------------------------------*/ #define StateFunc(state) ((state)->nl_function_eval) -#define StateBCPressure(state) ((state)->bc_pressure)//dok +#define StateBCPressure(state) ((state)->bc_pressure)//dok #define StateProblemData(state) ((state)->problem_data) #define StateOldDensity(state) ((state)->old_density) #define StateOldPressure(state) ((state)->old_pressure) @@ -77,14 +77,14 @@ typedef struct { #define StateTime(state) ((state)->time) #define StateJacEval(state) ((state)->richards_jacobian_eval) #define StateJac(state) ((state)->jacobian_matrix) -#define StateJacC(state) ((state)->jacobian_matrix_C)//dok -#define StateJacE(state) ((state)->jacobian_matrix_E)//dok -#define StateJacF(state) ((state)->jacobian_matrix_F)//dok +#define StateJacC(state) ((state)->jacobian_matrix_C)//dok +#define StateJacE(state) ((state)->jacobian_matrix_E)//dok +#define StateJacF(state) ((state)->jacobian_matrix_F)//dok #define StatePrecond(state) ((state)->precond) #define StateOutflow(state) ((state)->outflow) /*sk*/ #define StateEvapTrans(state) ((state)->evap_trans) /*sk*/ #define StateOvrlBcFlx(state) ((state)->ovrl_bc_flx) /*sk*/ -#define StateIceFraction(state) ((state)->ice_fraction) +#define StateIceImpedance(state) ((state)->ice_impedance) #define StateXvel(state) ((state)->x_velocity) //jjb #define StateYvel(state) ((state)->y_velocity) //jjb #define StateZvel(state) ((state)->z_velocity) //jjb From f0fd821675b4d00582bd658d516f80c8d3745528 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sat, 2 Sep 2023 12:57:15 +0200 Subject: [PATCH 09/19] Added warning log in richards_jacobian_eval.c --- pfsimulator/parflow_lib/richards_jacobian_eval.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index f709ffe12..3dd447848 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -255,6 +255,14 @@ void RichardsJacobianEval( double *z_mult_dat; //@RMM /* @RMM Flow Barrier / Boundary values */ +#ifdef HAVE_ECLM + amps_Printf("DEBUG: WARNING! eCLM ice impedance not yet implemented in richard_jacobian_eval.c\n") + //TODO: RichardsJacobianEval gets called when Solver.Nonlinear.UseJacobian=True. + // Ice impedance must also be implemented here. + //Vector *FBx = ice_impedance; + //Vector *FBy = ice_impedance; + //Vector *FBz = ice_impedance; +#endif Vector *FBx = ProblemDataFBx(problem_data); Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); From 67b4b33f0dfdfc532b0a1e2cd6fe709281677fbc Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sat, 2 Sep 2023 16:48:31 +0200 Subject: [PATCH 10/19] Added HAVE_ECLM debugging logs --- cmake/parflow_config.h.in | 3 +++ pfsimulator/parflow_lib/nl_function_eval.c | 3 ++- pfsimulator/parflow_lib/richards_jacobian_eval.c | 3 ++- pfsimulator/parflow_lib/solver_richards.c | 9 +++++++++ 4 files changed, 16 insertions(+), 2 deletions(-) diff --git a/cmake/parflow_config.h.in b/cmake/parflow_config.h.in index 6149881ee..2b2efda30 100644 --- a/cmake/parflow_config.h.in +++ b/cmake/parflow_config.h.in @@ -56,6 +56,9 @@ #cmakedefine PARFLOW_HAVE_OAS3 #cmakedefine HAVE_OAS3 +#cmakedefine PARFLOW_HAVE_ECLM +#cmakedefine HAVE_ECLM + #cmakedefine PARFLOW_HAVE_MALLOC_H #cmakedefine HAVE_MALLOC_H diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index a0450da6b..431967838 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -178,7 +178,8 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double *z_mult_dat; //@RMM #ifdef HAVE_ECLM - amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); + if (!amps_Rank(amps_CommWorld)) + amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); Vector *FBx = ice_impedance; Vector *FBy = ice_impedance; Vector *FBz = ice_impedance; diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index 3dd447848..e77780a4d 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -256,7 +256,8 @@ void RichardsJacobianEval( /* @RMM Flow Barrier / Boundary values */ #ifdef HAVE_ECLM - amps_Printf("DEBUG: WARNING! eCLM ice impedance not yet implemented in richard_jacobian_eval.c\n") + if (!amps_Rank(amps_CommWorld)) + amps_Printf("DEBUG: WARNING! eCLM ice impedance not yet implemented in richard_jacobian_eval.c\n"); //TODO: RichardsJacobianEval gets called when Solver.Nonlinear.UseJacobian=True. // Ice impedance must also be implemented here. //Vector *FBx = ice_impedance; diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index afe38aa8d..5d9b11231 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -4187,6 +4187,15 @@ TeardownRichards(PFModule * this_module) int k; log_file = OpenLogFile("SolverRichards"); +#if defined(HAVE_CLM) + fprintf(log_file, "Coupled land surface model: internal CLM\n\n"); +#elif defined(HAVE_ECLM) + fprintf(log_file, "Coupled land surface model: eCLM (experimental)\n\n"); +#elif defined(HAVE_OAS3) + fprintf(log_file, "Coupled land surface model: CLM3.5\n\n"); +#else + fprintf(log_file, "Coupled land surface model: none\n\n"); +#endif if (start_count >= 0) { From e8e5d3742c830a1d2273265b01d08e86f1a24a03 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Sun, 3 Sep 2023 12:30:15 +0200 Subject: [PATCH 11/19] Passed ice impedance to RichardsJacobianEval and KinsolPCInitInstanceXtra --- .../parflow_lib/kinsol_nonlin_solver.c | 20 ++++++-- pfsimulator/parflow_lib/kinsol_pc.c | 3 +- pfsimulator/parflow_lib/parflow_proto.h | 8 +-- .../parflow_lib/richards_jacobian_eval.c | 50 +++++++++++-------- 4 files changed, 49 insertions(+), 32 deletions(-) diff --git a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c index e3a7ed4f0..d2899fc4c 100644 --- a/pfsimulator/parflow_lib/kinsol_nonlin_solver.c +++ b/pfsimulator/parflow_lib/kinsol_nonlin_solver.c @@ -107,7 +107,9 @@ int KINSolInitPC( ProblemData *problem_data = StateProblemData(((State*)current_state)); Vector *saturation = StateSaturation(((State*)current_state)); Vector *density = StateDensity(((State*)current_state)); + Vector *ice_impedance = StateIceImpedance(((State*)current_state)); Vector *old_pressure = StateOldPressure(((State*)current_state)); + double dt = StateDt(((State*)current_state)); double time = StateTime(((State*)current_state)); @@ -125,7 +127,7 @@ int KINSolInitPC( * itself */ PFModuleReNewInstanceType(KinsolPCInitInstanceXtraInvoke, precond, (NULL, NULL, problem_data, NULL, - pressure, old_pressure, saturation, density, dt, time)); + pressure, old_pressure, saturation, density, ice_impedance, dt, time)); return(0); } @@ -353,7 +355,7 @@ PFModule *KinsolNonlinSolverInitInstanceXtra( instance_xtra->precond = PFModuleNewInstanceType(KinsolPCInitInstanceXtraInvoke, public_xtra->precond, (problem, grid, problem_data, temp_data, - NULL, NULL, NULL, NULL, 0, 0)); + NULL, NULL, NULL, NULL, NULL, 0, 0)); else instance_xtra->precond = NULL; @@ -361,11 +363,14 @@ PFModule *KinsolNonlinSolverInitInstanceXtra( PFModuleNewInstanceType(NlFunctionEvalInitInstanceXtraInvoke, public_xtra->nl_function_eval, (problem, grid, temp_data)); - if (public_xtra->richards_jacobian_eval != NULL) + if (public_xtra->richards_jacobian_eval != NULL) { /* Initialize instance for nonsymmetric matrix */ + if (!amps_Rank(amps_CommWorld)) {} + amps_Printf("DEBUG: richards_jacobian_eval is not NULL!!\n"); instance_xtra->richards_jacobian_eval = PFModuleNewInstanceType(RichardsJacobianEvalInitInstanceXtraInvoke, public_xtra->richards_jacobian_eval, (problem, grid, problem_data, temp_data, 0)); + } else instance_xtra->richards_jacobian_eval = NULL; } @@ -375,7 +380,7 @@ PFModule *KinsolNonlinSolverInitInstanceXtra( PFModuleReNewInstanceType(KinsolPCInitInstanceXtraInvoke, instance_xtra->precond, (problem, grid, problem_data, temp_data, - NULL, NULL, NULL, NULL, 0, 0)); + NULL, NULL, NULL, NULL, NULL, 0, 0)); PFModuleReNewInstanceType(NlFunctionEvalInitInstanceXtraInvoke, instance_xtra->nl_function_eval, (problem, grid, temp_data)); @@ -612,6 +617,8 @@ PFModule *KinsolNonlinSolverNewPublicXtra() InputError("Invalid switch value <%s> for key <%s>", switch_name, key); } } + if (!amps_Rank(amps_CommWorld)) + amps_Printf("DEBUG: Solver.Nonlinear.UseJacobian=%i\n", switch_value); NA_FreeNameArray(switch_na); sprintf(key, "Solver.Nonlinear.DerivativeEpsilon"); @@ -671,12 +678,15 @@ PFModule *KinsolNonlinSolverNewPublicXtra() public_xtra->neq = ((public_xtra->max_restarts) + 1) * (public_xtra->krylov_dimension); - if (public_xtra->matvec != NULL) + if (public_xtra->matvec != NULL) { + if (!amps_Rank(amps_CommWorld)) + amps_Printf("DEBUG: matvec is not NULL!!\n"); public_xtra->richards_jacobian_eval = PFModuleNewModuleType( RichardsJacobianEvalNewPublicXtraInvoke, RichardsJacobianEval, ("Solver.Nonlinear.Jacobian")); + } else public_xtra->richards_jacobian_eval = NULL; diff --git a/pfsimulator/parflow_lib/kinsol_pc.c b/pfsimulator/parflow_lib/kinsol_pc.c index def622e2a..3666c3ad5 100644 --- a/pfsimulator/parflow_lib/kinsol_pc.c +++ b/pfsimulator/parflow_lib/kinsol_pc.c @@ -91,6 +91,7 @@ PFModule *KinsolPCInitInstanceXtra( Vector * old_pressure, Vector * saturation, Vector * density, + Vector * ice_impedance, double dt, double time) { @@ -146,7 +147,7 @@ PFModule *KinsolPCInitInstanceXtra( else if (pressure != NULL) { PFModuleInvokeType(RichardsJacobianEvalInvoke, (instance_xtra->discretization), - (pressure, old_pressure, &PC, &JC, saturation, density, problem_data, dt, + (pressure, old_pressure, &PC, &JC, saturation, density, ice_impedance, problem_data, dt, time, pc_matrix_type)); PFModuleReNewInstanceType(PrecondInitInstanceXtraInvoke, (instance_xtra->precond), diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index 64e894676..1e9596317 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -308,12 +308,12 @@ void KinsolNonlinSolverFreePublicXtra(void); int KinsolNonlinSolverSizeOfTempData(void); typedef void (*KinsolPCInvoke) (Vector *rhs); -typedef PFModule * (*KinsolPCInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, Vector *pressure, Vector *old_pressure, Vector *saturation, Vector *density, double dt, double time); +typedef PFModule * (*KinsolPCInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, Vector *pressure, Vector *old_pressure, Vector *saturation, Vector *density, Vector *ice_impedance, double dt, double time); typedef PFModule *(*KinsolPCNewPublicXtraInvoke) (char *name, char *pc_name); /* kinsol_pc.c */ void KinsolPC(Vector *rhs); -PFModule *KinsolPCInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, Vector *pressure, Vector *old_pressure, Vector *saturation, Vector *density, double dt, double time); +PFModule *KinsolPCInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, Vector *pressure, Vector *old_pressure, Vector *saturation, Vector *density, Vector *ice_impedance, double dt, double time); void KinsolPCFreeInstanceXtra(void); PFModule *KinsolPCNewPublicXtra(char *name, char *pc_name); void KinsolPCFreePublicXtra(void); @@ -1050,12 +1050,12 @@ void DeleteSubregion(SubregionArray *sr_array, int index); void AppendSubregionArray(SubregionArray *sr_array_0, SubregionArray *sr_array_1); -typedef void (*RichardsJacobianEvalInvoke) (Vector *pressure, Vector *old_pressure, Matrix **ptr_to_J, Matrix **ptr_to_JC, Vector *saturation, Vector *density, ProblemData *problem_data, double dt, double time, int symm_part); +typedef void (*RichardsJacobianEvalInvoke) (Vector *pressure, Vector *old_pressure, Matrix **ptr_to_J, Matrix **ptr_to_JC, Vector *saturation, Vector *density, Vector *ice_impedance, ProblemData *problem_data, double dt, double time, int symm_part); typedef PFModule *(*RichardsJacobianEvalInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, int symmetric_jac); typedef PFModule *(*RichardsJacobianEvalNewPublicXtraInvoke) (char *name); /* richards_jacobian_eval.c */ int KINSolMatVec(void *current_state, N_Vector x, N_Vector y, int *recompute, N_Vector pressure); -void RichardsJacobianEval(Vector *pressure, Vector *old_pressure, Matrix **ptr_to_J, Matrix **ptr_to_JC, Vector *saturation, Vector *density, ProblemData *problem_data, double dt, double time, int symm_part); +void RichardsJacobianEval(Vector *pressure, Vector *old_pressure, Matrix **ptr_to_J, Matrix **ptr_to_JC, Vector *saturation, Vector *density, Vector *ice_impedance, ProblemData *problem_data, double dt, double time, int symm_part); PFModule *RichardsJacobianEvalInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data, int symmetric_jac); void RichardsJacobianEvalFreeInstanceXtra(void); PFModule *RichardsJacobianEvalNewPublicXtra(char *name); diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index e77780a4d..95fed7186 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -146,6 +146,7 @@ int KINSolMatVec( Vector *old_pressure = StateOldPressure(((State*)current_state)); Vector *saturation = StateSaturation(((State*)current_state)); Vector *density = StateDensity(((State*)current_state)); + Vector *ice_impedance = StateIceImpedance(((State*)current_state)); ProblemData *problem_data = StateProblemData(((State*)current_state)); double dt = StateDt(((State*)current_state)); double time = StateTime(((State*)current_state)); @@ -164,7 +165,7 @@ int KINSolMatVec( if (*recompute) { PFModuleInvokeType(RichardsJacobianEvalInvoke, richards_jacobian_eval, - (pressure, old_pressure, &J, &JC, saturation, density, problem_data, + (pressure, old_pressure, &J, &JC, saturation, density, ice_impedance, problem_data, dt, time, 0)); *recompute = 0; @@ -185,20 +186,21 @@ int KINSolMatVec( * pressure values. */ void RichardsJacobianEval( - Vector * pressure, /* Current pressure values */ - Vector * old_pressure, /* Pressure values at previous timestep */ - Matrix ** ptr_to_J, /* Pointer to the J pointer - this will be set - * to instance_xtra pointer at end */ - Matrix ** ptr_to_JC, /* Pointer to the JC pointer - this will be set - * to instance_xtra pointer at end */ - Vector * saturation, /* Saturation / work vector */ - Vector * density, /* Density vector */ - ProblemData *problem_data, /* Geometry data for problem */ - double dt, /* Time step size */ - double time, /* New time value */ - int symm_part) /* Specifies whether to compute just the - * symmetric part of the Jacobian (1), or the - * full Jacobian */ + Vector * pressure, /* Current pressure values */ + Vector * old_pressure, /* Pressure values at previous timestep */ + Matrix ** ptr_to_J, /* Pointer to the J pointer - this will be set + * to instance_xtra pointer at end */ + Matrix ** ptr_to_JC, /* Pointer to the JC pointer - this will be set + * to instance_xtra pointer at end */ + Vector * saturation, /* Saturation / work vector */ + Vector * density, /* Density vector */ + Vector * ice_impedance, /* Ice impedance vector */ + ProblemData *problem_data, /* Geometry data for problem */ + double dt, /* Time step size */ + double time, /* New time value */ + int symm_part) /* Specifies whether to compute just the + * symmetric part of the Jacobian (1), or the + * full Jacobian */ { PUSH_NVTX("RichardsJacobianEval",1) @@ -257,16 +259,20 @@ void RichardsJacobianEval( /* @RMM Flow Barrier / Boundary values */ #ifdef HAVE_ECLM if (!amps_Rank(amps_CommWorld)) - amps_Printf("DEBUG: WARNING! eCLM ice impedance not yet implemented in richard_jacobian_eval.c\n"); - //TODO: RichardsJacobianEval gets called when Solver.Nonlinear.UseJacobian=True. - // Ice impedance must also be implemented here. - //Vector *FBx = ice_impedance; - //Vector *FBy = ice_impedance; - //Vector *FBz = ice_impedance; -#endif + { + amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); + //*(int*)0 = 0; + } + + Vector *FBx = ice_impedance; + Vector *FBy = ice_impedance; + Vector *FBz = ice_impedance; +#else Vector *FBx = ProblemDataFBx(problem_data); Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); +#endif + Subvector *FBx_sub, *FBy_sub, *FBz_sub; //@RMM double *FBx_dat, *FBy_dat, *FBz_dat; //@RMM From 0a48120bc2394964e683064b7be04afe7372d124 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Wed, 21 Feb 2024 16:13:05 +0100 Subject: [PATCH 12/19] Fixed H2OSOIL_LIQ calculation Subsurface layer thickness should be dz*dz_mult --- pfsimulator/amps/oas3/oas3_coupler.h | 8 ++++---- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 9 +++++---- pfsimulator/amps/oas3/send_fld2_clm.F90 | 5 +++-- pfsimulator/parflow_lib/solver_richards.c | 5 +++-- 4 files changed, 15 insertions(+), 12 deletions(-) diff --git a/pfsimulator/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index dc1bbd3ae..0acfabc0d 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -22,10 +22,10 @@ void oas_pfl_rcv(int *kid, int *kstep, double *kdata, #define send_fld2_clm send_fld2_clm_ -#define CALL_send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep, porosity, dz) \ - send_fld2_clm(pressure, saturation, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep, porosity, dz) -void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, int*iy, - int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep, double *porosity, double *dz); +#define CALL_send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep, porosity, dz_mult, dz) \ + send_fld2_clm(pressure, saturation, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep, porosity, dz_mult, &dz) +void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, int *iy, + int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep, double *porosity, double *dz_mult, double *dz); #define receive_fld2_clm receive_fld2_clm_ #define CALL_receive_fld2_clm(evap_trans, ice_impedance, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 7cb9b34fb..0f112bdb7 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -123,7 +123,7 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, if (ierror /= 0) call oasis_abort (comp_id, 'oas_pfl_define', 'oasis_enddef failed') end subroutine oas_pfl_define - subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep, porosity, dz) bind(c, name='send_fld2_clm_') + subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep, porosity,dz_mult,dz) bind(c, name='send_fld2_clm_') !---------------------------------------------------------------------------- ! Sends pressure head and soil liquid water from ParFlow to eCLM ! pressure is converted from [m] to [mm] @@ -137,7 +137,8 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n saturation((nx+2)*(ny+2)*(nz+2)), & ! saturation [-] topo((nx+2)*(ny+2)*(nz+2)), & ! topography mask (0 for inactive, 1 for active) porosity((nx+2)*(ny+2)*(nz+2)), & ! porosity [m^3/m^3] - dz((nx+2)*(ny+2)*(nz+2)) ! subsurface layer thickness [m] + dz_mult((nx+2)*(ny+2)*(nz+2)), & ! dz scaling values [unitless] + dz ! size of gridcell in z-direction [m] ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the ! extra "+2" terms account for the ghost nodes/halo points ! Local variables @@ -162,8 +163,8 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n do k = 1, nlevgrnd ! eCLM: 1=topmost layer, nlevgrnd=deepest layer z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) - pressure_3d(i,j,k) = pressure(l)*1000.0 ! multiply these quantities by 1000 - h2osoi_liq_3d(i,j,k) = saturation(l)*porosity(l)*dz(l)*1000 ! to convert from [m] to [mm] + pressure_3d(i,j,k) = pressure(l)*1000.0 ! multiply these quantities by 1000 + h2osoi_liq_3d(i,j,k) = saturation(l)*porosity(l)*dz*dz_mult(l)*1000 ! to convert from [m] to [mm] end do end if end do diff --git a/pfsimulator/amps/oas3/send_fld2_clm.F90 b/pfsimulator/amps/oas3/send_fld2_clm.F90 index 2c6ecf729..753a838c5 100644 --- a/pfsimulator/amps/oas3/send_fld2_clm.F90 +++ b/pfsimulator/amps/oas3/send_fld2_clm.F90 @@ -1,4 +1,4 @@ -SUBROUTINE send_fld2_clm(pressure,saturation,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep,porosity,dz) +SUBROUTINE send_fld2_clm(pressure,saturation,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep,porosity,dz_mult,dz) !---------------------------------------------------------------------------- ! @@ -46,7 +46,8 @@ SUBROUTINE send_fld2_clm(pressure,saturation,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep saturation((nx+2)*(ny+2)*(nz+2)), &! saturation (-) topo((nx+2)*(ny+2)*(nz+2)), &! mask (0 for inactive, 1 for active) porosity((nx+2)*(ny+2)*(nz+2)), &! porosity [m^3/m^3] - dz((nx+2)*(ny+2)*(nz+2)) ! subsurface layer thickness [m] + dz_mult((nx+2)*(ny+2)*(nz+2)), &! dz scaling values [unitless] + dz &! size of gridcell in z-direction [m] ! All vecotrs from parflow on grid w/ ghost nodes for current proc !Local Variables diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 5d9b11231..71c059ce1 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -1757,6 +1757,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time ForSubgridI(is, GridSubgrids(grid)) { int ix, iy, nx, ny, nz, nx_f, ny_f; + double dz; subgrid = GridSubgrid(grid, is); @@ -1774,10 +1775,10 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time nx = SubgridNX(subgrid); ny = SubgridNY(subgrid); nz = SubgridNZ(subgrid); + dz = SubgridDZ(subgrid); nx_f = SubvectorNX(et_sub); ny_f = SubvectorNY(et_sub); - sp = SubvectorData(s_sub); pp = SubvectorData(p_sub); et = SubvectorData(et_sub); @@ -1787,7 +1788,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time dz_dat = SubvectorData(dz_sub); //CPS amps_Printf("Calling oasis send/receive for time %3.1f \n", t); CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, - t,po_dat,dz_dat); + t,po_dat,dz_dat,dz); amps_Sync(amps_CommWorld); CALL_receive_fld2_clm(et, ice_impedance, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); } From 612082caf8ce76e27466df88d050708589e1e19d Mon Sep 17 00:00:00 2001 From: kvrigor Date: Fri, 23 Feb 2024 15:35:02 +0100 Subject: [PATCH 13/19] Ice impedance bugfixes - Set initial ice impedance values to 1.0 - solver_richards.c: ice_impedance in wrong arg position when calling NonlinSolverInvoke - parflow_proto.h: Make eCLM-related fields easier to track by moving arg position of ice_impedance right next to evap_trans --- CMakeLists.txt | 4 ++-- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 2 +- pfsimulator/parflow_lib/nl_function_eval.c | 10 ++++------ pfsimulator/parflow_lib/parflow_proto.h | 8 ++++---- pfsimulator/parflow_lib/problem_phase_rel_perm.c | 3 ++- pfsimulator/parflow_lib/richards_jacobian_eval.c | 6 ------ pfsimulator/parflow_lib/solver_richards.c | 8 ++++---- 7 files changed, 17 insertions(+), 24 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index bb45fc071..915da014e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -499,10 +499,10 @@ if ( ${PARFLOW_HAVE_ECLM} ) endif (HAVE_CLM) if (NOT HAVE_OAS3) message(FATAL_ERROR "ERROR: OASIS is required when compiling with PARFLOW_HAVE_ECLM=ON (check OAS3_ROOT variable)") - else (HAVE_OAS3) + else() message(WARNING "Enabling components for ParFlow-eCLM coupled model (experimental feature!)") set (HAVE_ECLM ${PARFLOW_HAVE_ECLM}) - endif (HAVE_OAS3) + endif (NOT HAVE_OAS3) endif ( ${PARFLOW_HAVE_ECLM} ) # # Parflow specific configuration options diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 0f112bdb7..7426d3976 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -212,7 +212,7 @@ subroutine receive_fld2_clm(evap_trans, ice_impedance, topo, ix, iy, nx, ny, nz, ! Save ET fluxes to ParFlow evap_trans vector evap_trans = 0. - ice_impedance = 0. + ice_impedance = 1. top_z_level = get_top_z_level(nx, ny, nz, topo) do i = 1, nx do j = 1, ny diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index 9fcf4b810..b2bc14bce 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -85,8 +85,8 @@ void KINSolFunctionEval( double dt = StateDt(((State*)current_state)); double time = StateTime(((State*)current_state)); Vector *evap_trans = StateEvapTrans(((State*)current_state)); - Vector *ovrl_bc_flx = StateOvrlBcFlx(((State*)current_state)); Vector *ice_impedance = StateIceImpedance(((State*)current_state)); + Vector *ovrl_bc_flx = StateOvrlBcFlx(((State*)current_state)); /* velocity vectors jjb */ Vector *x_velocity = StateXvel(((State*)current_state)); @@ -98,7 +98,7 @@ void KINSolFunctionEval( PFModuleInvokeType(NlFunctionEvalInvoke, nl_function_eval, (pressure, fval, problem_data, saturation, old_saturation, density, old_density, dt, time, old_pressure, evap_trans, - ovrl_bc_flx, ice_impedance, x_velocity, y_velocity, z_velocity)); + ice_impedance, ovrl_bc_flx, x_velocity, y_velocity, z_velocity)); return; } @@ -119,8 +119,8 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double time, /* New time value */ Vector * old_pressure, Vector * evap_trans, /*sk sink term from land surface model*/ - Vector * ovrl_bc_flx, /*sk overland flow boundary fluxes*/ Vector * ice_impedance, /* soil ice impedance */ + Vector * ovrl_bc_flx, /*sk overland flow boundary fluxes*/ Vector * x_velocity, /* velocity vectors jjb */ Vector * y_velocity, Vector * z_velocity) @@ -178,8 +178,6 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double *z_mult_dat; //@RMM #ifdef HAVE_ECLM - if (!amps_Rank(amps_CommWorld)) - amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); Vector *FBx = ice_impedance; Vector *FBy = ice_impedance; Vector *FBz = ice_impedance; @@ -189,10 +187,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); #endif + Subvector *FBx_sub, *FBy_sub, *FBz_sub; //@RMM double *FBx_dat=NULL, *FBy_dat=NULL, *FBz_dat=NULL; //@RMM - double gravity = ProblemGravity(problem); double viscosity = ProblemPhaseViscosity(problem, 0); diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h index 1e9596317..43a3c5dc8 100644 --- a/pfsimulator/parflow_lib/parflow_proto.h +++ b/pfsimulator/parflow_lib/parflow_proto.h @@ -293,14 +293,14 @@ PFModule *InputRFNewPublicXtra(char *geom_name); void InputRFFreePublicXtra(void); int InputRFSizeOfTempData(void); -typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef int (*NonlinSolverInvoke) (Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ice_impedance, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NonlinSolverInitInstanceXtraInvoke) (Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); /* kinsol_nonlin_solver.c */ int KINSolInitPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtemp1, N_Vector vtemp2, void *nl_function, double uround, long int *nfePtr, void *current_state); int KINSolCallPC(int neq, N_Vector pressure, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vtem, N_Vector ftem, void *nl_function, double uround, long int *nfePtr, void *current_state); void PrintFinalStats(FILE *out_file, long int *integer_outputs_now, long int *integer_outputs_total); -int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +int KinsolNonlinSolver(Vector *pressure, Vector *density, Vector *old_density, Vector *saturation, Vector *old_saturation, double t, double dt, ProblemData *problem_data, Vector *old_pressure, Vector *evap_trans, Vector *ice_impedance, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *KinsolNonlinSolverInitInstanceXtra(Problem *problem, Grid *grid, ProblemData *problem_data, double *temp_data); void KinsolNonlinSolverFreeInstanceXtra(void); PFModule *KinsolNonlinSolverNewPublicXtra(void); @@ -433,12 +433,12 @@ void FreeTempVector(Vector *vector); /* new_endpts.c */ void NewEndpts(double *alpha, double *beta, double *pp, int *size_ptr, int n, double *a_ptr, double *b_ptr, double *cond_ptr, double ereps); -typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +typedef void (*NlFunctionEvalInvoke) (Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ice_impedance, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); typedef PFModule *(*NlFunctionEvalInitInstanceXtraInvoke) (Problem *problem, Grid *grid, double *temp_data); /* nl_function_eval.c */ void KINSolFunctionEval(int size, N_Vector pressure, N_Vector fval, void *current_state); -void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ovrl_bc_flx, Vector *ice_impedance, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); +void NlFunctionEval(Vector *pressure, Vector *fval, ProblemData *problem_data, Vector *saturation, Vector *old_saturation, Vector *density, Vector *old_density, double dt, double time, Vector *old_pressure, Vector *evap_trans, Vector *ice_impedance, Vector *ovrl_bc_flx, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity); PFModule *NlFunctionEvalInitInstanceXtra(Problem *problem, Grid *grid, double *temp_data); void NlFunctionEvalFreeInstanceXtra(void); PFModule *NlFunctionEvalNewPublicXtra(char *name); diff --git a/pfsimulator/parflow_lib/problem_phase_rel_perm.c b/pfsimulator/parflow_lib/problem_phase_rel_perm.c index 5716d8475..26ade3e43 100644 --- a/pfsimulator/parflow_lib/problem_phase_rel_perm.c +++ b/pfsimulator/parflow_lib/problem_phase_rel_perm.c @@ -273,6 +273,8 @@ static inline double VanGLookupSpline( int max = num_sample_points + 1; // This table goes from 0 to fabs(min_pressure_head) + //if (!amps_Rank(amps_CommWorld)) + // amps_Printf("DEBUG[VanGLookupSpline() @ problem_phase_rel_perm.c]: pressure_head = %3.6e \n",pressure_head); assert(pressure_head >= 0); // SGS TODO add warning in output? @@ -702,7 +704,6 @@ void PhaseRelPerm( else { double head = fabs(ppdat[ipp]) / (pddat[ipd] * gravity); - prdat[ipr] = VanGLookupSpline(head, dummy1->lookup_tables[ir], CALCFCN); diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index 95fed7186..f9174044e 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -258,12 +258,6 @@ void RichardsJacobianEval( /* @RMM Flow Barrier / Boundary values */ #ifdef HAVE_ECLM - if (!amps_Rank(amps_CommWorld)) - { - amps_Printf("DEBUG: Initializing Flow barrier vectors with ice impedance values from eCLM.\n"); - //*(int*)0 = 0; - } - Vector *FBx = ice_impedance; Vector *FBy = ice_impedance; Vector *FBz = ice_impedance; diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 71c059ce1..b585ca55d 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -242,6 +242,7 @@ typedef struct { Vector *mask; Vector *evap_trans; /* sk: Vector that contains the sink terms from the land surface model */ + Vector *ice_impedance; /* soil ice impedance from eCLM */ Vector *evap_trans_sum; /* running sum of evaporation and transpiration */ Vector *overland_sum; Vector *ovrl_bc_flx; /* vector containing outflow at the boundary */ @@ -249,7 +250,6 @@ typedef struct { Vector *x_velocity; /* vector containing x-velocity face values */ Vector *y_velocity; /* vector containing y-velocity face values */ Vector *z_velocity; /* vector containing z-velocity face values */ - Vector *ice_impedance; /* soil ice impedance from eCLM */ #ifdef HAVE_CLM /* RM: vars for pf printing of clm output */ Vector *eflx_lh_tot; /* total LH flux from canopy height to atmosphere [W/m^2] */ @@ -852,7 +852,7 @@ SetupRichards(PFModule * this_module) InitVectorAll(instance_xtra->evap_trans, 0.0); instance_xtra->ice_impedance = NewVectorType(grid, 1, 1, vector_cell_centered); - InitVectorAll(instance_xtra->ice_impedance, 0.0); + InitVectorAll(instance_xtra->ice_impedance, 1.0); if (public_xtra->evap_trans_file) { @@ -1795,7 +1795,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); - handle = InitVectorUpdate(instance_xtra->ice_impedance, VectorUpdateAll); + handle = InitVectorUpdate(instance_xtra->ice_impedance, VectorUpdateAll); FinalizeVectorUpdate(handle); #endif // end to HAVE_OAS3 CALL @@ -2888,8 +2888,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time problem_data, instance_xtra->old_pressure, evap_trans, - instance_xtra->ovrl_bc_flx, instance_xtra->ice_impedance, + instance_xtra->ovrl_bc_flx, instance_xtra->x_velocity, instance_xtra->y_velocity, instance_xtra->z_velocity)); From ff4c7b02afa44692863fd08305accb4126acfcc6 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Thu, 29 Feb 2024 15:56:57 +0100 Subject: [PATCH 14/19] Don't initialized ice impedance when not in coupled mode --- pfsimulator/parflow_lib/solver_richards.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index b585ca55d..258c056a0 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -851,8 +851,12 @@ SetupRichards(PFModule * this_module) instance_xtra->evap_trans = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->evap_trans, 0.0); +#ifdef HAVE_ECLM instance_xtra->ice_impedance = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->ice_impedance, 1.0); +#else + instance_xtra->ice_impedance = NULL; +#fi if (public_xtra->evap_trans_file) { From c12995ae470a8c8129b72b13063b936ca0114bed Mon Sep 17 00:00:00 2001 From: kvrigor Date: Thu, 29 Feb 2024 15:59:45 +0100 Subject: [PATCH 15/19] Fixed wrong preprocessor syntax --- pfsimulator/parflow_lib/solver_richards.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 697fd5c53..14d754113 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -856,7 +856,7 @@ SetupRichards(PFModule * this_module) InitVectorAll(instance_xtra->ice_impedance, 1.0); #else instance_xtra->ice_impedance = NULL; -#fi +#endif if (public_xtra->evap_trans_file) { From 6dbd5fbdd48ac5fcfefeda2cb18e7e3547b8c573 Mon Sep 17 00:00:00 2001 From: kvrigor Date: Fri, 1 Mar 2024 16:19:07 +0100 Subject: [PATCH 16/19] Fixed sytax error at send_fld2_clm.F90 --- pfsimulator/amps/oas3/send_fld2_clm.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pfsimulator/amps/oas3/send_fld2_clm.F90 b/pfsimulator/amps/oas3/send_fld2_clm.F90 index 753a838c5..079831fd1 100644 --- a/pfsimulator/amps/oas3/send_fld2_clm.F90 +++ b/pfsimulator/amps/oas3/send_fld2_clm.F90 @@ -47,7 +47,7 @@ SUBROUTINE send_fld2_clm(pressure,saturation,topo,ix,iy,nx,ny,nz,nx_f,ny_f,pstep topo((nx+2)*(ny+2)*(nz+2)), &! mask (0 for inactive, 1 for active) porosity((nx+2)*(ny+2)*(nz+2)), &! porosity [m^3/m^3] dz_mult((nx+2)*(ny+2)*(nz+2)), &! dz scaling values [unitless] - dz &! size of gridcell in z-direction [m] + dz ! size of gridcell in z-direction [m] ! All vecotrs from parflow on grid w/ ghost nodes for current proc !Local Variables From 7c37ce8827d2d10e5ac86ec5be2f4340a0ebf456 Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Mon, 21 Sep 2026 14:21:11 +0200 Subject: [PATCH 17/19] eCLM-ParFlow: Couple effective porosity from eCLM - Receive ECLM_EFFPOROSITY from eCLM; keep ParFlow's porosity at uncoupled grid points - Update the received porosity after each converged step, capped so saturation stays <= 1 and floored at 0.01 - Include in storage term with separate porosity at both time levels and a freeze term (ice_liq_delta) for the liquid eCLM phase change --- pfsimulator/amps/oas3/oas3_coupler.h | 6 +- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 37 +++++-- pfsimulator/parflow_lib/nl_function_eval.c | 25 +++++ pfsimulator/parflow_lib/problem.c | 9 ++ pfsimulator/parflow_lib/problem.h | 9 ++ pfsimulator/parflow_lib/solver_richards.c | 107 ++++++++++++++++++++- 6 files changed, 178 insertions(+), 15 deletions(-) diff --git a/pfsimulator/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index 59a3d59ce..54fd37636 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -28,7 +28,7 @@ void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep, double *porosity, double *dz); #define receive_fld2_clm receive_fld2_clm_ -#define CALL_receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ - receive_fld2_clm(evap_trans, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) -void receive_fld2_clm(double *evap_trans, double *topo, int *ix, int*iy, +#define CALL_receive_fld2_clm(evap_trans, porosity_recv, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ + receive_fld2_clm(evap_trans, porosity_recv, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) +void receive_fld2_clm(double *evap_trans, double *porosity_recv, double *topo, int *ix, int*iy, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep); diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index a4a779bec..161b94cfb 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -30,12 +30,15 @@ module oas_pfl_mod ! Local variables used by OASIS integer :: comp_id, rank, ierror - integer :: soilliq_id, psi_id, et_id ! coupling field IDs + integer :: soilliq_id, psi_id, et_id, eff_porosity_id ! coupling field IDs ! TODO: Get these values from eCLM instead of hardcoding it here integer :: nlevsoi = 20 ! Number of hydrologically active soil layers in eCLM integer :: nlevgrnd = 25 ! Total number of soil layers in eCLM + ! Marks a cell eCLM does not couple. Must match pfl_uncoupled in eCLM's lnd2atmMod.F90 + real(kind=8), parameter, public :: PFL_UNCOUPLED = -9999.d0 + contains subroutine oas_pfl_init(argc) bind(c, name='oas_pfl_init_') @@ -116,6 +119,8 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_PSI') call oasis_def_var (soilliq_id, "PFL_SOILLIQ", part_id, var_nodims, OASIS_Out, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_SOILLIQ') + call oasis_def_var (eff_porosity_id, "PFL_EFFPOROSITY", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) + if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_EFFPOROSITY') call oasis_enddef ( ierror ) if (ierror /= 0) call oasis_abort (comp_id, 'oas_pfl_define', 'oasis_enddef failed') @@ -176,16 +181,19 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n deallocate(pressure_3d) end subroutine send_fld2_clm - subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') - !----------------------------------------------- - ! Receives evapotranspiration fluxes from eCLM. - !----------------------------------------------- + subroutine receive_fld2_clm(evap_trans, porosity_recv, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') + !------------------------------------------------------------------------- + ! Receives evapotranspiration fluxes and effective porosity from eCLM. + ! + ! ParFlow ownly keep own porosity, where eCLM is not coupled (PFL_UNCOUPLED). + !------------------------------------------------------------------------- integer, intent(in) :: ix, iy, & ! starting coordinate of the subgrid (UNUSED) nx, ny, nz, & ! subgrid dimensions nx_f, ny_f ! dimensions of ET subvector real(kind=8), intent(in) :: pstep, & ! current model time in hours topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) - real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)) ! evapotranspiration [1/hrs] + real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)), & ! evapotranspiration [1/hrs] + porosity_recv((nx+2)*(ny+2)*(nz+2)) ! effective porosity from eCLM [m3/m3] ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the ! extra "+2" terms account for the ghost nodes/halo points @@ -196,14 +204,18 @@ subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, ps integer :: z ! subsurface level (z=nz topmost layer, z=1 deepest layer) integer :: top_z_level(nx,ny) ! topmost z level of active ParFlow cells real(kind=8), allocatable :: evap_trans_3d(:,:,:) ! Root ET fluxes received from eCLM [1/hrs] + real(kind=8), allocatable :: porosity_3d(:,:,:) ! Effective porosity received from eCLM [m3/m3] - ! Receive ET fluxes from eCLM + ! Receive ET fluxes and effective porosity from eCLM allocate(evap_trans_3d(nx,ny,nlevsoi)) + allocate(porosity_3d(nx,ny,nlevgrnd)) seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) + call oasis_get(eff_porosity_id, seconds_elapsed, porosity_3d, ierror) - ! Save ET fluxes to ParFlow evap_trans vector + ! Save the received fields to the ParFlow vectors evap_trans = 0. + porosity_recv = PFL_UNCOUPLED top_z_level = get_top_z_level(nx, ny, nz, topo) do i = 1, nx do j = 1, ny @@ -213,10 +225,17 @@ subroutine receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, ps l = flattened_array_index(i, j, z, nx_f, ny_f) evap_trans(l) = evap_trans_3d(i,j,k) end do + do k = 1, nlevgrnd + z = top_z_level(i,j) - (k-1) + if (z < 1) exit ! eCLM column reaches deeper than ParFlow + l = flattened_array_index(i, j, z, nx_f, ny_f) + porosity_recv(l) = porosity_3d(i,j,k) + end do end if end do end do + deallocate(porosity_3d) deallocate(evap_trans_3d) end subroutine receive_fld2_clm @@ -258,4 +277,4 @@ integer function flattened_array_index(i, j, k, nx, ny) flattened_array_index = (i+1) + (j*nx) + (k*nx*ny) end function flattened_array_index -end module oas_pfl_mod \ No newline at end of file +end module oas_pfl_mod diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index 5785e30f6..6b3d167a9 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -221,6 +221,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double *q_overlnd_x_ = NULL, *q_overlnd_y_ = NULL; Vector *porosity = ProblemDataPorosity(problem_data); +#ifdef HAVE_OAS3 + Vector *old_porosity = ProblemDataOldPorosity(problem_data); + Vector *ice_liq_delta = ProblemDataIceLiqDelta(problem_data); +#endif Vector *permeability_x = ProblemDataPermeabilityX(problem_data); Vector *permeability_y = ProblemDataPermeabilityY(problem_data); Vector *permeability_z = ProblemDataPermeabilityZ(problem_data); @@ -258,12 +262,18 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Subgrid *subgrid; Subvector *p_sub, *d_sub, *od_sub, *s_sub, *os_sub, *po_sub, *op_sub, *ss_sub, *et_sub; +#ifdef HAVE_OAS3 + Subvector *opo_sub, *ild_sub; +#endif Subvector *f_sub, *rp_sub, *permx_sub, *permy_sub, *permz_sub; Grid *grid = VectorGrid(pressure); Grid *grid2d = VectorGrid(x_sl); double *pp, *odp, *sp, *osp, *pop, *fp, *dp, *rpp, *opp, *ss, *et; +#ifdef HAVE_OAS3 + double *opop, *ildp; +#endif double *permxp, *permyp, *permzp; // Velocity data jjb @@ -338,6 +348,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ s_sub = VectorSubvector(saturation, is); os_sub = VectorSubvector(old_saturation, is); po_sub = VectorSubvector(porosity, is); +#ifdef HAVE_OAS3 + opo_sub = VectorSubvector(old_porosity, is); + ild_sub = VectorSubvector(ice_liq_delta, is); +#endif f_sub = VectorSubvector(fval, is); /* @RMM added to provide access to zmult */ @@ -389,6 +403,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ opp = SubvectorData(op_sub); osp = SubvectorData(os_sub); pop = SubvectorData(po_sub); +#ifdef HAVE_OAS3 + opop = SubvectorData(opo_sub); + ildp = SubvectorData(ild_sub); +#endif fp = SubvectorData(f_sub); GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, @@ -401,7 +419,14 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double del_x_slope = 1.0; double del_y_slope = 1.0; +#ifdef HAVE_OAS3 + /* Porosity varies over time, so storage terms cannot be factored. + * Ice_liq_delta is here to avoid double-counting on solver retries. */ + fp[ip] = (sp[ip] * dp[ip] * pop[ipo] - osp[ip] * odp[ip] * opop[ipo] + + ildp[ipo] * dp[ip]) * vol * del_x_slope * del_y_slope * z_mult_dat[ip]; +#else fp[ip] = (sp[ip] * dp[ip] - osp[ip] * odp[ip]) * pop[ipo] * vol * del_x_slope * del_y_slope * z_mult_dat[ip]; +#endif }); } diff --git a/pfsimulator/parflow_lib/problem.c b/pfsimulator/parflow_lib/problem.c index c8fe8c958..de8def67b 100644 --- a/pfsimulator/parflow_lib/problem.c +++ b/pfsimulator/parflow_lib/problem.c @@ -457,6 +457,11 @@ ProblemData *NewProblemData( ProblemDataPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); +#ifdef HAVE_OAS3 + ProblemDataOldPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataIceLiqDelta(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); +#endif + ProblemDataBCPressureData(problem_data) = NewBCPressureData(); ProblemDataWellData(problem_data) = NewWellData(); @@ -490,6 +495,10 @@ void FreeProblemData( FreeReservoirData(ProblemDataReservoirData(problem_data)); FreeBCPressureData(ProblemDataBCPressureData(problem_data)); FreeVector(ProblemDataPorosity(problem_data)); +#ifdef HAVE_OAS3 + FreeVector(ProblemDataOldPorosity(problem_data)); + FreeVector(ProblemDataIceLiqDelta(problem_data)); +#endif FreeVector(ProblemDataPermeabilityX(problem_data)); FreeVector(ProblemDataPermeabilityY(problem_data)); FreeVector(ProblemDataPermeabilityZ(problem_data)); diff --git a/pfsimulator/parflow_lib/problem.h b/pfsimulator/parflow_lib/problem.h index e331d6b6b..041a165be 100644 --- a/pfsimulator/parflow_lib/problem.h +++ b/pfsimulator/parflow_lib/problem.h @@ -155,6 +155,11 @@ typedef struct { Vector *porosity; +#ifdef HAVE_OAS3 + Vector *old_porosity; /* porosity used by previous timestep */ + Vector *ice_liq_delta; /* liquid-equivalent volume frozen [-] */ +#endif + Vector *specific_storage; //sk Vector *FBx; //RMM @@ -288,6 +293,10 @@ typedef struct { #define ProblemDataPermeabilityY(problem_data) ((problem_data)->permeability_y) #define ProblemDataPermeabilityZ(problem_data) ((problem_data)->permeability_z) #define ProblemDataPorosity(problem_data) ((problem_data)->porosity) +#ifdef HAVE_OAS3 +#define ProblemDataOldPorosity(problem_data) ((problem_data)->old_porosity) +#define ProblemDataIceLiqDelta(problem_data) ((problem_data)->ice_liq_delta) +#endif #define ProblemDataFBx(problem_data) ((problem_data)->FBx) //RMM #define ProblemDataFBy(problem_data) ((problem_data)->FBy) //RMM #define ProblemDataFBz(problem_data) ((problem_data)->FBz) //RMM diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 9e83629e8..08210a1d2 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -331,6 +331,10 @@ typedef struct { Vector *mask; Vector *evap_trans; /* sk: Vector that contains the sink terms from the land surface model */ +#ifdef HAVE_OAS3 + Vector *porosity_recv; /* effective porosity received from eCLM; 0 = uncoupled cell */ + Vector *old_porosity_recv; /* effective porosity at the previous coupling step */ +#endif Vector *evap_trans_sum; /* running sum of evaporation and transpiration */ Vector *overland_sum; Vector *ovrl_bc_flx; /* vector containing outflow at the boundary */ @@ -549,6 +553,10 @@ SetupRichards(PFModule * this_module) /* Do turning bands (and other stuff maybe) */ PFModuleInvokeType(SetProblemDataInvoke, set_problem_data, (problem_data)); +#ifdef HAVE_OAS3 + /* Init with no porosity change (first coupled step carries no freeze term) */ + PFVCopy(ProblemDataPorosity(problem_data), ProblemDataOldPorosity(problem_data)); +#endif ComputeTopAndBottom(problem, problem_data); if (public_xtra->print_top || public_xtra->write_silo_top) @@ -1080,6 +1088,15 @@ SetupRichards(PFModule * this_module) instance_xtra->evap_trans = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->evap_trans, 0.0); +#ifdef HAVE_OAS3 + instance_xtra->porosity_recv = NewVectorType(grid, 1, 1, vector_cell_centered); + instance_xtra->old_porosity_recv = NewVectorType(grid, 1, 1, vector_cell_centered); + /* Must match PFL_UNCOUPLED in oas_pfl_mod.F90 */ + InitVectorAll(instance_xtra->porosity_recv, -9999.0); + InitVectorAll(instance_xtra->old_porosity_recv, -9999.0); + InitVectorAll(ProblemDataIceLiqDelta(problem_data), 0.0); +#endif + if (public_xtra->evap_trans_file) { sprintf(filename, "%s", public_xtra->evap_trans_filename); @@ -1940,8 +1957,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time #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; + Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub, *prcv_sub; + double *pp, *sp, *et, *ms, *po_dat, *dz_dat, *prcv_dat; double sw_lat = .0; double sw_lon = .0; #endif @@ -2157,6 +2174,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time m_sub = VectorSubvector(instance_xtra->mask, is); po_sub = VectorSubvector(porosity, is); dz_sub = VectorSubvector(instance_xtra->dz_mult, is); + prcv_sub = VectorSubvector(instance_xtra->porosity_recv, is); ix = SubgridIX(subgrid); iy = SubgridIY(subgrid); @@ -2173,15 +2191,18 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); + prcv_dat = SubvectorData(prcv_sub); //CPS amps_Printf("Calling oasis send/receive for time %3.1f \n", t); CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t, po_dat, dz_dat); amps_Sync(amps_CommWorld); - CALL_receive_fld2_clm(et, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); + CALL_receive_fld2_clm(et, prcv_dat, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); } amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); + handle = InitVectorUpdate(instance_xtra->porosity_recv, VectorUpdateAll); + FinalizeVectorUpdate(handle); #endif // end to HAVE_OAS3 CALL // IMF: Added to include CLM dumps in file_number updating. @@ -2993,10 +3014,86 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time problem_data)); } +#ifdef HAVE_OAS3 + PFVCopy(porosity, ProblemDataOldPorosity(problem_data)); +#endif PFVCopy(instance_xtra->density, instance_xtra->old_density); PFVCopy(instance_xtra->saturation, instance_xtra->old_saturation); PFVCopy(instance_xtra->pressure, instance_xtra->old_pressure); + +#ifdef HAVE_OAS3 + /* Implement the effective porosity from eCLM. Cells eCLM does not couple keep + * the porosity ParFlow read from its own input. */ + { + GrGeomSolid *gr_domain = ProblemDataGrDomain(problem_data); + Vector *old_porosity = ProblemDataOldPorosity(problem_data); + Vector *ice_liq_delta = ProblemDataIceLiqDelta(problem_data); + Subvector *po_s, *opo_s, *ild_s, *os_s, *rcv_s, *orcv_s; + double *po_d, *opo_d, *ild_d, *os_d, *rcv_d, *orcv_d; + Subgrid *sg; + int i, j, k, r, is; + int ix, iy, iz, nx, ny, nz; + + /* Ice occupies h2osoi_ice/denice of pore volume but removes only + * h2osoi_ice/denh2o of liquid. */ + const double ice_water_ratio = 0.917; + const double porosity_min = 0.01; + /* Threshold for non-coupled values, needs to be out of physical range */ + const double coupled_min = -1000.0; + + ForSubgridI(is, GridSubgrids(grid)) + { + sg = GridSubgrid(grid, is); + + po_s = VectorSubvector(porosity, is); + opo_s = VectorSubvector(old_porosity, is); + ild_s = VectorSubvector(ice_liq_delta, is); + os_s = VectorSubvector(instance_xtra->old_saturation, is); + rcv_s = VectorSubvector(instance_xtra->porosity_recv, is); + orcv_s = VectorSubvector(instance_xtra->old_porosity_recv, is); + + po_d = SubvectorData(po_s); + opo_d = SubvectorData(opo_s); + ild_d = SubvectorData(ild_s); + os_d = SubvectorData(os_s); + rcv_d = SubvectorData(rcv_s); + orcv_d = SubvectorData(orcv_s); + + r = SubgridRX(sg); + ix = SubgridIX(sg); iy = SubgridIY(sg); iz = SubgridIZ(sg); + nx = SubgridNX(sg); ny = SubgridNY(sg); nz = SubgridNZ(sg); + + GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, + { + int ipo = SubvectorEltIndex(po_s, i, j, k); + int is_ = SubvectorEltIndex(os_s, i, j, k); + + if (rcv_d[ipo] > coupled_min) + { + /* No previous value means a fresh start or a restart: skip the freeze. */ + ild_d[ipo] = (orcv_d[ipo] > coupled_min) + ? (orcv_d[ipo] - rcv_d[ipo]) * ice_water_ratio + : 0.0; + + /* Cap the pore-space loss at the liquid left, so ice expansion cannot + * demand a saturation above one. Releases once the cell has drained. */ + double liq_target = opo_d[ipo] * os_d[is_] - ild_d[ipo]; + po_d[ipo] = pfmax(pfmax(rcv_d[ipo], liq_target), porosity_min); + } + else + { + ild_d[ipo] = 0.0; + } + }); + } + + PFVCopy(instance_xtra->porosity_recv, instance_xtra->old_porosity_recv); + + handle = InitVectorUpdate(porosity, VectorUpdateAll); + FinalizeVectorUpdate(handle); + } +#endif } else /* Not converged, so decrease time step */ { @@ -4904,6 +5001,10 @@ TeardownRichards(PFModule * this_module) FreeVector(instance_xtra->q_overlnd_y); FreeVector(instance_xtra->evap_trans); +#ifdef HAVE_OAS3 + FreeVector(instance_xtra->porosity_recv); + FreeVector(instance_xtra->old_porosity_recv); +#endif FreeVector(instance_xtra->evap_trans_sum); FreeVector(instance_xtra->overland_sum); From c5da35d12c19b5c5b393006265158a6793620bba Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 22 Sep 2026 09:43:19 +0200 Subject: [PATCH 18/19] eCLM-ParFlow: use HAVE_ECLM instead of HAVE_OAS3 for porosity coupling --- pfsimulator/parflow_lib/nl_function_eval.c | 12 ++++---- pfsimulator/parflow_lib/problem.c | 8 ++--- pfsimulator/parflow_lib/problem.h | 8 ++--- pfsimulator/parflow_lib/solver_richards.c | 36 ++++++++++++---------- 4 files changed, 29 insertions(+), 35 deletions(-) diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index 0274015d2..ec82af676 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -221,7 +221,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double *q_overlnd_x_ = NULL, *q_overlnd_y_ = NULL; Vector *porosity = ProblemDataPorosity(problem_data); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM Vector *old_porosity = ProblemDataOldPorosity(problem_data); Vector *ice_liq_delta = ProblemDataIceLiqDelta(problem_data); #endif @@ -269,7 +269,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Subgrid *subgrid; Subvector *p_sub, *d_sub, *od_sub, *s_sub, *os_sub, *po_sub, *op_sub, *ss_sub, *et_sub; -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM Subvector *opo_sub, *ild_sub; #endif Subvector *f_sub, *rp_sub, *permx_sub, *permy_sub, *permz_sub; @@ -278,7 +278,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ Grid *grid2d = VectorGrid(x_sl); double *pp, *odp, *sp, *osp, *pop, *fp, *dp, *rpp, *opp, *ss, *et; -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM double *opop, *ildp; #endif double *permxp, *permyp, *permzp; @@ -355,7 +355,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ s_sub = VectorSubvector(saturation, is); os_sub = VectorSubvector(old_saturation, is); po_sub = VectorSubvector(porosity, is); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM opo_sub = VectorSubvector(old_porosity, is); ild_sub = VectorSubvector(ice_liq_delta, is); #endif @@ -410,7 +410,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ opp = SubvectorData(op_sub); osp = SubvectorData(os_sub); pop = SubvectorData(po_sub); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM opop = SubvectorData(opo_sub); ildp = SubvectorData(ild_sub); #endif @@ -426,7 +426,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double del_x_slope = 1.0; double del_y_slope = 1.0; -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM /* Porosity varies over time, so storage terms cannot be factored. * Ice_liq_delta is here to avoid double-counting on solver retries. */ fp[ip] = (sp[ip] * dp[ip] * pop[ipo] - osp[ip] * odp[ip] * opop[ipo] diff --git a/pfsimulator/parflow_lib/problem.c b/pfsimulator/parflow_lib/problem.c index f140fa201..138952158 100644 --- a/pfsimulator/parflow_lib/problem.c +++ b/pfsimulator/parflow_lib/problem.c @@ -457,11 +457,9 @@ ProblemData *NewProblemData( ProblemDataPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM ProblemDataOldPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); ProblemDataIceLiqDelta(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); -#endif -#ifdef HAVE_ECLM ProblemDataIceImpedance(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); #endif @@ -498,11 +496,9 @@ void FreeProblemData( FreeReservoirData(ProblemDataReservoirData(problem_data)); FreeBCPressureData(ProblemDataBCPressureData(problem_data)); FreeVector(ProblemDataPorosity(problem_data)); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM FreeVector(ProblemDataOldPorosity(problem_data)); FreeVector(ProblemDataIceLiqDelta(problem_data)); -#endif -#ifdef HAVE_ECLM FreeVector(ProblemDataIceImpedance(problem_data)); #endif FreeVector(ProblemDataPermeabilityX(problem_data)); diff --git a/pfsimulator/parflow_lib/problem.h b/pfsimulator/parflow_lib/problem.h index 3dc522c68..0a2a5c07a 100644 --- a/pfsimulator/parflow_lib/problem.h +++ b/pfsimulator/parflow_lib/problem.h @@ -155,11 +155,9 @@ typedef struct { Vector *porosity; -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM Vector *old_porosity; /* porosity used by previous timestep */ Vector *ice_liq_delta; /* liquid-equivalent volume frozen [-] */ -#endif -#ifdef HAVE_ECLM Vector *ice_impedance; /* soil ice impedance from eCLM [-] */ #endif @@ -296,11 +294,9 @@ typedef struct { #define ProblemDataPermeabilityY(problem_data) ((problem_data)->permeability_y) #define ProblemDataPermeabilityZ(problem_data) ((problem_data)->permeability_z) #define ProblemDataPorosity(problem_data) ((problem_data)->porosity) -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM #define ProblemDataOldPorosity(problem_data) ((problem_data)->old_porosity) #define ProblemDataIceLiqDelta(problem_data) ((problem_data)->ice_liq_delta) -#endif -#ifdef HAVE_ECLM #define ProblemDataIceImpedance(problem_data) ((problem_data)->ice_impedance) #endif #define ProblemDataFBx(problem_data) ((problem_data)->FBx) //RMM diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index de4807aa0..2ccc4b62b 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -332,7 +332,7 @@ typedef struct { Vector *mask; Vector *evap_trans; /* sk: Vector that contains the sink terms from the land surface model */ -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM Vector *porosity_recv; /* effective porosity received from eCLM; 0 = uncoupled cell */ Vector *old_porosity_recv; /* effective porosity at the previous coupling step */ #endif @@ -554,7 +554,7 @@ SetupRichards(PFModule * this_module) /* Do turning bands (and other stuff maybe) */ PFModuleInvokeType(SetProblemDataInvoke, set_problem_data, (problem_data)); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM /* Init with no porosity change (first coupled step carries no freeze term) */ PFVCopy(ProblemDataPorosity(problem_data), ProblemDataOldPorosity(problem_data)); #endif @@ -1089,16 +1089,14 @@ SetupRichards(PFModule * this_module) instance_xtra->evap_trans = NewVectorType(grid, 1, 1, vector_cell_centered); InitVectorAll(instance_xtra->evap_trans, 0.0); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM instance_xtra->porosity_recv = NewVectorType(grid, 1, 1, vector_cell_centered); instance_xtra->old_porosity_recv = NewVectorType(grid, 1, 1, vector_cell_centered); /* Must match PFL_UNCOUPLED in oas_pfl_mod.F90 */ InitVectorAll(instance_xtra->porosity_recv, -9999.0); InitVectorAll(instance_xtra->old_porosity_recv, -9999.0); InitVectorAll(ProblemDataIceLiqDelta(problem_data), 0.0); -#ifdef HAVE_ECLM InitVectorAll(ProblemDataIceImpedance(problem_data), 1.0); -#endif #endif if (public_xtra->evap_trans_file) @@ -1961,8 +1959,12 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time #ifdef HAVE_OAS3 Grid *grid = (instance_xtra->grid); Subgrid *subgrid; - Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub, *prcv_sub, *iimp_sub; - double *pp, *sp, *et, *ms, *po_dat, *dz_dat, *prcv_dat, *iimp_dat; + Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; + double *pp, *sp, *et, *ms, *po_dat, *dz_dat; +#ifdef HAVE_ECLM + Subvector *prcv_sub, *iimp_sub; + double *prcv_dat, *iimp_dat; +#endif double sw_lat = .0; double sw_lon = .0; #endif @@ -2176,13 +2178,13 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time p_sub = VectorSubvector(instance_xtra->pressure, is); s_sub = VectorSubvector(instance_xtra->saturation, is); et_sub = VectorSubvector(evap_trans, is); -#ifdef HAVE_ECLM - iimp_sub = VectorSubvector(ProblemDataIceImpedance(problem_data), is); -#endif m_sub = VectorSubvector(instance_xtra->mask, is); po_sub = VectorSubvector(porosity, is); dz_sub = VectorSubvector(instance_xtra->dz_mult, is); +#ifdef HAVE_ECLM prcv_sub = VectorSubvector(instance_xtra->porosity_recv, is); + iimp_sub = VectorSubvector(ProblemDataIceImpedance(problem_data), is); +#endif ix = SubgridIX(subgrid); iy = SubgridIY(subgrid); @@ -2196,13 +2198,13 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time sp = SubvectorData(s_sub); pp = SubvectorData(p_sub); et = SubvectorData(et_sub); -#ifdef HAVE_ECLM - iimp_dat = SubvectorData(iimp_sub); -#endif ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); +#ifdef HAVE_ECLM prcv_dat = SubvectorData(prcv_sub); + iimp_dat = SubvectorData(iimp_sub); +#endif //CPS amps_Printf("Calling oasis send/receive for time %3.1f \n", t); CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t, po_dat, dz_dat, dz); @@ -2216,9 +2218,9 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); +#ifdef HAVE_ECLM handle = InitVectorUpdate(instance_xtra->porosity_recv, VectorUpdateAll); FinalizeVectorUpdate(handle); -#ifdef HAVE_ECLM handle = InitVectorUpdate(ProblemDataIceImpedance(problem_data), VectorUpdateAll); FinalizeVectorUpdate(handle); #endif @@ -3033,7 +3035,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time problem_data)); } -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM PFVCopy(porosity, ProblemDataOldPorosity(problem_data)); #endif PFVCopy(instance_xtra->density, instance_xtra->old_density); @@ -3041,7 +3043,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time instance_xtra->old_saturation); PFVCopy(instance_xtra->pressure, instance_xtra->old_pressure); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM /* Implement the effective porosity from eCLM. Cells eCLM does not couple keep * the porosity ParFlow read from its own input. */ { @@ -5031,7 +5033,7 @@ TeardownRichards(PFModule * this_module) FreeVector(instance_xtra->q_overlnd_y); FreeVector(instance_xtra->evap_trans); -#ifdef HAVE_OAS3 +#ifdef HAVE_ECLM FreeVector(instance_xtra->porosity_recv); FreeVector(instance_xtra->old_porosity_recv); #endif From 15015c3b0ca7b9a6c6111d4a520b43322b8d8f5a Mon Sep 17 00:00:00 2001 From: Stefan Poll Date: Tue, 22 Sep 2026 17:00:47 +0200 Subject: [PATCH 19/19] eCLM-ParFlow: compute ice impedance from the effective porosity - Remove PFL_ICE_IMPEDANCE; derive icefrac = 1 - porosity_recv/static_porosity - Impedance 10^(-e_ice*icefrac) per cell face in x, y and z; multiplies the flow barriers - New key Solver.ECLM.IceImpedanceFactor (default 6.0 = eCLM e_ice, 0 = off); add IceImpedance (z-face) as possible output --- pfsimulator/amps/oas3/oas3_coupler.h | 8 +- pfsimulator/amps/oas3/oas_pfl_mod.F90 | 24 ++-- pfsimulator/parflow_lib/nl_function_eval.c | 26 +++- pfsimulator/parflow_lib/problem.c | 10 +- pfsimulator/parflow_lib/problem.h | 10 +- .../parflow_lib/richards_jacobian_eval.c | 30 +++- pfsimulator/parflow_lib/solver_richards.c | 128 ++++++++++++++++-- 7 files changed, 189 insertions(+), 47 deletions(-) diff --git a/pfsimulator/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index 419e91c46..356e23c59 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -29,10 +29,10 @@ void send_fld2_clm(double *pressure, double *saturation, double *topo, int *ix, #define receive_fld2_clm receive_fld2_clm_ #ifdef HAVE_ECLM -/* eCLM also sends the effective porosity and the soil ice impedance */ -#define CALL_receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ - receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) -void receive_fld2_clm(double *evap_trans, double *porosity_recv, double *ice_impedance, double *topo, int *ix, int*iy, +/* eCLM also sends the effective porosity; ParFlow derives the ice impedance from it */ +#define CALL_receive_fld2_clm(evap_trans, porosity_recv, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ + receive_fld2_clm(evap_trans, porosity_recv, topo, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &pstep) +void receive_fld2_clm(double *evap_trans, double *porosity_recv, double *topo, int *ix, int*iy, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep); #else #define CALL_receive_fld2_clm(evap_trans, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) \ diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index 543c24e55..3760bcccb 100644 --- a/pfsimulator/amps/oas3/oas_pfl_mod.F90 +++ b/pfsimulator/amps/oas3/oas_pfl_mod.F90 @@ -30,7 +30,7 @@ module oas_pfl_mod ! Local variables used by OASIS integer :: comp_id, rank, ierror - integer :: soilliq_id, psi_id, et_id, eff_porosity_id, ice_impedance_id ! coupling field IDs + integer :: soilliq_id, psi_id, et_id, eff_porosity_id ! coupling field IDs ! TODO: Get these values from eCLM instead of hardcoding it here integer :: nlevsoi = 20 ! Number of hydrologically active soil layers in eCLM @@ -121,8 +121,6 @@ subroutine oas_pfl_define(nx, ny, pdx, pdy, ix, iy, sw_lon, sw_lat, nlon, nlat, if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_SOILLIQ') call oasis_def_var (eff_porosity_id, "PFL_EFFPOROSITY", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_EFFPOROSITY') - call oasis_def_var (ice_impedance_id, "PFL_ICE_IMPEDANCE", part_id, var_nodims, OASIS_In, OASIS_Real, ierror) - if (ierror /= 0) call oasis_abort(comp_id, 'oas_pfl_define', 'oasis_def_var failed for PFL_ICE_IMPEDANCE') call oasis_enddef ( ierror ) if (ierror /= 0) call oasis_abort (comp_id, 'oas_pfl_define', 'oasis_enddef failed') @@ -184,10 +182,9 @@ subroutine send_fld2_clm(pressure, saturation, topo, ix, iy, nx, ny, nz, nx_f, n deallocate(pressure_3d) end subroutine send_fld2_clm - subroutine receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') + subroutine receive_fld2_clm(evap_trans, porosity_recv, topo, ix, iy, nx, ny, nz, nx_f, ny_f, pstep) bind(c, name='receive_fld2_clm_') !------------------------------------------------------------------------- - ! Receives evapotranspiration fluxes, effective porosity and soil ice - ! impedance from eCLM. + ! Receives evapotranspiration fluxes and effective porosity from eCLM. ! ! ParFlow ownly keep own porosity, where eCLM is not coupled (PFL_UNCOUPLED). !------------------------------------------------------------------------- @@ -197,8 +194,7 @@ subroutine receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, ix, real(kind=8), intent(in) :: pstep, & ! current model time in hours topo((nx+2)*(ny+2)*(nz+2)) ! topography mask (0 for inactive, 1 for active) real(kind=8), intent(inout) :: evap_trans((nx+2)*(ny+2)*(nz+2)), & ! evapotranspiration [1/hrs] - porosity_recv((nx+2)*(ny+2)*(nz+2)), & ! effective porosity from eCLM [m3/m3] - ice_impedance((nx+2)*(ny+2)*(nz+2)) ! soil ice impedance [-] + porosity_recv((nx+2)*(ny+2)*(nz+2)) ! effective porosity from eCLM [m3/m3] ! (nx+2)*(ny+2)*(nz+2) = total number of subgrid cells; the ! extra "+2" terms account for the ghost nodes/halo points @@ -210,30 +206,25 @@ subroutine receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, ix, integer :: top_z_level(nx,ny) ! topmost z level of active ParFlow cells real(kind=8), allocatable :: evap_trans_3d(:,:,:) ! Root ET fluxes received from eCLM [1/hrs] real(kind=8), allocatable :: porosity_3d(:,:,:) ! Effective porosity received from eCLM [m3/m3] - real(kind=8), allocatable :: ice_impedance_3d(:,:,:) ! Soil ice impedance received from eCLM [unitless; values from 0 to 1] ! Receive ET fluxes and effective porosity from eCLM allocate(evap_trans_3d(nx,ny,nlevsoi)) allocate(porosity_3d(nx,ny,nlevgrnd)) - allocate(ice_impedance_3d(nx,ny,nlevgrnd)) seconds_elapsed = nint(pstep*3600.d0) call oasis_get(et_id, seconds_elapsed, evap_trans_3d, ierror) call oasis_get(eff_porosity_id, seconds_elapsed, porosity_3d, ierror) - call oasis_get(ice_impedance_id, seconds_elapsed, ice_impedance_3d, ierror) ! Save the received fields to the ParFlow vectors evap_trans = 0. porosity_recv = PFL_UNCOUPLED - ice_impedance = 1. top_z_level = get_top_z_level(nx, ny, nz, topo) do i = 1, nx do j = 1, ny if (top_z_level(i,j) > 0) then ! ***** Subsurface level indexing convention ***** - do k = 1, nlevgrnd ! eCLM: 1=topmost layer, nlevsoi=deepest layer + do k = 1, nlevsoi ! eCLM: 1=topmost layer, nlevsoi=deepest layer z = top_z_level(i,j) - (k-1) ! ParFlow: 1=deepest layer, nz=topmost layer l = flattened_array_index(i, j, z, nx_f, ny_f) - if (k <= nlevsoi) evap_trans(l) = evap_trans_3d(i,j,k) - ice_impedance(l) = ice_impedance_3d(i,j,k) + evap_trans(l) = evap_trans_3d(i,j,k) end do do k = 1, nlevgrnd z = top_z_level(i,j) - (k-1) @@ -245,7 +236,8 @@ subroutine receive_fld2_clm(evap_trans, porosity_recv, ice_impedance, topo, ix, end do end do - deallocate(evap_trans_3d, porosity_3d, ice_impedance_3d) + deallocate(porosity_3d) + deallocate(evap_trans_3d) end subroutine receive_fld2_clm function get_top_z_level(nx, ny, nz, topo) diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index ec82af676..f33e432ef 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -245,15 +245,16 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double *z_mult_dat; //@RMM /* @RMM Flow Barrier / Boundary values */ -#ifdef HAVE_ECLM - /* Soil ice impedance from eCLM replaces the flow barriers */ - Vector *FBx = ProblemDataIceImpedance(problem_data); - Vector *FBy = ProblemDataIceImpedance(problem_data); - Vector *FBz = ProblemDataIceImpedance(problem_data); -#else Vector *FBx = ProblemDataFBx(problem_data); Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); +#ifdef HAVE_ECLM + /* Soil ice impedance, applied on top of the flow barriers */ + Vector *IIx = ProblemDataIceImpedanceX(problem_data); + Vector *IIy = ProblemDataIceImpedanceY(problem_data); + Vector *IIz = ProblemDataIceImpedanceZ(problem_data); + Subvector *IIx_sub, *IIy_sub, *IIz_sub; + double *IIx_dat, *IIy_dat, *IIz_dat; #endif Subvector *FBx_sub, *FBy_sub, *FBz_sub; //@RMM double *FBx_dat = NULL, *FBy_dat = NULL, *FBz_dat = NULL; //@RMM @@ -558,6 +559,14 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ FBx_dat = SubvectorData(FBx_sub); FBy_dat = SubvectorData(FBy_sub); FBz_dat = SubvectorData(FBz_sub); +#ifdef HAVE_ECLM + IIx_sub = VectorSubvector(IIx, is); + IIy_sub = VectorSubvector(IIy, is); + IIz_sub = VectorSubvector(IIz, is); + IIx_dat = SubvectorData(IIx_sub); + IIy_dat = SubvectorData(IIy_sub); + IIz_dat = SubvectorData(IIz_sub); +#endif GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, { @@ -873,6 +882,11 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ u_right = u_right * FBx_dat[ip]; u_front = u_front * FBy_dat[ip]; u_upper = u_upper * FBz_dat[ip]; +#ifdef HAVE_ECLM + u_right = u_right * IIx_dat[ip]; + u_front = u_front * IIy_dat[ip]; + u_upper = u_upper * IIz_dat[ip]; +#endif /* velocity data jjb */ vx[vxi] = u_right / (ffx * z_mult_dat[ip] * del_y_slope); diff --git a/pfsimulator/parflow_lib/problem.c b/pfsimulator/parflow_lib/problem.c index 138952158..cd29f0107 100644 --- a/pfsimulator/parflow_lib/problem.c +++ b/pfsimulator/parflow_lib/problem.c @@ -460,7 +460,10 @@ ProblemData *NewProblemData( #ifdef HAVE_ECLM ProblemDataOldPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); ProblemDataIceLiqDelta(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); - ProblemDataIceImpedance(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataStaticPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataIceImpedanceX(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataIceImpedanceY(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataIceImpedanceZ(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); #endif ProblemDataBCPressureData(problem_data) = NewBCPressureData(); @@ -499,7 +502,10 @@ void FreeProblemData( #ifdef HAVE_ECLM FreeVector(ProblemDataOldPorosity(problem_data)); FreeVector(ProblemDataIceLiqDelta(problem_data)); - FreeVector(ProblemDataIceImpedance(problem_data)); + FreeVector(ProblemDataStaticPorosity(problem_data)); + FreeVector(ProblemDataIceImpedanceX(problem_data)); + FreeVector(ProblemDataIceImpedanceY(problem_data)); + FreeVector(ProblemDataIceImpedanceZ(problem_data)); #endif FreeVector(ProblemDataPermeabilityX(problem_data)); FreeVector(ProblemDataPermeabilityY(problem_data)); diff --git a/pfsimulator/parflow_lib/problem.h b/pfsimulator/parflow_lib/problem.h index 0a2a5c07a..38342cbf8 100644 --- a/pfsimulator/parflow_lib/problem.h +++ b/pfsimulator/parflow_lib/problem.h @@ -158,7 +158,10 @@ typedef struct { #ifdef HAVE_ECLM Vector *old_porosity; /* porosity used by previous timestep */ Vector *ice_liq_delta; /* liquid-equivalent volume frozen [-] */ - Vector *ice_impedance; /* soil ice impedance from eCLM [-] */ + Vector *static_porosity; /* porosity from the input (eCLM watsat) */ + Vector *ice_impedance_x; /* soil ice impedance on face ip|ip+1 [-] */ + Vector *ice_impedance_y; /* soil ice impedance on face ip|ip+sy [-] */ + Vector *ice_impedance_z; /* soil ice impedance on face ip|ip+sz [-] */ #endif Vector *specific_storage; //sk @@ -297,7 +300,10 @@ typedef struct { #ifdef HAVE_ECLM #define ProblemDataOldPorosity(problem_data) ((problem_data)->old_porosity) #define ProblemDataIceLiqDelta(problem_data) ((problem_data)->ice_liq_delta) -#define ProblemDataIceImpedance(problem_data) ((problem_data)->ice_impedance) +#define ProblemDataStaticPorosity(problem_data) ((problem_data)->static_porosity) +#define ProblemDataIceImpedanceX(problem_data) ((problem_data)->ice_impedance_x) +#define ProblemDataIceImpedanceY(problem_data) ((problem_data)->ice_impedance_y) +#define ProblemDataIceImpedanceZ(problem_data) ((problem_data)->ice_impedance_z) #endif #define ProblemDataFBx(problem_data) ((problem_data)->FBx) //RMM #define ProblemDataFBy(problem_data) ((problem_data)->FBy) //RMM diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index 7260e15ab..fe8a6c0b5 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -292,15 +292,16 @@ void RichardsJacobianEval( double *z_mult_dat; //@RMM /* RMM Flow Barrier / Boundary values */ -#ifdef HAVE_ECLM - /* Soil ice impedance from eCLM replaces the flow barriers */ - Vector *FBx = ProblemDataIceImpedance(problem_data); - Vector *FBy = ProblemDataIceImpedance(problem_data); - Vector *FBz = ProblemDataIceImpedance(problem_data); -#else Vector *FBx = ProblemDataFBx(problem_data); Vector *FBy = ProblemDataFBy(problem_data); Vector *FBz = ProblemDataFBz(problem_data); +#ifdef HAVE_ECLM + /* Soil ice impedance, applied on top of the flow barriers */ + Vector *IIx = ProblemDataIceImpedanceX(problem_data); + Vector *IIy = ProblemDataIceImpedanceY(problem_data); + Vector *IIz = ProblemDataIceImpedanceZ(problem_data); + Subvector *IIx_sub, *IIy_sub, *IIz_sub; + double *IIx_dat, *IIy_dat, *IIz_dat; #endif Subvector *FBx_sub, *FBy_sub, *FBz_sub; //@RMM double *FBx_dat, *FBy_dat, *FBz_dat; //@RMM @@ -597,6 +598,14 @@ void RichardsJacobianEval( FBx_dat = SubvectorData(FBx_sub); FBy_dat = SubvectorData(FBy_sub); FBz_dat = SubvectorData(FBz_sub); +#ifdef HAVE_ECLM + IIx_sub = VectorSubvector(IIx, is); + IIy_sub = VectorSubvector(IIy, is); + IIz_sub = VectorSubvector(IIz, is); + IIx_dat = SubvectorData(IIx_sub); + IIy_dat = SubvectorData(IIy_sub); + IIz_dat = SubvectorData(IIz_sub); +#endif GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, { @@ -669,6 +678,9 @@ void RichardsJacobianEval( double x_coeff = FBx_dat[ip] * dt * ffx * (1.0 / dx) * z_mult_dat[ip] * PMean(pp[ip], pp[ip + 1], permxp[ip], permxp[ip + 1]) / viscosity; +#ifdef HAVE_ECLM + x_coeff *= IIx_dat[ip]; +#endif double sym_west_temp = (-x_coeff @@ -699,6 +711,9 @@ void RichardsJacobianEval( double y_coeff = FBy_dat[ip] * dt * ffy * (1.0 / dy) * z_mult_dat[ip] * PMean(pp[ip], pp[ip + sy_v], permyp[ip], permyp[ip + sy_v]) / viscosity; +#ifdef HAVE_ECLM + y_coeff *= IIy_dat[ip]; +#endif double sym_south_temp = -y_coeff * RPMean(updir, 0.0, prod, prod_no) * y_dir_g_c; //RMM TFG contributions, SYMM @@ -735,6 +750,9 @@ void RichardsJacobianEval( double z_coeff = FBz_dat[ip] * dt * ffz * PMeanDZ(permzp[ip], permzp[ip + sz_v], z_mult_dat[ip], z_mult_dat[ip + sz_v]) / viscosity; +#ifdef HAVE_ECLM + z_coeff *= IIz_dat[ip]; +#endif double sym_lower_temp = -z_coeff * (1.0 / (dz * Mean(z_mult_dat[ip], z_mult_dat[ip + sz_v]))) * RPMean(lower_cond, upper_cond, prod, diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 2ccc4b62b..02dabef8c 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -69,6 +69,9 @@ typedef struct { double drop_tol; int max_iterations; int max_convergence_failures; /* maximum number of convergence failures that are allowed */ +#ifdef HAVE_ECLM + double ice_impedance_factor; /* e_ice in imped calculation; 0 switches impedance off */ +#endif int lsm; /* land surface model */ int terrain_following_grid; /* @RMM flag for terrain following grid in NL fn eval, sets sslopes=toposl */ int variable_dz; /* @RMM flag for variable dz-multipliers */ @@ -1096,7 +1099,44 @@ SetupRichards(PFModule * this_module) InitVectorAll(instance_xtra->porosity_recv, -9999.0); InitVectorAll(instance_xtra->old_porosity_recv, -9999.0); InitVectorAll(ProblemDataIceLiqDelta(problem_data), 0.0); - InitVectorAll(ProblemDataIceImpedance(problem_data), 1.0); + InitVectorAll(ProblemDataIceImpedanceX(problem_data), 1.0); + InitVectorAll(ProblemDataIceImpedanceY(problem_data), 1.0); + InitVectorAll(ProblemDataIceImpedanceZ(problem_data), 1.0); + + /* Keep the input porosity (eCLM watsat) to derive the ice fraction from the + * effective porosity. Zero outside the domain, thus also marks domain cells. */ + { + GrGeomSolid *gr_domain = ProblemDataGrDomain(problem_data); + Vector *porosity = ProblemDataPorosity(problem_data); + Vector *static_porosity = ProblemDataStaticPorosity(problem_data); + Subvector *po_s, *spo_s; + double *po_d, *spo_d; + Subgrid *sg; + int i, j, k, r, is; + int ix, iy, iz, nx, ny, nz; + + InitVectorAll(static_porosity, 0.0); + ForSubgridI(is, GridSubgrids(grid)) + { + sg = GridSubgrid(grid, is); + po_s = VectorSubvector(porosity, is); + spo_s = VectorSubvector(static_porosity, is); + po_d = SubvectorData(po_s); + spo_d = SubvectorData(spo_s); + + r = SubgridRX(sg); + ix = SubgridIX(sg); iy = SubgridIY(sg); iz = SubgridIZ(sg); + nx = SubgridNX(sg); ny = SubgridNY(sg); nz = SubgridNZ(sg); + + GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, + { + int ipo = SubvectorEltIndex(po_s, i, j, k); + spo_d[ipo] = po_d[ipo]; + }); + } + handle = InitVectorUpdate(static_porosity, VectorUpdateAll); + FinalizeVectorUpdate(handle); + } #endif if (public_xtra->evap_trans_file) @@ -1962,8 +2002,8 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; double *pp, *sp, *et, *ms, *po_dat, *dz_dat; #ifdef HAVE_ECLM - Subvector *prcv_sub, *iimp_sub; - double *prcv_dat, *iimp_dat; + Subvector *prcv_sub; + double *prcv_dat; #endif double sw_lat = .0; double sw_lon = .0; @@ -2183,7 +2223,6 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time dz_sub = VectorSubvector(instance_xtra->dz_mult, is); #ifdef HAVE_ECLM prcv_sub = VectorSubvector(instance_xtra->porosity_recv, is); - iimp_sub = VectorSubvector(ProblemDataIceImpedance(problem_data), is); #endif ix = SubgridIX(subgrid); @@ -2203,14 +2242,13 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time dz_dat = SubvectorData(dz_sub); #ifdef HAVE_ECLM prcv_dat = SubvectorData(prcv_sub); - iimp_dat = SubvectorData(iimp_sub); #endif //CPS amps_Printf("Calling oasis send/receive for time %3.1f \n", t); CALL_send_fld2_clm(pp, sp, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t, po_dat, dz_dat, dz); amps_Sync(amps_CommWorld); #ifdef HAVE_ECLM - CALL_receive_fld2_clm(et, prcv_dat, iimp_dat, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); + CALL_receive_fld2_clm(et, prcv_dat, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); #else CALL_receive_fld2_clm(et, ms, ix, iy, nx, ny, nz, nx_f, ny_f, t); #endif @@ -2221,8 +2259,6 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time #ifdef HAVE_ECLM handle = InitVectorUpdate(instance_xtra->porosity_recv, VectorUpdateAll); FinalizeVectorUpdate(handle); - handle = InitVectorUpdate(ProblemDataIceImpedance(problem_data), VectorUpdateAll); - FinalizeVectorUpdate(handle); #endif #endif // end to HAVE_OAS3 CALL @@ -3060,8 +3096,9 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time * h2osoi_ice/denh2o of liquid. */ const double ice_water_ratio = 0.917; const double porosity_min = 0.01; - /* Threshold for non-coupled values, needs to be out of physical range */ - const double coupled_min = -1000.0; + /* Threshold for non-coupled values. Needs to be <0 but rather close + * due to interpolation artefacts. */ + const double coupled_min = -0.01; ForSubgridI(is, GridSubgrids(grid)) { @@ -3109,6 +3146,69 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time }); } + /* Ice impedance per cell face, as in eCLM: imped = 10^(-e_ice*icefrac) with + * icefrac = theta_ice/watsat = 1 - phi_eff/watsat, averaged over the two cells + * of the face. Uses the received porosity, thus so it follows the + * ice eCLM have. Faces to cells outside the domain keep 1. */ + if (public_xtra->ice_impedance_factor != 0.0) + { + const double e_ice = public_xtra->ice_impedance_factor; + const double icefrac_min = 1.0e-10; + Vector *static_porosity = ProblemDataStaticPorosity(problem_data); + Subvector *spo_s, *iix_s, *iiy_s, *iiz_s; + double *spo_d, *iix_d, *iiy_d, *iiz_d; + int sy, sz; + +/* Ice fraction of cell n, -1 where eCLM does not couple. Fractions below + * icefrac_min are remap round-off, so an ice-free cell gets exactly 1. */ +#define ICE_FRACTION(n) ((rcv_d[n] > coupled_min && spo_d[n] > 0.0) \ + ? ((1.0 - rcv_d[n] / spo_d[n] < icefrac_min) ? 0.0 \ + : pfmin(1.0 - rcv_d[n] / spo_d[n], 1.0)) : -1.0) +/* With one coupled side, use it alone similar to eCLM column bottom treatment */ +#define FACE_ICE_FRACTION(fa, fb) (((fa) >= 0.0 && (fb) >= 0.0) \ + ? 0.5 * ((fa) + (fb)) : pfmax(pfmax((fa), (fb)), 0.0)) + + ForSubgridI(is, GridSubgrids(grid)) + { + sg = GridSubgrid(grid, is); + + rcv_s = VectorSubvector(instance_xtra->porosity_recv, is); + spo_s = VectorSubvector(static_porosity, is); + iix_s = VectorSubvector(ProblemDataIceImpedanceX(problem_data), is); + iiy_s = VectorSubvector(ProblemDataIceImpedanceY(problem_data), is); + iiz_s = VectorSubvector(ProblemDataIceImpedanceZ(problem_data), is); + + rcv_d = SubvectorData(rcv_s); + spo_d = SubvectorData(spo_s); + iix_d = SubvectorData(iix_s); + iiy_d = SubvectorData(iiy_s); + iiz_d = SubvectorData(iiz_s); + + sy = SubvectorNX(rcv_s); + sz = SubvectorNY(rcv_s) * sy; + + r = SubgridRX(sg); + ix = SubgridIX(sg); iy = SubgridIY(sg); iz = SubgridIZ(sg); + nx = SubgridNX(sg); ny = SubgridNY(sg); nz = SubgridNZ(sg); + + GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz, + { + int ip = SubvectorEltIndex(rcv_s, i, j, k); + double f_c = ICE_FRACTION(ip); + double f_x = ICE_FRACTION(ip + 1); + double f_y = ICE_FRACTION(ip + sy); + double f_z = ICE_FRACTION(ip + sz); + + /* Same face convention as FBx/FBy/FBz: [ip] is the face to ip+1, ip+sy, ip+sz */ + iix_d[ip] = (spo_d[ip + 1] > 0.0) ? pow(10.0, -e_ice * FACE_ICE_FRACTION(f_c, f_x)) : 1.0; + iiy_d[ip] = (spo_d[ip + sy] > 0.0) ? pow(10.0, -e_ice * FACE_ICE_FRACTION(f_c, f_y)) : 1.0; + iiz_d[ip] = (spo_d[ip + sz] > 0.0) ? pow(10.0, -e_ice * FACE_ICE_FRACTION(f_c, f_z)) : 1.0; + }); + } +#undef ICE_FRACTION +#undef FACE_ICE_FRACTION + } + PFVCopy(instance_xtra->porosity_recv, instance_xtra->old_porosity_recv); handle = InitVectorUpdate(porosity, VectorUpdateAll); @@ -4047,7 +4147,7 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time if (public_xtra->write_netcdf_ice_impedance) { sprintf(nc_postfix, "%05d", instance_xtra->file_number); - WritePFNC(file_prefix, nc_postfix, t, ProblemDataIceImpedance(problem_data), + WritePFNC(file_prefix, nc_postfix, t, ProblemDataIceImpedanceZ(problem_data), public_xtra->numVarTimeVariant, "ice_impedance", 3, false, public_xtra->numVarIni); any_file_dumped = 1; @@ -6382,6 +6482,12 @@ SolverRichardsNewPublicXtra(char *name) sprintf(key, "%s.MaxConvergenceFailures", name); public_xtra->max_convergence_failures = GetIntDefault(key, 3); +#ifdef HAVE_ECLM + /* Same default as e_ice in eCLM's clm_varcon.F90 */ + sprintf(key, "%s.ECLM.IceImpedanceFactor", name); + public_xtra->ice_impedance_factor = GetDoubleDefault(key, 6.0); +#endif + if (public_xtra->max_convergence_failures > 9) { amps_Printf("Warning: Input variable <%s> \n", key);