CRNS measurement operator with better skin depth evaluation

bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90

@@ -77,7 +77,6 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
 #endif
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
-       depth_obs_p, &
        sc_p, idx_obs_nc_p, &
 #endif
 #endif
@@ -92,10 +91,10 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
        clm_obs, &
        var_id_obs_nc, dim_nx, dim_ny, &
        clmobs_lon, clmobs_lat, clmobs_layer, clmobs_dr, clm_obserr, &
-       crns_flag, depth_obs, dampfac_state_time_dependent_in, dampfac_param_time_dependent_in
+       dampfac_state_time_dependent_in, dampfac_param_time_dependent_in
   use mod_tsmp, &
       only: idx_map_subvec2state_fortran, tag_model_parflow, enkf_subvecsize, &
-      nx_glob, ny_glob, nz_glob, &
+      nx_glob, ny_glob, nz_glob, crns_flag, &
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
       xcoord, ycoord, zcoord, xcoord_fortran, ycoord_fortran, &
@@ -148,7 +147,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   logical :: is_multi_observation_files
   character (len = 110) :: current_observation_filename
   integer,allocatable :: local_dis(:),local_dim(:)
-  integer :: k_count,nsc !hcp
+  integer :: k_count !,nsc !hcp
   real    :: sum_interp_weights
 
 #ifndef PARFLOW_STAND_ALONE
@@ -202,8 +201,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   call mpi_bcast(dim_obs, 1, MPI_INTEGER, 0, comm_filter, ierror)
   ! Switch for vector of observation errors
   call mpi_bcast(multierr, 1, MPI_INTEGER, 0, comm_filter, ierror)
-  ! broadcast crns_flag
-  call mpi_bcast(crns_flag, 1, MPI_INTEGER, 0, comm_filter, ierror)
+  !! broadcast crns_flag
+  !call mpi_bcast(crns_flag, 1, MPI_INTEGER, 0, comm_filter, ierror)
   ! broadcast dim_ny and dim_nx
   if(point_obs.eq.0) then
      call mpi_bcast(dim_nx, 1, MPI_INTEGER, 0, comm_filter, ierror)
@@ -279,11 +278,6 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
              allocate(pressure_obserr(dim_obs))
         endif
 
-        if (crns_flag.eq.1) then
-             if (allocated(depth_obs)) deallocate(depth_obs)
-             allocate(depth_obs(dim_obs))
-             depth_obs(:)=0.d0
-        endif
         if(allocated(idx_obs_nc)) deallocate(idx_obs_nc)
         allocate(idx_obs_nc(dim_obs))
         if(allocated(x_idx_obs_nc))deallocate(x_idx_obs_nc)
@@ -341,7 +335,6 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
       ! if exist ParFlow-type obs
      call mpi_bcast(pressure_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror)
      if(multierr.eq.1) call mpi_bcast(pressure_obserr, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror)
-     if(crns_flag.eq.1) call mpi_bcast(depth_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror)
      call mpi_bcast(idx_obs_nc, dim_obs, MPI_INTEGER, 0, comm_filter, ierror)
      ! broadcast xyz indices
      call mpi_bcast(x_idx_obs_nc, dim_obs, MPI_INTEGER, 0, comm_filter, ierror)
@@ -553,10 +546,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
         allocate(pressure_obserr_p(dim_obs_p))
      endif
      if(crns_flag.eq.1) then 
-        if (allocated(depth_obs_p)) deallocate(depth_obs_p)
-        allocate(depth_obs_p(dim_obs_p))
         if (allocated(sc_p)) deallocate(sc_p)
-        allocate(sc_p(dim_obs_p))
+        allocate(sc_p(nz_glob, dim_obs_p))
         if (allocated(idx_obs_nc_p)) deallocate(idx_obs_nc_p)
         allocate(idx_obs_nc_p(dim_obs_p))
      endif
@@ -618,6 +609,11 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
      end do
   else if (point_obs.eq.1) then
 
+     !hcp
+     if(crns_flag.eq.1) then
+         idx_obs_nc(:)=nx_glob*(y_idx_obs_nc(:)-1)+x_idx_obs_nc(:)
+     endif
+     !hcp fin
      count = 1
      do i = 1, dim_obs
         obs(i) = pressure_obs(i)  
@@ -631,9 +627,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
               obs_p(count) = pressure_obs(i)
               if(multierr.eq.1) pressure_obserr_p(count) = pressure_obserr(i)
               if(crns_flag.eq.1) then
-                  depth_obs_p(count) = depth_obs(i)
                   idx_obs_nc_p(count)=idx_obs_nc(i)
