eCLM-PDAF with local filters (#16)

* CLMSA: Use `cradius` in km in `init_dim_obs_l_pdaf`

* LESTKF/LETKF: no deallocation of observations lat/lon arrays

* init_dim_l_clm: dependent on domain_p

Flexibility needed for different dimensions of local domains

* local filter for eCLM: flexible index handling with index trafo

- STATE_LOC2CLM_C_P: Returns the CLM-column c for the local domain
domain_p (from the column, the gridcell can be derived)

* local filters: INIT_DIM_L_CLM, add hydrologically-active condition

Author: jjokella

interface/framework/g2l_state_pdaf.F90

@@ -52,14 +52,10 @@ 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
-#if defined CLMFIVE
-  USE decompMod, ONLY: get_proc_bounds
-#else  
-  USE decompMod, ONLY: get_proc_bounds_atm
-#endif
+  USE enkf_clm_mod, ONLY: g2l_state_clm
 #endif
 
-  USE iso_c_binding, ONLY: c_loc
+  ! USE iso_c_binding, ONLY: c_loc
 
   IMPLICIT NONE
 
@@ -94,17 +90,7 @@ 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
-  ! beg and end gridcell for atm
-#if defined CLMFIVE
-  call get_proc_bounds(begg, endg)
-#else  
-  call get_proc_bounds_atm(begg, endg)
-#endif
-  n_domain = endg - begg + 1
-  DO i = 0, dim_l-1
-     nshift_p = domain_p + i * n_domain
-     state_l(i+1) = state_p(nshift_p)
-  ENDDO
+  call g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l)
 #endif
 
 END SUBROUTINE g2l_state_pdaf

interface/framework/init_dim_l_pdaf.F90

@@ -48,10 +48,9 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
   USE mod_tsmp, ONLY: tag_model_parflow, &
        tag_model_clm, model
   USE mod_tsmp, &
-       ONLY: init_parf_l_size
+       ONLY: init_dim_l_pfl
 #ifdef CLMSA
-  USE enkf_clm_mod, &
-       ONLY: init_clm_l_size
+  USE enkf_clm_mod, ONLY: init_dim_l_clm
 #endif
   IMPLICIT NONE
 
@@ -69,22 +68,27 @@ SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l)
 ! ****************************************
 ! *** Initialize local state dimension ***
 ! ****************************************
-#if (defined PARFLOW_STAND_ALONE || defined COUP_OAS_PFL)
-  if (model.eq.tag_model_parflow) then
+#if defined PARFLOW_STAND_ALONE
+  ! Set the size of the local analysis domain
+  call init_dim_l_pfl(dim_l)
+#endif
+
+#if defined COUP_OAS_PFL
+  if (model == tag_model_parflow) then
      ! Set the size of the local analysis domain 
-     call init_parf_l_size(dim_l)
+     call init_dim_l_pfl(dim_l)
   end if
-#endif  
-
-#ifndef CLMSA
-  if (model.eq.tag_model_clm) then
+  if (model == tag_model_clm) then
      ! Set the size of the local analysis domain   
      dim_l = 1     
   end if   
-#else
+#endif
+
+#if defined CLMSA
   ! Set the size of the local analysis domain  
   ! for clm stand alone mode only
-  call init_clm_l_size(dim_l)
+
+  call init_dim_l_clm(domain_p, dim_l)
 #endif
 
 END SUBROUTINE init_dim_l_pdaf

interface/framework/init_dim_obs_l_pdaf.F90

@@ -71,6 +71,18 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
        point_obs, model
 #endif
 
+#if defined CLMSA
+  USE enkf_clm_mod, ONLY: state_loc2clm_c_p
+  use shr_kind_mod, only: r8 => shr_kind_r8
+
+#ifdef CLMFIVE
+  USE GridcellType, ONLY: grc
+  USE ColumnType, ONLY : col
+#else
+  USE clmtype, ONLY : clm3
+#endif
+#endif
+
   USE, INTRINSIC :: iso_c_binding, ONLY: C_F_POINTER
 
   IMPLICIT NONE
@@ -95,11 +107,26 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
   INTEGER  :: domain_p_coord   ! Current local analysis domain for coord arrays
 
   !kuw
-  integer :: dx,dy, max_var_id, ierror
+#if defined CLMSA
+  real :: dx,dy
+#else
+  integer :: dx,dy
+#endif
+  integer :: max_var_id, ierror
   integer :: obsind(dim_obs)
   real    :: obsdist(dim_obs)
   ! kuw end
 
+#if defined CLMSA
+  ! INTEGER :: dim_l
+  ! INTEGER :: ncellxy
+  ! INTEGER :: k
+  real(r8), pointer :: lon(:)
+  real(r8), pointer :: lat(:)
+  integer, pointer :: mycgridcell(:) !Pointer for CLM3.5/CLM5.0 col->gridcell index arrays
+  REAL :: yhalf
+#endif
+
   ! **********************************************
   ! *** Initialize local observation dimension ***
   ! **********************************************
