Snow-DA

bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90

@@ -135,6 +135,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_snow
   !hcp end
 #endif
 #endif
@@ -966,7 +967,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
                    end if
 #endif
                  else
-                   obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   if(clmupdate_snow.ne.0) then
+                     ! Snow-DA: no layer in state vector variables
+                     obs_index_p(cnt) = g-begg+1
+                   else
+                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   end if
                  end if
 
                  !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)

bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90

@@ -131,6 +131,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_snow
   !hcp end
 #endif
 #endif
@@ -959,7 +960,12 @@ 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))
+                   if(clmupdate_snow.ne.0) then
+                     ! Snow-DA: no layer in state vector variables
+                     obs_index_p(cnt) = g-begg+1
+                   else
+                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   end if
                  end if
 
                  !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)

bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90

@@ -65,6 +65,8 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmupdate_swc")     :: clmupdate_swc
   integer(c_int),bind(C,name="clmupdate_T")     :: clmupdate_T  ! by hcp
   integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
+  integer(c_int),bind(C,name="clmupdate_snow")    :: clmupdate_snow
+  integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning
   integer(c_int),bind(C,name="clmprint_swc")      :: clmprint_swc
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
@@ -147,6 +149,25 @@ subroutine define_clm_statevec(mype)
       error stop "Not implemented: clmupdate_texture.eq.2"
     endif
 
+    ! Snow assimilation
+    ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5
+    ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell
+    ! Case 2: Assimilation of snow water equivalent same as above
+    if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then 
+      error stop "Not implemented: clmupdate_snow.eq.1 or clmupdate_snow.eq.1"
+    endif
+    ! Case 3: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    if(clmupdate_snow.eq.3) then
+      error stop "Not implemented: clmupdate_snow.eq.3"
+    endif
+    ! Case 4: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    ! and dz
+    if(clmupdate_snow.eq.4) then
+      error stop "Not implemented: clmupdate_snow.eq.4"
+    endif
+
     !hcp LST DA
     if(clmupdate_T.eq.1) then
       clm_varsize      =  endg-begg+1 
@@ -287,6 +308,14 @@ subroutine set_clm_statevec(tstartcycle, mype)
       end do
     endif
 
+    ! Snow assimilation state vector
+    ! Case 1: Snow depth
+    ! Case 2: SWE
+    ! Case 3: Snow depth + SWE
+    if(clmupdate_snow.ne.0) then
+      error stop "Not implemented: clmupdate_snow.ne.0"
+    endif
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -534,6 +563,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       call clm_texture_to_parameters
     endif
 
+    ! Snow assimilation:
+    ! Case 1: Snow depth
+    ! Case 2: Snow water equivalent
+    ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector
+    ! Write updated snow variable back to CLM and then repartition snow and adjust related variables
+    if(clmupdate_snow.ne.0) then
+      error stop "Not implemented: clmupdate_snow.ne.0"
+    endif
+
   end subroutine update_clm
 
   subroutine clm_correct_texture()