GRACE-DA in eCLM-PDAF

docs/users_guide/introduction_to_tsmp_pdaf/README.md

@@ -40,6 +40,7 @@ In alphabetic order (to be extended):
 * Stefan Poll
 * Mukund Pondkule
 * Prabhakar Shrestha
+* Anne Springer
 * Lukas Strebel
 
 ## About this documentation

interface/framework/add_obs_error_pdaf.F90

@@ -50,8 +50,11 @@ SUBROUTINE add_obs_error_pdaf(step, dim_obs, C_p)
 ! !USES:
   USE mod_assimilation, &
        ONLY: rms_obs, obs_pdaf2nc
+  USE mod_assimilation, ONLY: obscov
 
   USE mod_read_obs, ONLY: multierr,clm_obserr, pressure_obserr
+  USE mod_read_obs, ONLY: vec_useObs_global
+  USE enkf_clm_mod, ONLY: clmupdate_tws
   USE mod_parallel_pdaf, ONLY: mype_world
   USE mod_parallel_pdaf, ONLY: abort_parallel
   USE mod_tsmp, ONLY: point_obs
@@ -72,7 +75,9 @@ SUBROUTINE add_obs_error_pdaf(step, dim_obs, C_p)
 
 ! *** local variables ***
   INTEGER :: i          ! index of observation component
+  INTEGER :: j          ! index of observation component
   REAL :: variance_obs  ! variance of observations
+  REAL :: clm_obserr_model(dim_obs) ! errors of observations in the model domain
 
 
 ! **********************
@@ -98,6 +103,8 @@ SUBROUTINE add_obs_error_pdaf(step, dim_obs, C_p)
 
   if(multierr==1) then
 
+    NOGRACE: if (clmupdate_tws/=1) then
+
     ! Check that point observations are used
     if (.not. point_obs == 1) then
       print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR(3) `point_obs.eq.1` needed for using obs_pdaf2nc."
@@ -111,6 +118,24 @@ SUBROUTINE add_obs_error_pdaf(step, dim_obs, C_p)
       C_p(i,i) = C_p(i,i) + pressure_obserr(obs_pdaf2nc(i))*pressure_obserr(obs_pdaf2nc(i))
 #endif
     enddo
+
+    else NOGRACE
+
+    clm_obserr_model = pack(clm_obserr,vec_useObs_global)
+    do i = 1,dim_obs
+      C_p(i,i) = C_p(i,i) + clm_obserr_model(i) ! we put already variances in GRACE observation files, so no square need
+    end do
+
+    end if NOGRACE
+
   endif
 
+  if(multierr==2) then
+    do i = 1, size(obscov,1)
+      do j = 1, size(obscov,2)
+        C_p(i,j) = C_p(i,j)+obscov(i,j)
+      end do
+    end do
+  end if
+
 END SUBROUTINE add_obs_error_pdaf

interface/framework/g2l_obs_pdaf.F90