@@ -241,10 +268,44 @@ SUBROUTINE init_dim_obs_l_pdaf(domain_p, step, dim_obs_f, dim_obs_l)
      end if
   else 
      if(model == tag_model_clm) THEN
+
+#ifdef CLMSA
+    ! Lon/Lat information from CLM
+#ifdef CLMFIVE
+    ! Obtain CLM lon/lat information
+    lon   => grc%londeg
+    lat   => grc%latdeg
+    ! Obtain CLM column-gridcell information
+    mycgridcell => col%gridcell
+#else
+    lon   => clm3%g%londeg
+    lat   => clm3%g%latdeg
+    mycgridcell => clm3%g%l%c%gridcell
+#endif
+#endif
+
      do i = 1,dim_obs
+#ifdef CLMSA
+        ! Units: lat/lon (degrees)
+        ! More doc on following lines: See `localize_covar_pdaf`
+
+       ! Compared to LOCALIZE_COVAR_PDAF: No OBS_PDAF2NC. This is in
+       ! order to have OBS_INDEX_L return a NC-ordered array, not
+       ! PDAF-ordered array.
+        dx = abs(clmobs_lon(i) - lon(mycgridcell(state_loc2clm_c_p(domain_p))))
+        dy = abs(clmobs_lat(i) - lat(mycgridcell(state_loc2clm_c_p(domain_p))))
+        IF (dx > 180.0) THEN
+          dx = 360.0 - dx
+        END IF
+        yhalf = ( clmobs_lat(i) + lat(mycgridcell(state_loc2clm_c_p(domain_p))) ) / 2.0
+        dx = dx * cos(yhalf * 3.14159265358979323846 / 180.0)
+        dist = 111.19492664455873 * sqrt(real(dx)**2 + real(dy)**2)
+#else
+        ! Units: Index numbering
         dx = abs(longxy_obs(i) - longxy(domain_p))
         dy = abs(latixy_obs(i) - latixy(domain_p))
         dist = sqrt(real(dx)**2 + real(dy)**2)
+#endif
         obsdist(i) = dist
         if (dist <= real(cradius)) then
            dim_obs_l = dim_obs_l + 1

interface/framework/init_n_domains_pdaf.F90

@@ -51,23 +51,16 @@ SUBROUTINE init_n_domains_pdaf(step, n_domains_p)
   USE mod_assimilation, &
       ONLY: dim_state_p
   USE mod_tsmp, &
-      ONLY: init_n_domains_size
+      ONLY: init_n_domains_pfl
 #if defined CLMSA
-#if defined CLMFIVE
-  USE decompMod, ONLY: get_proc_bounds
-#else
-  USE decompMod, ONLY: get_proc_bounds_atm
-#endif
+  USE enkf_clm_mod, ONLY: init_n_domains_clm
 #endif
 
   IMPLICIT NONE
 
 ! !ARGUMENTS:
   INTEGER, INTENT(in)  :: step        ! Current time step
   INTEGER, INTENT(out) :: n_domains_p ! PE-local number of analysis domains
-#if defined CLMSA
-  INTEGER :: begg, endg   ! per-proc gridcell ending gridcell indices
-#endif
 
 ! !CALLING SEQUENCE:
 ! Called by: PDAF_lseik_update   (as U_init_n_domains)
@@ -78,27 +71,22 @@ SUBROUTINE init_n_domains_pdaf(step, n_domains_p)
 ! ************************************
 ! *** Initialize number of domains ***
 ! ************************************
-#if (defined PARFLOW_STAND_ALONE || defined COUP_OAS_PFL)
-  if (model.eq.tag_model_parflow) then
-     ! Here simply the process-local state dimension
-     call init_n_domains_size(n_domains_p)
-  end if
+#if defined PARFLOW_STAND_ALONE
+  call init_n_domains_pfl(n_domains_p)
 #endif   
 
-#ifndef CLMSA
+#if defined COUP_OAS_PFL
+  if (model == tag_model_parflow) then
+     call init_n_domains_pfl(n_domains_p)
+  end if
   if (model == tag_model_clm) then
-     ! Here simply the process-local state dimension  
+     ! For CLM coupled with ParFlow, use the process-local state dimension
      n_domains_p = dim_state_p
-  end if   
-#else
-  ! beg and end gridcell for atm
-#if defined CLMFIVE
-  call get_proc_bounds(begg, endg)
-#else  
-  call get_proc_bounds_atm(begg, endg)
+  end if
 #endif
-  ! Here simply the process-local state dimension  
-  n_domains_p = endg - begg + 1
+
+#if defined CLMSA
+  call init_n_domains_clm(n_domains_p)
 #endif
 
 END SUBROUTINE init_n_domains_pdaf

interface/framework/init_obs_l_pdaf.F90

@@ -74,6 +74,7 @@ SUBROUTINE init_obs_l_pdaf(domain_p, step, dim_obs_l, observation_l)
   observation_l(:) = 0.0
 
 !#ifndef CLMSA