-                  Allocate(sc_p(count)%scol_obs_in(nz_glob-z_idx_obs_nc(i)+1))       
+                  !Allocate(sc_p(count)%scol_obs_in(nz_glob))       
               endif
               obs_nc2pdaf(local_dis(mype_filter+1)+count) = i
               count = count + 1
@@ -642,11 +637,10 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
      end do
      do i = 1, dim_obs_p
       if(crns_flag.eq.1) then 
-        nsc=size(sc_p(i)%scol_obs_in(:))
-        do k = 1, nsc
+        do k = 1, nz_glob
           k_count=idx_obs_nc_p(i)+(k-1)*nx_glob*ny_glob
           do j = 1, enkf_subvecsize
-             if (k_count .eq. idx_map_subvec2state_fortran(j)) sc_p(i)%scol_obs_in(nsc-k+1)=j
+             if (k_count .eq. idx_map_subvec2state_fortran(j)) sc_p(nz_glob-k+1,i)=j
           enddo
         enddo
       endif

bldsva/intf_DA/pdaf/framework/mod_assimilation.F90

@@ -72,12 +72,11 @@ MODULE mod_assimilation
   INTEGER, ALLOCATABLE :: local_dims_obs(:) ! Array for process-local observation dimensions
   INTEGER, ALLOCATABLE :: obs_nc2pdaf(:)  ! mapping ordering of obs between netcdf input and internal ordering in pdaf
   REAL, ALLOCATABLE :: pressure_obserr_p(:) ! Vector holding observation errors for paraflow run at each PE-local domain 
-  REAL, ALLOCATABLE :: depth_obs_p(:) ! Vector holding observation errors for paraflow run at each PE-local domain 
   !hcp
-  type :: scoltype
-       integer, dimension(:), allocatable :: scol_obs_in
-  endtype
-  type(scoltype), dimension(:), allocatable :: sc_p 
+  !type :: scoltype
+  !     integer, dimension(:), allocatable :: scol_obs_in
+  !endtype
+  integer, dimension(:,:), allocatable :: sc_p  !soil moisture of a soil column distributed over the procs of PE-local (i_z, i_obs)
   real, allocatable    :: idx_obs_nc_p(:)        
   INTEGER :: toffset      ! offset time step to shift all the assimilation steps
   !end hcp  

bldsva/intf_DA/pdaf/framework/mod_read_obs.F90

@@ -50,8 +50,7 @@ module mod_read_obs
   ! Flag: Use vector of observation errors in observation file
   integer :: multierr=0
   integer :: dim_nx, dim_ny
-  integer :: crns_flag=0   !hcp
-  real, allocatable :: depth_obs(:)   !hcp
+!  integer :: crns_flag=0   !hcp
   real, allocatable :: dampfac_state_time_dependent_in(:)
   real, allocatable :: dampfac_param_time_dependent_in(:)
 contains
@@ -73,7 +72,7 @@ subroutine read_obs_nc(current_observation_filename)
     !      only: comm_filter
     use mod_tsmp, &
         only: point_obs, obs_interp_switch, is_dampfac_state_time_dependent, &
-        is_dampfac_param_time_dependent
+        is_dampfac_param_time_dependent, crns_flag
     use netcdf
     implicit none
     integer :: ncid
@@ -252,16 +251,10 @@ subroutine read_obs_nc(current_observation_filename)
             end if
         endif
 
-        has_depth = nf90_inq_varid(ncid, depth_name, depth_varid)
-        if(has_depth == nf90_noerr) then
-            crns_flag = 1
-            if(allocated(depth_obs)) deallocate(depth_obs)
-            allocate(depth_obs(dim_obs))
-            call check(nf90_get_var(ncid, depth_varid, depth_obs))
-            if (screen > 2) then
-                print *, "TSMP-PDAF mype(w)=", mype_world, ": depth_obs=", depth_obs
-            end if
-        end if
+        !has_depth = nf90_inq_varid(ncid, depth_name, depth_varid)
+        !if(has_depth == nf90_noerr) then
+        !    crns_flag = 1
+        !end if
 
         ! Read the surface pressure and idxerature data from the file.
         ! Since we know the contents of the file we know that the data
@@ -290,6 +283,15 @@ subroutine read_obs_nc(current_observation_filename)
 
         call check( nf90_inq_varid(ncid, Z_IDX_NAME, z_idx_varid) )
         call check( nf90_get_var(ncid, z_idx_varid, z_idx_obs_nc) )
