diff --git a/cmake/parflow_config.h.in b/cmake/parflow_config.h.in index 4bee36e07..bfbfaea7b 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/amps/oas3/oas3_coupler.h b/pfsimulator/amps/oas3/oas3_coupler.h index 59a3d59ce..356e23c59 100644 --- a/pfsimulator/amps/oas3/oas3_coupler.h +++ b/pfsimulator/amps/oas3/oas3_coupler.h @@ -22,13 +22,21 @@ 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_ +#ifdef HAVE_ECLM +/* 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) \ 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, int *nx, int *ny, int *nz, int *nx_f, int *ny_f, double *pstep); +#endif diff --git a/pfsimulator/amps/oas3/oas_pfl_mod.F90 b/pfsimulator/amps/oas3/oas_pfl_mod.F90 index a4a779bec..3760bcccb 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,12 +119,14 @@ 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') 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] @@ -135,7 +140,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 @@ -160,8 +166,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 @@ -176,16 +182,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 +205,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 +226,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 +278,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/amps/oas3/send_fld2_clm.F90 b/pfsimulator/amps/oas3/send_fld2_clm.F90 index e4740a95f..aa75b2025 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 vectors from parflow on grid w/ ghost nodes for current proc !Local Variables diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index 5785e30f6..f33e432ef 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_ECLM + 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); @@ -244,6 +248,14 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ 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 @@ -258,12 +270,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_ECLM + 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_ECLM + double *opop, *ildp; +#endif double *permxp, *permyp, *permzp; // Velocity data jjb @@ -338,6 +356,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_ECLM + 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 +411,10 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ opp = SubvectorData(op_sub); osp = SubvectorData(os_sub); pop = SubvectorData(po_sub); +#ifdef HAVE_ECLM + 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 +427,14 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ double del_x_slope = 1.0; double del_y_slope = 1.0; +#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] + + 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 }); } @@ -526,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, { @@ -841,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 c8fe8c958..cd29f0107 100644 --- a/pfsimulator/parflow_lib/problem.c +++ b/pfsimulator/parflow_lib/problem.c @@ -457,6 +457,15 @@ ProblemData *NewProblemData( ProblemDataPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); +#ifdef HAVE_ECLM + ProblemDataOldPorosity(problem_data) = NewVectorType(grid, 1, 1, vector_cell_centered); + ProblemDataIceLiqDelta(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(); ProblemDataWellData(problem_data) = NewWellData(); @@ -490,6 +499,14 @@ void FreeProblemData( FreeReservoirData(ProblemDataReservoirData(problem_data)); FreeBCPressureData(ProblemDataBCPressureData(problem_data)); FreeVector(ProblemDataPorosity(problem_data)); +#ifdef HAVE_ECLM + FreeVector(ProblemDataOldPorosity(problem_data)); + FreeVector(ProblemDataIceLiqDelta(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)); FreeVector(ProblemDataPermeabilityZ(problem_data)); diff --git a/pfsimulator/parflow_lib/problem.h b/pfsimulator/parflow_lib/problem.h index e331d6b6b..38342cbf8 100644 --- a/pfsimulator/parflow_lib/problem.h +++ b/pfsimulator/parflow_lib/problem.h @@ -155,6 +155,15 @@ typedef struct { Vector *porosity; +#ifdef HAVE_ECLM + Vector *old_porosity; /* porosity used by previous timestep */ + Vector *ice_liq_delta; /* liquid-equivalent volume frozen [-] */ + 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 Vector *FBx; //RMM @@ -288,6 +297,14 @@ 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_ECLM +#define ProblemDataOldPorosity(problem_data) ((problem_data)->old_porosity) +#define ProblemDataIceLiqDelta(problem_data) ((problem_data)->ice_liq_delta) +#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 #define ProblemDataFBz(problem_data) ((problem_data)->FBz) //RMM diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c index 0fd566eb4..fe8a6c0b5 100644 --- a/pfsimulator/parflow_lib/richards_jacobian_eval.c +++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c @@ -295,6 +295,14 @@ void RichardsJacobianEval( 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 @@ -590,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, { @@ -662,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 @@ -692,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 @@ -728,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 9e83629e8..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 */ @@ -274,6 +277,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_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_qx_overland; /* write q_overlnd_x? */ @@ -331,6 +335,10 @@ typedef struct { Vector *mask; Vector *evap_trans; /* sk: Vector that contains the sink terms from the land surface model */ +#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 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 +557,10 @@ SetupRichards(PFModule * this_module) /* Do turning bands (and other stuff maybe) */ PFModuleInvokeType(SetProblemDataInvoke, set_problem_data, (problem_data)); +#ifdef HAVE_ECLM + /* 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 +1092,53 @@ 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->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); + 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) { sprintf(filename, "%s", public_xtra->evap_trans_filename); @@ -1942,6 +2001,10 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time Subgrid *subgrid; Subvector *p_sub, *s_sub, *et_sub, *m_sub, *po_sub, *dz_sub; double *pp, *sp, *et, *ms, *po_dat, *dz_dat; +#ifdef HAVE_ECLM + Subvector *prcv_sub; + double *prcv_dat; +#endif double sw_lat = .0; double sw_lon = .0; #endif @@ -2147,6 +2210,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); @@ -2157,31 +2221,45 @@ 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); +#ifdef HAVE_ECLM + prcv_sub = VectorSubvector(instance_xtra->porosity_recv, is); +#endif ix = SubgridIX(subgrid); iy = SubgridIY(subgrid); 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); ms = SubvectorData(m_sub); po_dat = SubvectorData(po_sub); dz_dat = SubvectorData(dz_sub); +#ifdef HAVE_ECLM + prcv_dat = SubvectorData(prcv_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); + t, po_dat, dz_dat, dz); amps_Sync(amps_CommWorld); +#ifdef HAVE_ECLM + 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 } amps_Sync(amps_CommWorld); handle = InitVectorUpdate(evap_trans, VectorUpdateAll); FinalizeVectorUpdate(handle); +#ifdef HAVE_ECLM + handle = InitVectorUpdate(instance_xtra->porosity_recv, VectorUpdateAll); + FinalizeVectorUpdate(handle); +#endif #endif // end to HAVE_OAS3 CALL // IMF: Added to include CLM dumps in file_number updating. @@ -2993,10 +3071,150 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time problem_data)); } +#ifdef HAVE_ECLM + 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_ECLM + /* 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 <0 but rather close + * due to interpolation artefacts. */ + const double coupled_min = -0.01; + + 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; + } + }); + } + + /* 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); + FinalizeVectorUpdate(handle); + } +#endif } else /* Not converged, so decrease time step */ { @@ -3711,6 +3929,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_impedance || public_xtra->write_netcdf_overland_sum || public_xtra->write_netcdf_overland_bc_flux) { @@ -3924,6 +4143,16 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time any_file_dumped = 1; } +#ifdef HAVE_ECLM + if (public_xtra->write_netcdf_ice_impedance) + { + sprintf(nc_postfix, "%05d", instance_xtra->file_number); + WritePFNC(file_prefix, nc_postfix, t, ProblemDataIceImpedanceZ(problem_data), + public_xtra->numVarTimeVariant, "ice_impedance", 3, + false, public_xtra->numVarIni); + any_file_dumped = 1; + } +#endif if (public_xtra->print_evaptrans_sum || public_xtra->write_silo_evaptrans_sum @@ -4904,6 +5133,10 @@ TeardownRichards(PFModule * this_module) FreeVector(instance_xtra->q_overlnd_y); FreeVector(instance_xtra->evap_trans); +#ifdef HAVE_ECLM + FreeVector(instance_xtra->porosity_recv); + FreeVector(instance_xtra->old_porosity_recv); +#endif FreeVector(instance_xtra->evap_trans_sum); FreeVector(instance_xtra->overland_sum); @@ -4984,6 +5217,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) { @@ -6240,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); @@ -6651,6 +6899,15 @@ SolverRichardsNewPublicXtra(char *name) } public_xtra->write_netcdf_evaptrans = switch_value; + 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_impedance = switch_value; + sprintf(key, "NetCDF.WriteEvapTransSum"); switch_name = GetStringDefault(key, "False"); switch_value = NA_NameToIndexExitOnError(switch_na, switch_name, key); @@ -6758,6 +7015,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_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 73fc5ad1d..418913010 100644 --- a/pfsimulator/parflow_lib/write_parflow_netcdf.c +++ b/pfsimulator/parflow_lib/write_parflow_netcdf.c @@ -961,6 +961,48 @@ int LookUpInventory(char * varName, varNCData **myVarNCData, int *netCDFIDs) return evaptransVarID; } + if (strcmp(varName, "ice_impedance") == 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 iceImpedanceVarID; + int res = nc_def_var(netCDFIDs[0], varName, (*myVarNCData)->ncType, (*myVarNCData)->dimSize, + (*myVarNCData)->dimIDs, &iceImpedanceVarID); + 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], iceImpedanceVarID, NC_CHUNKED, chunksize); + } + if (enable_netcdf_compression) { + nc_def_var_deflate(netCDFIDs[0],iceImpedanceVarID,0,1,compression_level); + } + + } + if (res == NC_ENAMEINUSE) + { + res = nc_inq_varid(netCDFIDs[0], varName, &iceImpedanceVarID); + } + return iceImpedanceVarID; + } + if (strcmp(varName, "evaptrans_sum") == 0) { *myVarNCData = (varNCData*)malloc(sizeof(varNCData));