+  ! Index array OBS_INDEX_L (returns nc-ordered index) set in subroutine INIT_DIM_OBS_L_PDAF
   do i=1,dim_obs_l
     observation_l(i) = obs(obs_index_l(i)) 
   end do

interface/framework/l2g_state_pdaf.F90

@@ -55,11 +55,7 @@ SUBROUTINE l2g_state_pdaf(step, domain_p, dim_l, state_l, dim_p, state_p)
   USE iso_c_binding, ONLY: c_loc
 
 #if defined CLMSA
-#if defined CLMFIVE
-  USE decompMod, ONLY: get_proc_bounds
-#else  
-  USE decompMod, ONLY: get_proc_bounds_atm
-#endif
+  use enkf_clm_mod, ONLY: l2g_state_clm
 #endif
 
   IMPLICIT NONE
@@ -95,17 +91,7 @@ SUBROUTINE l2g_state_pdaf(step, domain_p, dim_l, state_l, dim_p, state_p)
   end if   
   !call l2g_state(domain_p, c_loc(state_p), dim_l, c_loc(state_l))
 #else
-  ! beg and end gridcell for atm
-#if defined CLMFIVE
-  call get_proc_bounds(begg, endg)
-#else  
-  call get_proc_bounds_atm(begg, endg)
-#endif
-  n_domain = endg - begg + 1
-  DO i = 0, dim_l-1
-     nshift_p = domain_p + i * n_domain
-     state_p(nshift_p) = state_l(i+1)
-  ENDDO
+  call l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p)
 #endif
 
 END SUBROUTINE l2g_state_pdaf

interface/framework/localize_covar_pdaf.F90

@@ -38,7 +38,7 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
    USE shr_kind_mod , only : r8 => shr_kind_r8
    USE mod_read_obs, ONLY: clmobs_lon
    USE mod_read_obs, ONLY: clmobs_lat
-   USE enkf_clm_mod, ONLY: init_clm_l_size, clmupdate_T
+   USE enkf_clm_mod, ONLY: clmupdate_T
    USE enkf_clm_mod, ONLY: clm_begc
    USE enkf_clm_mod, ONLY: clm_endc
    USE enkf_clm_mod, ONLY: state_pdaf2clm_c_p
@@ -226,7 +226,7 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
          ! dx = abs(longxy_obs(j) - longxy(state_pdaf2clm_c_p(i)))
          ! dy = abs(latixy_obs(j) - latixy(state_pdaf2clm_c_p(i)))
 
-         ! Units: lat/lon
+         ! Units: lat/lon (degrees)
          dx = abs(clmobs_lon(obs_pdaf2nc(j)) - lon(mycgridcell(state_pdaf2clm_c_p(i))))
          dy = abs(clmobs_lat(obs_pdaf2nc(j)) - lat(mycgridcell(state_pdaf2clm_c_p(i))))
 
@@ -267,7 +267,7 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
          ! dx = abs(longxy_obs(j) - longxy_obs(i))
          ! dy = abs(latixy_obs(j) - latixy_obs(i))
 
-         ! Units: lat/lon
+         ! Units: lat/lon (degrees)
          dx = abs(clmobs_lon(obs_pdaf2nc(j)) - clmobs_lon(obs_pdaf2nc(i)))
          dy = abs(clmobs_lat(obs_pdaf2nc(j)) - clmobs_lat(obs_pdaf2nc(i)))
 

interface/framework/mod_read_obs.F90

@@ -517,8 +517,8 @@ subroutine clean_obs_nc()
     !if(allocated(y_idx_obs_nc))deallocate(y_idx_obs_nc)
     !if(allocated(z_idx_obs_nc))deallocate(z_idx_obs_nc)
     !kuw: clean clm observations
-    IF (.NOT. filtertype == 8) THEN
-      ! For LEnKF lat/lon are used
+    IF (.NOT. filtertype == 5 .AND. .NOT. filtertype == 7 .AND. .NOT. filtertype == 8) THEN
+      ! For LETKF, LESTKF, LEnKF lat/lon are used
       if(allocated(clmobs_lon))deallocate(clmobs_lon)
       if(allocated(clmobs_lat))deallocate(clmobs_lat)
     END IF

interface/framework/mod_tsmp.F90

@@ -91,19 +91,19 @@ end subroutine update_tsmp
     end interface
 
      interface
-        subroutine init_n_domains_size(n_domains_p) bind(c)
+        subroutine init_n_domains_pfl(n_domains_p) bind(c)
             use iso_c_binding
             import
             INTEGER(c_int) :: n_domains_p ! PE-local number of analysis domains
-        end subroutine init_n_domains_size
+        end subroutine init_n_domains_pfl
     end interface
 
      interface
-        subroutine init_parf_l_size(dim_l) bind(c)
+        subroutine init_dim_l_pfl(dim_l) bind(c)
             use iso_c_binding
             import
               INTEGER(c_int) :: dim_l ! Local state dimension
-        end subroutine init_parf_l_size
+        end subroutine init_dim_l_pfl
     end interface
 
 !!$    interface