eCLM frozen soil coupling

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
 

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

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
+end module oas_pfl_mod

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 

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);

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));

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