diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 23864426..759ad83d 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/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 diff --git a/bldsva/intf_DA/pdaf/framework/mod_assimilation.F90 b/bldsva/intf_DA/pdaf/framework/mod_assimilation.F90 index 98710e9d..64f4730c 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_assimilation.F90 +++ b/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 diff --git a/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 b/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 index b2afcb22..1e9ffda5 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 +++ b/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) diff --git a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 index c3b6b158..f4ba8f4b 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 +++ b/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 diff --git a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index e3d56a96..53182e7e 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/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,9 +133,8 @@ 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 @@ -143,23 +142,45 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) 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) diff --git a/bldsva/intf_DA/pdaf/model/common/enkf.h b/bldsva/intf_DA/pdaf/model/common/enkf.h index 2dbccb4b..72810ba6 100755 --- a/bldsva/intf_DA/pdaf/model/common/enkf.h +++ b/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; diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c index f0c35bed..eae2a82f 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/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); diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index ea36b667..5e7c340b 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/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`.