CRNS method 2

bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90

@@ -49,7 +49,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 !
 ! !USES:
    USE mod_assimilation, &
-        ONLY: obs_index_p, &
+        ONLY: obs_index_p, obs_p, &
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
         sc_p, &
@@ -70,7 +70,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 
 #if defined CLMSA
    USE enkf_clm_mod, & 
-        ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T
+        ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd
 #endif
   IMPLICIT NONE
 
@@ -80,7 +80,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
   INTEGER, INTENT(in) :: dim_obs_p          ! Dimension of observed state
   REAL, INTENT(in)    :: state_p(dim_p)     ! PE-local model state
   REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PE-local observed state
-  integer :: i, j, k
+  integer :: i, j, k, z, n
   integer :: icorner
   logical :: lpointobs       !If true: no special observation; use point observation
 ! !CALLING SEQUENCE:
@@ -100,6 +100,21 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 integer :: nsc
 ! end of hcp 
 
+#ifndef PARFLOW_STAND_ALONE
+#ifndef OBS_ONLY_PARFLOW
+! Variables used in crns version 2
+REAL :: weights_r1(920), weights_r2(920), weights_r3(920)
+Real :: weights_layer(8)
+Integer :: nweights(8)
+Real  :: d86_r1, d86_r2, d86_r3
+REAL :: r1 = 1.0
+REAL :: r2 = 20.0
+REAL :: r3 = 85.0
+REAL :: bd, y
+REAL :: sum_r1, sum_r2, sum_r3, totw
+#endif
+#endif
+
 
 ! *********************************************
 ! *** Perform application of measurement    ***
@@ -194,6 +209,94 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 #endif
 #endif
 
+#ifndef PARFLOW_STAND_ALONE
+#ifndef OBS_ONLY_PARFLOW
+ if (crns_flag.EQ.2) then
+   lpointobs = .false.
+   ! CRNS implementation based on Schrön et al. 2017 using
+   ! d86 for 3 different radius values
+   DO i = 1, dim_obs_p
+     ! Bulk density average for the 8 considered layers
+     IF (clmcrns_bd > 0.0) THEN
+       bd = clmcrns_bd
+     ELSE
+       bd = 0.0
+       DO j = 1, 8
+         bd = bd + soilstate_inst%bd_col(obs_index_p(i),j) ! bulk density
+       END DO
+         bd = bd / 8.0 * 0.001 ! average and convert from kg/m^3 to g/cm^3
+     ENDIF
+     ! CRNS observed value
+     y = obs_p(i) ! CRNS observation
+     ! Penetration depth calculations D86(bd, r, y)
+     d86_r1 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r1))*(20+y)/(0.0429+y)))
+     d86_r2 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r2))*(20+y)/(0.0429+y)))
+     d86_r3 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r3))*(20+y)/(0.0429+y)))
+     ! Then calculate the weights for thin (1mm) slices of the layers for 85cm
+     sum_r1 = 0.0
+     sum_r2 = 0.0
+     sum_r3 = 0.0
+     DO j = 1, 920 ! depth in mm but in calculation used in cm:
+       weights_r1(j) = exp(-2*(j/10)/d86_r1)
+       weights_r2(j) = exp(-2*(j/10)/d86_r2)
+       weights_r3(j) = exp(-2*(j/10)/d86_r3)
+
+       sum_r1 = sum_r1 + weights_r1(j)
+       sum_r2 = sum_r2 + weights_r2(j)
+       sum_r3 = sum_r3 + weights_r3(j)
+     END DO
+     ! Normalize the weights:
+     weights_r1(:) = weights_r1(:) / sum_r1
+     weights_r2(:) = weights_r2(:) / sum_r2
+     weights_r3(:) = weights_r3(:) / sum_r3
+     ! assign average weights to each layer
+     nweights(:) = 0
+     weights_layer(:) = 0.0
+     ! z index for different layers, manually here, could be done better with model layer depth
+     DO j = 1, 920
+       IF (j > 680) then
+         z = 8
+       ELSEIF (j < 680 .and. j > 480) then
+         z = 7
+       ELSEIF (j < 480 .and. j > 320) then
+         z = 6
+       ELSEIF (j < 320 .and. j > 200) then
+         z = 5
+       ELSEIF (j < 200 .and. j > 120) then
+         z = 4
+       ELSEIF (j < 120 .and. j > 60) then
+         z = 3
+       ELSEIF (j < 60 .and. j > 20) then
+         z = 2
+       ELSEIF (j < 20) then
+         z = 1
+       ENDIF
+       weights_layer(z) = weights_layer(z) + weights_r1(j) + weights_r2(j) + weights_r3(j)
+       nweights(z) = nweights(z) + 1
+     END DO
+     ! Normalize the weights
+     totw = 0.0
+     DO j = 1, 8
+        weights_layer(j) = weights_layer(j) / (nweights(j) * 3)
+        totw = totw + weights_layer(j)
+     END DO
+     weights_layer(:) = weights_layer(:) / totw
+     ! Finally use the weights to calculate the weighted average of the state variable
+     avesm = 0.0
+     DO j = 1, 8
+       avesm = avesm + weights_layer(j) * state_p(obs_index_p(i) + (j-1))
+       ! This assumes that obs_index_p(i) for obs i is the index of 
+       ! the first layer of the gridcell where obs i is
+     END DO
+     ! Assign new average as the state variable
+     m_state_p(i) = avesm
+     ! end loop over observations
+   END DO
+ end if
+#endif
+#endif
+
+
  if(obs_interp_switch == 1) then
 
       lpointobs = .false.

bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90

@@ -58,6 +58,7 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmstatevec_max_layer")  :: clmstatevec_max_layer
   integer(c_int),bind(C,name="clmt_printensemble")       :: clmt_printensemble
   integer(c_int),bind(C,name="clmwatmin_switch")         :: clmwatmin_switch
+  real(c_double),bind(C,name="clmcrns_bd")      :: clmcrns_bd
 
   integer  :: nstep     ! time step index
   real(r8) :: dtime     ! time step increment (sec)

bldsva/intf_DA/pdaf/model/common/enkf.h

@@ -134,3 +134,4 @@ GLOBAL double pf_dampfac_state;
 GLOBAL double dampfac_state_time_dependent;
 GLOBAL double dampfac_param_time_dependent;
 GLOBAL double da_crns_depth_tol;
+GLOBAL double clmcrns_bd;

bldsva/intf_DA/pdaf/model/common/read_enkfpar.c

@@ -102,6 +102,7 @@ void read_enkfpar(char *parname)
   obs_interp_switch     = iniparser_getint(pardict,"DA:obs_interp_switch",0);
   crns_flag             = iniparser_getint(pardict,"DA:crns_flag",0);
   da_crns_depth_tol     = iniparser_getdouble(pardict,"DA:da_crns_depth_tol",0.01);
+  clmcrns_bd            = iniparser_getdouble(pardict, "DA:crns_bd", -1.0);
   da_print_obs_index    = iniparser_getint(pardict,"DA:print_obs_index",0);
   total_steps = (int) (t_sim/da_interval);
   tstartcycle = (int) (t_start/da_interval);

bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90

@@ -58,6 +58,7 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmstatevec_max_layer")  :: clmstatevec_max_layer
   integer(c_int),bind(C,name="clmt_printensemble")       :: clmt_printensemble
   integer(c_int),bind(C,name="clmwatmin_switch")         :: clmwatmin_switch
+  real(c_double),bind(C,name="clmcrns_bd")      :: clmcrns_bd
 
   integer  :: nstep     ! time step index
   real(r8) :: dtime     ! time step increment (sec)

doc/content/setup_tsmp/input_enkfpf.md

@@ -701,16 +701,35 @@ Effect of `obs_interp_switch=1`:
 
 ### DA:crns_flag ###
 
-`DA:crns_flag`: (int) Set to 1 will read the parflow soil moisture data 
-as averaged soil moisture with the conventional weighting profile proposed 
-in Schroen etal HESS 2017 to model CRNS observation. Default is zero.
+`DA:crns_flag`: (int) Flag for using CRNS-observations
+
+- `1`: ParFlow-based CRNS-measurement operator. ParFlow soil moisture
+data is averaged and used in conventional weighting profile proposed
+in Schroen etal HESS 2017 for modelling CRNS observations.
+
+- `2`: CLM-based CRNS-measurement operator. CLM-soil moisture are
+    weighted and used in the operator. CRNS implementation based on
+    Schrön et al. 2017 using d86 for 3 different radius values. Note
+    that the penetration depths are computed using the observations y.
+
+Default: `0` (no CRNS-observations).
+
+### DA:crns_bd ###
+
+`DA:crns_bd`: (int) Only used for `DA:crns_flag=2`. Prescribed bulk
+density average for the eight soil layers considered in
+CRNS-measurement operator.
+
+Default: `-1.0` (bulk density computed from CLM's
+`soilstate_inst%bd_col`).
 
 ### DA:da_crns_depth_tol ###
 
-`DA:da_crns_depth_tol`: (double) Convergence criteria for weighting
-procedure that "generate[s] a weighted average of [soil-moisture]
-point measurements that can be compared with [...] [a] CRNS product"
-(Schrön et al, 2017 http://dx.doi.org/10.5194/hess-21-5009-2017)
+`DA:da_crns_depth_tol`: (double) Only used for
+`DA:crns_flag=1`. Convergence criteria for weighting procedure that
+"generate[s] a weighted average of [soil-moisture] point measurements
+that can be compared with [...] [a] CRNS product" (Schrön et al, 2017
+http://dx.doi.org/10.5194/hess-21-5009-2017)
 
 Only used if CRNS-observations are turned on by adding parameter
 `depth` in observation input.