@@ -55,6 +55,7 @@ SUBROUTINE g2l_obs_pdaf(domain_p, step, dim_obs_f, dim_obs_l, mstate_f, &
 ! !USES:
   USE mod_assimilation, &
        ONLY: obs_index_l, obs_nc2pdaf
+  USE enkf_clm_mod, ONLY: clmupdate_tws
 
   IMPLICIT NONE
 
@@ -89,8 +90,12 @@ SUBROUTINE g2l_obs_pdaf(domain_p, step, dim_obs_f, dim_obs_l, mstate_f, &
   ! Index array OBS_NC2PDAF set in subroutine INIT_DIM_OBS_F_PDAF
   ! Index array OBS_INDEX_L (returns nc-ordered index) set in subroutine INIT_DIM_OBS_L_PDAF
   do i=1,dim_obs_l
+   GRACE: if (clmupdate_tws==1) then
+    mstate_l(i) = mstate_f(obs_index_l(i))
+   else GRACE
     !mstate_l(i) = mstate_f(obs_index_l(i))
     mstate_l(i) = mstate_f(obs_nc2pdaf(obs_index_l(i)))
+   end if GRACE
   end do
 
 END SUBROUTINE g2l_obs_pdaf

interface/framework/g2l_state_pdaf.F90

@@ -52,6 +52,8 @@ SUBROUTINE g2l_state_pdaf(step, domain_p, dim_p, state_p, dim_l, state_l)
   USE mod_tsmp, &
        ONLY: nx_local, ny_local
 #if defined CLMSA
+  use enkf_clm_mod, only: hactiveg_levels, num_layer, state_setup, num_hactiveg_patch, hactiveg_patch, clm_varsize_tws
+  USE enkf_clm_mod, ONLY: clmupdate_tws
   USE enkf_clm_mod, ONLY: g2l_state_clm
 #endif
 
@@ -68,6 +70,7 @@ SUBROUTINE g2l_state_pdaf(step, domain_p, dim_p, state_p, dim_l, state_l)
   REAL, TARGET, INTENT(out)   :: state_l(dim_l) ! State vector on local analysis domain
 
   INTEGER :: i, n_domain, nshift_p
+  INTEGER :: sub, j, g
   INTEGER :: begg, endg   ! per-proc gridcell ending gridcell indices
 ! !CALLING SEQUENCE:
 ! Called by: PDAF_lseik_update    (as U_g2l_state)
@@ -90,7 +93,157 @@ SUBROUTINE g2l_state_pdaf(step, domain_p, dim_p, state_p, dim_l, state_l)
   end if
   !call g2l_state(domain_p, c_loc(state_p), dim_l, c_loc(state_l))
 #else
+  if (clmupdate_tws/=1) then
   call g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l)
+  end if
+
+  if (clmupdate_tws==1) then
+
+
+    ! first depth dependent variables --> liq and ice (together in statevector or not)
+    ! dim_l is number of layers of gridcell + 3 (3 for other compartments that are added to the statevector)
+    ! lets loop over known layers
+
+    ! first do canopy water to check how many variables are for this gridcell available next to snow and surface water as well as the layers
+    ! compartments to subtract from dim_l to get layers variables:
+
+    if (clm_varsize_tws(5)/=0) then
+      sub=3
+    else
+      sub=2
+    end if
+
+
+    select case (state_setup)
+    case(0) ! liq and ice seperated
+      g = hactiveg_levels(domain_p,1)
+      do i = 1, (dim_l-sub)/2 ! two entries for liq and ice seperated
+        do j = 1, num_layer(i) ! i is the layer that we are in right now
+          if (g==hactiveg_levels(j,i)) then ! if the counter is the gridcell of the local domain, we know the position in the statevector
+
+            if (i == 1) then ! if first layer
+              state_l(i) = state_p(j)     ! first liquid water as it is first in the statevector
+              state_l(i+(dim_l-3)/2) = state_p(j+clm_varsize_tws(1))
+            else
+              state_l(i) = state_p(j + sum(num_layer(1:i-1)))
+              state_l(i+(dim_l-3)/2) = state_p(j + sum(num_layer(1:i-1)) + clm_varsize_tws(1))
+            end if
+
+          end if
+        end do
+      end do
+
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          if (sub==3) then
+            state_l(dim_l-2) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+            state_l(dim_l-1) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3))
+            state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3) + clm_varsize_tws(4))
+          else
+            state_l(dim_l-1) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+            state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3))
+          end if
+        end if
+      end do
+
+    case(1)
+
+      g = hactiveg_levels(domain_p,1)
+      do i = 1, dim_l-sub ! liq and ice added up
+        do j = 1, num_layer(i) ! i is the layer that we are in right now
+          if (g==hactiveg_levels(j,i)) then ! if the counter is the gridcell of the local domain, we know the position in the statevector
+
+            if (i == 1) then ! if first layer
+              state_l(i) = state_p(j)     ! first liquid water as it is first in the statevector
+            else
+              state_l(i) = state_p(j + sum(num_layer(1:i-1)))
+            end if
+
+          end if
+        end do
+      end do
+
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          if (sub==3) then
+            state_l(dim_l-2) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+            state_l(dim_l-1) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3))
+            state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3) + clm_varsize_tws(4))
+          else
+            state_l(dim_l-1) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+            state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3))
+          end if
+        end if
+      end do
+
+    case(2) ! only tws in statevector
+      g = hactiveg_levels(domain_p,1)
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          state_l(1) = state_p(j)
+        end if
+      end do
+
+    case(3)
+
+      g = hactiveg_levels(domain_p,1)
+
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          state_l(1) = state_p(j)
+          state_l(2) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+        end if
+      end do
+
+    case(4)
+
+      g = hactiveg_levels(domain_p,1)
+
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          state_l(1) = state_p(j)
+          state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+        end if
+      end do
+
+      if (dim_l==3) then
+        do j = 1, num_layer(8)
+          if (g==hactiveg_levels(j,8)) then
+            state_l(2) = state_p(j + clm_varsize_tws(1))
+          end if
+        end do
+      end if
+
+    case(5)
+
+      g = hactiveg_levels(domain_p,1)
+
+      do j = 1, num_layer(1)
+        if (g==hactiveg_levels(j,1)) then
+          state_l(1) = state_p(j)
+          state_l(dim_l) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2) + clm_varsize_tws(3))
+        end if
+      end do
+
+      if (dim_l>=3) then
+        do j = 1, num_layer(4)
+          if (g==hactiveg_levels(j,4)) then
+            state_l(2) = state_p(j + clm_varsize_tws(1))
+          end if
+        end do
+      end if
+
+      if (dim_l>=4) then
+        do j = 1, num_layer(13)
+          if (g==hactiveg_levels(j,13)) then
+            state_l(3) = state_p(j + clm_varsize_tws(1) + clm_varsize_tws(2))
+          end if
+        end do
+      end if
+
+    end select
+
+  end if
 #endif
 
 END SUBROUTINE g2l_state_pdaf

