LAI-DA for eCLM-PDAF

interface/framework/obs_op_pdaf.F90

@@ -71,6 +71,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 #if defined CLMSA
    USE enkf_clm_mod, &
         ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd
+   USE enkf_clm_mod, &
+        ONLY : clmupdate_lai, clm_begp, clm_endp, clm_patch2gc, clm_patchwt
 #ifdef CLMFIVE
    USE clm_instMod, &
      ONLY : soilstate_inst
@@ -128,6 +130,36 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 ! If no special observation operator is compiled, use point observations
 lpointobs = .true.
 
+! LAI assimilation assuming statevec contains leafc, slatop, dsladlai for each patch
+! two component op to get to LAI obs
+! 1 calculate tlai from statevec per patch
+! state_p contains leafc of varsize then slatop of varsize and then dsadlai of varsize
+! 2 average patches according to weight per gridcell
+if (clmupdate_lai==2) then
+  lpointobs = .false. ! disable general obs_op
+  ! Note: This fails if there ever is an observation pointing to a gridcell with 0 patches
+  ! this should not happen, but would cause div by zero.
+
+  do i = 1, dim_obs_p ! loop over all observations
+    avesm = 0.0 ! use avesm as grid cell average collecter
+    do j = 1, clm_endp-clm_begp ! loop over all patches
+      write(*,*) 'DEBUG LAI : obs_index_p(i), clm_patch2gc(j)', obs_index_p(i), clm_patch2gc(j)
+      if (obs_index_p(i)==clm_patch2gc(j)) then
+        write(*,*) 'DEBUG LAI : state dsladlai', state_p(j+2*clm_varsize)
+        if (state_p(j+2*clm_varsize)>0.0) then
+          avesm = avesm + clm_patchwt(j) * ((state_p(j+1*clm_varsize)*(exp(state_p(j)*state_p(j+2*clm_varsize)) - 1.0))/state_p(j+2*clm_varsize)) ! formula for tlai from leafc
+        else
+          avesm = avesm +  clm_patchwt(j) *(state_p(j+1*clm_varsize)*state_p(j)) ! 2nd formula
+        endif ! dsladlai decider
+      endif ! this patch is in gridcell of observation
+      write(*,*) 'DEBUG LAI : average ', avesm
+    enddo ! all patches for obs checked
+    m_state_p(i) = avesm ! == lai of gridcell where the observation is
+  enddo ! all observations
+endif ! clmupdate_lai == 2 : m_state_p now contains lai for each gridcell with an observation.
+
+
+
 #if defined CLMSA
 if (clmupdate_T==1) then
 

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_lai;
+GLOBAL int clmupdate_lai_params;
 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_lai         = iniparser_getint(pardict,"CLM:update_lai",0);
+  clmupdate_lai_params  = iniparser_getint(pardict,"CLM:update_lai_params",0);
   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);