PDAF: Soil moisture state vector with gridcell averages (#268)

* CLM:statevec_colmean

- introduce hydrologically active grid cells only
- save prior column means in `clm_statevec_orig`, used for computing update increment
- update column values with increment-factor from mean

* define_clm_statevec: re-arranging state-vector definition

* 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]
```

Author: jjokella

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)

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

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;