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 fafa5a261..235eb6f96 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 @@ -959,7 +959,16 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if #endif else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) +#ifdef CLMFIVE + if(clmstatevec_only_active.eq.1) then + obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) + else + +#endif + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) +#ifdef CLMFIVE + end if +#endif end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 59eac4b35..da359d14f 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -41,6 +41,7 @@ module enkf_clm_mod integer :: clm_begc,clm_endc integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) + real(r8),allocatable :: clm_statevec_orig(:) integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model @@ -54,6 +55,7 @@ module enkf_clm_mod #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol + integer(c_int),bind(C,name="clmstatevec_colmean") :: clmstatevec_colmean integer(c_int),bind(C,name="clmstatevec_only_active") :: clmstatevec_only_active integer(c_int),bind(C,name="clmstatevec_max_layer") :: clmstatevec_max_layer integer(c_int),bind(C,name="clmt_printensemble") :: clmt_printensemble @@ -98,6 +100,7 @@ subroutine define_clm_statevec(mype) integer :: jj integer :: c integer :: g + integer :: cg integer :: cc integer :: cccheck @@ -122,140 +125,160 @@ subroutine define_clm_statevec(mype) clm_begp = begp clm_endp = endp + ! Soil Moisture DA: State vector index arrays if(clmupdate_swc.eq.1) then + + ! 1) COL/GRC: CLM->PDAF + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + do i=1,nlevsoi + do c=clm_begc,clm_endc + ! Default: inactive + state_clm2pdaf_p = ispval + end do + end do + + ! All column variables in state vector if(clmstatevec_allcol.eq.1) then + ! Only hydrologically active columns if(clmstatevec_only_active .eq. 1) then - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) - cc = 0 do i=1,nlevsoi - do c=clm_begc,clm_endc - ! Only take into account layers above input maximum layer - if(i<=clmstatevec_max_layer) then + ! Only take into account layers above input maximum layer + if(i<=clmstatevec_max_layer) then + + do c=clm_begc,clm_endc ! Only take into account hydrologically active columns ! and layers above bedrock if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then cc = cc + 1 state_clm2pdaf_p(c,i) = cc - else - state_clm2pdaf_p(c,i) = ispval end if - else - state_clm2pdaf_p(c,i) = ispval - end if - end do - end do + end do - ! Set `clm_varsize`, even though it is currently not used - ! for `clmupdate_swc.eq.1` - clm_varsize = cc - clm_statevecsize = cc - - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) - - cc = 0 - - do i=1,nlevsoi - do c=clm_begc,clm_endc - ! Only take into account layers above input maximum layer - if(i<=clmstatevec_max_layer) then - ! Only take into account hydrologically active columns - ! and layers above bedrock - if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then - cc = cc + 1 - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i - end if - end if - end do + end if end do + ! All column variables in state vector simplifying the indexing else - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) - do i=1,nlevsoi do c=clm_begc,clm_endc state_clm2pdaf_p(c,i) = (c - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) end do end do - ! #cols values per grid-cell - clm_varsize = (endc-begc+1) * nlevsoi - clm_statevecsize = (endc-begc+1) * nlevsoi + end if - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Gridcell values or averages in state vector + else + + ! Only hydrologically active columns + if(clmstatevec_only_active.eq.1) then cc = 0 + do i=1,nlevsoi + ! Only layers above max_layer + if(i<=clmstatevec_max_layer) then + + do g=clm_begg,clm_endg + + newgridcell = .true. + + do c=clm_begc,clm_endc + if(col%gridcell(c) == g) then + ! All hydrologically active columns above + ! bedrock in a gridcell point to the state + ! vector index of the gridcell + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + if(newgridcell) then + ! Update the index if first col found for grc, + ! otherwise reproduce previous index + cc = cc + 1 + newgridcell = .false. + end if + state_clm2pdaf_p(c,i) = cc + end if + end if + end do + + end do + end if + end do + else do i=1,nlevsoi do c=clm_begc,clm_endc - cc = cc + 1 - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i + ! All columns in a gridcell are assigned the updated + ! gridcell-SWC + state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) end do end do - end if - else - - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + end if - do i=1,nlevsoi - do c=clm_begc,clm_endc - ! All columns in a gridcell are assigned the updated - ! gridcell-SWC - state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) - end do - end do + ! 2) COL/GRC: STATEVECSIZE + if(clmstatevec_only_active.eq.1) then + ! Use iterator cc for setting state vector size. + ! + ! Set `clm_varsize`, even though it is currently not used + ! for `clmupdate_swc.eq.1` + clm_varsize = cc + clm_statevecsize = cc + else + if(clmstatevec_allcol.eq.1) then + ! #cols * #levels + clm_varsize = (endc-begc+1) * nlevsoi + clm_statevecsize = (endc-begc+1) * nlevsoi + else + ! #grcs * #levels + clm_varsize = (endg-begg+1) * nlevsoi + clm_statevecsize = (endg-begg+1) * nlevsoi + end if + end if - ! One value per grid-cell - clm_varsize = (endg-begg+1) * nlevsoi - clm_statevecsize = (endg-begg+1) * nlevsoi + ! 3) COL/GRC: PDAF->CLM + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Defaults + do cc=1,clm_statevecsize + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do - cc = 0 + do cc=1,clm_statevecsize - do i=1,nlevsoi - do j=clm_begg,clm_endg - - ! SWC from the first column of each gridcell - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - cc = cc + 1 - ! Possibliy: Add state_pdaf2clm_g_p - state_pdaf2clm_c_p(cc) = jj - state_pdaf2clm_j_p(cc) = i - end if - end if - end do + lay: do i=1,nlevsoi + do c=clm_begc,clm_endc + if (state_clm2pdaf_p(c,i) == cc) then + ! Set column index and then exit loop + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + exit lay + end if end do - end do - + end do lay +#ifdef PDAF_DEBUG + ! Check that all state vectors have been assigned c, i + if(state_pdaf2clm_c_p(cc) == ispval) then + write(*,*) 'cc: ', cc + error stop "state_pdaf2clm_c_p not set at cc" + end if + if(state_pdaf2clm_j_p(cc) == ispval) then + write(*,*) 'cc: ', cc + error stop "state_pdaf2clm_j_p not set at cc" + end if +#endif + end do - end if endif if(clmupdate_swc.eq.2) then @@ -288,6 +311,14 @@ subroutine define_clm_statevec(mype) allocate(clm_statevec(clm_statevecsize)) end if + ! Allocate statevector-duplicate for saving original column mean + ! values used in computing increments during updating the state + ! vector in column-mean-mode. + IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) + if (clmupdate_swc.ne.0 .and. clmstatevec_colmean.ne.0) then + allocate(clm_statevec_orig(clm_statevecsize)) + end if + !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp if ((clmupdate_T.ne.0)) then !hcp @@ -322,7 +353,8 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,cc=0,offset=0 + integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -352,9 +384,51 @@ subroutine set_clm_statevec(tstartcycle, mype) if(clmupdate_swc.ne.0) then ! write swc values to state vector - do cc = 1, clm_statevecsize - clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) - end do + if (clmstatevec_colmean.eq.1) then + + do cc = 1, clm_statevecsize + + clm_statevec(cc) = 0.0 + n_c = 0 + + ! Get gridcell and layer + g = col%gridcell(state_pdaf2clm_c_p(cc)) + j = state_pdaf2clm_j_p(cc) + + ! Loop over all columns + do c=clm_begc,clm_endc + ! Select columns in gridcell g + if(col%gridcell(c).eq.g) then + ! Select hydrologically active columns + if(col%hydrologically_active(c)) then + ! Add active column to swc-sum + clm_statevec(cc) = clm_statevec(cc) + swc(c,j) + n_c = n_c + 1 + end if + end if + end do + + if(n_c == 0) then + write(*,*) "WARNING: Gridcell g=", g + write(*,*) "WARNING: Layer j=", j + write(*,*) "Grid cell g at layer j without hydrologically active column! Setting SWC as in gridcell mode." + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + else + ! Normalize sum to average + clm_statevec(cc) = clm_statevec(cc) / real(n_c, r8) + end if + + ! Save prior column mean state vector for computing + ! increment in updating the state vector + clm_statevec_orig(cc) = clm_statevec(cc) + + end do + + else + do cc = 1, clm_statevecsize + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + end if endif !hcp LAI @@ -534,7 +608,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if - swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + if (clmstatevec_colmean.eq.1) then + ! Update SWC column value with the increment-factor + ! of the state vector update (state vector updates + ! are means of cols in grc) + swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) + else + ! Update SWC with updated state vector + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + end if if(swc_update.le.watmin_check) then swc(j,i) = watmin_set diff --git a/bldsva/intf_DA/pdaf/model/common/enkf.h b/bldsva/intf_DA/pdaf/model/common/enkf.h index fba54ce56..ce5dd918f 100755 --- a/bldsva/intf_DA/pdaf/model/common/enkf.h +++ b/bldsva/intf_DA/pdaf/model/common/enkf.h @@ -90,6 +90,7 @@ GLOBAL int clmupdate_texture; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; +GLOBAL int clmstatevec_colmean; GLOBAL int clmstatevec_only_active; GLOBAL int clmstatevec_max_layer; GLOBAL int clmt_printensemble; diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c index dc8046a55..e391f916c 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c @@ -81,6 +81,7 @@ void read_enkfpar(char *parname) clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); + clmstatevec_colmean = iniparser_getint(pardict,"CLM:statevec_colmean",0); clmstatevec_only_active = iniparser_getint(pardict,"CLM:statevec_only_active",0); clmstatevec_max_layer = iniparser_getint(pardict,"CLM:statevec_max_layer",25); clmt_printensemble = iniparser_getint(pardict,"CLM:t_printensemble",-1); @@ -142,6 +143,14 @@ void read_enkfpar(char *parname) exit(1); } + /* Check: Consistency of `statevec_allcol` and `statevec_colmean` */ + if (clmstatevec_allcol != 0 && clmstatevec_colmean != 0){ + printf("clmstatevec_allcol=%d\n", clmstatevec_allcol); + printf("clmstatevec_colmean=%d\n", clmstatevec_colmean); + printf("Error: Either clmstatevec_allcol or clmstatevec_colmean must turned off, i.e. equal to 0.\n"); + exit(1); + } + /* Assign model specifier (0=clm, 1=parflow, 2=cosmo) */ /* Order: ParFlow, CLM, COSMO */ if (mype_model < nprocpf) { diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 59eac4b35..da359d14f 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -41,6 +41,7 @@ module enkf_clm_mod integer :: clm_begc,clm_endc integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) + real(r8),allocatable :: clm_statevec_orig(:) integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model @@ -54,6 +55,7 @@ module enkf_clm_mod #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol + integer(c_int),bind(C,name="clmstatevec_colmean") :: clmstatevec_colmean integer(c_int),bind(C,name="clmstatevec_only_active") :: clmstatevec_only_active integer(c_int),bind(C,name="clmstatevec_max_layer") :: clmstatevec_max_layer integer(c_int),bind(C,name="clmt_printensemble") :: clmt_printensemble @@ -98,6 +100,7 @@ subroutine define_clm_statevec(mype) integer :: jj integer :: c integer :: g + integer :: cg integer :: cc integer :: cccheck @@ -122,140 +125,160 @@ subroutine define_clm_statevec(mype) clm_begp = begp clm_endp = endp + ! Soil Moisture DA: State vector index arrays if(clmupdate_swc.eq.1) then + + ! 1) COL/GRC: CLM->PDAF + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + do i=1,nlevsoi + do c=clm_begc,clm_endc + ! Default: inactive + state_clm2pdaf_p = ispval + end do + end do + + ! All column variables in state vector if(clmstatevec_allcol.eq.1) then + ! Only hydrologically active columns if(clmstatevec_only_active .eq. 1) then - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) - cc = 0 do i=1,nlevsoi - do c=clm_begc,clm_endc - ! Only take into account layers above input maximum layer - if(i<=clmstatevec_max_layer) then + ! Only take into account layers above input maximum layer + if(i<=clmstatevec_max_layer) then + + do c=clm_begc,clm_endc ! Only take into account hydrologically active columns ! and layers above bedrock if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then cc = cc + 1 state_clm2pdaf_p(c,i) = cc - else - state_clm2pdaf_p(c,i) = ispval end if - else - state_clm2pdaf_p(c,i) = ispval - end if - end do - end do + end do - ! Set `clm_varsize`, even though it is currently not used - ! for `clmupdate_swc.eq.1` - clm_varsize = cc - clm_statevecsize = cc - - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) - - cc = 0 - - do i=1,nlevsoi - do c=clm_begc,clm_endc - ! Only take into account layers above input maximum layer - if(i<=clmstatevec_max_layer) then - ! Only take into account hydrologically active columns - ! and layers above bedrock - if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then - cc = cc + 1 - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i - end if - end if - end do + end if end do + ! All column variables in state vector simplifying the indexing else - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) - do i=1,nlevsoi do c=clm_begc,clm_endc state_clm2pdaf_p(c,i) = (c - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) end do end do - ! #cols values per grid-cell - clm_varsize = (endc-begc+1) * nlevsoi - clm_statevecsize = (endc-begc+1) * nlevsoi + end if - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Gridcell values or averages in state vector + else + + ! Only hydrologically active columns + if(clmstatevec_only_active.eq.1) then cc = 0 + do i=1,nlevsoi + ! Only layers above max_layer + if(i<=clmstatevec_max_layer) then + + do g=clm_begg,clm_endg + + newgridcell = .true. + + do c=clm_begc,clm_endc + if(col%gridcell(c) == g) then + ! All hydrologically active columns above + ! bedrock in a gridcell point to the state + ! vector index of the gridcell + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + if(newgridcell) then + ! Update the index if first col found for grc, + ! otherwise reproduce previous index + cc = cc + 1 + newgridcell = .false. + end if + state_clm2pdaf_p(c,i) = cc + end if + end if + end do + + end do + end if + end do + else do i=1,nlevsoi do c=clm_begc,clm_endc - cc = cc + 1 - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i + ! All columns in a gridcell are assigned the updated + ! gridcell-SWC + state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) end do end do - end if - else - - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + end if - do i=1,nlevsoi - do c=clm_begc,clm_endc - ! All columns in a gridcell are assigned the updated - ! gridcell-SWC - state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) - end do - end do + ! 2) COL/GRC: STATEVECSIZE + if(clmstatevec_only_active.eq.1) then + ! Use iterator cc for setting state vector size. + ! + ! Set `clm_varsize`, even though it is currently not used + ! for `clmupdate_swc.eq.1` + clm_varsize = cc + clm_statevecsize = cc + else + if(clmstatevec_allcol.eq.1) then + ! #cols * #levels + clm_varsize = (endc-begc+1) * nlevsoi + clm_statevecsize = (endc-begc+1) * nlevsoi + else + ! #grcs * #levels + clm_varsize = (endg-begg+1) * nlevsoi + clm_statevecsize = (endg-begg+1) * nlevsoi + end if + end if - ! One value per grid-cell - clm_varsize = (endg-begg+1) * nlevsoi - clm_statevecsize = (endg-begg+1) * nlevsoi + ! 3) COL/GRC: PDAF->CLM + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Defaults + do cc=1,clm_statevecsize + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do - cc = 0 + do cc=1,clm_statevecsize - do i=1,nlevsoi - do j=clm_begg,clm_endg - - ! SWC from the first column of each gridcell - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - cc = cc + 1 - ! Possibliy: Add state_pdaf2clm_g_p - state_pdaf2clm_c_p(cc) = jj - state_pdaf2clm_j_p(cc) = i - end if - end if - end do + lay: do i=1,nlevsoi + do c=clm_begc,clm_endc + if (state_clm2pdaf_p(c,i) == cc) then + ! Set column index and then exit loop + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + exit lay + end if end do - end do - + end do lay +#ifdef PDAF_DEBUG + ! Check that all state vectors have been assigned c, i + if(state_pdaf2clm_c_p(cc) == ispval) then + write(*,*) 'cc: ', cc + error stop "state_pdaf2clm_c_p not set at cc" + end if + if(state_pdaf2clm_j_p(cc) == ispval) then + write(*,*) 'cc: ', cc + error stop "state_pdaf2clm_j_p not set at cc" + end if +#endif + end do - end if endif if(clmupdate_swc.eq.2) then @@ -288,6 +311,14 @@ subroutine define_clm_statevec(mype) allocate(clm_statevec(clm_statevecsize)) end if + ! Allocate statevector-duplicate for saving original column mean + ! values used in computing increments during updating the state + ! vector in column-mean-mode. + IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) + if (clmupdate_swc.ne.0 .and. clmstatevec_colmean.ne.0) then + allocate(clm_statevec_orig(clm_statevecsize)) + end if + !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp if ((clmupdate_T.ne.0)) then !hcp @@ -322,7 +353,8 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,cc=0,offset=0 + integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -352,9 +384,51 @@ subroutine set_clm_statevec(tstartcycle, mype) if(clmupdate_swc.ne.0) then ! write swc values to state vector - do cc = 1, clm_statevecsize - clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) - end do + if (clmstatevec_colmean.eq.1) then + + do cc = 1, clm_statevecsize + + clm_statevec(cc) = 0.0 + n_c = 0 + + ! Get gridcell and layer + g = col%gridcell(state_pdaf2clm_c_p(cc)) + j = state_pdaf2clm_j_p(cc) + + ! Loop over all columns + do c=clm_begc,clm_endc + ! Select columns in gridcell g + if(col%gridcell(c).eq.g) then + ! Select hydrologically active columns + if(col%hydrologically_active(c)) then + ! Add active column to swc-sum + clm_statevec(cc) = clm_statevec(cc) + swc(c,j) + n_c = n_c + 1 + end if + end if + end do + + if(n_c == 0) then + write(*,*) "WARNING: Gridcell g=", g + write(*,*) "WARNING: Layer j=", j + write(*,*) "Grid cell g at layer j without hydrologically active column! Setting SWC as in gridcell mode." + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + else + ! Normalize sum to average + clm_statevec(cc) = clm_statevec(cc) / real(n_c, r8) + end if + + ! Save prior column mean state vector for computing + ! increment in updating the state vector + clm_statevec_orig(cc) = clm_statevec(cc) + + end do + + else + do cc = 1, clm_statevecsize + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + end if endif !hcp LAI @@ -534,7 +608,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if - swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + if (clmstatevec_colmean.eq.1) then + ! Update SWC column value with the increment-factor + ! of the state vector update (state vector updates + ! are means of cols in grc) + swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) + else + ! Update SWC with updated state vector + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + end if if(swc_update.le.watmin_check) then swc(j,i) = watmin_set diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index 4d1082a1a..987b15dea 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -56,6 +56,7 @@ update_T = print_swc = print_et = statevec_allcol = +statevec_colmean = statevec_only_active = statevec_max_layer = t_printensemble = @@ -495,7 +496,7 @@ further information, see source code. Default: `0`. ### CLM:statevec_allcol ### `CLM:statevec_allcol`: (integer) Switch for using all SWC columns of a -CLM5 gridcell in the state vector. +eCLM gridcell in the state vector. If `0` (default): Only one SWC value per grid cell is saved in the state vector. @@ -503,15 +504,31 @@ state vector. If `1`: `#columns` SWC values per grid cell are saved in the state vector. -### CLM:statevec_only_active ### +Remark: `CLM:statevec_allcol` should not be switched on together with +`CLM:statevec_colmean`. -**Not yet in main branch** +### CLM:statevec_colmean ### + +`CLM:statevec_colmean`: (integer) Switch for using gridcell averages +of SWC (which is a column-variable) in the state vector. + +If `0` (default): Only the SWC value from the **first** column is +saved per grid cell in the state vector. + +If `1`: The gridcell mean over all (possibly hydrologically active) +SWC values is saved per grid cell in the state vector. + +Remark: `CLM:statevec_colmean` should not be switched on together with +`CLM:statevec_allcol`. + +### CLM:statevec_only_active ### `CLM:statevec_only_active`: (integer) Switch for using in the state vector only (1) hydrologically active columns of a CLM5 gridcell and (2) only layers until bedrock. -Only used, when `CLM:statevec_allcol` is switched on. +Only used for `CLM:update_swc`. Therein, only used, when either +`CLM:statevec_colmean` or `CLM:statevec_allcol` is used. If `0` (default): Use all columns and all layers. @@ -786,6 +803,7 @@ Default: 0, output turned off. | | `print_swc` | 0 | | | `print_et` | 0 | | | `statevec_allcol` | 0 | + | | `statevec_colmean` | 0 | | | `statevec_only_active` | 0 | | | `statevec_max_layer` | 25 | | | `t_printensemble` | -1 |