From 8f3afb750d494058b02524b4c2b66ce18764ab3c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 25 Apr 2025 15:15:40 +0200 Subject: [PATCH 01/10] first draft: gridcell mean over columns --- bldsva/intf_DA/pdaf/model/common/enkf.h | 1 + .../intf_DA/pdaf/model/common/read_enkfpar.c | 9 ++++ .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 46 +++++++++++++++++-- 3 files changed, 52 insertions(+), 4 deletions(-) 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 45b849411..311969658 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 @@ -54,6 +54,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 @@ -322,7 +323,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 + real :: 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 +354,45 @@ 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.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 + ! Add columns in gridcell g + if(col%gridcell(c).eq.g) then + ! Only hydrologically active columns above bedrock + if(col%hydrologically_active(c) .and. j<=col%nbedrock(c)) then + clm_statevec(cc) = clm_statevec(cc) + swc(c,j) + n_c = n_c + 1.0 + end if + end if + end do + + if(n_c == 0.0) then + write(*,*) "WARNING: Gridcell g=", g + write(*,*) "Grid cell g without hydrologically active column! Setting SWC as before from first column." + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + else + ! Compute final average + clm_statevec(cc) = clm_statevec(cc) / n_c + end if + + 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 From 7808a83f768127e7d1854a5096e57dbbfd30289c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 6 May 2025 17:01:26 +0200 Subject: [PATCH 02/10] colmean update - introduce hydrologically active grid cells only - save prior column means in `clm_statevec_orig` - update column values with increment-factor from mean --- .../pdaf/framework/init_dim_obs_pdaf.F90 | 11 +- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 137 ++++++++++++++---- 2 files changed, 118 insertions(+), 30 deletions(-) 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/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 311969658..b8969612b 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 @@ -99,6 +100,7 @@ subroutine define_clm_statevec(mype) integer :: jj integer :: c integer :: g + integer :: cg integer :: cc integer :: cccheck @@ -215,17 +217,56 @@ subroutine define_clm_statevec(mype) 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 - ! 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) + if(clmstatevec_only_active.eq.1) then + cc = 0 + do i=1,nlevsoi + do g=clm_begg,clm_endg + + newgridcell = .true. + + do c=clm_begc,clm_endc + ! All hydrologically active columns above max layer + ! and bedrock in a gridcell are assigned the updated + ! gridcell-SWC + if(i<=clmstatevec_max_layer) then + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + if(col%gridcell(c) == g) then + if(newgridcell) then + ! Update the index if first col for grc, + ! otherwise reproduce previous index + cc = cc + 1 + newgridcell = .false. + end if + state_clm2pdaf_p(c,i) = cc + end if + else + state_clm2pdaf_p(c,i) = ispval + end if + else + state_clm2pdaf_p(c,i) = ispval + end if + end do + + end do end do - end do + else + 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 + end if ! One value per grid-cell - clm_varsize = (endg-begg+1) * nlevsoi - clm_statevecsize = (endg-begg+1) * nlevsoi + if(clmstatevec_only_active.eq.1) then + clm_varsize = cc + clm_statevecsize = cc + else + clm_varsize = (endg-begg+1) * nlevsoi + clm_statevecsize = (endg-begg+1) * nlevsoi + end if IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) allocate(state_pdaf2clm_c_p(clm_statevecsize)) @@ -234,27 +275,50 @@ subroutine define_clm_statevec(mype) cc = 0 - 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 + if(clmstatevec_only_active .eq. 1) then + do i=1,nlevsoi + do g=clm_begg,clm_endg + + if(i<=clmstatevec_max_layer) then + ! SWC from the first column of each gridcell + newgridcell = .true. + do c=clm_begc,clm_endc + cg = col%gridcell(c) + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + if (cg .eq. g) then + if (newgridcell) then + newgridcell = .false. + cc = cc + 1 + ! Possibly: Add state_pdaf2clm_g_p + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end if + end if + end if + end do end if end do end do - end do - - + else + do i=1,nlevsoi + do g=clm_begg,clm_endg + ! SWC from the first column of each gridcell + newgridcell = .true. + do c=clm_begc,clm_endc + cg = col%gridcell(c) + if (cg .eq. g) then + if (newgridcell) then + newgridcell = .false. + cc = cc + 1 + ! Possibly: Add state_pdaf2clm_g_p + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end if + end if + end do + end do + end do + end if end if endif @@ -289,6 +353,11 @@ subroutine define_clm_statevec(mype) allocate(clm_statevec(clm_statevecsize)) end if + 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 @@ -369,8 +438,9 @@ subroutine set_clm_statevec(tstartcycle, mype) do c=clm_begc,clm_endc ! Add columns in gridcell g if(col%gridcell(c).eq.g) then - ! Only hydrologically active columns above bedrock - if(col%hydrologically_active(c) .and. j<=col%nbedrock(c)) then + ! If there are hydrologically inactive columns in the + ! same gridcell with active columns + if(col%hydrologically_active(c)) then clm_statevec(cc) = clm_statevec(cc) + swc(c,j) n_c = n_c + 1.0 end if @@ -379,13 +449,17 @@ subroutine set_clm_statevec(tstartcycle, mype) if(n_c == 0.0) then write(*,*) "WARNING: Gridcell g=", g - write(*,*) "Grid cell g without hydrologically active column! Setting SWC as before from first column." + write(*,*) "WARNING: Layer j=", j + write(*,*) "Grid cell g at layer j without hydrologically active column! Setting SWC as before from first column." clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) else ! Compute final average clm_statevec(cc) = clm_statevec(cc) / n_c end if + ! Save prior state vector for increment + clm_statevec_orig(cc) = clm_statevec(cc) + end do else @@ -574,6 +648,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + if (clmstatevec_colmean.eq.1) then + ! Update with the factor of the update for the column mean + swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) + end if + if(swc_update.le.watmin_check) then swc(j,i) = watmin_set else if(swc_update.ge.watsat(j,i)) then From 8ce706bc052fc649767b2a8211dc0660e1b56a61 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 7 May 2025 13:41:10 +0200 Subject: [PATCH 03/10] if-condition-fix --- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 b8969612b..c6b9fa1c7 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 @@ -354,7 +354,7 @@ subroutine define_clm_statevec(mype) end if IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) - if (clmupdate_swc.ne.0 .and clmstatevec_colmean.ne.0) then + if (clmupdate_swc.ne.0 .and. clmstatevec_colmean.ne.0) then allocate(clm_statevec_orig(clm_statevecsize)) end if From bf9d69afd880f1d67bcbf48549ed92aa2a45b100 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 23 May 2025 07:56:21 +0200 Subject: [PATCH 04/10] doc: input documentation for `CLM:statevec_colmean` --- doc/content/setup_tsmp/input_enkfpf.md | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) 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 | From f5ddad02d799eabf42d0d0b6c502fc5fc78b9cc0 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 10:32:14 +0200 Subject: [PATCH 05/10] define_clm_statevec: re-arrange statevector definition --- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 222 +++++++----------- 1 file changed, 90 insertions(+), 132 deletions(-) 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 748e7c829..bb7ef8db5 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 @@ -125,129 +125,89 @@ 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 - - ! 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 + end do - 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 - - 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 - cc = cc + 1 - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i - end do - end do - end if + ! Gridcell values or averages in state vector else - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) - + ! Only hydrologically active columns if(clmstatevec_only_active.eq.1) then + cc = 0 + do i=1,nlevsoi - do g=clm_begg,clm_endg + ! Only layers above max_layer + if(i<=clmstatevec_max_layer) then - newgridcell = .true. + do g=clm_begg,clm_endg - do c=clm_begc,clm_endc - ! All hydrologically active columns above max layer - ! and bedrock in a gridcell are assigned the updated - ! gridcell-SWC - if(i<=clmstatevec_max_layer) then - if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then - if(col%gridcell(c) == g) then + 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 for grc, + ! 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 - else - state_clm2pdaf_p(c,i) = ispval end if - else - state_clm2pdaf_p(c,i) = ispval - end if - end do + end do - end do + end do + end if end do else do i=1,nlevsoi @@ -259,68 +219,66 @@ subroutine define_clm_statevec(mype) end do end if - ! One value per grid-cell - if(clmstatevec_only_active.eq.1) then - clm_varsize = cc - clm_statevecsize = cc + end if + + ! 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 - 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)) + ! 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)) - cc = 0 + ! Defaults + do cc=1,clm_statevecsize + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do - if(clmstatevec_only_active .eq. 1) then - do i=1,nlevsoi - do g=clm_begg,clm_endg + do cc=1,clm_statevecsize - if(i<=clmstatevec_max_layer) then - ! SWC from the first column of each gridcell - newgridcell = .true. - do c=clm_begc,clm_endc - cg = col%gridcell(c) - if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then - if (cg .eq. g) then - if (newgridcell) then - newgridcell = .false. - cc = cc + 1 - ! Possibly: Add state_pdaf2clm_g_p - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i - end if - end if - end if - end do - end if - end do - end do - else - do i=1,nlevsoi - do g=clm_begg,clm_endg - ! SWC from the first column of each gridcell - newgridcell = .true. - do c=clm_begc,clm_endc - cg = col%gridcell(c) - if (cg .eq. g) then - if (newgridcell) then - newgridcell = .false. - cc = cc + 1 - ! Possibly: Add state_pdaf2clm_g_p - state_pdaf2clm_c_p(cc) = c - state_pdaf2clm_j_p(cc) = i - end if - end if - end do - end do - end do + lay: do i=1,nlevsoi + col: 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 col + 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 From 7ce2aa8e32c19788e934b26dd3b5def568499b34 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 10:47:42 +0200 Subject: [PATCH 06/10] doc updates --- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) 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 bb7ef8db5..c7f932da9 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 @@ -311,6 +311,9 @@ 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)) @@ -394,11 +397,11 @@ subroutine set_clm_statevec(tstartcycle, mype) ! Loop over all columns do c=clm_begc,clm_endc - ! Add columns in gridcell g + ! Select columns in gridcell g if(col%gridcell(c).eq.g) then - ! If there are hydrologically inactive columns in the - ! same gridcell with active columns + ! 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.0 end if @@ -408,14 +411,15 @@ subroutine set_clm_statevec(tstartcycle, mype) if(n_c == 0.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 before from first column." + 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 ! Compute final average clm_statevec(cc) = clm_statevec(cc) / n_c end if - ! Save prior state vector for increment + ! Save prior column mean state vector for computing + ! increment in updating the state vector clm_statevec_orig(cc) = clm_statevec(cc) end do From ba76793b87500373dd0271baf515833c02519ace Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 10:48:37 +0200 Subject: [PATCH 07/10] make n_c integer --- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) 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 c7f932da9..2c600b409 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 @@ -354,7 +354,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) integer :: i,j,jj,g,c,cc=0,offset=0 - real :: n_c + 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 @@ -389,7 +389,7 @@ subroutine set_clm_statevec(tstartcycle, mype) do cc = 1, clm_statevecsize clm_statevec(cc) = 0.0 - n_c = 0.0 + n_c = 0 ! Get gridcell and layer g = col%gridcell(state_pdaf2clm_c_p(cc)) @@ -403,19 +403,19 @@ subroutine set_clm_statevec(tstartcycle, mype) 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.0 + n_c = n_c + 1 end if end if end do - if(n_c == 0.0) then + 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 - ! Compute final average - clm_statevec(cc) = clm_statevec(cc) / n_c + ! Normalize sum to average + clm_statevec(cc) = clm_statevec(cc) / real(n_c, r8) end if ! Save prior column mean state vector for computing From 202d3b7e4b697629c53f6d242d549f9d63d7a7c2 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 10:48:51 +0200 Subject: [PATCH 08/10] refactor update --- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) 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 2c600b409..c320027c9 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 @@ -608,11 +608,14 @@ 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 with the factor of the update for the column mean + ! 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 From a3e8c02fffc787e79f65432b4d7df81b070bfa67 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 10:50:16 +0200 Subject: [PATCH 09/10] update clm5_0 module --- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 276 ++++++++++++------ 1 file changed, 179 insertions(+), 97 deletions(-) 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..c320027c9 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)) - - 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 + end if - ! One value per grid-cell - clm_varsize = (endg-begg+1) * nlevsoi - clm_statevecsize = (endg-begg+1) * nlevsoi + ! 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 - 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)) + ! 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)) - cc = 0 + ! Defaults + do cc=1,clm_statevecsize + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do - 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 - end do - end do + do cc=1,clm_statevecsize + lay: do i=1,nlevsoi + col: 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 col + 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 From 421056d7b21c35b1a5de81c64e09fa44171aa120 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Jun 2025 12:59:01 +0200 Subject: [PATCH 10/10] bugfix: remove loop-name `col` Error message: ``` eclm/enkf_clm_mod_5.F90(264): error #6401: The attributes of this name conflict with those made accessible by a USE statement. [COL] ``` --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 4 ++-- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) 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 c320027c9..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 @@ -256,14 +256,14 @@ subroutine define_clm_statevec(mype) do cc=1,clm_statevecsize lay: do i=1,nlevsoi - col: do c=clm_begc,clm_endc + 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 col + end do end do lay #ifdef PDAF_DEBUG 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 c320027c9..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 @@ -256,14 +256,14 @@ subroutine define_clm_statevec(mype) do cc=1,clm_statevecsize lay: do i=1,nlevsoi - col: do c=clm_begc,clm_endc + 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 col + end do end do lay #ifdef PDAF_DEBUG