+        !hcp     
+        if  (crns_flag .EQ. 1) then
+            z_idx_obs_nc(:)=1
+            !if ((maxval(z_idx_obs_nc).NE.1) .OR. (minval(z_idx_obs_nc).NE.1)) then
+            !   write(*,*) 'For crns average mode parflow obs layer iz must be 1'
+            !   stop
+            !endif 
+        endif
+        !end hcp
         if (screen > 2) then
             print *, "TSMP-PDAF mype(w)=", mype_world, ": z_idx_obs_nc=", z_idx_obs_nc
         end if
@@ -501,7 +503,6 @@ subroutine clean_obs_nc()
     !if(allocated(x_idx_obs_nc))deallocate(x_idx_obs_nc)
     !if(allocated(y_idx_obs_nc))deallocate(y_idx_obs_nc)
     !if(allocated(z_idx_obs_nc))deallocate(z_idx_obs_nc)
-    if(allocated(depth_obs))deallocate(depth_obs)
     !kuw: clean clm observations
     if(allocated(clmobs_lon))deallocate(clmobs_lon)
     if(allocated(clmobs_lat))deallocate(clmobs_lat)

bldsva/intf_DA/pdaf/framework/mod_tsmp.F90

@@ -34,10 +34,12 @@ module mod_tsmp
     integer(c_int), bind(c)  :: tag_model_clm = 0
     integer(c_int), bind(c)  :: tag_model_parflow = 1
     integer(c_int), bind(c)  :: tag_model_cosmo   = 2
+    integer(c_int), bind(c)  :: crns_flag
     type(c_ptr), bind(c)     :: pf_statevec
     type(c_ptr), bind(c)     :: xcoord, ycoord, zcoord
     real(c_double), pointer  :: xcoord_fortran(:), ycoord_fortran(:), zcoord_fortran(:)
     real(c_double), pointer  :: pf_statevec_fortran(:)
+    real(c_double), bind(c)  :: da_crns_depth_tol
     type(c_ptr), bind(c)     :: idx_map_subvec2state
     integer(c_int), pointer  :: idx_map_subvec2state_fortran(:)
     type(c_ptr), bind(c)     :: soilay

bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90

@@ -52,18 +52,18 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
         ONLY: obs_index_p, &
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
-        depth_obs_p, & !obs_p
         sc_p, &
 #endif
 #endif
         obs_interp_indices_p, &
         obs_interp_weights_p
-   USE mod_read_obs, ONLY: crns_flag !clm_obs
    use mod_tsmp, &
        only: obs_interp_switch, &
        soilay, &
        soilay_fortran, &
-       nz_glob
+       nz_glob, &
+       da_crns_depth_tol, &
+       crns_flag
 !      tcycle
 
    USE, INTRINSIC :: iso_c_binding
@@ -96,7 +96,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! soide=(/0.d0,  0.02d0,  0.05d0,  0.1d0,  0.17d0, 0.3d0,  0.5d0, &
 !                0.8d0,   1.3d0,   2.d0,  3.d0, 5.d0,  12.d0/) !soil depth
 
-real(8) :: tot, avesm
+real(8) :: tot, avesm, avesm_temp, Dp
 integer :: nsc
 ! end of hcp 
 
@@ -133,33 +133,54 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
  if (crns_flag.EQ.1) then
-
+    !Schroen et al HESS 2017 modelled CRNS averaging
     lpointobs = .false.
-
      call C_F_POINTER(soilay,soilay_fortran,[nz_glob])
      Allocate(soide(0:nz_glob))
      soide(0)=0.d0
      do i=1,nz_glob
        soide(i)=soide(i-1)+soilay_fortran(nz_glob-i+1) 
      enddo
      do i = 1, dim_obs_p
-       nsc= size(sc_p(i)%scol_obs_in(:))
+
+       !Initial average soil moisture for 1st iteration
        avesm=0.d0
-       do j=1, nsc-1
-           avesm=avesm+(1.d0-0.5d0*(soide(j)+soide(j-1))/depth_obs_p(i))*(soide(j)-soide(j-1)) &
-                 *state_p(sc_p(i)%scol_obs_in(j))/depth_obs_p(i)
-       enddo
-       avesm=avesm+(1.d0-0.5d0*(depth_obs_p(i)+soide(nsc-1))/depth_obs_p(i))*(depth_obs_p(i)-soide(nsc-1)) &
-             *state_p(sc_p(i)%scol_obs_in(nsc))/depth_obs_p(i)
-       tot=0.d0
-       do j=1, nsc-1
-           tot=tot+(1.d0-0.5d0*(soide(j)+soide(j-1))/depth_obs_p(i))*(soide(j)-soide(j-1)) &
-               /depth_obs_p(i)
+       do j=1,nz_glob
+            avesm=avesm+(soide(j)-soide(j-1))*state_p(sc_p(j,i))/soide(nz_glob)
        enddo
