Snow-DA for eCLM-PDAF

interface/framework/init_dim_obs_f_pdaf.F90

@@ -145,6 +145,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
@@ -986,7 +987,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
                    if(clmstatevec_only_active==1) then
                      obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i))
                    else
-                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                     if(clmupdate_snow/=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
 #else
                    obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))

interface/framework/init_dim_obs_pdaf.F90

@@ -141,6 +141,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
@@ -979,7 +980,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
                    if(clmstatevec_only_active==1) then
                      obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i))
                    else
-                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                     if(clmupdate_snow/=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
 #else
                    obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))

interface/model/clm3_5/enkf_clm_mod.F90

@@ -66,6 +66,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
@@ -148,6 +150,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 
@@ -288,6 +309,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
@@ -535,6 +564,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()

interface/model/common/enkf.h

@@ -88,6 +88,8 @@ GLOBAL int nx_local,ny_local,nz_local;
 GLOBAL int clmupdate_swc;
 GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
+GLOBAL int clmupdate_snow;
+GLOBAL int clmupdate_snow_repartitioning;
 GLOBAL int clmprint_swc;
 GLOBAL int clmprint_et;
 GLOBAL int clmstatevec_allcol;

interface/model/common/read_enkfpar.c

@@ -78,6 +78,8 @@ void read_enkfpar(char *parname)
   clmupdate_swc         = iniparser_getint(pardict,"CLM:update_swc",1);
   clmupdate_T           = iniparser_getint(pardict,"CLM:update_T",0);
   clmupdate_texture     = iniparser_getint(pardict,"CLM:update_texture",0);
+  clmupdate_snow        = iniparser_getint(pardict,"CLM:update_snow",0);
+  clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",3);
   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);