interface/framework/init_dim_l_pdaf.F90

@@ -49,7 +49,11 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
        tag_model_clm, model
   USE mod_tsmp, &
        ONLY: init_dim_l_pfl
+  use clm_varpar   , only : nlevsoi
 #ifdef CLMSA
+  USE enkf_clm_mod, &
+       ONLY: hactiveg_levels, num_layer, state_setup, num_hactiveg_patch, hactiveg_patch, clm_varsize_tws
+  USE enkf_clm_mod, ONLY: clmupdate_tws
   USE enkf_clm_mod, ONLY: init_dim_l_clm
 #endif
   IMPLICIT NONE
@@ -59,6 +63,8 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
   INTEGER, INTENT(in)  :: domain_p ! Current local analysis domain
   INTEGER, INTENT(out) :: dim_l    ! Local state dimension
 
+  integer :: g, i, count
+
 ! !CALLING SEQUENCE:
 ! Called by: PDAF_lseik_update   (as U_init_dim_l)
 ! Called by: PDAF_lestkf_update  (as U_init_dim_l)
@@ -68,6 +74,7 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
 ! ****************************************
 ! *** Initialize local state dimension ***
 ! ****************************************
+
 #if defined PARFLOW_STAND_ALONE
   ! Set the size of the local analysis domain
   call init_dim_l_pfl(dim_l)
@@ -88,7 +95,91 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
   ! Set the size of the local analysis domain
   ! for clm stand alone mode only
 
+  if (clmupdate_tws/=1) then
   call init_dim_l_clm(domain_p, dim_l)
+  end if
+
+  ! Set the size of the local analysis domain
+  ! for clm stand alone mode only
+  if (clmupdate_tws==1) then
+
+    ! initialize local state dimension --> different for each domain as a domain is a gridcell, number of layers in statevector differs per gridcell
+    ! --> check with begg and endg as well as domain_p and hactiveg_levels what is going on
+    ! hactiveg_levels in level 1 at position domain_p gives gridcell index of domain_p --> go through rest of hactiveg_levels and check for number of layers
+    dim_l = 0
+    g = hactiveg_levels(domain_p,1)
+
+    select case(state_setup)
+    case(0)
+      do i = 1,nlevsoi
+        do count = 1, num_layer(i)
+          if (g==hactiveg_levels(count,i)) then
+            dim_l = dim_l+2
+          end if
+        end do
+      end do
+
+      ! snow and surface water
+      dim_l = dim_l+2
+
+      ! canopy water
+      if (clm_varsize_tws(5)/=0) then
+        dim_l = dim_l+1
+      end if
+
+    case(1)
+      do i = 1,nlevsoi
+        do count = 1, num_layer(i)
+          if (g==hactiveg_levels(count,i)) then
+            dim_l = dim_l+1
+          end if
+        end do
+      end do
+
+      ! snow and surface water
+      dim_l = dim_l+2
+
+      if (clm_varsize_tws(5)/=0) then
+        dim_l = dim_l+1
+      end if
+
+    case(2) ! only TWS in statevector
+
+      dim_l=1
+
+    case(3) ! sum over all soil layers and snow in statevector
+
+      dim_l=2
+
+    case(4)
+
+      dim_l=2
+
+      do count = 1, num_layer(8)
+        if (g==hactiveg_levels(count,8)) then
+          dim_l = dim_l+1
+        end if
+      end do
+
+    case(5)
+
+      dim_l=2
+
+      do count = 1, num_layer(4)
+        if (g==hactiveg_levels(count,4)) then
+          dim_l = dim_l+1
+        end if
+      end do
+
+      do count = 1, num_layer(13)
+        if (g==hactiveg_levels(count,13)) then
+          dim_l = dim_l+1
+        end if
+      end do
+
+    end select
+
+  end if
 #endif
 
 END SUBROUTINE init_dim_l_pdaf