-       tot=tot+(1.d0-0.5d0*(depth_obs_p(i)+soide(nsc-1))/depth_obs_p(i))*(depth_obs_p(i)-soide(nsc-1)) &
-          /depth_obs_p(i)
+       avesm_temp=0.d0
+
+       !iteration
+       do while (abs(avesm-avesm_temp)/avesm .GE. da_crns_depth_tol)
+          !Averaging, conventional profile, Schroen et al HESS 2017 Eq. (3)
+          avesm_temp=avesm
+          Dp=0.058d0/(avesm+0.0829d0)
+
+          !Sum weight*soil_moisture
+          avesm=0.d0; nsc=nz_glob
+          do j=1,nz_glob
+             if ((soide(j-1).LT.Dp).AND.(Dp.LE.soide(j))) then
+               nsc=j
+             endif
+          enddo
+          do j=1, nsc-1
+              avesm=avesm+(1.d0-0.5d0*(soide(j)+soide(j-1))/Dp)*(soide(j)-soide(j-1)) &
+                    *state_p(sc_p(j,i))/Dp
+          enddo
+          avesm=avesm+(1.d0-0.5d0*(Dp+soide(nsc-1))/Dp)*(Dp-soide(nsc-1)) &
+             *state_p(sc_p(nsc,i))/Dp
+
+          !Sum weight
+          tot=0.d0
+          do j=1, nsc-1
+              tot =   tot+(1.d0-0.5d0*(soide(j)+soide(j-1))/Dp)*(soide(j)-soide(j-1)) &
+                                                    /Dp
+          enddo
+          tot  =  tot+(1.d0-0.5d0*(Dp+soide(nsc-1))/Dp)*(Dp-soide(nsc-1)) &
+                                               /Dp
 
-       avesm=avesm/tot
+          avesm=avesm/tot
+       enddo
        m_state_p(i)=avesm
      enddo
      deallocate(soide)

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

@@ -102,6 +102,7 @@ GLOBAL int pf_aniso_use_parflow;
 GLOBAL int is_dampfac_state_time_dependent;
 GLOBAL int is_dampfac_param_time_dependent;
 GLOBAL int pf_dampswitch_sm;
+GLOBAL int crns_flag;
 extern int model;
 extern int mype_model;
 extern int npes_model;
@@ -128,3 +129,4 @@ GLOBAL double pf_dampfac_param;
 GLOBAL double pf_dampfac_state;
 GLOBAL double dampfac_state_time_dependent;
 GLOBAL double dampfac_param_time_dependent;
+GLOBAL double da_crns_depth_tol;

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

@@ -94,6 +94,8 @@ void read_enkfpar(char *parname)
   screen_wrapper        = iniparser_getint(pardict,"DA:screen_wrapper",1);
   point_obs             = iniparser_getint(pardict,"DA:point_obs",1);
   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);
   total_steps = (int) (t_sim/da_interval);
   tstartcycle = (int) (t_start/da_interval);
 

doc/content/setup_tsmp/input_enkfpf.md

@@ -65,6 +65,8 @@ startreal =
 da_interval       =
 stat_dumpoffset   =
 point_obs =
+crns_flag =
+da_crns_depth_tol =
 ```
 
 In the following the individual entries of `enkfpf.par` are described:
@@ -559,6 +561,29 @@ Effect of `obs_interp_switch=1`:
   of the simulated measurements from the value at the surrounding grid
   points in the x-y-plane.
 
+### 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: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)
+
+Only used if CRNS-observations are turned on by adding parameter
+`depth` in observation input.
+
+Relevant to assimilating PF averaged soil moisture data only for the
+modelled CRNS data in Schroen etal in HESS 2017. 
+
+This is a convergence criteria (between 0 and 1) for the averaging
+skin depth (compare the weighting procedure described in Schrön et al,
+2017, section 2.3).
+
 ## Parameter Summary ##
 
  | section   | parameter               | default value |
@@ -603,6 +628,8 @@ Effect of `obs_interp_switch=1`:
  |           | `screen_wrapper`        | 1             |
  |           | `point_obs`             | 1             |
  |           | `obs_interp_switch`     | 0             |
+ |           | `crns_flag`             | 1             |
+ |           | `da_crns_depth_tol`     | 0.01          |
 
 Default values for parameter file `enkfpf.par`.