diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 3161901e6..dba27c54a 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -6,5 +6,6 @@ Wolfgang Kurtz Stefan Poll Mukund Pondkule Prabhakar Shrestha +Anne Springer Lukas Strebel diff --git a/Makefile b/Makefile index 06ed26693..9478c096c 100644 --- a/Makefile +++ b/Makefile @@ -138,7 +138,8 @@ SRC_PDAF_GEN = PDAF_analysis_utils.F90 \ PDAFlocalomi_put_state_nondiagR.F90 \ PDAFlocalomi_put_state_nondiagR_si.F90 \ PDAFlocalomi_put_state_si.F90 \ - PDAF_correlation_function.F90 + PDAF_correlation_function.F90 \ + PDAF_reset_dim_p.F90 # Specific PDAF-routines for SEIK SRC_SEIK = PDAF_seik_init.F90 \ diff --git a/docs/_toc.yml b/docs/_toc.yml index 19e344131..3d53c18d1 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -28,6 +28,7 @@ parts: - file: users_guide/running_tsmp_pdaf/input_cmd - file: users_guide/running_tsmp_pdaf/input_obs - file: users_guide/running_tsmp_pdaf/input_enkfpf + - file: users_guide/running_tsmp_pdaf/tsmp_pdaf_omi - file: users_guide/debugging_tsmp_pdaf/README title: Debugging TSMP-PDAF diff --git a/docs/users_guide/running_tsmp_pdaf/input_cmd.md b/docs/users_guide/running_tsmp_pdaf/input_cmd.md index 17862f935..9940e0a92 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_cmd.md +++ b/docs/users_guide/running_tsmp_pdaf/input_cmd.md @@ -51,6 +51,15 @@ Details: `subtype` (integer) Parameter subtype, different options for each filter. See [](cmd:command-line-examples). +## use_omi ## + +`use_omi` (logical) Controls whether to use OMI interface. + +- `.true.`: OMI interface is used +- `.false.`: OMI interface is not used + +See [](omi:tsmp-pdaf-with-pdaf-omi). + ## obs_filename ## `obs_filename` (string) Prefix for observation files. diff --git a/docs/users_guide/running_tsmp_pdaf/tsmp_pdaf_omi.md b/docs/users_guide/running_tsmp_pdaf/tsmp_pdaf_omi.md new file mode 100644 index 000000000..f7dbfb3b5 --- /dev/null +++ b/docs/users_guide/running_tsmp_pdaf/tsmp_pdaf_omi.md @@ -0,0 +1,14 @@ +(omi:tsmp-pdaf-with-pdaf-omi)= +# TSMP-PDAF with PDAF-OMI + +PDAF-OMI for multivariate data assimiliaton. +It is design to handle different observation types (currently soil moisture and TWS) automatically. +Additional to current observation files, the observation type has to be included (SM, GRACE). + +See also create observation script for details: https://icg4geo.icg.kfa-juelich.de/ExternalRepos/tsmp-pdaf/tsmp-pdaf-observation-scripts/-/tree/main/omi_obs_data?ref_type=heads + +Both global and local filters can be used. To enable multi-scale data assimilation, different localization radii for different observation types can be passed. Note that the localization radius for SM is currently in km and for GRACE in #gridcells. + +The framework generates a state vector for each type individually before the assimilation, some things would need to be adapted when mutliple observation types are assimilated at the same timestep. Currently, one observation file only consists of one observation type. As SM observations are usually assimilated at noon and GRACE observations are assimilated at the end of the month at midnight, this should not provide any problems. + +If questions arise contact ewerdwalbesloh@geod.uni-bonn.de diff --git a/interface/framework/Makefile b/interface/framework/Makefile index ea27c0e5c..35719909e 100644 --- a/interface/framework/Makefile +++ b/interface/framework/Makefile @@ -75,10 +75,17 @@ MOD_ASSIM = mod_parallel_pdaf.o \ parser_mpi.o \ mod_read_obs.o +# Routines of observation handling (PDAF-OMI) +MOD_USER_PDAFOMI = obs_GRACE_pdafomi.o \ + obs_SM_pdafomi.o + # Model routines used with PDAF OBJ_MODEL_PDAF =pdaf_terrsysmp.o\ integrate_pdaf.o +# Routines of observation handling (PDAF-OMI) +OBJ_USER_PDAFOMI = callback_obs_pdafomi.o + # Interface to PDAF - model sided OBJ_PDAF_INT = init_parallel_pdaf.o \ finalize_pdaf.o \ @@ -121,7 +128,7 @@ OBJ_USER_LOCAL = init_n_domains_pdaf.o \ localize_covar_pdaf.o # Full list of user-supplied routines for online modes -OBJ_PDAF_USER = $(OBJ_USER_GEN) $(OBJ_USER_SEIK) $(OBJ_USER_ENKF) $(OBJ_USER_LOCAL) +OBJ_PDAF_USER = $(OBJ_USER_GEN) $(OBJ_USER_SEIK) $(OBJ_USER_ENKF) $(OBJ_USER_LOCAL) $(OBJ_USER_PDAFOMI) ###################################################### @@ -143,10 +150,12 @@ info: $(PROG) : $(LIBMODEL) libpdaf-d.a \ - $(MODULES) $(MOD_ASSIM) $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER) + $(MODULES) $(MOD_ASSIM) $(MOD_USER_PDAFOMI) \ + $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER) $(PREP_C) $(LD) $(OPT_LNK) -o $@ \ - $(MODULES) $(MOD_ASSIM) $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER) \ - -L$(BASEDIR)/lib -lpdaf-d $(LIBMODEL) $(LIBS) $(LINK_LIBS) + $(MODULES) $(MOD_ASSIM) $(MOD_USER_PDAFOMI) \ + $(OBJ_MODEL_PDAF) $(OBJ_PDAF_INT) $(OBJ_PDAF_USER) \ + -L$(BASEDIR)/lib $(LIBMODEL) -lpdaf-d $(LIBS) $(LINK_LIBS) ###################################################### diff --git a/interface/framework/assimilate_pdaf.F90 b/interface/framework/assimilate_pdaf.F90 index 5e754b7df..5492ef2f3 100644 --- a/interface/framework/assimilate_pdaf.F90 +++ b/interface/framework/assimilate_pdaf.F90 @@ -45,7 +45,14 @@ SUBROUTINE assimilate_pdaf() ONLY: abort_parallel, mype_world USE mod_assimilation, & ! Variables for assimilation ONLY: filtertype - ! USE PDAF_interfaces_module ! Check consistency of PDAF calls + USE mod_assimilation, ONLY: use_omi +#ifdef CLMFIVE + USE PDAF_interfaces_module, & ! Check consistency of PDAF calls + ONLY: PDAFomi_assimilate_local, PDAFomi_assimilate_global, & + PDAFomi_assimilate_lenkf, PDAF_get_localfilter, & + PDAFomi_assimilate_enkf_nondiagR, & + PDAFomi_assimilate_global_nondiagR, PDAFomi_assimilate_local_nondiagR +#endif IMPLICIT NONE @@ -57,6 +64,7 @@ SUBROUTINE assimilate_pdaf() ! Local variables INTEGER :: status_pdaf ! PDAF status flag + INTEGER :: localfilter ! Flag for domain-localized filter (1=true) ! ! External subroutines @@ -102,6 +110,19 @@ SUBROUTINE assimilate_pdaf() ! EXTERNAL :: likelihood_hyb_l_pdaf, & ! Compute local likelihood awith hybrid weight for an ensemble member ! prodRinvA_hyb_l_pdaf ! Provide product R^-1 A for some matrix A including hybrid weight + ! Interface to PDAF-OMI for local and global filters + EXTERNAL :: init_dim_obs_pdafomi, & ! Get dimension of full obs. vector for PE-local domain + obs_op_pdafomi, & ! Obs. operator for full obs. vector for PE-local domain + init_dim_obs_l_pdafomi, & ! Get dimension of obs. vector for local analysis domain + localize_covar_pdafomi, & ! Apply localization to covariance matrix in LEnKF + add_obs_err_pdafomi, & ! Add obs. error covariance R to HPH in EnKF + init_obscovar_pdafomi, & ! Initialize obs error covar R in EnKF + prodRinvA_pdafomi, & ! Provide product R^-1 A for some matrix A + prodRinvA_l_pdafomi ! Provide product R^-1 A for some local matrix A + + + + ! *** Switch on debug output *** ! *** for main process *** #ifdef PDAF_DEBUG @@ -112,6 +133,43 @@ SUBROUTINE assimilate_pdaf() ! *** Call assimilation routine *** ! ********************************* + OMI: IF (use_omi) THEN +#ifdef CLMFIVE + CALL PDAF_get_localfilter(localfilter) + + IF (localfilter == 1) THEN + + CALL PDAFomi_assimilate_local_nondiagR(collect_state_pdaf, distribute_state_pdaf, & + init_dim_obs_pdafomi, obs_op_pdafomi, prepoststep_ens_pdaf, init_n_domains_pdaf, & + init_dim_l_pdaf, init_dim_obs_l_pdafomi, prodRinvA_l_pdafomi, g2l_state_pdaf, & + l2g_state_pdaf, next_observation_pdaf, status_pdaf) + + + ELSE + + IF (filtertype == 8) THEN + ! LEnKF has its own OMI interface routine + CALL PDAFomi_assimilate_lenkf(collect_state_pdaf, distribute_state_pdaf, & + init_dim_obs_pdafomi, obs_op_pdafomi, prepoststep_ens_pdaf, & + localize_covar_pdafomi, next_observation_pdaf, status_pdaf) + + ELSEIF (filtertype == 2) then ! non diagonal R for EnKF has its own callback routine + CALL PDAFomi_assimilate_enkf_nondiagR(collect_state_pdaf, distribute_state_pdaf, & + init_dim_obs_pdafomi, obs_op_pdafomi, add_obs_err_pdafomi, init_obscovar_pdafomi, & + prepoststep_ens_pdaf, next_observation_pdaf, status_pdaf) + + ELSE + + CALL PDAFomi_assimilate_global_nondiagR(collect_state_pdaf, distribute_state_pdaf, & + init_dim_obs_pdafomi, obs_op_pdafomi, prodRinvA_pdafomi, & + prepoststep_ens_pdaf, next_observation_pdaf, status_pdaf) + + ENDIF + + ENDIF +#endif + ELSE OMI + ! IF (filtertype == 1) THEN ! CALL PDAF_assimilate_seik(collect_state_pdaf, distribute_state_pdaf, & ! init_dim_obs_pdaf, obs_op_pdaf, init_obs_pdaf, prepoststep_ens_pdaf, & @@ -182,6 +240,8 @@ SUBROUTINE assimilate_pdaf() ! likelihood_pdaf, next_observation_pdaf, status_pdaf) END IF + END IF OMI + ! Check for errors during execution of PDAF IF (status_pdaf /= 0) THEN diff --git a/interface/framework/callback_obs_pdafomi.F90 b/interface/framework/callback_obs_pdafomi.F90 new file mode 100644 index 000000000..6ac260d9e --- /dev/null +++ b/interface/framework/callback_obs_pdafomi.F90 @@ -0,0 +1,319 @@ +!> callback_obs_pdafomi +!! +!! This file provides interface routines between the call-back routines +!! of PDAF and the observation-specific routines in PDAF-OMI. This structure +!! collects all calls to observation-specific routines in this single file +!! to make it easier to find the routines that need to be adapted. +!! +!! The routines here are mainly pure pass-through routines. Thus they +!! simply call one of the routines from PDAF-OMI. Partly some addtional +!! variable is required, e.g. to specify the offset of an observation +!! in the observation vector containing all observation types. These +!! cases are described in the routines. +!! +!! **Adding an observation type:** +!! When adding an observation type, one has to add one module +!! obs_OBSTYPE_pdafomi (based on the template obs_OBSTYPE_pdafomi_TEMPLATE.F90). +!! In addition one has to add a call to the different routines include +!! in this file. It is recommended to keep the order of the calls +!! consistent over all files. +!! +!! __Revision history:__ +!! * 2019-12 - Lars Nerger - Initial code +!! * Later revisions - see repository log +!! +!------------------------------------------------------------------------------- + +!> Call-back routine for init_dim_obs +!! +!! This routine calls the observation-specific +!! routines init_dim_obs_TYPE. +!! + +! Author: Yorck Ewerdwalbesloh + +#ifdef CLMFIVE + SUBROUTINE init_dim_obs_pdafomi(step, dim_obs) + + use enkf_clm_mod, only: clmupdate_swc, clmupdate_tws + + ! Include functions for different observations + USE obs_GRACE_pdafomi, ONLY: assim_GRACE, init_dim_obs_GRACE + USE obs_SM_pdafomi, ONLY: assim_SM, init_dim_obs_SM + !USE obs_ST_pdafomi, ONLY: assim_C, init_dim_obs_C + + use mod_assimilation, only: screen + USE mod_parallel_pdaf, only: mype_world + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(out) :: dim_obs !< Dimension of full observation vector + + ! *** Local variables *** + INTEGER :: dim_obs_GRACE ! Observation dimensions + INTEGER :: dim_obs_SM ! Observation dimensions + !INTEGER :: dim_obs_C ! Observation dimensions + + + ! ********************************************* + ! *** Initialize full observation dimension *** + ! ********************************************* + + ! Initialize number of observations + dim_obs_GRACE = 0 + dim_obs_SM = 0 + !dim_obs_C = 0 + + + assim_SM = .true. + assim_GRACE = .true. + + + + ! Call observation-specific routines + ! The routines are independent, so it is not relevant + ! in which order they are called + + if (mype_world==0 .and. screen > 2) then + write(*,*)'Call dimension initialization' + end if + + IF (assim_GRACE) CALL init_dim_obs_GRACE(step, dim_obs_GRACE) + IF (assim_SM) CALL init_dim_obs_SM(step, dim_obs_SM) + !IF (assim_C) CALL init_dim_obs_C(step, dim_obs_C) + + dim_obs = dim_obs_GRACE + dim_obs_SM! + dim_obs_C + + END SUBROUTINE init_dim_obs_pdafomi + + + + !------------------------------------------------------------------------------- + !> Call-back routine for obs_op + !! + !! This routine calls the observation-specific + !! routines obs_op_TYPE. + !! + SUBROUTINE obs_op_pdafomi(step, dim_p, dim_obs, state_p, ostate) + + ! Include functions for different observations + USE obs_GRACE_pdafomi, ONLY: obs_op_GRACE + USE obs_SM_pdafomi, ONLY: obs_op_SM + !USE obs_C_pdafomi, ONLY: obs_op_C + + use mod_assimilation, only: screen + USE mod_parallel_pdaf, only: mype_world + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< Dimension of full observed state + REAL, INTENT(in) :: state_p(dim_p) !< PE-local model state + REAL, INTENT(inout) :: ostate(dim_obs) !< PE-local full observed state + + + ! ****************************************************** + ! *** Apply observation operator H on a state vector *** + ! ****************************************************** + + ! The order of these calls is not relevant as the setup + ! of the overall observation vector is defined by the + ! order of the calls in init_dim_obs_pdafomi + + if (mype_world==0 .and. screen > 2) then + write(*,*)'Call observation operators' + end if + + CALL obs_op_GRACE(dim_p, dim_obs, state_p, ostate) + CALL obs_op_SM(dim_p, dim_obs, state_p, ostate) + !CALL obs_op_C(dim_p, dim_obs, state_p, ostate) + + END SUBROUTINE obs_op_pdafomi + + + + !------------------------------------------------------------------------------- + !> Call-back routine for init_dim_obs_l + !! + !! This routine calls the routine PDAFomi_init_dim_obs_l + !! for each observation type + !! + SUBROUTINE init_dim_obs_l_pdafomi(domain_p, step, dim_obs, dim_obs_l) + + ! Include functions for different observations + USE obs_GRACE_pdafomi, ONLY: init_dim_obs_l_GRACE + USE obs_SM_pdafomi, ONLY: init_dim_obs_l_SM + !USE obs_C_pdafomi, ONLY: init_dim_obs_l_C + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: domain_p !< Index of current local analysis domain + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(in) :: dim_obs !< Full dimension of observation vector + INTEGER, INTENT(out) :: dim_obs_l !< Local dimension of observation vector + + + ! ********************************************** + ! *** Initialize local observation dimension *** + ! ********************************************** + + ! Call init_dim_obs_l specific for each observation + + CALL init_dim_obs_l_GRACE(domain_p, step, dim_obs, dim_obs_l) + CALL init_dim_obs_l_SM(domain_p, step, dim_obs, dim_obs_l) + !CALL init_dim_obs_l_C(domain_p, step, dim_obs, dim_obs_l) + + END SUBROUTINE init_dim_obs_l_pdafomi + + + + !------------------------------------------------------------------------------- + !> Call-back routine for localize_covar + !! + !! This routine calls the routine PDAFomi_localize_covar + !! for each observation type to apply covariance + !! localization in the LEnKF. + !! + SUBROUTINE localize_covar_pdafomi(dim_p, dim_obs, HP_p, HPH) + + ! Include functions for different observations + USE obs_GRACE_pdafomi, ONLY: localize_covar_GRACE + USE obs_SM_pdafomi, ONLY: localize_covar_SM + !USE obs_C_pdafomi, ONLY: localize_covar_C + + ! Include information on model grid + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< number of observations + REAL, INTENT(inout) :: HP_p(dim_obs, dim_p) !< PE local part of matrix HP + REAL, INTENT(inout) :: HPH(dim_obs, dim_obs) !< Matrix HPH + + ! *** local variables *** + INTEGER :: i, j, cnt ! Counters + REAL, ALLOCATABLE :: coords_p(:,:) ! Coordinates of PE-local state vector entries + + + ! ********************** + ! *** INITIALIZATION *** + ! ********************** + + ! Initialize coordinate array + + ALLOCATE(coords_p(2, dim_p)) + + + ! ************************************* + ! *** Apply covariance localization *** + ! ************************************* + + ! Call localize_covar specific for each observation + + CALL localize_covar_GRACE(dim_p, dim_obs, HP_p, HPH, coords_p) + CALL localize_covar_SM(dim_p, dim_obs, HP_p, HPH, coords_p) + !CALL localize_covar_C(dim_p, dim_obs, HP_p, HPH, coords_p) + + + ! **************** + ! *** Clean up *** + ! **************** + + DEALLOCATE(coords_p) + + END SUBROUTINE localize_covar_pdafomi + + + SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C) + + USE obs_GRACE_pdafomi, ONLY: add_obs_err_GRACE + USE obs_SM_pdafomi, ONLY: add_obs_err_SM + + IMPLICIT NONE + + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of obs. vector + REAL, INTENT(inout) :: C(dim_obs, dim_obs) ! Matrix to that the observation + ! error covariance matrix is added + + + INTEGER :: i ! index of observation component + REAL :: variance_obs ! variance of observations + CALL add_obs_err_GRACE(step, dim_obs, C) + CALL add_obs_err_SM(step, dim_obs, C) + !CALL add_obs_err_C(step, dim_obs, C) + + END SUBROUTINE add_obs_err_pdafomi + + + SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + + USE obs_GRACE_pdafomi, ONLY: init_obscovar_GRACE + USE obs_SM_pdafomi, ONLY: init_obscovar_SM + IMPLICIT NONE + + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector + REAL, INTENT(out) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix + REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PE-local observation vector + LOGICAL, INTENT(out) :: isdiag ! Whether the observation error covar. matrix is diagonal + + integer :: i + REAL :: variance_obs + + + CALL init_obscovar_GRACE(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + CALL init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + !CALL init_obscovar_C(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + + END SUBROUTINE init_obscovar_pdafomi + + + SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p) + + use obs_GRACE_pdafomi, ONLY: prodRinvA_GRACE + use obs_SM_pdafomi, ONLY: prodRinvA_SM + + implicit none + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_p(dim_obs_p) ! PE-local vector of observations + REAL, INTENT(in) :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine + REAL, INTENT(out) :: C_p(dim_obs_p,rank) ! Output matrix + + CALL prodRinvA_GRACE(step, dim_obs_p, rank, obs_p, A_p, C_p) + CALL prodRinvA_SM(step, dim_obs_p, rank, obs_p, A_p, C_p) + !CALL prodRinvA_C(step, dim_obs_p, rank, obs_p, A_p, C_p) + + END SUBROUTINE prodRinvA_pdafomi + + + SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + + use obs_GRACE_pdafomi, ONLY: prodRinvA_l_GRACE + use obs_SM_pdafomi, ONLY: prodRinvA_l_SM + + implicit none + + INTEGER, INTENT(in) :: domain_p ! Current local analysis domain + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs_l ! Dimension of local observation vector + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local vector of observations + REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine + REAL, INTENT(out) :: C_l(dim_obs_l, rank) ! Output matrix + + CALL prodRinvA_l_GRACE(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + CALL prodRinvA_l_SM(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + !CALL prodRinvA_l_C(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) + + END SUBROUTINE prodRinvA_l_pdafomi + +#endif diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index a0280f9c6..53ad99e7e 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -966,7 +966,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) ! Error if observation deeper than clmstatevec_max_layer if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." print *, "i=", i print *, "c=", c print *, "clmobs_layer(i)=", clmobs_layer(i) diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 6ce590e63..3631a9352 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -959,7 +959,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) ! Error if observation deeper than clmstatevec_max_layer if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." print *, "i=", i print *, "c=", c print *, "clmobs_layer(i)=", clmobs_layer(i) diff --git a/interface/framework/init_pdaf.F90 b/interface/framework/init_pdaf.F90 index 20fc01e84..53b52e238 100644 --- a/interface/framework/init_pdaf.F90 +++ b/interface/framework/init_pdaf.F90 @@ -80,6 +80,7 @@ SUBROUTINE init_pdaf() type_winf, limit_winf, & type_hyb, hyb_gamma, hyb_kappa, & pf_res_type, pf_noise_type, pf_noise_amp + USE mod_assimilation, ONLY: use_omi USE mod_tsmp, & ONLY: pf_statevecsize, nprocpf, tag_model_parflow, tag_model_clm, nprocclm, pf_statevec, pf_statevec_fortran, & idx_map_subvec2state, idx_map_subvec2state_fortran, model @@ -206,7 +207,9 @@ SUBROUTINE init_pdaf() #ifdef PDAF_DEBUG ! Debug output: local state dimension array - if (mype_model == 0) WRITE(*, '(a,x,a,i5,x,a,x)', advance="no") "TSMP-PDAF-debug", "mype(w)=", mype_world, "init_pdaf: dim_state_p_count in modified:" + if (mype_model == 0) WRITE(*, '(a,x,a,i5,x,a,x)', advance="no") & + "TSMP-PDAF-debug", "mype(w)=", mype_world, & + "init_pdaf: dim_state_p_count in modified:" if (mype_model == 0) WRITE(*, *) dim_state_p_count #endif @@ -251,6 +254,8 @@ SUBROUTINE init_pdaf() type_sqrt = 0 ! SEIK/LSEIK/ESTKF/LESTKF: Type of transform matrix square-root incremental = 0 ! SEIK/LSEIK: (1) to perform incremental updating + use_omi = .false. ! Default: Do not use OMI interface + !EnKF rank_analysis_enkf = 0 ! EnKF: rank to be considered for inversion of HPH in analysis step diff --git a/interface/framework/init_pdaf_parse.F90 b/interface/framework/init_pdaf_parse.F90 index 92c2dab4a..af8869f48 100644 --- a/interface/framework/init_pdaf_parse.F90 +++ b/interface/framework/init_pdaf_parse.F90 @@ -52,7 +52,17 @@ SUBROUTINE init_pdaf_parse() rms_obs, model_error, model_err_amp, incremental, type_forget, & forget, rank_analysis_enkf, locweight, cradius, & sradius, filename, type_trans, dim_obs, & - type_sqrt, obs_filename, dim_lag + type_sqrt, obs_filename, dim_lag, temp_mean_filename + USE mod_assimilation, ONLY: use_omi + + use mod_assimilation,& + only: cradius_GRACE, sradius_GRACE, & + cradius_SM, sradius_SM + +#ifdef CLMFIVE + use obs_GRACE_pdafomi, only: rms_obs_GRACE + use obs_SM_pdafomi, only: rms_obs_SM +#endif IMPLICIT NONE @@ -82,6 +92,16 @@ SUBROUTINE init_pdaf_parse() CALL parse(handle, toffset) handle = 'rms_obs' ! Assumed uniform RMS error of the observations CALL parse(handle, rms_obs) + +#ifdef CLMFIVE + rms_obs_GRACE = rms_obs ! backward compatibility + handle = 'rms_obs_GRACE' ! RMS error for GRACE observations + CALL parse(handle, rms_obs_GRACE) + rms_obs_SM = rms_obs ! backward compatibility + handle = 'rms_obs_SM' ! RMS error for SM observations + CALL parse(handle, rms_obs_SM) +#endif + handle = 'dim_obs' ! Number of observations CALL parse(handle, dim_obs) @@ -96,6 +116,8 @@ SUBROUTINE init_pdaf_parse() CALL parse(handle, subtype) handle = 'incremental' ! Set whether to use incremental updating CALL parse(handle, incremental) + handle = 'use_omi' ! Set whether to use OMI interface + CALL parse(handle, use_omi) ! Filter-specific settings handle = 'type_trans' ! Type of ensemble transformation in SEIK/ETKF/LSEIK/LETKF @@ -123,6 +145,20 @@ SUBROUTINE init_pdaf_parse() ! for 5th-order polynomial or radius for 1/e in exponential weighting CALL parse(handle, sradius) + ! Settings for different observation types + cradius_GRACE = cradius ! For backward compatibility + handle = 'cradius_GRACE' ! Set cut-off radius for GRACE observations + call parse(handle, cradius_GRACE) + sradius_GRACE = sradius ! For backward compatibility + handle = 'sradius_GRACE' ! Set support radius for GRACE observations + call parse(handle, sradius_GRACE) + cradius_SM = cradius ! For backward compatibility + handle = 'cradius_SM' ! Set cut-off radius for SM observations + call parse(handle, cradius_SM) + sradius_SM = sradius ! For backward compatibility + handle = 'sradius_SM' ! Set support radius for SM observations + call parse(handle, sradius_SM) + ! Setting for file output handle = 'filename' ! Set name of output file CALL parse(handle, filename) @@ -131,6 +167,10 @@ SUBROUTINE init_pdaf_parse() handle = 'obs_filename' call parse(handle, obs_filename) + ! *** Yorck: user defined filename for temporal mean of TWS to be subtracted in observation operator *** ! + handle = 'temp_mean_filename' + call parse(handle, temp_mean_filename) + !kuw: add smoother support handle = 'smoother_lag' call parse(handle, dim_lag) diff --git a/interface/framework/mod_assimilation.F90 b/interface/framework/mod_assimilation.F90 index 170b26175..ae1d957d2 100755 --- a/interface/framework/mod_assimilation.F90 +++ b/interface/framework/mod_assimilation.F90 @@ -65,12 +65,15 @@ MODULE mod_assimilation ! *** Variables specific for TSMP-PDAF *** ! gw - INTEGER, ALLOCATABLE :: dim_state_p_count(:) !Vector holding local state vector dimensions for processors of a single model communicator + ! Vector holding local state vector dimensions for processors of a single model communicator + INTEGER, ALLOCATABLE :: dim_state_p_count(:) ! gw end REAL, ALLOCATABLE :: obs(:) ! Vector holding all observations for Global domain INTEGER, ALLOCATABLE :: obs_index_l(:) ! Vector holding local state-vector indices of observations - INTEGER, ALLOCATABLE :: obs_interp_indices_p(:,:) ! Vector holding state-vector indices of grid cells surrounding interpolation for PE-local domain - INTEGER, ALLOCATABLE :: obs_interp_weights_p(:,:) ! Vector holding weights of grid cells surrounding observation for PE-local domain + ! Vector holding state-vector indices of grid cells surrounding interpolation for PE-local domain + INTEGER, ALLOCATABLE :: obs_interp_indices_p(:,:) + ! Vector holding weights of grid cells surrounding observation for PE-local domain + INTEGER, ALLOCATABLE :: obs_interp_weights_p(:,:) INTEGER, ALLOCATABLE :: local_dims_obs(:) ! Array for process-local observation dimensions INTEGER, ALLOCATABLE :: local_disp_obs(:) ! Observation displacement array for gathering. Displacement: #obs before current PE ! pdaf-ordered index: determined by domain-decomposition @@ -88,14 +91,29 @@ MODULE mod_assimilation REAL, ALLOCATABLE :: clm_obserr_p(:) ! Vector holding observation errors for CLM run at each PE-local domain REAL, ALLOCATABLE :: distance(:) ! Localization distance INTEGER, ALLOCATABLE :: global_to_local(:) ! Vector to map global index to local domain index - INTEGER, ALLOCATABLE :: longxy(:), latixy(:), longxy_obs(:), latixy_obs(:) ! longitude and latitude of grid cells and observation cells - INTEGER, ALLOCATABLE :: longxy_obs_floor(:), latixy_obs_floor(:) ! indices of grid cells with smaller lon/lat than observation location + ! longitude and latitude of grid cells and observation cells + INTEGER, ALLOCATABLE :: longxy(:), latixy(:), longxy_obs(:), latixy_obs(:) + ! indices of grid cells with smaller lon/lat than observation location + INTEGER, ALLOCATABLE :: longxy_obs_floor(:), latixy_obs_floor(:) INTEGER, ALLOCATABLE :: var_id_obs(:) ! for remote sensing data the variable identifier to group ! variables distributed over a grid surface area !kuw INTEGER, ALLOCATABLE :: obs_nc2pdaf(:) ! index for mapping mstate to local domain !kuw end + + ! Yorck + + ! interval until next observation, used by next_observation_pdaf.F90, + ! better solution for next assimilation time step + REAL :: da_interval_variable + ! has to be read from observation file --> no empty observation files have to be written + REAL, ALLOCATABLE :: obscov(:,:) ! observation covariance matrix + REAL, ALLOCATABLE :: obscov_inv(:,:) ! inverse of the observation covariance matrix + character (len = 110) :: temp_mean_filename ! User defined filename of temporal mean + + ! END Yorck + ! Multi-scale DA ! store the maximum and minimum limits for remote sensing data with @@ -134,6 +152,10 @@ MODULE mod_assimilation ! ! Settings for observations - available as command line options INTEGER :: delt_obs ! time step interval between assimilation steps REAL :: rms_obs ! RMS error size for observation generation + + REAL :: rms_obs_GRACE + REAL :: rms_obs_SM + INTEGER :: dim_obs ! Number of observations ! ! General control of PDAF - available as command line options @@ -198,6 +220,7 @@ MODULE mod_assimilation ! (6) hybrid 3D-Var using LESTKF for ensemble update ! (7) hybrid 3D-Var using ESTKF for ensemble update INTEGER :: incremental ! Perform incremental updating in LSEIK + LOGICAL :: use_omi ! Set whether OMI interface is used INTEGER :: dim_lag ! Number of time instances for smoother ! ! Filter settings - available as command line options @@ -251,6 +274,12 @@ MODULE mod_assimilation ! (2) 5th-order polynomial weight function REAL :: sradius ! Support radius for 5th order polynomial ! or radius for 1/e for exponential weighting + + + REAL :: cradius_GRACE + REAL :: sradius_GRACE + REAL :: cradius_SM + REAL :: sradius_SM ! ! SEIK-subtype4/LSEIK-subtype4/ESTKF/LESTKF INTEGER :: type_sqrt ! Type of the transform matrix square-root ! (0) symmetric square root diff --git a/interface/framework/mod_read_obs.F90 b/interface/framework/mod_read_obs.F90 index 9c202acda..fd0824203 100755 --- a/interface/framework/mod_read_obs.F90 +++ b/interface/framework/mod_read_obs.F90 @@ -68,6 +68,172 @@ module mod_read_obs real, allocatable :: dampfac_param_time_dependent_in(:) contains + + !> @author Yorck Ewerdwalbesloh + !> @date 17.03.2025 + !> @brief Read NetCDF observation file for different observation types + !! to read only those observations that should be read in for the specified type + !> @param[in] Name of observation file, Name of observation type + !> @param[inout] Full observation dimension, full observation vector, uncertainty information, coordinates (lon and lat) + !> @details + !> This subroutine reads the observation file and return the data + + subroutine read_obs_nc_type(current_observation_filename, & + current_observation_type, dim_obs_g, obs_g, & + lon_obs_g, lat_obs_g, layer_obs_g, & + dr_obs_g, obserr_g, obscov_g) + use netcdf, only: nf90_max_name + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_noerr + use netcdf, only: nf90_get_var + use mod_assimilation, only: screen + implicit none + + integer :: ncid + integer :: dimid, status + integer :: haserr + + character (len = *), intent(in) :: current_observation_filename + character (len = *), intent(in) :: current_observation_type + INTEGER, INTENT(inout) :: dim_obs_g !< Dimension of full observation vector + REAL, allocatable, INTENT(inout) :: obs_g(:) + REAL, allocatable, INTENT(inout) :: lon_obs_g(:) + REAL, allocatable, INTENT(inout) :: lat_obs_g(:) + INTEGER, allocatable, INTENT(inout) :: layer_obs_g(:) + REAL, allocatable, INTENT(inout) :: dr_obs_g(:) + REAL, allocatable, INTENT(inout) :: obserr_g(:) + REAL, allocatable, INTENT(inout) :: obscov_g(:,:) + + integer :: clmobs_varid, dr_varid, clmobs_lon_varid, clmobs_lat_varid, & + clmobs_layer_varid, clmobserr_varid, clmobscov_varid, obstype_varid + + character (len = *), parameter :: dim_name = "dim_obs" + character (len = *), parameter :: obs_name = "obs_clm" + character (len = *), parameter :: dr_name = "dr" + character (len = *), parameter :: lon_name = "lon" + character (len = *), parameter :: lat_name = "lat" + character (len = *), parameter :: layer_name = "layer" + character (len = *), parameter :: obserr_name = "obserr_clm" + character (len = *), parameter :: obscov_name = "obscov_clm" + character (len = *), parameter :: type_name = "type_clm" + character(len = nf90_max_name) :: RecordDimName + + character(len=20), allocatable :: obs_type(:) + integer, allocatable :: indices(:) + + integer :: has_obs_clm, dim_obs + + + + call check( nf90_open(current_observation_filename, nf90_nowrite, ncid) ) + call check(nf90_inq_dimid(ncid, dim_name, dimid)) + call check(nf90_inquire_dimension(ncid, dimid, recorddimname, dim_obs)) + + has_obs_clm = nf90_inq_varid(ncid, obs_name, clmobs_varid) + + if(has_obs_clm == nf90_noerr) then + + if(allocated(obs_type)) deallocate(obs_type) + allocate(obs_type(dim_obs)) + + call check(nf90_inq_varid(ncid, type_name, obstype_varid)) + call check(nf90_get_var(ncid, obstype_varid, obs_type)) + + if (trim(obs_type(1)) /= trim(current_observation_type)) then + + dim_obs_g = 0 + + if(allocated(obs_g)) deallocate(obs_g) + allocate(obs_g(dim_obs_g)) + + if(allocated(lon_obs_g)) deallocate(lon_obs_g) + allocate(lon_obs_g(dim_obs_g)) + + if(allocated(lat_obs_g)) deallocate(lat_obs_g) + allocate(lat_obs_g(dim_obs_g)) + + if(allocated(layer_obs_g)) deallocate(layer_obs_g) + allocate(layer_obs_g(dim_obs_g)) + + if(allocated(dr_obs_g)) deallocate(dr_obs_g) + allocate(dr_obs_g(dim_obs_g)) + + else + + if(allocated(obs_g)) deallocate(obs_g) + allocate(obs_g(dim_obs)) + + if(allocated(lon_obs_g)) deallocate(lon_obs_g) + allocate(lon_obs_g(dim_obs)) + + if(allocated(lat_obs_g)) deallocate(lat_obs_g) + allocate(lat_obs_g(dim_obs)) + + if(allocated(layer_obs_g)) deallocate(layer_obs_g) + allocate(layer_obs_g(dim_obs)) + + if(allocated(dr_obs_g)) deallocate(dr_obs_g) + allocate(dr_obs_g(dim_obs)) + + call check(nf90_get_var(ncid, clmobs_varid, obs_g)) + + call check( nf90_inq_varid(ncid, lon_name, clmobs_lon_varid) ) + call check( nf90_get_var(ncid, clmobs_lon_varid, lon_obs_g) ) + + call check( nf90_inq_varid(ncid, lat_name, clmobs_lat_varid) ) + call check( nf90_get_var(ncid, clmobs_lat_varid, lat_obs_g) ) + + haserr = nf90_inq_varid(ncid, layer_name, clmobs_layer_varid) + if (haserr == nf90_noerr) then + call check( nf90_get_var(ncid, clmobs_layer_varid, layer_obs_g) ) + else + layer_obs_g(:) = 1 + end if + + call check( nf90_inq_varid(ncid, dr_name, dr_varid) ) + call check( nf90_get_var(ncid, dr_varid, dr_obs_g) ) + + haserr = nf90_inq_varid(ncid, obserr_name, clmobserr_varid) + + if(haserr == nf90_noerr) then + + multierr = 1 + + if(allocated(obserr_g)) deallocate(obserr_g) + allocate(obserr_g(dim_obs)) + + call check(nf90_get_var(ncid, clmobserr_varid, obserr_g)) + + end if + + haserr = nf90_inq_varid(ncid, obscov_name, clmobscov_varid) + if(haserr == nf90_noerr) then + + multierr = 2 + + if(allocated(obscov_g)) deallocate(obscov_g) + allocate(obscov_g(dim_obs,dim_obs)) + + call check(nf90_get_var(ncid, clmobscov_varid, obscov_g)) + + end if + + dim_obs_g = dim_obs + + end if + + end if + + end subroutine read_obs_nc_type + + + + + !> @author Wolfgang Kurtz, Guowei He, Mukund Pondkule !> @date 03.03.2023 !> @brief Read NetCDF observation file @@ -604,7 +770,7 @@ subroutine check_n_observationfile(fn,nn) end subroutine check_n_observationfile - !> @author Yorck Ewerdwalbesloh, Johannes Keller + !> @author Anne Springer, Yorck Ewerdwalbesloh, Johannes Keller !> @date 11.09.2023 !> @brief Return data assimilation interval from file !> @param[in] fn Filename of the observation file @@ -621,6 +787,7 @@ subroutine check_n_observationfile_da_interval(fn,aa) use shr_kind_mod, only: r8 => shr_kind_r8 use netcdf, only: nf90_max_name, nf90_open, nf90_nowrite, & nf90_inq_varid, nf90_get_var, nf90_close, nf90_noerr + use mod_assimilation, only: use_omi implicit none @@ -629,9 +796,12 @@ subroutine check_n_observationfile_da_interval(fn,aa) integer :: ncid, varid, status !,dimid - character (len = *), parameter :: varname = "da_interval" + character (len = nf90_max_name) :: varname real(r8) :: dtime ! land model time step (sec) + varname = "da_interval " + + !character (len = *), parameter :: dim_name = "dim_obs" !character(len = nf90_max_name) :: recorddimname @@ -665,4 +835,350 @@ subroutine check(status) end subroutine check +#ifdef CLMFIVE + !> @author Anne Springer, adaptation for TSMP2 by Yorck Ewerdwalbesloh + !> @date 04.12.2023 + !> @brief Return set zero interval for running mean of model variables from file + !> @param[in] fn Filename of the observation file + !> @param[out] nn number of hours until setting zero + !> @details + !> Reads the content of the variable name `set_zero` from NetCDF + !> file `fn` using subroutines from the NetCDF module. + !> The result is returned in `nn`. + !> + !> The result is used to reset the running average of state variables. + subroutine check_n_observationfile_set_zero(fn,nn) + use shr_kind_mod, only: r8 => shr_kind_r8 + use netcdf, only: nf90_max_name, nf90_open, nf90_nowrite, & + nf90_inq_varid, nf90_get_var, nf90_close, nf90_noerr + use clm_varcon, only: ispval + use clm_time_manager, only : get_step_size + + implicit none + + character(len=*),intent(in) :: fn + integer, intent(out) :: nn + + integer :: ncid, varid, status !,dimid + character (len = *), parameter :: varname = "set_zero" + real(r8) :: dtime ! land model time step (sec) + + !character (len = *), parameter :: dim_name = "dim_obs" + !character(len = nf90_max_name) :: recorddimname + + dtime = get_step_size() + + call check(nf90_open(fn, nf90_nowrite, ncid)) + !call check(nf90_inq_dimid(ncid, dim_name, dimid)) + !call check(nf90_inquire_dimension(ncid, dimid, recorddimname, nn)) + status = nf90_inq_varid(ncid, varname, varid) + if (status == nf90_noerr) then + call check(nf90_inq_varid(ncid, varname, varid)) + call check( nf90_get_var(ncid, varid, nn) ) + call check(nf90_close(ncid)) + ! at this point: half hourly time steps, this is adjusted here. In the GRACE files, set_zero is set up as hours + ! --> is adjusted using information from inside CLM + if (nn/=ispval) then + nn = nn*INT(3600/dtime) + end if + else + nn = ispval + end if + + end subroutine check_n_observationfile_set_zero +#endif + + !> @author Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Return next observation type from next observation file + !> @param[in] fn Filename of the observation file + !> @param[out] obs_type_str next observation type + !> @details + !> Reads the content of the variable name `type_clm` from NetCDF + !> file `fn` using subroutines from the NetCDF module. + !> The first entry of the vector is returned in `obs_type_str`. + !> + !> The result is used to reset the observation type and to initialize the next assimilation cycle. + subroutine check_n_observationfile_next_type(fn, obs_type_str) + use netcdf, only: nf90_max_name + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_noerr + use netcdf, only: nf90_get_var + use netcdf, only: nf90_close + use mod_assimilation, only: screen + implicit none + + character(len=*), intent(in) :: fn + character(len=*), intent(out) :: obs_type_str + + integer :: ncid, status, obstype_varid, dimid + integer :: dim_obs + character (len = *), parameter :: dim_name = "dim_obs" + character (len = *), parameter :: type_name = "type_clm" + character(len=20), allocatable :: obs_type_lok(:) + character(len = nf90_max_name) :: RecordDimName + + + call check( nf90_open(fn, nf90_nowrite, ncid) ) + call check(nf90_inq_dimid(ncid, dim_name, dimid)) + call check(nf90_inquire_dimension(ncid, dimid, recorddimname, dim_obs)) + + if(allocated(obs_type_lok)) deallocate(obs_type_lok) + allocate(obs_type_lok(dim_obs)) + + obs_type_str = '' + + status = nf90_inq_varid(ncid, "type_clm", obstype_varid) + if (status == nf90_noerr) then + call check(nf90_inq_varid(ncid, "type_clm", obstype_varid)) + call check(nf90_get_var(ncid, obstype_varid, obs_type_lok)) + obs_type_str = trim(obs_type_lok(1)) + end if + call check(nf90_close(ncid)) + + if(allocated(obs_type_lok)) deallocate(obs_type_lok) + + end subroutine check_n_observationfile_next_type + + +#ifdef CLMFIVE + !> @author Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Update observation type for next assimilation cycle + !> @param[in] obs_type_st next observation type + !> @details + !> Updates the observation type for the next assimilation cycle when using the OMI interface + subroutine update_obs_type(obs_type_str) + use enkf_clm_mod, only: clmupdate_tws, clmupdate_swc, clmupdate_T, clmupdate_texture + use mod_parallel_pdaf, only: abort_parallel + implicit none + + character(len=*), intent(in) :: obs_type_str + + select case (trim(adjustl(obs_type_str))) + case ('GRACE') + clmupdate_tws = 1 + clmupdate_swc = 0 + clmupdate_T = 0 + clmupdate_texture = 0 + + case ('SM') + clmupdate_tws = 0 + clmupdate_swc = 1 + clmupdate_T = 0 + clmupdate_texture = 0 + + case default + write(*,*) 'ERROR: Unknown obs_type_str in update_obs_type:', trim(obs_type_str) + call abort_parallel() + end select + end subroutine update_obs_type + + + !> @author Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Defines an index domain for GRACE assimilation + !> @param[in] lon_clmobs longitude of the observations + !> @param[in] lat_clmobs latitude of the observations + !> @param[in] dim_obs number of observations + !> @param[out] longxy local x-indices of the gridcells + !> @param[out] latixy local y-indices of the gridcells + !> @param[out] longxy_obs local x-indices of the observations + !> @param[out] latixy_obs local y-indices of the observations + !> @details + !> Generates an index grid instead of the lon/lat grid + subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, & + longxy, latixy, longxy_obs, latixy_obs) + + use mpi, only: MPI_INTEGER + use mpi, only: MPI_DOUBLE_PRECISION + use mpi, only: MPI_IN_PLACE + use mpi, only: MPI_SUM + use mpi, only: MPI_2INTEGER + use mpi, only: MPI_MINLOC + + use spmdMod, only : npes, iam + use domainMod, only : ldomain, lon1d, lat1d + use decompMod, only : get_proc_total, get_proc_bounds, ldecomp + use GridcellType, only: grc + use shr_kind_mod, only: r8 => shr_kind_r8 + use enkf_clm_mod, only: hactiveg_levels, num_hactiveg + !USE mod_parallel_pdaf, & + ! ONLY: mpi_2integer, mpi_minloc + USE mod_parallel_pdaf, & + ONLY: comm_filter, npes_filter, abort_parallel, & + mype_world, mype_filter + real, intent(in) :: lon_clmobs(:) + real, intent(in) :: lat_clmobs(:) + integer, intent(in) :: dim_obs + integer, allocatable, intent(inout) :: longxy(:) + integer, allocatable, intent(inout) :: latixy(:) + integer, allocatable, intent(inout) :: longxy_obs(:) + integer, allocatable, intent(inout) :: latixy_obs(:) + integer :: ni, nj, ii, jj, kk, cid, ier, ncells, nlunits, & + ncols, npatches, ncohorts, counter, i, g, ll + real :: minlon, minlat, maxlon, maxlat + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + real(r8), allocatable :: longxy_obs_lokal(:), latixy_obs_lokal(:) + + INTEGER, allocatable :: in_mpi_(:,:), out_mpi_(:,:) + + integer :: begg, endg ! per-proc gridcell ending gridcell indices + + real(r8) :: lat1, lon1, lat2, lon2, a, c, R, pi + + real(r8) :: dist + real(r8), allocatable :: min_dist(:) + integer, allocatable :: min_g(:) + + integer :: ierror + + integer :: lok_lon, lok_lat + + + + lon => grc%londeg + lat => grc%latdeg + + + ni = ldomain%ni + nj = ldomain%nj + + ! get total number of gridcells, landunits, + ! columns, patches and cohorts on processor + + call get_proc_total(iam, ncells, nlunits, ncols, npatches, ncohorts) + + ! beg and end gridcell + call get_proc_bounds(begg=begg, endg=endg) + + if (allocated(longxy)) deallocate(longxy) + if (allocated(latixy)) deallocate(latixy) + allocate(longxy(num_hactiveg), stat=ier) + allocate(latixy(num_hactiveg), stat=ier) + + + longxy(:) = 0 + latixy(:) = 0 + + + counter = 1 + do ii = 1, nj + do jj = 1, ni + cid = (ii-1)*ni + jj + do ll = 1, num_hactiveg + kk = hactiveg_levels(ll,1) + if(cid == ldecomp%gdc2glo(kk)) then + latixy(counter) = ii + longxy(counter) = jj + counter = counter + 1 + end if + end do + end do + end do + + if (allocated(min_dist)) deallocate(min_dist) + allocate(min_dist(dim_obs)) + min_dist(:) = huge(min_dist) + + if (allocated(min_g)) deallocate(min_g) + allocate(min_g(dim_obs)) + + R = 6371.0 + pi = 3.14159265358979323846 + do i = 1, dim_obs + do g = begg, endg + + ! check distance from each grid point to observation location --> take the coordinate in local system that equals + ! the one of the closest coordinate + lat1 = lat(g) * pi / 180.0 + lon1 = lon(g) * pi / 180.0 + lat2 = lat_clmobs(i) * pi / 180.0 + lon2 = lon_clmobs(i) * pi / 180.0 + + a = sin((lat2 - lat1) / 2)**2 + cos(lat1) * cos(lat2) * sin((lon2 - lon1) / 2)**2 + c = 2 * atan2(sqrt(a), sqrt(1 - a)) + dist = R * c + + if (dist < min_dist(i)) then + min_dist(i) = dist + min_g(i) = g + end if + end do + end do + + + IF (ALLOCATED(in_mpi_)) DEALLOCATE(in_mpi_) + ALLOCATE(in_mpi_(2,dim_obs)) + IF (ALLOCATED(out_mpi_)) DEALLOCATE(out_mpi_) + ALLOCATE(out_mpi_(2,dim_obs)) + + in_mpi_(1,:) = int(ceiling(min_dist)) + in_mpi_(2,:) = min_g + + if (allocated(longxy_obs_lokal)) deallocate(longxy_obs_lokal) + if (allocated(latixy_obs_lokal)) deallocate(latixy_obs_lokal) + + allocate(longxy_obs_lokal(dim_obs)) + allocate(latixy_obs_lokal(dim_obs)) + + do i =1, dim_obs + outer: do ii = 1, nj + do jj = 1, ni + cid = (ii-1)*ni + jj + do kk = begg, endg + if (kk == in_mpi_(2,i)) then + if(cid == ldecomp%gdc2glo(kk)) then + if (min_dist(i)<30) then + latixy_obs_lokal(i) = ii + longxy_obs_lokal(i) = jj + else + longxy_obs_lokal(i) = -9999 + latixy_obs_lokal(i) = -9999 + end if + exit outer + end if + end if + end do + end do + end do outer + end do + + + if (allocated(longxy_obs)) deallocate(longxy_obs) + if (allocated(latixy_obs)) deallocate(latixy_obs) + allocate(longxy_obs(dim_obs), stat=ier) + allocate(latixy_obs(dim_obs), stat=ier) + + in_mpi_(2,:) = longxy_obs_lokal + call mpi_allreduce(in_mpi_,out_mpi_, dim_obs, mpi_2integer, mpi_minloc, comm_filter, ierror) + longxy_obs(:) = out_mpi_(2,:) + + in_mpi_(2,:) = latixy_obs_lokal + call mpi_allreduce(in_mpi_,out_mpi_, dim_obs, mpi_2integer, mpi_minloc, comm_filter, ierror) + latixy_obs(:) = out_mpi_(2,:) + + deallocate(longxy_obs_lokal) + deallocate(latixy_obs_lokal) + deallocate(in_mpi_) + deallocate(out_mpi_) + deallocate(min_dist) + deallocate(min_g) + + + if (mype_filter == 0) then + print*, "longxy_obs = ", longxy_obs + print*, "latixy_obs = ", latixy_obs + end if + + end subroutine domain_def_clm +#endif + + end module mod_read_obs diff --git a/interface/framework/next_observation_pdaf.F90 b/interface/framework/next_observation_pdaf.F90 index 37b78f8a4..f11c64d91 100644 --- a/interface/framework/next_observation_pdaf.F90 +++ b/interface/framework/next_observation_pdaf.F90 @@ -52,7 +52,7 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! ! !USES: USE mod_assimilation, & - ONLY: delt_obs, toffset, screen + ONLY: delt_obs, toffset, screen, da_interval_variable USE mod_parallel_pdaf, & ONLY: mype_world USE mod_tsmp, & @@ -62,9 +62,17 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) USE mod_tsmp, ONLY: flexible_da_interval USE mod_assimilation, & ONLY: obs_filename + USE mod_assimilation, ONLY: use_omi use mod_read_obs, & - only: check_n_observationfile - use mod_read_obs, ONLY: check_n_observationfile_da_interval + only: check_n_observationfile, check_n_observationfile_da_interval, check_n_observationfile_set_zero, & + check_n_observationfile_next_type +#ifdef CLMFIVE + use mod_read_obs, only: update_obs_type + use clm_time_manager, & + only: get_nstep + use clm_varcon, only: set_averaging_to_zero, ispval + use enkf_clm_mod, only: clmupdate_tws +#endif IMPLICIT NONE ! !ARGUMENTS: @@ -80,7 +88,10 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) !kuw: local variables integer :: counter integer :: no_obs + integer :: nstep character (len = 110) :: fn + character(len=32) :: obs_type_str + logical :: file_exists !kuw end REAL :: da_interval_new @@ -185,7 +196,6 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) - ! IF (stepnow + nsteps <= total_steps) THEN ! if (2<1) then ! ! *** During the assimilation process *** @@ -213,6 +223,52 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! doexit = ?? !print *, "next_observation_pdaf finished" +#ifdef CLMSA +#ifdef CLMFIVE + OMI:if (use_omi) then + if (clmupdate_tws/=0) then ! only update set_zero when GRACE is assimilated at the current time step + nstep = get_nstep() + if (stepnow/=toffset) then + write(fn, '(a, i5.5)') trim(obs_filename)//'.', stepnow + call check_n_observationfile_set_zero(fn, set_averaging_to_zero) + if (set_averaging_to_zero/=ispval) then + set_averaging_to_zero = set_averaging_to_zero+nstep + end if + + if (mype_world==0 .and. screen > 2) then + write(*,*) 'set_averaging_to_zero (in next_observation_pdaf):',set_averaging_to_zero + end if + end if + end if + + ! update observation type with next file + write(fn, '(a, i5.5)') trim(obs_filename)//'.', stepnow + delt_obs + if (mype_world==0 .and. screen > 2) then + write(*,*)'next_observation_pdaf: fn = ', fn + write(*,*)'Call check_n_observationfile_next_type' + end if + + inquire(file=fn, exist=file_exists) + if (.not. file_exists) then + if (mype_world == 0 .and. screen > 2) then + write(*,*) 'next_observation_pdaf: skipping setting next observation type as no next file available' + end if + else + call check_n_observationfile_next_type(fn, obs_type_str) + if (trim(obs_type_str) /= '') then + call update_obs_type(obs_type_str) + end if + + if (mype_world==0 .and. screen > 2) then + write(*,*)'next_type (in next_observation_pdaf):',trim(obs_type_str) + end if + end if + + end if OMI +#endif +#endif + END SUBROUTINE next_observation_pdaf + diff --git a/interface/framework/obs_GRACE_pdafomi.F90 b/interface/framework/obs_GRACE_pdafomi.F90 new file mode 100644 index 000000000..f8318a1e7 --- /dev/null +++ b/interface/framework/obs_GRACE_pdafomi.F90 @@ -0,0 +1,1347 @@ +!> PDAF-OMI observation module for type GRACE observations +!! +!! This module handles operations for one data type (called 'module-type' below): +!! OBSTYPE = GRACE +!! +!! __Observation type GRACE:__ +!! The observation type GRACE are TWSA observations (gridded) +!! +!! The subroutines in this module are for the particular handling of +!! a single observation type. +!! The routines are called by the different call-back routines of PDAF +!! usually by callback_obs_pdafomi.F90 +!! Most of the routines are generic so that in practice only 2 routines +!! need to be adapted for a particular data type. These are the routines +!! for the initialization of the observation information (init_dim_obs) +!! and for the observation operator (obs_op). +!! +!! The module and the routines are named according to the observation type. +!! This allows to distinguish the observation type and the routines in this +!! module from other observation types. +!! +!! The module uses two derived data types (obs_f and obs_l), which contain +!! all information about the full and local observations. Only variables +!! of the type obs_f need to be initialized in this module. The variables +!! in the type obs_l are initilized by the generic routines from PDAFomi. +!! +!! +!! These 2 routines need to be adapted for the particular observation type: +!! * init_dim_obs_OBSTYPE \n +!! Count number of process-local and full observations; +!! initialize vector of observations and their inverse variances; +!! initialize coordinate array and index array for indices of +!! observed elements of the state vector. +!! * obs_op_OBSTYPE \n +!! observation operator to get full observation vector of this type. Here +!! one has to choose a proper observation operator or implement one. +!! +!! In addition, there are two optional routines, which are required if filters +!! with localization are used: +!! * init_dim_obs_l_OBSTYPE \n +!! Only required if domain-localized filters (e.g. LESTKF, LETKF) are used: +!! Count number of local observations of module-type according to +!! their coordinates (distance from local analysis domain). Initialize +!! module-internal distances and index arrays. +!! * localize_covar_OBSTYPE \n +!! Only required if the localized EnKF is used: +!! Apply covariance localization in the LEnKF. +!! +!! __Revision history:__ +!! * 2019-06 - Lars Nerger - Initial code +!! * Later revisions - see repository log +!! +#ifdef CLMFIVE +MODULE obs_GRACE_pdafomi + + USE mod_parallel_pdaf, & + ONLY: mype_filter ! Rank of filter process + USE PDAFomi, & + ONLY: obs_f, obs_l ! Declaration of observation data types + + IMPLICIT NONE + SAVE + + PUBLIC + + ! Variables which are inputs to the module (usually set in init_pdaf) + LOGICAL :: assim_GRACE !< Whether to assimilate this data type + REAL :: rms_obs_GRACE !< Observation error standard deviation (for constant errors) + logical, allocatable :: vec_useObs(:) + ! vector of number of points for each GRACE observation, same dimension as observation vector + integer, allocatable :: vec_numPoints_global(:) + ! vector that tells if an observation of used (1) or not (0), same dimension as observation vector, global + logical, allocatable :: vec_useObs_global(:) + + real, allocatable :: tws_temp_mean(:,:) ! temporal mean for TWS + real, allocatable :: tws_temp_mean_d(:) ! temporal mean for TWS, vectorized with the same bounds as local process + real, allocatable :: lon_temp_mean(:,:) ! corresponding longitude + real, allocatable :: lat_temp_mean(:,:) ! corresponding latitude + + ! longitude and latitude of grid cells and observation cells + INTEGER, ALLOCATABLE :: longxy(:), latixy(:), longxy_obs(:), latixy_obs(:) + + ! One can declare further variables, e.g. for file names which can + ! be use-included in init_pdaf() and initialized there. + + + ! ********************************************************* + ! *** Data type obs_f defines the full observations by *** + ! *** internally shared variables of the module *** + ! ********************************************************* + + ! Relevant variables that can be modified by the user: + ! TYPE obs_f + ! ---- Mandatory variables to be set in INIT_DIM_OBS ---- + ! INTEGER :: doassim ! Whether to assimilate this observation type + ! INTEGER :: disttype ! Type of distance computation to use for localization + ! ! (0) Cartesian, (1) Cartesian periodic + ! ! (2) simplified geographic, (3) geographic haversine function + ! INTEGER :: ncoord ! Number of coordinates use for distance computation + ! INTEGER, ALLOCATABLE :: id_obs_p(:,:) ! Indices of observed field in state vector (process-local) + ! + ! ---- Optional variables - they can be set in INIT_DIM_OBS ---- + ! REAL, ALLOCATABLE :: icoeff_p(:,:) ! Interpolation coefficients for obs. operator + ! REAL, ALLOCATABLE :: domainsize(:) ! Size of domain for periodicity (<=0 for no periodicity) + ! + ! ---- Variables with predefined values - they can be changed in INIT_DIM_OBS ---- + ! INTEGER :: obs_err_type=0 ! Type of observation error: (0) Gauss, (1) Laplace + ! INTEGER :: use_global_obs=1 ! Whether to use (1) global full obs. + ! ! or (0) obs. restricted to those relevant for a process domain + ! REAL :: inno_omit=0.0 ! Omit obs. if squared innovation larger this factor times + ! ! observation variance + ! REAL :: inno_omit_ivar=1.0e-12 ! Value of inverse variance to omit observation + ! END TYPE obs_f + + ! Data type obs_l defines the local observations by internally shared variables of the module + + ! *********************************************************************** + + ! Declare instances of observation data types used here + ! We use generic names here, but one could rename the variables + TYPE(obs_f), TARGET, PUBLIC :: thisobs ! full observation + TYPE(obs_l), TARGET, PUBLIC :: thisobs_l ! local observation + + !$OMP THREADPRIVATE(thisobs_l) + + + !------------------------------------------------------------------------------- + + CONTAINS + + !> Initialize information on the module-type observation + !! + !! The routine is called by each filter process. + !! at the beginning of the analysis step before + !! the loop through all local analysis domains. + !! + !! It has to count the number of observations of the + !! observation type handled in this module according + !! to the current time step for all observations + !! required for the analyses in the loop over all local + !! analysis domains on the PE-local state domain. + !! + !! The following four variables have to be initialized in this routine + !! * thisobs\%doassim - Whether to assimilate this type of observations + !! * thisobs\%disttype - type of distance computation for localization with this observaton + !! * thisobs\%ncoord - number of coordinates used for distance computation + !! * thisobs\%id_obs_p - array with indices of module-type observation in process-local state vector + !! + !! Optional is the use of + !! * thisobs\%icoeff_p - Interpolation coefficients for obs. operator (only if interpolation is used) + !! * thisobs\%domainsize - Size of domain for periodicity for disttype=1 (<0 for no periodicity) + !! * thisobs\%obs_err_type - Type of observation errors for particle filter and NETF (default: 0=Gaussian) + !! * thisobs\%use_global obs - Whether to use global observations or restrict the observations to the relevant ones + !! (default: 1=use global full observations) + !! * thisobs\%inno_omit - Omit obs. if squared innovation larger this factor times observation variance + !! (default: 0.0, omission is active if >0) + !! * thisobs\%inno_omit_ivar - Value of inverse variance to omit observation + !! (default: 1.0e-12, change this if this value is not small compared to actual obs. error) + !! + !! Further variables are set when the routine PDAFomi_gather_obs is called. + !! + + + !> @author Anne Springer, Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Initialized the observation dimension and related arrays for GRACE observations + !> @param[in] step current time step + !> @param[in,out] dim_obs dimension of full observation vector + SUBROUTINE init_dim_obs_GRACE(step, dim_obs) + + USE mpi, ONLY: MPI_INTEGER + USE mpi, ONLY: MPI_SUM + USE mpi, ONLY: MPI_2INTEGER + USE mpi, ONLY: MPI_MAXLOC + USE mpi, ONLY: MPI_IN_PLACE + + USE mod_parallel_pdaf, & + ONLY: mype_filter, comm_filter, npes_filter, abort_parallel, & + mype_world + + USE PDAFomi, & + ONLY: PDAFomi_gather_obs + USE mod_assimilation, & + ONLY: filtertype, cradius_GRACE, obs_filename, temp_mean_filename, screen, obscov, obscov_inv + + use mod_assimilation, only: obs_nc2pdaf, obs_pdaf2nc, & + local_dims_obs, & + local_disp_obs + + use mod_read_obs, only: read_obs_nc_type, domain_def_clm, multierr + + use enkf_clm_mod, only: num_layer, hactiveg_levels + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use GridcellType, only: grc + + use clm_varcon, only: spval + + USE mod_parallel_pdaf, & + ONLY: mype_world + + use decompMod , only : get_proc_bounds + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(inout) :: dim_obs !< Dimension of full observation vector + + ! *** Local variables *** + INTEGER :: i, j, count_points, c, l, k ! Counters + INTEGER :: cnt, cnt0 ! Counters + INTEGER :: dim_obs_p ! Number of process-local observations + REAL, ALLOCATABLE :: obs_field(:,:) ! Observation field read from file + REAL, ALLOCATABLE :: obs_p(:) ! PE-local observation vector + REAL, ALLOCATABLE :: obs_g(:) ! Global observation vector + REAL, ALLOCATABLE :: ivar_obs_p(:) ! PE-local inverse observation error variance + REAL, ALLOCATABLE :: ocoord_p(:,:) ! PE-local observation coordinates + REAL, ALLOCATABLE :: clm_obscov(:,:) ! full observation error covariance matrix before removing observations that cannot be seen by enough gridcells + CHARACTER(len=2) :: stepstr ! String for time step + character (len = 110) :: current_observation_filename + + + character(len=20) :: obs_type_name ! name of observation type (e.g. GRACE, SM, ST, ...) + + integer :: numPoints ! minimum number of points so that the GRACE observation is used + integer, allocatable :: vec_numPoints(:) ! number of model grid cells that are in a radius of dr around the GRACE observation + INTEGER, allocatable :: in_mpi(:,:), out_mpi(:,:) + REAL, ALLOCATABLE :: lon_obs(:) + REAL, ALLOCATABLE :: lat_obs(:) + INTEGER, ALLOCATABLE :: layer_obs(:) + REAL, ALLOCATABLE :: dr_obs(:) + REAL, ALLOCATABLE :: obserr(:) + + integer :: begp, endp ! per-proc beginning and ending pft indices + integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begl, endl ! per-proc beginning and ending landunit indices + integer :: begg, endg ! per-proc gridcell ending gridcell indices + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + real(r8) :: pi + + integer :: ierror, cnt_p + + real :: deltax, deltay + + INTEGER, ALLOCATABLE :: ipiv(:) + real(r8), ALLOCATABLE :: work(:) + + real(r8) :: work_query + integer :: lwork + + integer :: countR, countC, info + + + ! ********************************************* + ! *** Initialize full observation dimension *** + ! ********************************************* + + !IF (mype_filter==0) & + IF (mype_filter==0) & + WRITE (*,*) 'Assimilate observations - obs type GRACE' + + ! Store whether to assimilate this observation type (used in routines below) + + IF (assim_GRACE) thisobs%doassim = 1 + ! Specify type of distance computation + thisobs%disttype = 0 ! 0=Cartesian + + ! Number of coordinates used for distance computation + ! The distance compution starts from the first row + thisobs%ncoord = 2 + + + ! ********************************** + ! *** Read PE-local observations *** + ! ********************************** + + ! read observations from nc file --> call function in mod_read_obs. + ! Idea: when you have multiple observation types in one file, also pass the + ! observation type, here 'GRACE' to the function. You have to give each + ! observation in the file an information which type it is. This way, the + ! output in this function is only the GRACE observation (or soil moisture, + ! ...; dependent on what you want to implement) + + + obs_type_name = 'GRACE' + + ! now call function to get observations + + if (mype_filter==0 .and. screen > 2) then + write(*,*)'load observations from type GRACE' + end if + write(current_observation_filename, '(a, i5.5)') trim(obs_filename)//'.', step + call read_obs_nc_type(current_observation_filename, obs_type_name, & + dim_obs, obs_g, lon_obs, lat_obs, layer_obs, & + dr_obs, obserr, clm_obscov) + if (mype_filter==0 .and. screen > 2) then + write(*,*)'Done: load observations from type GRACE' + end if + + if (dim_obs == 0) then + if (mype_filter==0 .and. screen > 2) then + write(*,*)'TSMP-PDAF mype(w) =', mype_world, & + ': No observations of type GRACE found in file ', & + trim(current_observation_filename) + end if + dim_obs_p = 0 + ALLOCATE(obs_p(1)) + ALLOCATE(ivar_obs_p(1)) + ALLOCATE(ocoord_p(2, 1)) + ALLOCATE(thisobs%id_obs_p(1, 1)) + thisobs%infile=0 + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_GRACE, dim_obs) + return + end if + thisobs%infile=1 + + ! *********************************************************** + ! *** Count available observations for the process domain *** + ! *** and initialize index and coordinate arrays. *** + ! *********************************************************** + + call domain_def_clm(lon_obs, lat_obs, dim_obs, longxy, latixy, longxy_obs, latixy_obs) + + IF (ALLOCATED(vec_useObs)) DEALLOCATE(vec_useObs) + ALLOCATE(vec_useObs(dim_obs)) + IF (ALLOCATED(vec_numPoints)) DEALLOCATE(vec_numPoints) + ALLOCATE(vec_numPoints(dim_obs)) + vec_numPoints=0 + IF (ALLOCATED(vec_numPoints_global)) DEALLOCATE(vec_numPoints_global) + ALLOCATE(vec_numPoints_global(dim_obs)) + IF (ALLOCATED(vec_useObs_global)) DEALLOCATE(vec_useObs_global) + ALLOCATE(vec_useObs_global(dim_obs)) + vec_useObs_global = .true. + IF (ALLOCATED(in_mpi)) DEALLOCATE(in_mpi) + ALLOCATE(in_mpi(2,dim_obs)) + IF (ALLOCATED(out_mpi)) DEALLOCATE(out_mpi) + ALLOCATE(out_mpi(2,dim_obs)) + + + ! additions for GRACE assimilation, it can be the case that not enough + ! CLM gridpoints lie in the neighborhood of a GRACE observation + ! if this is the case, the GRACE observations cannot be reproduced in a + ! satisfactory manner and is not used in the assimilation + ! count grdicells that are in a certain radius. This effect is especially + ! present when the applied GRACE resolution is high or for + ! observation lying directly at the coast + + lon => grc%londeg + lat => grc%latdeg + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + + + do i = 1, dim_obs + + count_points = 0 + + do c = 1, num_layer(1) + + j = hactiveg_levels(c,1) + deltax = abs(longxy(c)-longxy_obs(i)) + deltay = abs(latixy(c)-latixy_obs(i)) + + if((sqrt(real(deltax)**2 + real(deltay)**2)) < dr_obs(1)/0.11) then + count_points = count_points+1 + end if + + end do + + vec_numPoints(i) = count_points + + end do + + ! get vec_numPoints from all processes and add them up together via mpi_allreduce + ! then we have for each observation the number of all points that lie within their neighborhood, not only from one process + + call mpi_allreduce(vec_numPoints,vec_numPoints_global, dim_obs, mpi_integer, mpi_sum, comm_filter, ierror) + + ! only observations should be used that "see" enough gridcells + pi = 3.14159265358979323846 + numPoints = int(ceiling((pi*(dr_obs(1)/0.11)**2)/2)) + if (screen > 2) then + if (mype_filter==0) then + print *, "Minimum number of points for using one observation is ", numPoints + end if + end if + + vec_useObs_global = merge(vec_useObs_global,.false.,vec_numPoints_global>=numPoints) + vec_useObs = vec_useObs_global + + if (screen > 2) then + if (mype_filter==0) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_f_pdaf: vec_useObs_global=", vec_useObs_global + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_f_pdaf: vec_numPoints_global=", vec_numPoints_global + end if + end if + + ! The observation should be reproduced by the process that has the most grid points inside the observation radius + ! However, in the observation operator, all points will be used by using mpi + + in_mpi(1,:) = vec_numPoints + in_mpi(2,:) = mype_filter + + call mpi_allreduce(in_mpi,out_mpi, dim_obs, mpi_2integer, mpi_maxloc, comm_filter, ierror) + + vec_useObs = merge(vec_useObs,.false.,out_mpi(2,:)==mype_filter) + + IF (ALLOCATED(in_mpi)) DEALLOCATE(in_mpi) + IF (ALLOCATED(out_mpi)) DEALLOCATE(out_mpi) + + dim_obs_p = count(vec_useObs) + + if(allocated(thisobs%id_obs_p)) deallocate(thisobs%id_obs_p) + ALLOCATE(thisobs%id_obs_p(1, begg:endg)) + thisobs%id_obs_p(1, :) = 0 + + do i = 1, dim_obs + + if (vec_useObs_global(i)) then + + do c = 1, num_layer(1) + + j = hactiveg_levels(c,1) + deltax = abs(longxy(c)-longxy_obs(i)) + deltay = abs(latixy(c)-latixy_obs(i)) + + if((sqrt(real(deltax)**2 + real(deltay)**2))<=dr_obs(1)/0.11) then + thisobs%id_obs_p(1, j) = i + end if + + end do + + end if + + end do + + + ! now obtain information about which observation is on which PE --> necessary for full VCV matrix later on + ! Allocate array of PE-local observation dimensions + IF (ALLOCATED(local_dims_obs)) DEALLOCATE(local_dims_obs) + ALLOCATE(local_dims_obs(npes_filter)) + + ! Gather array of PE-local observation dimensions + call mpi_allgather(dim_obs_p, 1, MPI_INTEGER, local_dims_obs, 1, MPI_INTEGER, & + comm_filter, ierror) + + ! Allocate observation displacement array local_disp_obs + IF (ALLOCATED(local_disp_obs)) DEALLOCATE(local_disp_obs) + ALLOCATE(local_disp_obs(npes_filter)) + + ! Set observation displacement array local_disp_obs + local_disp_obs(1) = 0 + do i = 2, npes_filter + local_disp_obs(i) = local_disp_obs(i-1) + local_dims_obs(i-1) + end do + + if (mype_filter==0 .and. screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_pdaf: local_disp_obs=", local_disp_obs + end if + + if (allocated(obs_nc2pdaf)) deallocate(obs_nc2pdaf) + allocate(obs_nc2pdaf(count(vec_useObs_global))) + obs_nc2pdaf = 0 + + if (allocated(obs_pdaf2nc)) deallocate(obs_pdaf2nc) + allocate(obs_pdaf2nc(count(vec_useObs_global))) + obs_pdaf2nc = 0 + + + cnt = 0 + j = 0 + do i = 1, dim_obs + if (vec_useObs_global(i)) then + j = j + 1 + end if + if (vec_useObs(i)) then + cnt = cnt + 1 + obs_nc2pdaf(j) = local_disp_obs(mype_filter+1) + cnt + obs_pdaf2nc(local_disp_obs(mype_filter+1) + cnt) = j + end if + end do + + ! collect values from all PEs, by adding all PE-local arrays (works + ! since only the subsection belonging to a specific PE is non-zero) + call mpi_allreduce(MPI_IN_PLACE,obs_pdaf2nc,count(vec_useObs_global),MPI_INTEGER,MPI_SUM,comm_filter,ierror) + call mpi_allreduce(MPI_IN_PLACE,obs_nc2pdaf,count(vec_useObs_global),MPI_INTEGER,MPI_SUM,comm_filter,ierror) + + + + IF (ALLOCATED(obs_p)) DEALLOCATE(obs_p) + ALLOCATE(obs_p(dim_obs_p)) + + IF (ALLOCATED(ivar_obs_p)) DEALLOCATE(ivar_obs_p) + ALLOCATE(ivar_obs_p(dim_obs_p)) + + IF (ALLOCATED(ocoord_p)) DEALLOCATE(ocoord_p) + ALLOCATE(ocoord_p(2, dim_obs_p)) + + if (multierr==0) then + cnt_p = 1 + do i = 1, dim_obs + if (vec_useObs(i)) then + ivar_obs_p(cnt_p) = 1.0/(rms_obs_GRACE*rms_obs_GRACE) + cnt_p = cnt_p + 1 + end if + end do + end if + + if (multierr==1) ivar_obs_p = pack(1/obserr, vec_useObs) + + if (multierr==2) then + + if (allocated(obscov)) deallocate(obscov) + allocate(obscov(count(vec_useObs_global), count(vec_useObs_global))) + countR = 1 + countC = 1 + do i = 1, dim_obs + if (vec_useObs_global(i)) then + do j = 1, dim_obs + if (vec_useObs_global(j)) then + obscov(countR, countC) = clm_obscov(i,j) + countC = countC + 1 + end if + end do + countC = 1 + countR = countR + 1 + end if + end do + + cnt_p = 1 + countC = 1 + do i = 1, dim_obs + if (vec_useObs(i)) then + ivar_obs_p(cnt_p) = 1.0/obscov(countC,countC) + cnt_p = cnt_p + 1 + end if + if (vec_useObs_global(i)) then + countC = countC + 1 + end if + end do + + end if + + cnt_p = 1 + do i = 1, dim_obs + if (vec_useObs(i)) then + obs_p(cnt_p) = obs_g(i) + ocoord_p(1,cnt_p) = longxy_obs(i) + ocoord_p(2,cnt_p) = latixy_obs(i) + cnt_p = cnt_p+1 + end if + end do + + dim_obs = count(vec_useObs_global) + + if (multierr==2) then ! compute inverse of covariance matrix for prodRinvA, has to be before PDAFomi_gather_obs because the routine changes dim_obs + + if (allocated(obscov_inv)) deallocate(obscov_inv) + allocate(obscov_inv(dim_obs, dim_obs)) + + obscov_inv = obscov + + if (allocated(ipiv)) deallocate(ipiv) + ALLOCATE(ipiv(dim_obs)) + + call dgetrf(dim_obs, dim_obs, obscov_inv, dim_obs, ipiv, info) + if (info /= 0) then + print *, "Error in dgetrf, info =", info + stop + end if + + lwork = -1 + call dgetri(dim_obs, obscov_inv, dim_obs, ipiv, work_query, lwork, info) + lwork = int(work_query) + if (allocated(work)) deallocate(work) + allocate(work(lwork)) + call dgetri(dim_obs, obscov_inv, dim_obs, ipiv, work, lwork, info) + if (info /= 0) then + print *, "Error in dgetri, info =", info + stop + end if + + end if + + + + ! **************************************** + ! *** Gather global observation arrays *** + ! **************************************** + + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_GRACE, dim_obs) + + ! ******************** + ! *** Finishing up *** + ! ******************** + + ! load temporal mean of TWS and vectorize for process from begg to endg + ! --> only has to be done once as the process does not change bounds + + if (.not. allocated(tws_temp_mean_d)) then + call read_temp_mean_model(temp_mean_filename) + + if (allocated(tws_temp_mean_d)) DEALLOCATE(tws_temp_mean_d) + ALLOCATE(tws_temp_mean_d(begg:endg)) + tws_temp_mean_d(:) = spval + + do j = begg,endg + outer: do l = 1,size(lon_temp_mean,1) + do k=1,size(lon_temp_mean,2) + if (lon_temp_mean(l,k)==lon(j) .and. lat_temp_mean(l,k)==lat(j)) then + tws_temp_mean_d(j) = tws_temp_mean(l,k) + exit outer + end if + end do + end do outer + + if (lon(j)/=lon_temp_mean(l,k) .or. lat(j)/=lat_temp_mean(l,k)) then + print *, "Attention: distributing model mean to clumps does not work properly" + print *, "idx_lon= ",l, "idx_lat= ",k + print *, "lon(j)= ", lon(j),"lon_temp_mean(idx_lon)= ",lon_temp_mean(l,k) + print *, "lat(j)= ", lat(j),"lat_temp_mean(idx_lat)= ",lat_temp_mean(l,k) + stop + end if + end do + + deallocate(tws_temp_mean) + deallocate(lon_temp_mean) + deallocate(lat_temp_mean) + + end if + + ! Deallocate all local arrays + DEALLOCATE(obs_g) + DEALLOCATE(obs_p, ocoord_p, ivar_obs_p) + + END SUBROUTINE init_dim_obs_GRACE + + + + !------------------------------------------------------------------------------- + !> Implementation of observation operator + !! + !! This routine applies the full observation operator + !! for the type of observations handled in this module. + !! + !! One can choose a proper observation operator from + !! PDAFOMI_OBS_OP or add one to that module or + !! implement another observation operator here. + !! + !! The routine is called by all filter processes. + !! + + !> @author Yorck Ewerdwalbesloh, Anne Springer + !> @date 29.10.2025 + !> @brief Observation operator for GRACE observations + !> @param[in] dim_p PE-local state dimension + !> @param[in] dim_obs Dimension of full observation vector + !> @param[in] state_p PE-local model state + !> @param[in,out] ostate Full observed state + SUBROUTINE obs_op_GRACE(dim_p, dim_obs, state_p, ostate) + + use mpi, only: MPI_DOUBLE_PRECISION + use mpi, only: MPI_SUM + + use enkf_clm_mod, & + only: clm_varsize_tws, state_setup, num_layer, hactiveg_levels + + use mod_parallel_pdaf, & + only: comm_filter + + use clm_varpar , only : nlevsoi + + use decompMod , only : get_proc_bounds + + use clm_varcon, only: spval + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use mod_assimilation, only: screen + + use PDAFomi_obs_f, only: PDAFomi_gather_obsstate + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< Dimension of full observed state (all observed fields) + REAL, INTENT(in) :: state_p(dim_p) !< PE-local model state + REAL, INTENT(inout) :: ostate(dim_obs) !< Full observed state + + REAL, allocatable :: tws_from_statevector(:) + real, allocatable :: ostate_p(:) + + integer :: begp, endp ! per-proc beginning and ending pft indices + integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begl, endl ! per-proc beginning and ending landunit indices + integer :: begg, endg ! per-proc gridcell ending gridcell indices + integer :: ierror + integer :: obs_point ! which observation is seen by which point? + + REAL:: m_state_sum(size(vec_useObs_global)) + REAL:: m_state_sum_global(size(vec_useObs_global)) + + integer :: count, j, g + + ! ****************************************************** + ! *** Apply observation operator H on a state vector *** + ! ****************************************************** + + IF (thisobs%dim_obs_p>0) THEN + if (allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(thisobs%dim_obs_p)) + ELSE + if (allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(1)) + + END IF + + if (thisobs%infile == 1) then ! as I need also tasks with dim_obs_p==0 to reproduce GRACE observations + ! I do not check whether dim_obs_p>0 but if I want to assimilate GRACE. A check still has to be included, + ! else there will be segmentation faults + + m_state_sum(:) = 0 + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + + if (allocated(tws_from_statevector)) deallocate(tws_from_statevector) + allocate(tws_from_statevector(begg:endg)) + + tws_from_statevector(begg:endg) = spval + + select case(state_setup) + case(0) ! all compartments included in state vector + + do j = 1,nlevsoi + do count = 1, num_layer(j) + g = hactiveg_levels(count,j) + if (j==1) then + tws_from_statevector(g) = 0._r8 + end if + if (j==1) then + ! liq + ice + tws_from_statevector(g) = tws_from_statevector(g) + state_p(count) + ! snow + tws_from_statevector(g) = tws_from_statevector(g) + & + state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2)) + !surface water + tws_from_statevector(g) = tws_from_statevector(g) + & + state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2)+ & + clm_varsize_tws(3)) + else + ! liq + ice + tws_from_statevector(g) = tws_from_statevector(g) + state_p(count+sum(num_layer(1:j-1))) + end if + end do + end do + + ! do count = 1, num_hactiveg_patch + ! g = hactiveg_patch(count) + ! ! canopy water + ! tws_from_statevector(g) = tws_from_statevector(g) + & + ! state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2)+ & + ! clm_varsize_tws(3)+ clm_varsize_tws(4)) + ! end do + + + case(1) ! TWS in statevector + + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + tws_from_statevector(g) = state_p(count) + end do + + case(2) ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + + do count = 1,num_layer(1) + g = hactiveg_levels(count,1) + tws_from_statevector(g) = state_p(count) + ! snow added for first layer + tws_from_statevector(g) = tws_from_statevector(g) + & + state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2) + & + clm_varsize_tws(3)) + end do + + do count = 1,num_layer(4) + g = hactiveg_levels(count,4) + tws_from_statevector(g) = tws_from_statevector(g) + state_p(count + clm_varsize_tws(1)) + end do + + do count = 1,num_layer(13) + g = hactiveg_levels(count,13) + tws_from_statevector(g) = tws_from_statevector(g) + state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2)) + end do + + end select + + + + ! subtract mean TWS to obtain TWSA model values + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + obs_point = thisobs%id_obs_p(1, g) + if (obs_point /= 0) then + if (tws_temp_mean_d(g)/=spval .and. tws_from_statevector(g)/=spval) then + m_state_sum(obs_point) = m_state_sum(obs_point) + tws_from_statevector(g)-tws_temp_mean_d(g) + else if (tws_temp_mean_d(g)==spval .and. .not. tws_from_statevector(g)==spval) then + print*, "error, tws temporal mean is spval and reproduced values is not spval for g = ", g + print*, "reproduced = ", tws_from_statevector(g) + stop + else if (.not. tws_temp_mean_d(g)==spval .and. tws_from_statevector(g)==spval) then + print*, "error, tws temporal mean is not spval and reproduced values is spval for g = ", g + print*, "temp_mean = ", tws_temp_mean_d(g) + stop + end if + end if + end do + + ! now get the sum of m_state_sum (sum over all TWSA values for the gridcells for one process) over all processes + call mpi_allreduce(m_state_sum, m_state_sum_global, & + size(vec_useObs_global), mpi_double_precision, & + mpi_sum, comm_filter, ierror) + m_state_sum_global = m_state_sum_global/vec_numPoints_global + + ostate_p = pack(m_state_sum_global, vec_useObs) + + end if + + ! *** Global: Gather full observed state vector + CALL PDAFomi_gather_obsstate(thisobs, ostate_p, ostate) + deallocate(ostate_p) + if (screen>2 .and. mype_filter==0 .and. thisobs%dim_obs_f>0) then + write(*,*)'m_state_sum_global = ', m_state_sum_global + end if + + END SUBROUTINE obs_op_GRACE + + + + !------------------------------------------------------------------------------- + !> Initialize local information on the module-type observation + !! + !! The routine is called during the loop over all local + !! analysis domains. It has to initialize the information + !! about local observations of the module type. It returns + !! number of local observations of the module type for the + !! current local analysis domain in DIM_OBS_L and the full + !! and local offsets of the observation in the overall + !! observation vector. + !! + !! This routine calls the routine PDAFomi_init_dim_obs_l + !! for each observation type. The call allows to specify a + !! different localization radius and localization functions + !! for each observation type and local analysis domain. + !! + + !> @author Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Initialize local observation dimension for GRACE observations + !> @param[in] domain_p Index of current local analysis domain + !> @param[in] step Current time step + !> @param[in] dim_obs Full dimension of observation vector + !> @param[in,out] dim_obs_l Local dimension of observation vector + SUBROUTINE init_dim_obs_l_GRACE(domain_p, step, dim_obs, dim_obs_l) + + ! Include PDAFomi function + USE PDAFomi, ONLY: PDAFomi_init_dim_obs_l + + ! Include localization radius and local coordinates + USE mod_assimilation, & + ONLY: cradius_GRACE, locweight, sradius_GRACE + + use clm_varcon, only: spval + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: domain_p !< Index of current local analysis domain + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(in) :: dim_obs !< Full dimension of observation vector + INTEGER, INTENT(inout) :: dim_obs_l !< Local dimension of observation vector + + REAL :: coords_l(2) ! Coordinates of local analysis domain + + + ! ********************************************** + ! *** Initialize local observation dimension *** + ! ********************************************** + ! count observations within a radius + + if (thisobs%infile==1) then + + ! get coords_l --> coordinates of local analysis domain, I can just set this to the rotated CLM coordinates + coords_l(1) = longxy(domain_p) + coords_l(2) = latixy(domain_p) + + else + coords_l(1) = spval + coords_l(2) = spval + end if + + CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, & + locweight, cradius_GRACE, sradius_GRACE, dim_obs_l) + + END SUBROUTINE init_dim_obs_l_GRACE + + + + !------------------------------------------------------------------------------- + !> Perform covariance localization for local EnKF on the module-type observation + !! + !! The routine is called in the analysis step of the localized + !! EnKF. It has to apply localization to the two matrices + !! HP and HPH of the analysis step for the module-type + !! observation. + !! + !! This routine calls the routine PDAFomi_localize_covar + !! for each observation type. The call allows to specify a + !! different localization radius and localization functions + !! for each observation type. + !! + + !> @author Yorck Ewerdwalbesloh + !> @date 29.10.2025 + !> @brief Covariance localization for GRACE observations, called only for the EnKF + !> @param[in] dim_p PE-local state dimension + !> @param[in] dim_obs Dimension of full observation vector + !> @param[in,out] HP_p PE-local part of matrix HP + !> @param[in,out] HPH Matrix HPH + !> @param[in,out] coords_p Coordinates of state vector elements + SUBROUTINE localize_covar_GRACE(dim_p, dim_obs, HP_p, HPH, coords_p) + + ! Include PDAFomi function + USE PDAFomi, ONLY: PDAFomi_localize_covar + + ! Include localization radius and local coordinates + USE mod_assimilation, & + ONLY: cradius_GRACE, locweight, sradius_GRACE + + use enkf_clm_mod, only: gridcell_state + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< Dimension of observation vector + REAL, INTENT(inout) :: HP_p(dim_obs, dim_p) !< PE local part of matrix HP + REAL, INTENT(inout) :: HPH(dim_obs, dim_obs) !< Matrix HPH + REAL, INTENT(inout) :: coords_p(:,:) !< Coordinates of state vector elements + + integer :: i + + + ! ************************************* + ! *** Apply covariance localization *** + ! ************************************* + + do i = 1,dim_p + coords_p(1,i) = longxy(gridcell_state(i)) + coords_p(2,i) = latixy(gridcell_state(i)) + end do + + + + CALL PDAFomi_localize_covar(thisobs, dim_p, locweight, cradius_GRACE, sradius_GRACE, & + coords_p, HP_p, HPH) + + END SUBROUTINE localize_covar_GRACE + + + subroutine add_obs_err_GRACE(step, dim_obs, C) + + use mod_assimilation, only: obscov, obs_pdaf2nc + use mod_read_obs, only: multierr + USE mod_parallel_pdaf, & + ONLY: npes_filter + + use PDAFomi, only: obsdims + + implicit none + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector + REAL, INTENT(inout) :: C(dim_obs,dim_obs) ! Matrix to that + ! observation covariance R is added + integer :: i, pe, cnt, j + INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector + INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + select case (multierr) + case(0,1) + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + C(i,i) = C(i,i) + 1.0/thisobs%ivar_obs_f(cnt) + cnt = cnt + 1 + end do + end do + case(2) + + do i=1, thisobs%dim_obs_f + do j=1, thisobs%dim_obs_f + C(i,j) = C(i,j) + obscov(obs_pdaf2nc(i),obs_pdaf2nc(j)) + end do + end do + end select + + + DEALLOCATE(id_start, id_end) + + end subroutine add_obs_err_GRACE + + + subroutine init_obscovar_GRACE(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + + use mod_read_obs, only: multierr + + USE mod_parallel_pdaf, & + ONLY: npes_filter + + use PDAFomi, only: obsdims, map_obs_id + + use mod_assimilation, only: obs_pdaf2nc, obscov + + implicit none + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector + REAL, INTENT(inout) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix + REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PE-local observation vector + LOGICAL, INTENT(inout) :: isdiag ! Whether the observation error covar. matrix is diagonal + + integer :: i, pe, cnt, j + INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector + INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + ! Initialize indices + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + ! Initialize mapping vector (to be used in PDAF_enkf_obs_ensemble) + cnt = 1 + IF (thisobs%obsid-1 > 0) cnt = cnt+ SUM(obsdims(:,1:thisobs%obsid-1)) + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + map_obs_id(i) = cnt + cnt = cnt + 1 + END DO + END DO + + select case(multierr) + case(0,1) + + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + covar(i, i) = covar(i, i) + 1.0/thisobs%ivar_obs_f(cnt) + cnt = cnt + 1 + ENDDO + ENDDO + + ! The matrix is diagonal + ! This setting avoids the computation of the SVD of COVAR + ! in PDAF_enkf_obs_ensemble + isdiag = .TRUE. + case(2) + + do i=1, thisobs%dim_obs_f + do j=1, thisobs%dim_obs_f + covar(i, i) = obscov(obs_pdaf2nc(i),obs_pdaf2nc(j)) + end do + end do + + isdiag = .FALSE. + + end select + + + DEALLOCATE(id_start, id_end) + + + end subroutine init_obscovar_GRACE + + + subroutine prodRinvA_GRACE(step, dim_obs_p, rank, obs_p, A_p, C_p) + + use mod_read_obs, only: multierr + use mod_assimilation, only: obscov_inv, obs_pdaf2nc + use shr_kind_mod, only: r8 => shr_kind_r8 + + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_p(dim_obs_p) ! PE-local vector of observations + REAL, INTENT(in) :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine + REAL, INTENT(inout) :: C_p(dim_obs_p,rank) ! Output matrix + + INTEGER :: i, j ! index of observation component + INTEGER :: off ! row offset in A_l and C_l + + real(r8) :: obscov_inv_l(thisobs%dim_obs_f,thisobs%dim_obs_f) ! errors of observations in the model domain + + off = thisobs%off_obs_f ! account for offset if multiple observation types are assimilated at once + + select case (multierr) + case(0,1) + do j=1, rank + do i=1, thisobs%dim_obs_f + C_p(i+off, j) = thisobs%ivar_obs_f(i) * A_p(i+off, j) + END DO + end do + + case(2) + do i =1, thisobs%dim_obs_f + do j = 1, thisobs%dim_obs_f + obscov_inv_l(i,j) = obscov_inv(obs_pdaf2nc(i),obs_pdaf2nc(j)) + end do + end do + C_p(off+1:off+thisobs%dim_obs_f,:) = matmul(obscov_inv_l,A_p(off+1:off+thisobs%dim_obs_f,:)) + end select + + + end subroutine prodRinvA_GRACE + + + subroutine prodRinvA_l_GRACE(domain_p, step, dim_obs, rank, obs_l, A_l, C_l) + + use shr_kind_mod, only: r8 => shr_kind_r8 + use mod_assimilation, only: obscov, obs_pdaf2nc, cradius_GRACE, locweight + use mod_read_obs, only: multierr + use PDAFomi, only: PDAFomi_observation_localization_weights + + implicit none + + INTEGER, INTENT(in) :: domain_p ! Current local analysis domain + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of local observation vector, multiple observation types possible, then we have to access with thisobs_l%dim_obs_l + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_l(dim_obs) ! Local vector of observations + REAL, INTENT(inout) :: A_l(dim_obs, rank) ! Input matrix from analysis routine + REAL, INTENT(out) :: C_l(dim_obs, rank) ! Output matrix + + INTEGER :: verbose ! Verbosity flag + INTEGER :: verbose_w ! Verbosity flag for weight computation + INTEGER, SAVE :: domain_save = -1 ! Save previous domain index + INTEGER :: wtype ! Type of weight function + INTEGER :: rtype ! Type of weight regulation + REAL, ALLOCATABLE :: weight(:) ! Localization weights + REAL, ALLOCATABLE :: A_obs(:,:) ! Array for a single row of A_l + REAL :: var_obs ! Variance of observation error + + INTEGER :: i, j + + INTEGER :: off ! row offset in A_l and C_l + INTEGER :: idummy ! Dummy to access nobs_all + + real(r8) :: ivariance_obs + + REAL(r8) :: obscov_l(thisobs_l%dim_obs_l,thisobs_l%dim_obs_l) ! local observation covariance matrix + REAL(r8) :: obscov_inv_l(thisobs_l%dim_obs_l,thisobs_l%dim_obs_l) ! inverse of local observation covariance matrix + + INTEGER, ALLOCATABLE :: ipiv(:) + real(r8), ALLOCATABLE :: work(:) + INTEGER :: info, lwork + real(r8) :: work_query + + real(r8) :: maxdiff + + + off = thisobs_l%off_obs_l + idummy = dim_obs + + IF ((domain_p <= domain_save .OR. domain_save < 0) .AND. mype_filter==0) THEN + verbose = 1 + ELSE + verbose = 0 + END IF + domain_save = domain_p + + ! Screen output + IF (verbose == 1) THEN + WRITE (*, '(8x, a, f12.3)') & + '--- Use global rms for observations of ', rms_obs_GRACE + WRITE (*, '(8x, a, 1x)') & + '--- Domain localization' + WRITE (*, '(12x, a, 1x, f12.2)') & + '--- Local influence radius', cradius_GRACE + + IF (locweight > 0) THEN + WRITE (*, '(12x, a)') & + '--- Use distance-dependent weight for observation errors' + + IF (locweight == 3) THEN + write (*, '(12x, a)') & + '--- Use regulated weight with mean error variance' + ELSE IF (locweight == 4) THEN + write (*, '(12x, a)') & + '--- Use regulated weight with single-point error variance' + END IF + END IF + ENDIF + + ALLOCATE(weight(thisobs_l%dim_obs_l)) + call PDAFomi_observation_localization_weights(thisobs_l, thisobs, rank, A_l, & + weight, verbose) + + select case(multierr) + case(0,1) + do j=1,rank + do i=1,thisobs_l%dim_obs_l + C_l(i+off,j) = thisobs_l%ivar_obs_l(i) * weight(i) * A_l(i+off, j) + end do + end do + + case(2) + + obscov_l = 0.0_r8 + + do i=1, thisobs_l%dim_obs_l ! fill local observation covariance matrix, invert it and apply it to A_l + do j=1, thisobs_l%dim_obs_l + obscov_l(i,j) = obscov(obs_pdaf2nc(thisobs_l%id_obs_l(i)),obs_pdaf2nc(thisobs_l%id_obs_l(j))) + end do + end do + + obscov_inv_l = obscov_l + + if (allocated(ipiv)) deallocate(ipiv) + ALLOCATE(ipiv(thisobs_l%dim_obs_l)) + + call dgetrf(thisobs_l%dim_obs_l, thisobs_l%dim_obs_l, obscov_inv_l, thisobs_l%dim_obs_l, ipiv, info) + if (info /= 0) then + print *, "Error in dgetrf, info =", info + stop + end if + + lwork = -1 + call dgetri(thisobs_l%dim_obs_l, obscov_inv_l, thisobs_l%dim_obs_l, ipiv, work_query, lwork, info) + lwork = int(work_query) + if (allocated(work)) deallocate(work) + allocate(work(lwork)) + call dgetri(thisobs_l%dim_obs_l, obscov_inv_l, thisobs_l%dim_obs_l, ipiv, work, lwork, info) + if (info /= 0) then + print *, "Error in dgetri, info =", info + stop + end if + + do j = 1,rank + do i = 1,thisobs_l%dim_obs_l + A_l(i+off,j) = weight(i)*A_l(i+off,j) + end do + end do + + C_l(off+1:off+thisobs_l%dim_obs_l,:) = matmul(obscov_inv_l,A_l(off+1:off+thisobs_l%dim_obs_l,:)) + + end select + + deallocate(weight) + + end subroutine prodRinvA_l_GRACE + + + !> @author Anne Springer + !> @date 29.10.2025 + !> @brief Read temporal mean of TWS from NetCDF file + !> @param[in] temp_mean_filename filename of NetCDF file containing temporal mean of TWS + subroutine read_temp_mean_model(temp_mean_filename) + + use netcdf, only: nf90_max_name + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_get_var + use netcdf, only: nf90_close + use mod_read_obs, only: check + + implicit none + integer :: ncid, dim_lon, dim_lat, lon_varid, lat_varid, tws_varid + character (len = *), parameter :: dim_lon_name = "lsmlon" + character (len = *), parameter :: dim_lat_name = "lsmlat" + character (len = *), parameter :: lon_name = "longitude" + character (len = *), parameter :: lat_name = "latitude" + character (len = *), parameter :: tws_name = "TWS" + character(len = nf90_max_name) :: RecordDimName + integer :: dimid_lon, dimid_lat, status + integer :: haserr + character (len = *), intent(in) :: temp_mean_filename + + call check(nf90_open(temp_mean_filename, nf90_nowrite, ncid)) + call check(nf90_inq_dimid(ncid, dim_lon_name, dimid_lon)) + call check(nf90_inq_dimid(ncid, dim_lat_name, dimid_lat)) + call check(nf90_inquire_dimension(ncid, dimid_lon, recorddimname, dim_lon)) + call check(nf90_inquire_dimension(ncid, dimid_lat, recorddimname, dim_lat)) + + if(allocated(lon_temp_mean))deallocate(lon_temp_mean) + if(allocated(lat_temp_mean))deallocate(lat_temp_mean) + if(allocated(tws_temp_mean))deallocate(tws_temp_mean) + + allocate(tws_temp_mean(dim_lon,dim_lat)) + allocate(lon_temp_mean(dim_lon,dim_lat)) + allocate(lat_temp_mean(dim_lon,dim_lat)) + + call check( nf90_inq_varid(ncid, lon_name, lon_varid)) + call check(nf90_get_var(ncid, lon_varid, lon_temp_mean)) + + call check( nf90_inq_varid(ncid, lat_name, lat_varid)) + call check(nf90_get_var(ncid, lat_varid, lat_temp_mean)) + + call check( nf90_inq_varid(ncid, tws_name, tws_varid)) + call check(nf90_get_var(ncid, tws_varid, tws_temp_mean)) + + call check( nf90_close(ncid) ) + + end subroutine read_temp_mean_model + + END MODULE obs_GRACE_pdafomi +#endif + + + diff --git a/interface/framework/obs_SM_pdafomi.F90 b/interface/framework/obs_SM_pdafomi.F90 new file mode 100644 index 000000000..882558277 --- /dev/null +++ b/interface/framework/obs_SM_pdafomi.F90 @@ -0,0 +1,1317 @@ +!> PDAF-OMI observation module for type SM observations +!! +!! This module handles operations for one data type (called 'module-type' below): +!! OBSTYPE = SM +!! +!! __Observation type SM:__ +!! The observation type SM are soil moisture observations +!! +!! The subroutines in this module are for the particular handling of +!! a single observation type. +!! The routines are called by the different call-back routines of PDAF +!! usually by callback_obs_pdafomi.F90 +!! Most of the routines are generic so that in practice only 2 routines +!! need to be adapted for a particular data type. These are the routines +!! for the initialization of the observation information (init_dim_obs) +!! and for the observation operator (obs_op). +!! +!! The module and the routines are named according to the observation type. +!! This allows to distinguish the observation type and the routines in this +!! module from other observation types. +!! +!! The module uses two derived data types (obs_f and obs_l), which contain +!! all information about the full and local observations. Only variables +!! of the type obs_f need to be initialized in this module. The variables +!! in the type obs_l are initilized by the generic routines from PDAFomi. +!! +!! +!! These 2 routines need to be adapted for the particular observation type: +!! * init_dim_obs_OBSTYPE \n +!! Count number of process-local and full observations; +!! initialize vector of observations and their inverse variances; +!! initialize coordinate array and index array for indices of +!! observed elements of the state vector. +!! * obs_op_OBSTYPE \n +!! observation operator to get full observation vector of this type. Here +!! one has to choose a proper observation operator or implement one. +!! +!! In addition, there are two optional routines, which are required if filters +!! with localization are used: +!! * init_dim_obs_l_OBSTYPE \n +!! Only required if domain-localized filters (e.g. LESTKF, LETKF) are used: +!! Count number of local observations of module-type according to +!! their coordinates (distance from local analysis domain). Initialize +!! module-internal distances and index arrays. +!! * localize_covar_OBSTYPE \n +!! Only required if the localized EnKF is used: +!! Apply covariance localization in the LEnKF. +!! +!! __Revision history:__ +!! * 2019-06 - Lars Nerger - Initial code +!! * Later revisions - see repository log +!! + +! Author: Yorck Ewerdwalbesloh, adaptations of original implementations of TSMP2-PDAF interface for OMI framework + + +#ifdef CLMFIVE +MODULE obs_SM_pdafomi + + USE mod_parallel_pdaf, & + ONLY: mype_filter ! Rank of filter process + USE PDAFomi, & + ONLY: obs_f, obs_l ! Declaration of observation data types + + IMPLICIT NONE + SAVE + + PUBLIC + + ! Variables which are inputs to the module (usually set in init_pdaf) + LOGICAL :: assim_SM !< Whether to assimilate this data type + REAL :: rms_obs_SM !< Observation error standard deviation (for constant errors) + + ! longitude and latitude of grid cells and observation cells + INTEGER, ALLOCATABLE :: longxy(:), latixy(:), longxy_obs(:), latixy_obs(:) + + ! One can declare further variables, e.g. for file names which can + ! be use-included in init_pdaf() and initialized there. + + + ! ********************************************************* + ! *** Data type obs_f defines the full observations by *** + ! *** internally shared variables of the module *** + ! ********************************************************* + + ! Relevant variables that can be modified by the user: + ! TYPE obs_f + ! ---- Mandatory variables to be set in INIT_DIM_OBS ---- + ! INTEGER :: doassim ! Whether to assimilate this observation type + ! INTEGER :: disttype ! Type of distance computation to use for localization + ! ! (0) Cartesian, (1) Cartesian periodic + ! ! (2) simplified geographic, (3) geographic haversine function + ! INTEGER :: ncoord ! Number of coordinates use for distance computation + ! INTEGER, ALLOCATABLE :: id_obs_p(:,:) ! Indices of observed field in state vector (process-local) + ! + ! ---- Optional variables - they can be set in INIT_DIM_OBS ---- + ! REAL, ALLOCATABLE :: icoeff_p(:,:) ! Interpolation coefficients for obs. operator + ! REAL, ALLOCATABLE :: domainsize(:) ! Size of domain for periodicity (<=0 for no periodicity) + ! + ! ---- Variables with predefined values - they can be changed in INIT_DIM_OBS ---- + ! INTEGER :: obs_err_type=0 ! Type of observation error: (0) Gauss, (1) Laplace + ! INTEGER :: use_global_obs=1 ! Whether to use (1) global full obs. + ! ! or (0) obs. restricted to those relevant for a process domain + ! REAL :: inno_omit=0.0 ! Omit obs. if squared innovation larger this factor times + ! ! observation variance + ! REAL :: inno_omit_ivar=1.0e-12 ! Value of inverse variance to omit observation + ! END TYPE obs_f + + ! Data type obs_l defines the local observations by internally shared variables of the module + + ! *********************************************************************** + + ! Declare instances of observation data types used here + ! We use generic names here, but one could rename the variables + TYPE(obs_f), TARGET, PUBLIC :: thisobs ! full observation + TYPE(obs_l), TARGET, PUBLIC :: thisobs_l ! local observation + + !$OMP THREADPRIVATE(thisobs_l) + + + !------------------------------------------------------------------------------- + + CONTAINS + + !> Initialize information on the module-type observation + !! + !! The routine is called by each filter process. + !! at the beginning of the analysis step before + !! the loop through all local analysis domains. + !! + !! It has to count the number of observations of the + !! observation type handled in this module according + !! to the current time step for all observations + !! required for the analyses in the loop over all local + !! analysis domains on the PE-local state domain. + !! + !! The following four variables have to be initialized in this routine + !! * thisobs\%doassim - Whether to assimilate this type of observations + !! * thisobs\%disttype - type of distance computation for localization with this observaton + !! * thisobs\%ncoord - number of coordinates used for distance computation + !! * thisobs\%id_obs_p - array with indices of module-type observation in process-local state vector + !! + !! Optional is the use of + !! * thisobs\%icoeff_p - Interpolation coefficients for obs. operator (only if interpolation is used) + !! * thisobs\%domainsize - Size of domain for periodicity for disttype=1 (<0 for no periodicity) + !! * thisobs\%obs_err_type - Type of observation errors for particle filter and NETF (default: 0=Gaussian) + !! * thisobs\%use_global obs - Whether to use global observations or restrict the observations to the relevant ones + !! (default: 1=use global full observations) + !! * thisobs\%inno_omit - Omit obs. if squared innovation larger this factor times observation variance + !! (default: 0.0, omission is active if >0) + !! * thisobs\%inno_omit_ivar - Value of inverse variance to omit observation + !! (default: 1.0e-12, change this if this value is not small compared to actual obs. error) + !! + !! Further variables are set when the routine PDAFomi_gather_obs is called. + !! + SUBROUTINE init_dim_obs_SM(step, dim_obs) + + USE mpi, ONLY: MPI_INTEGER + USE mpi, ONLY: MPI_DOUBLE_PRECISION + USE mpi, ONLY: MPI_IN_PLACE + USE mpi, ONLY: MPI_SUM + + USE mod_parallel_pdaf, & + ONLY: mype_filter, comm_filter, npes_filter, abort_parallel, & + mype_world + + USE mod_assimilation, & + ONLY: obs_index_p, obs_filename, & + obs_interp_indices_p, & + obs_interp_weights_p, & + obs_pdaf2nc, obs_nc2pdaf, & + local_dims_obs, & + local_disp_obs, & + ! dim_obs_p, & + longxy_obs_floor, latixy_obs_floor, & + screen, cradius_SM + + + USE PDAFomi, & + ONLY: PDAFomi_gather_obs, pi + USE mod_assimilation, & + ONLY: obs_filename, screen + + Use mod_read_obs, & + only: multierr, read_obs_nc_type + + use enkf_clm_mod, only: clmstatevec_allcol, clmstatevec_only_active, clmstatevec_max_layer, state_clm2pdaf_p + + use enkf_clm_mod, only: domain_def_clm + + use mod_parallel_pdaf, & + only: abort_parallel + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use GridcellType, only: grc + + use clm_varcon, only: spval + use clm_varcon, only: ispval + + USE mod_parallel_pdaf, & + ONLY: mype_world + + use decompMod , only : get_proc_bounds, get_proc_global + + use ColumnType , only : col + + use mod_tsmp, only: obs_interp_switch + + use enkf_clm_mod, only: get_interp_idx + + use mod_tsmp, only: point_obs + + use mod_tsmp, only: da_print_obs_index + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(inout) :: dim_obs !< Dimension of full observation vector + + ! *** Local variables *** + INTEGER :: i, j, c, g, cg ! Counters + INTEGER :: cnt ! Counters + INTEGER :: cnt_interp + INTEGER :: dim_obs_p ! Number of process-local observations + REAL, ALLOCATABLE :: obs_p(:) ! PE-local observation vector + REAL, ALLOCATABLE :: obs_g(:) ! Global observation vector + REAL, ALLOCATABLE :: ivar_obs_p(:) ! PE-local inverse observation error variance + REAL, ALLOCATABLE :: ocoord_p(:,:) ! PE-local observation coordinates + character (len = 110) :: current_observation_filename + + + character(len=20) :: obs_type_name ! name of observation type (e.g. GRACE, SM, ST, ...) + + REAL, ALLOCATABLE :: lon_obs(:) + REAL, ALLOCATABLE :: lat_obs(:) + INTEGER, ALLOCATABLE :: layer_obs(:) + REAL, ALLOCATABLE :: dr_obs(:) + REAL, ALLOCATABLE :: obserr(:) + REAL, ALLOCATABLE :: obscov(:,:) + + integer :: begp, endp ! per-proc beginning and ending pft indices + integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begl, endl ! per-proc beginning and ending landunit indices + integer :: begg, endg ! per-proc gridcell ending gridcell indices + + integer :: numg ! total number of gridcells across all processors + integer :: numl ! total number of landunits across all processors + integer :: numc ! total number of columns across all processors + integer :: nump ! total number of pfts across all processors + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + integer, pointer :: mycgridcell(:) + + integer :: ierror, cnt_p + + real :: deltax, deltay + + logical :: is_use_dr + logical :: obs_snapped !Switch for checking multiple observation counts + logical :: newgridcell + + INTEGER :: sum_dim_obs_p + + real :: sum_interp_weights + + character (len = 27) :: fn !TSMP-PDAF: function name for obs_index_p output + + + + ! ********************************************* + ! *** Initialize full observation dimension *** + ! ********************************************* + + IF (mype_filter==0) & + WRITE (*,*) 'Assimilate observations - obs type soil moisture' + + ! Store whether to assimilate this observation type (used in routines below) + + IF (assim_SM) thisobs%doassim = 1 + ! Specify type of distance computation + + thisobs%disttype = 3 ! 0=Cartesian, 3=geographic, distance computed with haversine formula + + ! Number of coordinates used for distance computation + ! The distance compution starts from the first row + thisobs%ncoord = 2 + + + ! ********************************** + ! *** Read PE-local observations *** + ! ********************************** + + + obs_type_name = 'SM' + + ! now call function to get observations + + if (mype_filter==0 .and. screen > 2) then + write(*,*)'load observations from type SM' + end if + write(current_observation_filename, '(a, i5.5)') trim(obs_filename)//'.', step + + + if (mype_filter == 0) then + call read_obs_nc_type(current_observation_filename, obs_type_name, & + dim_obs, obs_g, lon_obs, lat_obs, layer_obs, & + dr_obs, obserr, obscov) + end if + + call mpi_bcast(dim_obs, 1, MPI_INTEGER, 0, comm_filter, ierror) + + ! check if file contains observations from type SM + + if (dim_obs == 0) then + if (mype_filter==0 .and. screen > 2) then + write(*,*)'TSMP-PDAF mype(w) =', mype_world, & + ': No observations of type SM found in file ', & + trim(current_observation_filename) + end if + dim_obs_p = 0 + ALLOCATE(obs_p(1)) + ALLOCATE(ivar_obs_p(1)) + ALLOCATE(ocoord_p(2, 1)) + ALLOCATE(thisobs%id_obs_p(1, 1)) + thisobs%infile=0 + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_SM, dim_obs) + return + end if + + call mpi_bcast(multierr, 1, MPI_INTEGER, 0, comm_filter, ierror) + + ! Allocate observation arrays for non-root procs + ! ---------------------------------------------- + if (mype_filter /= 0) then ! for all non-master proc + if(allocated(obs_g)) deallocate(obs_g) + allocate(obs_g(dim_obs)) + if(allocated(lon_obs)) deallocate(lon_obs) + allocate(lon_obs(dim_obs)) + if(allocated(lat_obs)) deallocate(lat_obs) + allocate(lat_obs(dim_obs)) + if(allocated(dr_obs)) deallocate(dr_obs) + allocate(dr_obs(dim_obs)) + if(allocated(layer_obs)) deallocate(layer_obs) + allocate(layer_obs(dim_obs)) + if(multierr==1) then + if(allocated(obserr)) deallocate(obserr) + allocate(obserr(dim_obs)) + end if + end if + + call mpi_bcast(obs_g, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + if(multierr==1) call mpi_bcast(obserr, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(lon_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(lat_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(dr_obs, dim_obs, MPI_DOUBLE_PRECISION, 0, comm_filter, ierror) + call mpi_bcast(layer_obs, dim_obs, MPI_INTEGER, 0, comm_filter, ierror) + + + + if (mype_filter==0 .and. screen > 2) then + write(*,*)'Done: load observations from type SM' + end if + + + thisobs%infile = 1 + call domain_def_clm(lon_obs, lat_obs, dim_obs, longxy, latixy, longxy_obs, latixy_obs) + + ! Interpolation of measured states: Save the indices of the + ! nearest grid points + if (obs_interp_switch == 1) then + ! Get the floored values for latitudes and longitudes + call get_interp_idx(lon_obs, lat_obs, dim_obs, longxy_obs_floor, latixy_obs_floor) + end if + +#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 + + + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + call get_proc_global(numg, numl, numc, nump) + + dim_obs_p = 0 + + is_use_dr = .true. + + if(allocated(thisobs%id_obs_p)) deallocate(thisobs%id_obs_p) + allocate(thisobs%id_obs_p(1,endg-begg+1)) + thisobs%id_obs_p(1, :) = 0 + + cnt = 1 + + do i = 1, dim_obs + obs_snapped = .false. + do g = begg, endg + + newgridcell = .true. + + do c = begc,endc + + cg = mycgridcell(c) + + if(cg == g) then + + if(newgridcell) then + + if(is_use_dr) then + if (lon(g)>180) then + deltax = abs(lon(g)-lon_obs(i)-360) + else + deltax = abs(lon(g)-lon_obs(i)) + end if + deltay = abs(lat(g)-lat_obs(i)) + end if + + ! Assigning observations to grid cells according to + ! snapping distance or index arrays + if(((is_use_dr).and.(deltax<=dr_obs(1)).and.(deltay<=dr_obs(1))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(cnt)) .and. & + (latixy_obs(i) == latixy(cnt)))) then + + dim_obs_p = dim_obs_p + 1 + ! Use index array for setting the correct state vector index in `obs_id_p` + thisobs%id_obs_p(1,state_clm2pdaf_p(c,layer_obs(i))) = i + + if (obs_snapped) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR Observation snapped at multiple grid cells." + print *, "i=", i + if (is_use_dr) then + print *, "lon_obs(i)=", lon_obs(i) + print *, "lat_obs(i)=", lat_obs(i) + end if + call abort_parallel() + end if + + ! Set observation as counted + obs_snapped = .true. + newgridcell = .false. + + cnt=cnt+1 + end if + end if + end if + end do + end do + end do + + if (screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_pdaf: dim_obs_p=", dim_obs_p + end if + + + ! initalize OMI arrays + + IF (ALLOCATED(ivar_obs_p)) DEALLOCATE(ivar_obs_p) + ALLOCATE(ivar_obs_p(dim_obs_p)) + + IF (ALLOCATED(ocoord_p)) DEALLOCATE(ocoord_p) + ALLOCATE(ocoord_p(2, dim_obs_p)) + + + ! Dimension of full observation vector + ! ------------------------------------ + + ! add and broadcast size of PE-local observation dimensions using mpi_allreduce + call mpi_allreduce(dim_obs_p, sum_dim_obs_p, 1, MPI_INTEGER, MPI_SUM, & + comm_filter, ierror) + + ! Check sum of dimensions of PE-local observation vectors against + ! dimension of full observation vector + if (.not. sum_dim_obs_p == dim_obs) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR Sum of PE-local observation dimensions" + print *, "sum_dim_obs_p=", sum_dim_obs_p + print *, "dim_obs=", dim_obs + call abort_parallel() + end if + + ! Gather PE-local observation dimensions and displacements in arrays + ! ---------------------------------------------------------------- + + ! Allocate array of PE-local observation dimensions + IF (ALLOCATED(local_dims_obs)) DEALLOCATE(local_dims_obs) + ALLOCATE(local_dims_obs(npes_filter)) + + ! Gather array of PE-local observation dimensions + call mpi_allgather(dim_obs_p, 1, MPI_INTEGER, local_dims_obs, 1, MPI_INTEGER, & + comm_filter, ierror) + + ! Allocate observation displacement array local_disp_obs + IF (ALLOCATED(local_disp_obs)) DEALLOCATE(local_disp_obs) + ALLOCATE(local_disp_obs(npes_filter)) + + ! Set observation displacement array local_disp_obs + local_disp_obs(1) = 0 + do i = 2, npes_filter + local_disp_obs(i) = local_disp_obs(i-1) + local_dims_obs(i-1) + end do + + if (mype_filter==0 .and. screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_pdaf: local_disp_obs=", local_disp_obs + end if + + + ! Write index mapping array NetCDF->PDAF + ! -------------------------------------- + ! Set index mapping `obs_pdaf2nc` between observation order in + ! NetCDF input and observation order in pdaf as determined by domain + ! decomposition. + + ! Use-case: Correct index order in loops over NetCDF-observation + ! file input arrays. + + ! Trivial example: The order in the NetCDF file corresponds exactly + ! to the order in the domain decomposition in PDAF, e.g. for a + ! single PE per component model run. + + ! Non-trivial example: The first observation in the NetCDF file is + ! not located in the domain/subgrid of the first PE. Rather, the + ! second observation in the NetCDF file (`i=2`) is the only + ! observation (`cnt = 1`) in the subgrid of the first PE + ! (`mype_filter = 0`). This leads to a non-trivial index mapping, + ! e.g. `obs_pdaf2nc(1)==2`: + ! + ! i = 2 + ! cnt = 1 + ! mype_filter = 0 + ! + ! obs_pdaf2nc(local_disp_obs(mype_filter+1)+cnt) = i + !-> obs_pdaf2nc(local_disp_obs(1)+1) = 2 + !-> obs_pdaf2nc(1) = 2 + + + if (allocated(obs_pdaf2nc)) deallocate(obs_pdaf2nc) + allocate(obs_pdaf2nc(dim_obs)) + obs_pdaf2nc = 0 + if (allocated(obs_nc2pdaf)) deallocate(obs_nc2pdaf) + allocate(obs_nc2pdaf(dim_obs)) + obs_nc2pdaf = 0 + + + if (point_obs==1) then + + cnt = 1 + do i = 1, dim_obs + ! Many processes may not contain the observation / do not need + ! to snap it, so default true + obs_snapped = .true. + + do g = begg,endg + newgridcell = .true. + do c = begc,endc + cg = mycgridcell(c) + if(cg == g) then + if(newgridcell) then + + if(is_use_dr) then + if (lon(g)>180) then + deltax = abs(lon(g)-lon_obs(i)-360) + else + deltax = abs(lon(g)-lon_obs(i)) + end if + deltay = abs(lat(g)-lat_obs(i)) + end if + + if(((is_use_dr).and.(deltax<=dr_obs(1)).and.(deltay<=dr_obs(1))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. & + (latixy_obs(i) == latixy(g-begg+1)))) then +#ifdef CLMFIVE + if(state_clm2pdaf_p(c,1)==ispval) then + ! `ispval`: column not in state vector, most likely + ! because it is hydrologically inactive + + ! Observation not snapped, even though location is + ! right! + obs_snapped = .false. + + ! Do not use this column for snapping an + ! observation, instead cycle to next column + cycle + end if +#endif + obs_pdaf2nc(local_disp_obs(mype_filter+1)+cnt) = i + obs_nc2pdaf(i) = local_disp_obs(mype_filter+1)+cnt + cnt = cnt + 1 + + ! Observation snapped at location (possibly + ! overwriting a false from inactive column before) + obs_snapped = .true. + end if + + newgridcell = .false. + + end if + end if + end do + end do + + ! Warning, when an observation has not been snapped to any + ! active gridcell. + if(.not. obs_snapped) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observations exist at non-active gridcells." + print *, "Consider removing the following observation from the observation files." + print *, "Observation-index in NetCDF-file: i=", i + call abort_parallel() + end if + end do + + end if + + ! collect values from all PEs, by adding all PE-local arrays (works + ! since only the subsection belonging to a specific PE is non-zero) + call mpi_allreduce(MPI_IN_PLACE,obs_pdaf2nc,dim_obs,MPI_INTEGER,MPI_SUM,comm_filter,ierror) + call mpi_allreduce(MPI_IN_PLACE,obs_nc2pdaf,dim_obs,MPI_INTEGER,MPI_SUM,comm_filter,ierror) + + if (mype_filter==0 .and. screen > 2) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": init_dim_obs_pdaf: obs_pdaf2nc=", obs_pdaf2nc + end if + + ! Write process-local observation arrays + ! -------------------------------------- + IF (ALLOCATED(obs_p)) DEALLOCATE(obs_p) + ALLOCATE(obs_p(dim_obs_p)) + IF (ALLOCATED(obs_index_p)) DEALLOCATE(obs_index_p) + ALLOCATE(obs_index_p(dim_obs_p)) + if(obs_interp_switch == 1) then + ! Array for storing indices from states that are interpolated to observation locations + IF (ALLOCATED(obs_interp_indices_p)) DEALLOCATE(obs_interp_indices_p) + ALLOCATE(obs_interp_indices_p(dim_obs_p, 4)) ! Later 8 for 3D / ParFlow + IF (ALLOCATED(obs_interp_weights_p)) DEALLOCATE(obs_interp_weights_p) + ALLOCATE(obs_interp_weights_p(dim_obs_p, 4)) ! Later 8 for 3D / ParFlow + end if + + + cnt = 1 + + do i = 1, dim_obs + + do g = begg,endg + newgridcell = .true. + + do c = begc,endc + + cg = mycgridcell(c) + + if(cg == g) then + + if(newgridcell) then + + if(is_use_dr) then + if (lon(g)>180) then + deltax = abs(lon(g)-lon_obs(i)-360) + else + deltax = abs(lon(g)-lon_obs(i)) + end if + deltay = abs(lat(g)-lat_obs(i)) + end if + + if(((is_use_dr).and.(deltax<=dr_obs(1)).and.(deltay<=dr_obs(1))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. & + (latixy_obs(i) == latixy(g-begg+1)))) then + ! if haversine formula in distance calculation, the coordinates have to be converted to radians + if (thisobs%disttype/=3) then + ocoord_p(1,cnt) = lon_obs(i) + ocoord_p(2,cnt) = lat_obs(i) + else + ocoord_p(1,cnt) = lon_obs(i) * pi / 180.0 + ocoord_p(2,cnt) = lat_obs(i) * pi / 180.0 + end if + + ! Different settings of observation-location-index in + ! state vector depending on the method of state + ! vector assembling. + if(clmstatevec_allcol==1) then +#ifdef CLMFIVE + if(clmstatevec_only_active==1) then + + ! Error if observation deeper than clmstatevec_max_layer + if(layer_obs(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then + print *, "TSMP-PDAF mype(w)=", mype_world, & + ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "i=", i + print *, "c=", c + print *, "layer_obs(i)=", layer_obs(i) + print *, "col%nbedrock(c)=", col%nbedrock(c) + print *, "clmstatevec_max_layer=", clmstatevec_max_layer + call abort_parallel() + end if + obs_index_p(cnt) = state_clm2pdaf_p(c,layer_obs(i)) + else + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (layer_obs(i)-1)) + end if +#else + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) +#endif + else +#ifdef CLMFIVE + if(clmstatevec_only_active==1) then + obs_index_p(cnt) = state_clm2pdaf_p(c,layer_obs(i)) + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (layer_obs(i)-1)) + end if +#else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) +#endif + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = obs_g(i) + if(multierr==1) ivar_obs_p(cnt) = 1/(obserr(i)*obserr(i)) + if(multierr==0) ivar_obs_p(cnt) = 1/(rms_obs_SM*rms_obs_SM) + cnt = cnt + 1 + end if + + newgridcell = .false. + + end if + + end if + + end do + end do + + end do + + if(obs_interp_switch==1) then + ! loop over all obs and save the indices of the nearest grid + ! points to array obs_interp_indices_p and save the distance + ! weights to array obs_interp_weights_p (later normalized) + cnt = 1 + do i = 1, dim_obs + cnt_interp = 0 + do g = begg,endg + ! First: latitude and longitude smaller than observation location + if((longxy_obs_floor(i) == longxy(g-begg+1)) .and. (latixy_obs_floor(i) == latixy(g-begg+1))) then + + obs_interp_indices_p(cnt, 1) = g-begg+1 + ((endg-begg+1) * (layer_obs(i)-1)) + obs_interp_weights_p(cnt, 1) = sqrt(abs(lon(g)-lon_obs(i)) * & + abs(lon(g)-lon_obs(i)) + & + abs(lat(g)-lat_obs(i)) * & + abs(lat(g)-lat_obs(i))) + cnt_interp = cnt_interp + 1 + end if + ! Second: latitude larger than observation location, longitude smaller than observation location + if((longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs_floor(i) == latixy(g-begg+1))) then + obs_interp_indices_p(cnt, 2) = g-begg+1 + ((endg-begg+1) * (layer_obs(i)-1)) + obs_interp_weights_p(cnt, 2) = sqrt(abs(lon(g)-lon_obs(i)) * & + abs(lon(g)-lon_obs(i)) + & + abs(lat(g)-lat_obs(i)) * & + abs(lat(g)-lat_obs(i))) + cnt_interp = cnt_interp + 1 + end if + ! Third: latitude smaller than observation location, longitude larger than observation location + if((longxy_obs_floor(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1))) then + obs_interp_indices_p(cnt, 3) = g-begg+1 + ((endg-begg+1) * (layer_obs(i)-1)) + obs_interp_weights_p(cnt, 3) = sqrt(abs(lon(g)-lon_obs(i)) * & + abs(lon(g)-lon_obs(i)) + & + abs(lat(g)-lat_obs(i)) * & + abs(lat(g)-lat_obs(i))) + cnt_interp = cnt_interp + 1 + end if + ! Fourth: latitude and longitude larger than observation location + if((longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1))) then + obs_interp_indices_p(cnt, 4) = g-begg+1 + ((endg-begg+1) * (layer_obs(i)-1)) + obs_interp_weights_p(cnt, 4) = sqrt(abs(lon(g)-lon_obs(i)) * & + abs(lon(g)-lon_obs(i)) + & + abs(lat(g)-lat_obs(i)) * & + abs(lat(g)-lat_obs(i))) + cnt_interp = cnt_interp + 1 + end if + ! Check if all four corners are found + if(cnt_interp == 4) then + cnt = cnt + 1 + ! exit + end if + end do + end do + + do i = 1, dim_obs + + ! Sum of distance weights + sum_interp_weights = sum(obs_interp_weights_p(i, :)) + + do j = 1, 4 + ! Normalize distance weights + obs_interp_weights_p(i, j) = obs_interp_weights_p(i, j) / sum_interp_weights + end do + end do + + end if + +#ifdef PDAF_DEBUG + IF (da_print_obs_index > 0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn, "(a,i5.5,a,i5.5,a)") "obs_index_p_", mype_world, ".", step, ".txt" + OPEN(unit=71, file=fn, action="write") + DO i = 1, dim_obs_p + WRITE (71,"(i10)") obs_index_p(i) + END DO + CLOSE(71) + END IF +#endif + + + ! **************************************** + ! *** Gather global observation arrays *** + ! **************************************** + + CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & + thisobs%ncoord, cradius_SM, dim_obs) + + ! ******************** + ! *** Finishing up *** + ! ******************** + + ! Deallocate all local arrays + DEALLOCATE(obs_g) + DEALLOCATE(obs_p, ocoord_p, ivar_obs_p) + + END SUBROUTINE init_dim_obs_SM + + + + !------------------------------------------------------------------------------- + !> Implementation of observation operator + !! + !! This routine applies the full observation operator + !! for the type of observations handled in this module. + !! + !! One can choose a proper observation operator from + !! PDAFOMI_OBS_OP or add one to that module or + !! implement another observation operator here. + !! + !! The routine is called by all filter processes. + !! + SUBROUTINE obs_op_SM(dim_p, dim_obs, state_p, ostate) + + use mod_assimilation, only: obs_index_p + + use PDAFomi_obs_f, only: PDAFomi_gather_obsstate + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< Dimension of full observed state (all observed fields) + REAL, INTENT(in) :: state_p(dim_p) !< PE-local model state + REAL, INTENT(inout) :: ostate(dim_obs) !< Full observed state + + real, allocatable :: ostate_p(:) + + integer :: i + + + + ! ****************************************************** + ! *** Apply observation operator H on a state vector *** + ! ****************************************************** + + IF (thisobs%dim_obs_p>0) THEN + if (allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(thisobs%dim_obs_p)) + ELSE + if (allocated(ostate_p)) deallocate(ostate_p) + ALLOCATE(ostate_p(1)) + END IF + + + DO i = 1, thisobs%dim_obs_p + ostate_p(i) = state_p(obs_index_p(i)) + END DO + + ! *** Global: Gather full observed state vector + CALL PDAFomi_gather_obsstate(thisobs, ostate_p, ostate) + + deallocate(ostate_p) + + + END SUBROUTINE obs_op_SM + + + + !------------------------------------------------------------------------------- + !> Initialize local information on the module-type observation + !! + !! The routine is called during the loop over all local + !! analysis domains. It has to initialize the information + !! about local observations of the module type. It returns + !! number of local observations of the module type for the + !! current local analysis domain in DIM_OBS_L and the full + !! and local offsets of the observation in the overall + !! observation vector. + !! + !! This routine calls the routine PDAFomi_init_dim_obs_l + !! for each observation type. The call allows to specify a + !! different localization radius and localization functions + !! for each observation type and local analysis domain. + !! + SUBROUTINE init_dim_obs_l_SM(domain_p, step, dim_obs, dim_obs_l) + + ! Include PDAFomi function + USE PDAFomi, ONLY: PDAFomi_init_dim_obs_l, pi + + ! Include localization radius and local coordinates + USE mod_assimilation, & + ONLY: cradius_SM, locweight, sradius_SM, screen + + USE enkf_clm_mod, ONLY: state_loc2clm_c_p, clmstatevec_allcol, clmstatevec_only_active + + use shr_kind_mod, only: r8 => shr_kind_r8 + + use decompMod , only : get_proc_bounds + + +#ifdef CLMFIVE + USE GridcellType, ONLY: grc + USE ColumnType, ONLY : col +#else + USE clmtype, ONLY : clm3 +#endif + use clm_varcon, only: spval + + use mod_parallel_pdaf, & + ONLY: mype_world + + + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: domain_p !< Index of current local analysis domain + INTEGER, INTENT(in) :: step !< Current time step + INTEGER, INTENT(in) :: dim_obs !< Full dimension of observation vector + INTEGER, INTENT(inout) :: dim_obs_l !< Local dimension of observation vector + + REAL :: coords_l(2) ! Coordinates of local analysis domain + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + integer, pointer :: mycgridcell(:) !Pointer for CLM3.5/CLM5.0 col->gridcell index arrays + +#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 + + + ! ********************************************** + ! *** Initialize local observation dimension *** + ! ********************************************** + ! count observations within a radius + + if (thisobs%infile==1) then + + + if (clmstatevec_allcol==0 .and. clmstatevec_only_active==0) then + if (lon(state_loc2clm_c_p(domain_p))>180) then + coords_l(1) = lon(state_loc2clm_c_p(domain_p)) - 360.0 + else + coords_l(1) = lon(state_loc2clm_c_p(domain_p)) + end if + + else + + ! get coords_l --> coordinates of local analysis domain + if (lon(mycgridcell(state_loc2clm_c_p(domain_p)))>180) then + ! if SM should be assimilated. Else, state_loc2clm_c_p is not allocated + coords_l(1) = lon(mycgridcell(state_loc2clm_c_p(domain_p))) - 360.0 + else + coords_l(1) = lon(mycgridcell(state_loc2clm_c_p(domain_p))) + end if + coords_l(2) = lat(mycgridcell(state_loc2clm_c_p(domain_p))) + + end if + + if (thisobs%disttype==3) then ! if haversine formula in distance calculation, the coordinates have to be converted to radians + coords_l(1) = coords_l(1) * pi / 180.0 + coords_l(2) = coords_l(2) * pi / 180.0 + end if + + else + + coords_l(1) = spval + coords_l(2) = spval + + end if + + ! for disttype=3, the cradius and sradius have to passed in meters, + ! so I multiply by 1000 to be able to put it in km in the input file + + if (thisobs%disttype==3) then + CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, & + locweight, cradius_SM*1000.0, sradius_SM*1000.0, dim_obs_l) + else + CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, & + locweight, cradius_SM, sradius_SM, dim_obs_l) + end if + + + END SUBROUTINE init_dim_obs_l_SM + + + + !------------------------------------------------------------------------------- + !> Perform covariance localization for local EnKF on the module-type observation + !! + !! The routine is called in the analysis step of the localized + !! EnKF. It has to apply localization to the two matrices + !! HP and HPH of the analysis step for the module-type + !! observation. + !! + !! This routine calls the routine PDAFomi_localize_covar + !! for each observation type. The call allows to specify a + !! different localization radius and localization functions + !! for each observation type. + !! + SUBROUTINE localize_covar_SM(dim_p, dim_obs, HP_p, HPH, coords_p) + + ! Include PDAFomi function + USE PDAFomi, ONLY: PDAFomi_localize_covar + + ! Include localization radius and local coordinates + USE mod_assimilation, & + ONLY: cradius_SM, locweight, sradius_SM + + use enkf_clm_mod, only: state_pdaf2clm_c_p + + use shr_kind_mod, only: r8 => shr_kind_r8 + + USE GridcellType, ONLY: grc + USE ColumnType, ONLY : col + + IMPLICIT NONE + + ! *** Arguments *** + INTEGER, INTENT(in) :: dim_p !< PE-local state dimension + INTEGER, INTENT(in) :: dim_obs !< Dimension of observation vector + REAL, INTENT(inout) :: HP_p(dim_obs, dim_p) !< PE local part of matrix HP + REAL, INTENT(inout) :: HPH(dim_obs, dim_obs) !< Matrix HPH + REAL, INTENT(inout) :: coords_p(:,:) !< Coordinates of state vector elements + + integer :: i + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + integer, pointer :: mycgridcell(:) !Pointer for CLM3.5/CLM5.0 col->gridcell index arrays + +#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 + + + ! ************************************* + ! *** Apply covariance localization *** + ! ************************************* + + do i = 1,dim_p + if (lon(mycgridcell(state_pdaf2clm_c_p(i)))>180) then + coords_p(1,i) = lon(mycgridcell(state_pdaf2clm_c_p(i))) - 360.0 + else + coords_p(1,i) = lon(mycgridcell(state_pdaf2clm_c_p(i))) + end if + coords_p(2,i) = lat(mycgridcell(state_pdaf2clm_c_p(i))) + end do + + + + CALL PDAFomi_localize_covar(thisobs, dim_p, locweight, cradius_SM, sradius_SM, & + coords_p, HP_p, HPH) + + END SUBROUTINE localize_covar_SM + + + subroutine add_obs_err_SM(step, dim_obs, C) + + USE mod_parallel_pdaf, & + ONLY: npes_filter + + use PDAFomi, only: obsdims + + implicit none + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector + REAL, INTENT(inout) :: C(dim_obs,dim_obs) ! Matrix to that + ! observation covariance R is added + integer :: i, pe, cnt + INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector + INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + ! Initialize indices --> we only have information about local obs. dims per PE, so we get the global indices, more generalizable than using + ! the arrays initiliazed in init_dim_obs_SM as we can also consider different observation types in one observation file. Arrays from init_dim_obs_pdaf + ! (e.g. obs_nc2pdaf) may not be necessary anymore, @ Johannes, please have a check here., see also in PDAFomi_obs_f.F90, there the same code is used + ! addition: I also use now the obs_pdaf2nc for reordering the observation covariance matrix to the PDAF internal order + ! So for an obs type where correlations should be accounted for, this should not be removed! + + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + C(i,i) = C(i,i) + 1.0/thisobs%ivar_obs_f(cnt) + cnt = cnt + 1 + end do + end do + + DEALLOCATE(id_start, id_end) + + end subroutine add_obs_err_SM + + + subroutine init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) + + USE mod_parallel_pdaf, & + ONLY: npes_filter + + use PDAFomi, only: obsdims, map_obs_id + + implicit none + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector + REAL, INTENT(inout) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix + REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PE-local observation vector + LOGICAL, INTENT(inout) :: isdiag ! Whether the observation error covar. matrix is diagonal + + integer :: i, pe, cnt + INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector + INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector + + ALLOCATE(id_start(npes_filter), id_end(npes_filter)) + + ! Initialize indices --> we only have information about local obs. dims per PE, so we use the same logic as in add_obs_err_SM + pe = 1 + id_start(1) = 1 + IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) + id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 + DO pe = 2, npes_filter + id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) + IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) + id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 + END DO + + ! Initialize mapping vector (to be used in PDAF_enkf_obs_ensemble) --> has to be initialized here, else there will be errors! + cnt = 1 + IF (thisobs%obsid-1 > 0) cnt = cnt+ SUM(obsdims(:,1:thisobs%obsid-1)) + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + map_obs_id(i) = cnt + cnt = cnt + 1 + END DO + END DO + + cnt = 1 + DO pe = 1, npes_filter + DO i = id_start(pe), id_end(pe) + covar(i, i) = covar(i, i) + 1.0/thisobs%ivar_obs_f(cnt) ! the inverse of the observation variance is saved for each observation, so we do not need any other + ! array here. As we initiliazed the indices for each process, we also can just take index cnt instead of complicated mapping between nc and pdaf indices + cnt = cnt + 1 + ENDDO + ENDDO + + ! The matrix is diagonal + ! This setting avoids the computation of the SVD of COVAR + ! in PDAF_enkf_obs_ensemble + isdiag = .TRUE. + + DEALLOCATE(id_start, id_end) + + + end subroutine init_obscovar_SM + + + subroutine prodRinvA_SM(step, dim_obs_p, rank, obs_p, A_p, C_p) + + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_p(dim_obs_p) ! PE-local vector of observations + REAL, INTENT(in) :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine + REAL, INTENT(inout) :: C_p(dim_obs_p,rank) ! Output matrix + + INTEGER :: i, j ! index of observation component + INTEGER :: off ! row offset in A_l and C_l + + off = thisobs%off_obs_f + + do j=1, rank + do i=1, thisobs%dim_obs_f + C_p(i+off, j) = thisobs%ivar_obs_f(i) * A_p(i+off, j) + END DO + end do + + end subroutine prodRinvA_SM + + + subroutine prodRinvA_l_SM(domain_p, step, dim_obs, rank, obs_l, A_l, C_l) + + use shr_kind_mod, only: r8 => shr_kind_r8 + USE mod_assimilation, & + ONLY: cradius_SM, locweight, sradius_SM + use pdafomi, only: PDAFomi_observation_localization_weights + + implicit none + + INTEGER, INTENT(in) :: domain_p ! Current local analysis domain + INTEGER, INTENT(in) :: step ! Current time step + INTEGER, INTENT(in) :: dim_obs ! Dimension of local observation vector, multiple observation types possible, then we have to access with thisobs_l%dim_obs_l + INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix + REAL, INTENT(in) :: obs_l(dim_obs) ! Local vector of observations + REAL, INTENT(inout) :: A_l(dim_obs, rank) ! Input matrix from analysis routine + REAL, INTENT(out) :: C_l(dim_obs, rank) ! Output matrix + + INTEGER :: verbose ! Verbosity flag + INTEGER :: verbose_w ! Verbosity flag for weight computation + INTEGER, SAVE :: domain_save = -1 ! Save previous domain index + INTEGER :: wtype ! Type of weight function + INTEGER :: rtype ! Type of weight regulation + REAL, ALLOCATABLE :: weight(:) ! Localization weights + REAL, ALLOCATABLE :: A_obs(:,:) ! Array for a single row of A_l + REAL :: var_obs ! Variance of observation error + + INTEGER :: i, j + + INTEGER :: off ! row offset in A_l and C_l + INTEGER :: idummy ! Dummy to access nobs_all + + real(r8) :: ivariance_obs + + + off = thisobs_l%off_obs_l + idummy = dim_obs + + IF ((domain_p <= domain_save .OR. domain_save < 0) .AND. mype_filter==0) THEN + verbose = 1 + ELSE + verbose = 0 + END IF + domain_save = domain_p + + ! Screen output + IF (verbose == 1) THEN + WRITE (*, '(8x, a, f12.3)') & + '--- Use global rms for observations of ', rms_obs_SM + WRITE (*, '(8x, a, 1x)') & + '--- Domain localization' + WRITE (*, '(12x, a, 1x, f12.2)') & + '--- Local influence radius', cradius_SM + + IF (locweight > 0) THEN + WRITE (*, '(12x, a)') & + '--- Use distance-dependent weight for observation errors' + + IF (locweight == 3) THEN + write (*, '(12x, a)') & + '--- Use regulated weight with mean error variance' + ELSE IF (locweight == 4) THEN + write (*, '(12x, a)') & + '--- Use regulated weight with single-point error variance' + END IF + END IF + ENDIF + + ALLOCATE(weight(thisobs_l%dim_obs_l)) + call PDAFomi_observation_localization_weights(thisobs_l, thisobs, rank, A_l, & + weight, verbose) + + do j=1,rank + do i=1,thisobs_l%dim_obs_l + C_l(i+off,j) = thisobs_l%ivar_obs_l(i) * weight(i) * A_l(i+off, j) + end do + end do + + deallocate(weight) + + end subroutine prodRinvA_l_SM + + + END MODULE obs_SM_pdafomi +#endif diff --git a/interface/model/Makefile b/interface/model/Makefile index 185b906aa..580b8dc99 100644 --- a/interface/model/Makefile +++ b/interface/model/Makefile @@ -72,7 +72,7 @@ $(info $$OBJ is [${OBJ}]) all: libmodel.a libmodel.a: $(OBJ) - ar rcs $@ *.o + ar rcs $@ $(OBJ) ranlib $@ mv $@ $(TSMPPDAFLIBDIR) @echo "library compilation suceeded" @@ -93,4 +93,4 @@ wrapper_tsmp.o: $(PREP_C) $(FC) $(OPT) $(FCPP_FLAGS) $(FFLAGS) $(FINCS) -c $< -o $@ clean: - rm -f *.o *.mod $(PROG) libmodel.a + rm -f *.o *.mod $(PROG) libmodel.a diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index f943ad903..407f3a7d9 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -43,6 +43,7 @@ extern void clm_advance(int *ntstep, int *tstartcycle, int *mype); extern void update_clm(int *tstartcycle, int *mype); #if defined CLMSA extern void print_update_clm(int *ts, int *ttot); +extern void print_inc_clm(); #endif extern void write_clm_statistics(int *ts, int *ttot); extern void clm_finalize(); @@ -59,6 +60,7 @@ GLOBAL char pfoutfile_stat[500]; GLOBAL char pfproblemname[100]; GLOBAL char clminfile[100*2]; GLOBAL char outdir[100]; +GLOBAL char mean_filename[100]; /* integers */ GLOBAL int nprocpf; @@ -88,6 +90,7 @@ GLOBAL int nx_local,ny_local,nz_local; GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; +GLOBAL int clmupdate_tws; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; @@ -108,6 +111,12 @@ GLOBAL int pf_aniso_use_parflow; GLOBAL int is_dampfac_state_time_dependent; GLOBAL int is_dampfac_param_time_dependent; GLOBAL int pf_dampswitch_sm; +GLOBAL int TWS_smoother; +GLOBAL int state_setup; +GLOBAL int update_snow; +GLOBAL int remove_mean; +GLOBAL int exclude_greenland; +GLOBAL int set_zero_start; GLOBAL int crns_flag; GLOBAL int da_print_obs_index; extern int model; @@ -139,3 +148,4 @@ GLOBAL double dampfac_state_time_dependent; GLOBAL double dampfac_param_time_dependent; GLOBAL double da_crns_depth_tol; GLOBAL double clmcrns_bd; +GLOBAL double max_inc; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index e654b4190..c174cdfd5 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -88,6 +88,7 @@ void read_enkfpar(char *parname) clmt_printensemble = iniparser_getint(pardict,"CLM:t_printensemble",-2); clmwatmin_switch = iniparser_getint(pardict,"CLM:watmin_switch",0); clmswc_mask_snow = iniparser_getint(pardict,"CLM:swc_mask_snow",0); + clmupdate_tws = iniparser_getint(pardict,"CLM:update_tws",0); /* get settings for COSMO */ nproccosmo = iniparser_getint(pardict,"COSMO:nprocs",0); @@ -105,6 +106,15 @@ void read_enkfpar(char *parname) screen_wrapper = iniparser_getint(pardict,"DA:screen_wrapper",1); point_obs = iniparser_getint(pardict,"DA:point_obs",1); obs_interp_switch = iniparser_getint(pardict,"DA:obs_interp_switch",0); + + max_inc = iniparser_getdouble(pardict,"DA:max_inc",1.0); + TWS_smoother = iniparser_getint(pardict,"DA:TWS_smoother",0); + state_setup = iniparser_getint(pardict,"DA:state_setup",0); + update_snow = iniparser_getint(pardict,"DA:update_snow",0); + remove_mean = iniparser_getint(pardict,"DA:remove_mean",0); + exclude_greenland = iniparser_getint(pardict,"DA:exclude_greenland",0); + set_zero_start = iniparser_getint(pardict,"DA:set_zero_start",0); + crns_flag = iniparser_getint(pardict,"DA:crns_flag",0); da_crns_depth_tol = iniparser_getdouble(pardict,"DA:da_crns_depth_tol",0.01); clmcrns_bd = iniparser_getdouble(pardict, "DA:crns_bd", -1.0); diff --git a/interface/model/eclm/enkf_clm_5.F90 b/interface/model/eclm/enkf_clm_5.F90 index 873c9722a..59127e019 100644 --- a/interface/model/eclm/enkf_clm_5.F90 +++ b/interface/model/eclm/enkf_clm_5.F90 @@ -72,6 +72,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") #if defined CLMSA use enkf_clm_mod, only: define_clm_statevec #endif + use clm_varcon, only: averaging_var !!<< TSMP PDAF addition end implicit none @@ -185,6 +186,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") callcount=0) #if defined CLMSA + averaging_var=0 call define_clm_statevec(mype) #endif @@ -201,6 +203,8 @@ end subroutine clm_init !-------------------------------------------------------------------------- subroutine clm_advance(ntstep, tstartcycle, mype) bind(C,name="clm_advance") use cime_comp_mod, only : cime_run + use enkf_clm_mod, only : cleanup_clm_statevec + use enkf_clm_mod, only : define_clm_statevec use enkf_clm_mod, only : set_clm_statevec use iso_C_binding, only : c_int @@ -216,6 +220,11 @@ subroutine clm_advance(ntstep, tstartcycle, mype) bind(C,name="clm_advance") call cime_run(ntstep) #if defined CLMSA + ! TODO: Get the use_omi information here as IF-condition + call cleanup_clm_statevec() ! cleanup before defining statevec + call define_clm_statevec(mype) ! call define statevec not in the beginning + ! but here as we can define the statevec for each obs type + ! Calling PDAF Function to set state vector before assimiliation call set_clm_statevec(tstartcycle, mype) #endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 81ce1fd2f..58d203e06 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -56,6 +56,43 @@ module enkf_clm_mod 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="clmprint_swc") :: clmprint_swc + + ! Yorck + integer(c_int),bind(C,name="clmupdate_tws") :: clmupdate_tws + integer(c_int),bind(C,name="exclude_greenland") :: exclude_greenland + integer, dimension(1:5) :: clm_varsize_tws + real(r8),bind(C,name="max_inc") :: max_inc + integer(c_int),bind(C,name="TWS_smoother") :: TWS_smoother + integer(c_int),bind(C,name="state_setup") :: state_setup + integer(c_int),bind(C,name="set_zero_start") :: set_zero_start + integer, allocatable :: num_layer(:) + integer, allocatable :: num_layer_columns(:) + integer :: num_hactiveg, num_hactivec + + integer, allocatable :: hactiveg_levels(:,:) ! hydrolocial active filter for all levels (gridcell) + integer, allocatable :: hactivec_levels(:,:) ! hydrolocial active filter for all levels (column) + integer, allocatable :: gridcell_state(:) + + logical :: first_cycle = .TRUE. + + + ! OMI --> I want to update the observation type after each observation comes in. + ! problem: the observation type is updated before the update of the assimilation + ! this causes the wrong observation type to be used in the update + + ! idea: the state vector is newly initilized for each assimilation time step for the current variable + ! use a new variable, e.g. obs_type_update_tws, and make it equal to clmupdate_tws at the beginning of the initialization + ! this variable is then used in the update of the state vector + + integer :: obs_type_update_swc = 0 + integer :: obs_type_update_tws = 0 + integer :: obs_type_update_T = 0 + integer :: obs_type_update_texture = 0 + + + + ! end Yorck + #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol @@ -92,6 +129,8 @@ module enkf_clm_mod subroutine define_clm_statevec(mype) use decompMod , only : get_proc_bounds use clm_varpar , only : nlevsoi + use clm_varcon, only: set_averaging_to_zero + use PDAF_interfaces_module, only: PDAF_reset_dim_p implicit none @@ -118,6 +157,18 @@ subroutine define_clm_statevec(mype) clm_begp = begp clm_endp = endp + clm_statevecsize = 0 + clm_varsize = 0 + + ! check which observation type should be assimilated and design the state vector accordingly + ! I will introduce functions for each observation type so that other types can easily be included + + ! make update variable equal to the clmupdate variable + obs_type_update_swc = clmupdate_swc + obs_type_update_tws = clmupdate_tws + obs_type_update_T = clmupdate_T + obs_type_update_texture = clmupdate_texture + ! soil water content observations - case 1 if(clmupdate_swc==1) then call define_clm_statevec_swc(mype) @@ -144,13 +195,41 @@ subroutine define_clm_statevec(mype) end if !end hcp + ! TWS observations + if (clmupdate_tws==1) then + call define_clm_statevec_tws(mype) + end if + + ! + + ! Include your own state vector definitions here for different variables (ET, LAI, etc.) + + ! + + ! reset PDAF dimensions for multivariate assimilation, only not for first call as PDAF did not initalize yet + if (.not. first_cycle) then + call PDAF_reset_dim_p(clm_statevecsize,ierror) + end if + + if (first_cycle) then + ! possibility to assimilate GRACE not in the first month --> enkfpf.par file has information set_zero_start where the running average should be resetted + ! This is usually one month prior to the first GRACE observation. If it is not included in the file, it is resetted when the first GRACE observation + ! is assimilated. Afterwards, the normal set_zero information inside the observation file is used (see next_observation_pdaf for details). + if (set_zero_start/=0) then + set_averaging_to_zero = set_zero_start + end if + end if + + first_cycle = .FALSE. + + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, "define_clm_statevec: clm_statevecsize=", clm_statevecsize #endif IF (allocated(clm_statevec)) deallocate(clm_statevec) - if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0)) then + if ((clmupdate_swc/=0) .or. (clmupdate_T/=0) .or. (clmupdate_texture/=0) .or. (clmupdate_tws/=0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) end if @@ -160,7 +239,7 @@ subroutine define_clm_statevec(mype) ! values used in computing increments during updating the state ! vector in column-mean-mode. IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) - if (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) then + if ( (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_tws/=0 ) then allocate(clm_statevec_orig(clm_statevecsize)) end if @@ -170,6 +249,10 @@ subroutine define_clm_statevec(mype) error stop "Not implemented clmupdate_T.NE.0" end if + if (allocated(gridcell_state)) deallocate(gridcell_state) + allocate(gridcell_state(clm_statevecsize)) + + end subroutine define_clm_statevec @@ -216,11 +299,11 @@ subroutine define_clm_statevec_swc(mype) ! 1) COL/GRC: CLM->PDAF IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + allocate(state_clm2pdaf_p(clm_begc:clm_endc,nlevsoi)) do i=1,nlevsoi do c=clm_begc,clm_endc ! Default: inactive - state_clm2pdaf_p = ispval + state_clm2pdaf_p(c,i) = ispval end do end do @@ -367,6 +450,232 @@ subroutine define_clm_statevec_swc(mype) end subroutine define_clm_statevec_swc + !> @author Yorck Ewerdwalbesloh, Anne Springer + !> @date 29.10.2025 + !> @brief define the state vector for TWS assimilation + !> @param[in] mype MPI rank + subroutine define_clm_statevec_tws(mype) + use shr_kind_mod, only: r8 => shr_kind_r8 + use decompMod , only : get_proc_bounds + use clm_varpar , only : nlevsoi + use ColumnType , only : col + use PatchType, only: patch + use GridcellType, only: grc + + implicit none + + integer,intent(in) :: mype + + integer :: i + integer :: j + integer :: c + integer :: g + integer :: cc + + integer :: fa, fg + + integer :: begp, endp ! per-proc beginning and ending pft indices + integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begl, endl ! per-proc beginning and ending landunit indices + integer :: begg, endg ! per-proc gridcell ending gridcell indices + + logical, allocatable :: found(:) + + real(r8), pointer :: lon(:) + real(r8), pointer :: lat(:) + + lon => grc%londeg + lat => grc%latdeg + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + +#ifdef PDAF_DEBUG + WRITE(*,"(a,i5,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a)") & + "TSMP-PDAF mype(w)=", mype, " define_clm_statevec, CLM5-bounds (g,l,c,p)----",& + begg,",",endg,",",begl,",",endl,",",begc,",",endc,",",begp,",",endp," -------" +#endif + + clm_begg = begg + clm_endg = endg + clm_begc = begc + clm_endc = endc + clm_begp = begp + clm_endp = endp + + + ! first we build a filter to determine which columns are active / are not active + ! we also build a gridcell filter for gridcell averges + + num_hactiveg = 0 + num_hactivec = 0 + + if (allocated(found)) deallocate(found) + allocate(found(clm_begg:clm_endg)) + found(clm_begg:clm_endg) = .false. + + if (allocated(num_layer)) deallocate(num_layer) + allocate(num_layer(1:nlevsoi)) + num_layer(1:nlevsoi) = 0 + + if (allocated(num_layer_columns)) deallocate(num_layer_columns) + allocate(num_layer_columns(1:nlevsoi)) + num_layer_columns(1:nlevsoi) = 0 + + do c = clm_begc, clm_endc ! find hydrological active cells + + g = col%gridcell(c) ! gridcell of column + + ! greenland can be excluded from the statevector + if ((exclude_greenland==0) .or. & + (.not.(lon(g)<330 .and. lon(g)>180 .and. lat(g)>55))) then + + ! if the column is hydrologically active, add it, if the corresponding + ! gridcell is not found before, add also the gridcell + if (col%hydrologically_active(c)) then + + if (.not. found(g)) then ! if the gridcell is not found before + + found(g) = .true. + + do j = 1,nlevsoi ! add gridcell for each layer until bedrock + ! get number in layers + + if (j<=col%nbedrock(c)) then + num_layer(j) = num_layer(j) + 1 + end if + + end do + + num_hactiveg = num_hactiveg + 1 + + end if + + do j = 1,nlevsoi ! add column for each layer until bedrock + ! get number in layers + + if (j<=col%nbedrock(c)) then + num_layer_columns(j) = num_layer_columns(j) + 1 + end if + + end do + + num_hactivec = num_hactivec + 1 + + end if + end if + end do + + + ! allocate matrices that will hold the index of each hydrologically active component (gridcell, column, patch) in each depth + if (allocated(hactiveg_levels)) deallocate(hactiveg_levels) + if (allocated(hactivec_levels)) deallocate(hactivec_levels) + allocate(hactiveg_levels(1:num_hactiveg,1:nlevsoi)) + allocate(hactivec_levels(1:num_hactivec,1:nlevsoi)) + + ! now we fill these things with the columns and gridcells so that we can access all active things later on + + do j = 1,nlevsoi + + ! has to be inside the for lopp, else, the hactiveg_levels is only filled for the first level + found(clm_begg:clm_endg) = .false. + fa = 0 + fg = 0 + + do c = clm_begc, clm_endc + g = col%gridcell(c) ! gridcell of column + + if ((exclude_greenland==0) .or. (.not.(lon(g)<330 .and. lon(g)>180 .and. lat(g)>55))) then + + if (col%hydrologically_active(c)) then + + if (.not. found(g)) then ! if the gridcell is not found before + + found(g) = .true. + + if (j<=col%nbedrock(c)) then + fg = fg+1 + hactiveg_levels(fg,j) = g + end if + + end if + + if (j<=col%nbedrock(c)) then + fa = fa + 1 + hactivec_levels(fa,j) = c + end if + + end if + + end if + + end do + end do + + + if (allocated(found)) deallocate(found) + + ! now we have an array for the columns and gridcells of interest that we can + ! use when we fill the statevector and distribute the update + ! now lets find out the dimension of the state vector + + clm_varsize_tws(:) = 0 + + clm_statevecsize = 0 + + select case (state_setup) + case(0) + ! all compartments in state vector, liq and ice added up to not run into balancing errors due to different partioning + ! in different ensemble members in these variables even if the total amount of water is similar + + do j = 1,nlevsoi + clm_varsize_tws(1) = clm_varsize_tws(1) + num_layer(j) + clm_statevecsize = clm_statevecsize + num_layer(j) + + clm_varsize_tws(2) = 0 ! as liq + ice is added up, we do not need another variable for ice + clm_statevecsize = clm_statevecsize + 0 + end do + + ! snow + clm_varsize_tws(3) = num_layer(1) + clm_statevecsize = clm_statevecsize + num_layer(1) + + ! surface water + clm_varsize_tws(4) = num_layer(1) + clm_statevecsize = clm_statevecsize + num_layer(1) + + ! canopy water + clm_varsize_tws(5) = 0 + clm_statevecsize = clm_statevecsize + 0 + + case(1) ! TWS in statevector + + clm_varsize_tws(1) = num_layer(1) + clm_statevecsize = clm_statevecsize + num_layer(1) + clm_varsize_tws(2) = 0 + clm_varsize_tws(3) = 0 + clm_varsize_tws(4) = 0 + clm_varsize_tws(5) = 0 + + case(2) ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + + clm_varsize_tws(1) = num_layer(1) + clm_statevecsize = clm_statevecsize + num_layer(1) + + clm_varsize_tws(2) = num_layer(4) + clm_statevecsize = clm_statevecsize + num_layer(4) + + clm_varsize_tws(3) = num_layer(13) + clm_statevecsize = clm_statevecsize + num_layer(13) + + ! one variable for snow + clm_varsize_tws(4) = num_layer(1) + clm_statevecsize = clm_statevecsize + num_layer(1) + + end select + + end subroutine define_clm_statevec_tws + + subroutine cleanup_clm_statevec() implicit none @@ -457,6 +766,10 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + if (clmupdate_tws==1) then + call set_clm_statevec_tws + end if + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -536,6 +849,277 @@ subroutine set_clm_statevec_swc() end subroutine set_clm_statevec_swc + !> @author Yorck Ewerdwalbesloh, Anne Springer + !> @date 29.10.2025 + !> @brief set the state vector for TWS assimilation + subroutine set_clm_statevec_tws() + use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_varpar , only : nlevsoi + use ColumnType , only : col + use shr_kind_mod, only: r8 => shr_kind_r8 + use PatchType, only: patch + use GridcellType, only: grc + implicit none + + integer :: j,g,cc,count,c,count_c + integer :: n_c + real(r8) :: avg_sum + + real(r8), pointer :: liq(:,:) ! liquid water (kg/m2) + real(r8), pointer :: ice(:,:) ! ice lens (kg/m2) + real(r8), pointer :: snow(:) ! snow water (mm) + real(r8), pointer :: sfc(:) ! surface water + real(r8), pointer :: can(:) ! canopy water + real(r8), pointer :: TWS(:) ! TWS + + real(r8), pointer :: tws_state(:) + real(r8), pointer :: liq_state(:,:) + real(r8), pointer :: ice_state(:,:) + real(r8), pointer :: snow_state(:) + + cc = 0 + + tws_state => waterstate_inst%tws_state_before + liq_state => waterstate_inst%h2osoi_liq_state_before + ice_state => waterstate_inst%h2osoi_ice_state_before + snow_state => waterstate_inst%h2osno_state_before + + select case (TWS_smoother) + + case(0) ! instantaneous values + liq => waterstate_inst%h2osoi_liq_col + ice => waterstate_inst%h2osoi_ice_col + snow => waterstate_inst%h2osno_col + sfc => waterstate_inst%h2osfc_col + can => waterstate_inst%h2ocan_patch + TWS => waterstate_inst%tws_hactive + case(1) ! monthly means + liq => waterstate_inst%h2osoi_liq_col_mean + ice => waterstate_inst%h2osoi_ice_col_mean + snow => waterstate_inst%h2osno_col_mean + sfc => waterstate_inst%h2osfc_col_mean + can => waterstate_inst%h2ocan_patch_mean + TWS => waterstate_inst%tws_hactive_mean + end select + + select case (state_setup) + case(0) ! all compartments in state vector + cc = 1 + do j = 1,nlevsoi + do count = 1, num_layer(j) + g = hactiveg_levels(count,j) + + clm_statevec(cc) = 0.0 + n_c = 0 + + do count_c = 1,num_layer_columns(j) ! get average over liq+ice for all columns in gridcell g + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + clm_statevec(cc) = clm_statevec(cc) + liq(c,j) + ice(c,j) + n_c = n_c + 1 + end if + end do + + clm_statevec(cc) = clm_statevec(cc) / real(n_c, r8) + clm_statevec_orig(cc) = clm_statevec(cc) + + liq_state(g,j) = clm_statevec(cc) + + gridcell_state(cc) = g + + if (j==1) then ! for the first layer, we can fill the other compartments + + ! snow + clm_statevec(cc+sum(clm_varsize_tws(1:2))) = 0.0 + n_c = 0 + + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + clm_statevec(cc+sum(clm_varsize_tws(1:2))) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + snow(c) + n_c = n_c + 1 + end if + end do + + clm_statevec(cc+sum(clm_varsize_tws(1:2))) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) / real(n_c, r8) + clm_statevec_orig(cc+sum(clm_varsize_tws(1:2))) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + + snow_state(g) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + gridcell_state(cc+sum(clm_varsize_tws(1:2))) = g + + ! surface water + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = 0.0 + n_c = 0 + + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) + sfc(c) + n_c = n_c + 1 + end if + + end do + + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) / real(n_c, r8) + clm_statevec_orig(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) + gridcell_state(cc+sum(clm_varsize_tws(1:3))) = g + end if + + cc = cc+1 + end do + end do + + case(1) ! TWS in state vector + + cc = 1 + + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + clm_statevec(cc) = TWS(g) + clm_statevec_orig(cc) = TWS(g) + tws_state(g) = clm_statevec(cc) + gridcell_state(cc) = g + cc = cc+1 + end do + + case(2) ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + + cc = 1 + + ! surface soil moisture + do count = 1, num_layer(1) + clm_statevec(cc) = 0 + g = hactiveg_levels(count,1) + + do j = 1, 3 ! surface SM + avg_sum = 0 + n_c = 0 + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + avg_sum = avg_sum + liq(c,j) + ice(c,j) + n_c = n_c+1 + end if + + end do + + if (n_c/=0) then + + clm_statevec(cc) = clm_statevec(cc) + avg_sum / real(n_c, r8) + + end if + + end do + + clm_statevec_orig(cc) = clm_statevec(cc) + liq_state(g,1) = clm_statevec(cc) + + gridcell_state(cc) = g + + cc = cc + 1 + + end do + + cc = 1 + + ! root zone soil moisture + do count = 1, num_layer(4) + clm_statevec(cc+clm_varsize_tws(1)) = 0 + g = hactiveg_levels(count,4) + do j = 4, 12 + avg_sum = 0 + n_c = 0 + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + avg_sum = avg_sum + liq(c,j) + ice(c,j) + n_c = n_c+1 + end if + end do + + if (n_c/=0) then + + clm_statevec(cc+clm_varsize_tws(1)) = clm_statevec(cc+clm_varsize_tws(1)) + avg_sum / real(n_c, r8) + + end if + + end do + + clm_statevec_orig(cc+clm_varsize_tws(1)) = clm_statevec(cc+clm_varsize_tws(1)) + liq_state(g,2) = clm_statevec(cc+clm_varsize_tws(1)) + gridcell_state(cc+clm_varsize_tws(1)) = g + + cc = cc + 1 + + end do + + cc = 1 + + ! deep soil moisture + do count = 1, num_layer(13) + clm_statevec(cc+sum(clm_varsize_tws(1:2))) = 0 + g = hactiveg_levels(count,13) + do j = 13, nlevsoi + avg_sum = 0 + n_c = 0 + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (g==col%gridcell(c)) then + avg_sum = avg_sum + liq(c,j) + ice(c,j) + n_c = n_c+1 + end if + end do + + if (n_c/=0) then + + clm_statevec(cc+sum(clm_varsize_tws(1:2))) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + avg_sum / real(n_c, r8) + + end if + + end do + + clm_statevec_orig(cc+sum(clm_varsize_tws(1:2))) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + liq_state(g,3) = clm_statevec(cc+sum(clm_varsize_tws(1:2))) + gridcell_state(cc+sum(clm_varsize_tws(1:2))) = g + + cc = cc + 1 + + end do + + cc = 1 + + ! snow + do count = 1, num_layer(1) + + g = hactiveg_levels(count,1) + + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = 0 + n_c = 0 + do count_c = 1,num_layer_columns(1) + c = hactivec_levels(count_c,1) + if (g==col%gridcell(c)) then + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) + snow(c) + n_c = n_c+1 + end if + end do + + clm_statevec(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) / real(n_c, r8) + clm_statevec_orig(cc+sum(clm_varsize_tws(1:3))) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) + + snow_state(g) = clm_statevec(cc+sum(clm_varsize_tws(1:3))) + gridcell_state(cc+sum(clm_varsize_tws(1:3))) = g + + cc = cc+1 + + end do + + end select + + + end subroutine set_clm_statevec_tws + + + subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_time_manager , only : update_DA_nstep @@ -577,7 +1161,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN - IF(clmupdate_swc/=0) THEN + IF(obs_type_update_swc/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn5, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn5, action="write") @@ -595,7 +1179,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") #endif ! calculate shift when CRP data are assimilated - if(clmupdate_swc==2) then + if(obs_type_update_swc==2) then error stop "Not implemented: clmupdate_swc.eq.2" endif @@ -605,21 +1189,27 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call update_DA_nstep() ! write updated swc back to CLM - if(clmupdate_swc/=0) then + if(obs_type_update_swc/=0) then call update_clm_swc(tstartcycle, mype) endif !hcp: TG, TV - if(clmupdate_T==1) then + if(obs_type_update_T==1) then error stop "Not implemented: clmupdate_T.eq.1" endif ! end hcp TG, TV ! write updated texture back to CLM - if(clmupdate_texture/=0) then + if(obs_type_update_texture/=0) then call update_clm_texture(tstartcycle, mype) endif + if (obs_type_update_tws==1) then + call clm_update_tws + end if + + + end subroutine update_clm @@ -644,6 +1234,7 @@ subroutine update_clm_swc(tstartcycle, mype) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) real(r8), pointer :: snow_depth(:) + real(r8), pointer :: liq_inc(:,:), ice_inc(:,:), snow_inc(:) real(r8) :: rliq,rice real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) @@ -667,6 +1258,11 @@ subroutine update_clm_swc(tstartcycle, mype) snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + liq_inc => waterstate_inst%h2osoi_liq_col_inc + ice_inc => waterstate_inst%h2osoi_ice_col_inc + snow_inc => waterstate_inst%h2osno_col_inc + + ! Set minimum soil moisture for checking the state vector and ! for setting minimum swc for CLM if(clmwatmin_switch==3) then @@ -683,6 +1279,19 @@ subroutine update_clm_swc(tstartcycle, mype) watmin_set = 0.0 end if + do i = 1,nlevsoi + do j = clm_begc,clm_endc + + liq_inc(j,i) = h2osoi_liq(j,i) + ice_inc(j,i) = h2osoi_ice(j,i) + + if (i==1) then + snow_inc(j) = waterstate_inst%h2osno_col(j) + end if + + end do + end do + ! cc = 0 do i=1,nlevsoi ! CLM3.5: iterate over grid cells @@ -797,6 +1406,19 @@ subroutine update_clm_swc(tstartcycle, mype) END IF #endif + do i = 1,nlevsoi + do j=clm_begc,clm_endc + + liq_inc(j,i) = h2osoi_liq(j,i)-liq_inc(j,i) + ice_inc(j,i) = h2osoi_ice(j,i)-ice_inc(j,i) + + if (i==1) then + snow_inc(j) = waterstate_inst%h2osno_col(j)-snow_inc(j) + end if + + end do + end do + end subroutine update_clm_swc @@ -841,6 +1463,395 @@ subroutine update_clm_texture(tstartcycle, mype) end subroutine update_clm_texture + !> @author Yorck Ewerdwalbesloh, Anne Springer + !> @date 29.10.2025 + !> @brief update water storages from TWS assimilation + subroutine clm_update_tws() + + use clm_varpar , only : nlevsoi, nlevsno + use shr_kind_mod, only: r8 => shr_kind_r8 + use clm_varcon, only: watmin, denh2o, denice, averaging_var + use ColumnType, only : col + use clm_instMod, only : soilstate_inst, waterstate_inst + + implicit none + + integer :: c, j + + real(r8), pointer :: liq(:,:), ice(:,:), swc(:,:) + real(r8), pointer :: dz(:,:), zi(:,:), z(:,:) + real(r8), pointer :: watsat(:,:), snow(:), snow_depth(:) + integer, pointer :: snl(:) + + real(r8), pointer :: liq_mean(:,:), ice_mean(:,:), snow_mean(:) + real(r8), pointer :: liq_inc(:,:), ice_inc(:,:), snow_inc(:) + real(r8), pointer :: TWS(:) + + ! Pointer declarations + call assign_pointers() + + ! set averaging factor to zero + averaging_var = 0 + + ! Initialize increments + call initialize_increments() + + select case (state_setup) + case(0) ! all compartments in state vector + call update_state_0() + case(1) ! TWS in state vector + call update_state_1() + case(2) ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + call update_state_2() + end select + + call finalize_increments() + + contains + + subroutine assign_pointers() + + liq => waterstate_inst%h2osoi_liq_col + ice => waterstate_inst%h2osoi_ice_col + swc => waterstate_inst%h2osoi_vol_col + dz => col%dz + zi => col%zi + z => col%z + watsat => soilstate_inst%watsat_col + snow => waterstate_inst%h2osno_col + snow_depth => waterstate_inst%snow_depth_col + snl => col%snl + + select case (TWS_smoother) + case(0) + liq_mean => waterstate_inst%h2osoi_liq_col + ice_mean => waterstate_inst%h2osoi_ice_col + snow_mean => waterstate_inst%h2osno_col + case default + liq_mean => waterstate_inst%h2osoi_liq_col_mean + ice_mean => waterstate_inst%h2osoi_ice_col_mean + snow_mean => waterstate_inst%h2osno_col_mean + end select + + liq_inc => waterstate_inst%h2osoi_liq_col_inc + ice_inc => waterstate_inst%h2osoi_ice_col_inc + snow_inc => waterstate_inst%h2osno_col_inc + + TWS => waterstate_inst%tws_hactive + + end subroutine assign_pointers + + subroutine initialize_increments() + integer :: j, count, c + ! set increment values for inc output file to old values, then they can be updated with new values + liq_inc(:,:) = 0.0 + ice_inc(:,:) = 0.0 + snow_inc(:) = 0.0 + do j = 1,nlevsoi + do count = 1,num_layer_columns(j) + c = hactivec_levels(count,j) + + liq_inc(c,j) = liq(c,j) + ice_inc(c,j) = ice(c,j) + + if (j==1) then + snow_inc(c) = snow(c) + end if + end do + end do + end subroutine initialize_increments + + subroutine finalize_increments() + integer :: j, count, c + do j = 1,nlevsoi + do count = 1,num_layer_columns(j) + c = hactivec_levels(count,j) + + liq_inc(c,j) = liq(c,j)-liq_inc(c,j) + ice_inc(c,j) = ice(c,j)-ice_inc(c,j) + + if (j==1) then + snow_inc(c) = snow(c)-snow_inc(c) + end if + end do + end do + end subroutine finalize_increments + + + + + subroutine update_state_0() + real(r8) :: inc + integer :: c, j, count, g, cc, count_c + + ! Update soil + cc = 1 + do j = 1,nlevsoi + do count = 1, num_layer(j) + g = hactiveg_levels(count,j) + inc = clm_statevec(cc)-clm_statevec_orig(cc) + + if (abs(inc)>1.e-10_r8) then + do count_c = 1, num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (col%gridcell(c)==g) then + + call update_soil_layer(c, j, inc, liq_mean(c,j), ice_mean(c,j), clm_statevec_orig(cc)) + + end if + end do + end if + cc = cc+1 + end do + end do + + ! update snow + cc = 1 + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + inc = clm_statevec(cc+sum(clm_varsize_tws(1:2))) - clm_statevec_orig(cc+sum(clm_varsize_tws(1:2))) + if (abs(inc)>1.e-10_r8) then + do count_c = 1, num_layer_columns(1) + c = hactivec_levels(count_c,1) + if (col%gridcell(c)==g) then + + if (snl(c) < 0) then + call scale_snow(c, inc*snow_mean(c)/clm_statevec_orig(cc+sum(clm_varsize_tws(1:2)))) + end if + + if (snow(c) < 0._r8) then ! close snow layers if snow is negative + call zero_snow_layers(c) + end if + + end if + end do + end if + cc = cc+1 + end do + + end subroutine update_state_0 + + + + subroutine update_state_1() + real(r8) :: inc + integer :: c, j, count, g, cc, count_c + + cc = 1 + + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + inc = clm_statevec(cc)-clm_statevec_orig(cc) + + if (abs(inc)>1.e-10_r8) then + ! update soil water + do j = 1,nlevsoi + + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (col%gridcell(c)==g) then + call update_soil_layer(c, j, inc, liq_mean(c,j), ice_mean(c,j), clm_statevec_orig(cc)) + end if + end do + + end do + + ! update snow + + do count_c = 1,num_layer_columns(1) + c = hactivec_levels(count_c,1) + if (col%gridcell(c)==g) then + if (snl(c) < 0) then ! snow layers in the column + call scale_snow(c, inc*snow_mean(c)/clm_statevec_orig(cc)) + end if + if (snow(c) < 0._r8) then ! close snow layers if snow is negative + call zero_snow_layers(c) + end if + end if + end do + end if + + cc = cc+1 + + end do + + end subroutine update_state_1 + + + subroutine update_state_2() + real(r8) :: inc + integer :: c, j, count, g, cc, count_c + + ! surface soil moisture + cc = 1 + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + inc = clm_statevec(cc)-clm_statevec_orig(cc) + if (abs(inc)>1.e-10_r8) then + do j = 1,3 + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (col%gridcell(c)==g) then + + call update_soil_layer(c, j, inc, liq_mean(c,j), ice_mean(c,j), clm_statevec_orig(cc)) + + end if + end do + end do + end if + cc = cc+1 + end do + + ! root zone soil moisture + cc = 1 + do count = 1, num_layer(4) + g = hactiveg_levels(count,4) + inc = clm_statevec(cc+clm_varsize_tws(1))-clm_statevec_orig(cc+clm_varsize_tws(1)) + if (abs(inc)>1.e-10_r8) then + do j = 4,12 + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (col%gridcell(c)==g) then + + call update_soil_layer(c, j, inc, liq_mean(c,j), ice_mean(c,j), clm_statevec_orig(cc+clm_varsize_tws(1))) + + end if + end do + end do + end if + cc = cc+1 + end do + + ! deep soil moisture + cc = 1 + do count = 1, num_layer(13) + g = hactiveg_levels(count,13) + inc = clm_statevec(cc+sum(clm_varsize_tws(1:2)))-clm_statevec_orig(cc+sum(clm_varsize_tws(1:2))) + if (abs(inc)>1.e-10_r8) then + do j = 13,nlevsoi + do count_c = 1,num_layer_columns(j) + c = hactivec_levels(count_c,j) + if (col%gridcell(c)==g) then + + call update_soil_layer(c, j, inc, liq_mean(c,j), ice_mean(c,j), clm_statevec_orig(cc+sum(clm_varsize_tws(1:2)))) + + end if + end do + end do + end if + cc = cc+1 + end do + + ! update snow + cc = 1 + do count = 1, num_layer(1) + g = hactiveg_levels(count,1) + inc = clm_statevec(cc+sum(clm_varsize_tws(1:3))) - clm_statevec_orig(cc+sum(clm_varsize_tws(1:3))) + if (abs(inc)>1.e-10_r8) then + do count_c = 1, num_layer_columns(1) + c = hactivec_levels(count_c,1) + if (col%gridcell(c)==g) then + + if (snl(c) < 0) then + call scale_snow(c, inc*snow_mean(c)/clm_statevec_orig(cc+sum(clm_varsize_tws(1:3)))) + end if + + if (snow(c) < 0._r8) then ! close snow layers if snow is negative + call zero_snow_layers(c) + end if + + end if + end do + end if + cc = cc+1 + end do + + end subroutine update_state_2 + + + + subroutine update_soil_layer(c, j, inc, liq_mean_cj, ice_mean_cj, vec_orig) + real(r8), intent(in) :: inc, liq_mean_cj, ice_mean_cj, vec_orig + integer, intent(in) :: c, j + real(r8) :: inc_col, var_temp + + inc_col = inc * liq_mean_cj / vec_orig + if (inc_col /= inc_col) inc_col = 0.0_r8 + if (abs(inc_col) > max_inc * liq(c,j)) inc_col = sign(max_inc * liq(c,j), inc_col) + liq(c,j) = max(watmin, liq(c,j) + inc_col) + + inc_col = inc * ice_mean_cj / vec_orig + if (inc_col /= inc_col) inc_col = 0.0_r8 + if (abs(inc_col) > max_inc * ice(c,j)) inc_col = sign(max_inc * ice(c,j), inc_col) + ice(c,j) = max(0._r8, ice(c,j) + inc_col) + + swc(c,j) = liq(c,j)/(dz(c,j)*denh2o) + ice(c,j)/(dz(c,j)*denice) + + if (j > 1 .and. swc(c,j) - watsat(c,j) > 0.00001_r8) then + var_temp = watsat(c,j) / swc(c,j) + liq(c,j) = max(watmin, liq(c,j) * var_temp) + ice(c,j) = max(0._r8, watsat(c,j)*dz(c,j)*denice - liq(c,j)*denice/denh2o) + swc(c,j) = watsat(c,j) + if (abs(ice(c,j)) < 1.e-10_r8) ice(c,j) = 0._r8 + end if + end subroutine update_soil_layer + + + subroutine scale_snow(c, inc) + integer, intent(in) :: c + real(r8), intent(in) :: inc + real(r8) :: scale, inc_col + integer :: j + + inc_col=inc + if (inc_col /= inc_col) inc_col = 0.0 + + if (abs(inc_col)>max_inc*snow(c)) inc_col = sign(max_inc*snow(c),inc_col) + inc_col = snow(c) + inc_col + scale = inc_col / snow(c) + snow(c) = inc_col + do j = 0, snl(c)+1, -1 + liq(c,j) = liq(c,j)*scale + ice(c,j) = ice(c,j)*scale + dz(c,j) = dz(c,j)*scale + zi(c,j) = zi(c,j)*scale + z(c,j) = z(c,j)*scale + end do + + zi(c,snl(c)) = zi(c,snl(c))*scale + snow_depth(c) = snow_depth(c)*scale + end subroutine scale_snow + + + subroutine zero_snow_layers(c) + integer, intent(in) :: c + integer :: j + + if (snl(c) < 0) then + do j = 0, snl(c)+1, -1 + liq(c,j) = 0.0_r8 + ice(c,j) = 1.e-8_r8 + dz(c,j) = 1.e-8_r8 + zi(c,j-1) = sum(dz(c,j:0)) * -1._r8 + if (j == 0) then + z(c,j) = zi(c,j-1) / 2._r8 + else + z(c,j) = sum(zi(c,j-1:j)) / 2._r8 + end if + end do + else + liq(c,0) = 0.0_r8 + ice(c,0) = 1.e-8_r8 + dz(c,0) = 1.e-8_r8 + zi(c,-1) = dz(c,0) * -1._r8 + z(c,0) = zi(c,-1) / 2._r8 + end if + + snow_depth(c) = sum(dz(c,-nlevsno+1:0)) + snow(c) = sum(ice(c,-nlevsno+1:0)) + end subroutine zero_snow_layers + + end subroutine clm_update_tws subroutine clm_correct_texture() @@ -1239,7 +2250,7 @@ subroutine get_interp_idx(lon_clmobs, lat_clmobs, dim_obs, longxy_obs_floor, lat end subroutine get_interp_idx #if defined CLMSA - !> @author Johannes Keller + !> @author Johannes Keller, adaptations by Yorck Ewerdwalbesloh !> @date 24.04.2025 !> @brief Set number of local analysis domains N_DOMAINS_P !> @details @@ -1275,12 +2286,14 @@ subroutine init_n_domains_clm(n_domains_p) ! -> DIM_L: number of layers in gridcell n_domains_p = endg - begg + 1 end if - else + elseif (clmupdate_tws/=1) then ! Process-local number of gridcells Default, possibly not tested ! for other updates except SWC n_domains_p = endg - begg + 1 end if + NOGRACE: if (clmupdate_tws/=1) then + ! If only_active: Use clm2pdaf to check which columsn/gridcells ! are inside. Possibly: number of columns/gridcells reduced by ! hydrologically inactive columns/gridcells. @@ -1363,10 +2376,16 @@ subroutine init_n_domains_clm(n_domains_p) ! Possibly: Warning when final n_domains_p actually excludes ! hydrologically inactive gridcells + else NOGRACE + + n_domains_p = num_hactiveg + + end if NOGRACE + end subroutine init_n_domains_clm - !> @author Wolfgang Kurtz, Johannes Keller + !> @author Wolfgang Kurtz, Johannes Keller, Yorck Ewerdwalbesloh !> @date 20.11.2017 !> @brief Set local state vector dimension DIM_L local PDAF filters !> @details @@ -1381,6 +2400,8 @@ subroutine init_dim_l_clm(domain_p, dim_l) integer, intent(out) :: dim_l integer :: nshift + integer :: g, i, count + if(clmupdate_swc==1) then if(clmstatevec_only_active == 1) then ! Compare nlevsoi to clmstatevec_max_layer and bedrock if @@ -1407,9 +2428,53 @@ subroutine init_dim_l_clm(domain_p, dim_l) dim_l = 3*nlevsoi + nshift endif + if (clmupdate_tws==1) then + dim_l = 0 + g = hactiveg_levels(domain_p,1) + + select case(state_setup) + case(0) ! subdivided setup + do i = 1,nlevsoi + do count = 1, num_layer(i) + if (g==hactiveg_levels(count,i)) then + ! I could also check with col%nbedrock but then I would need the column index and not the gridcell index + 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(1) ! only TWS in statevector + + dim_l=1 + + case(2) ! aggregated setup + + 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 + endif + end subroutine init_dim_l_clm - !> @author Wolfgang Kurtz, Johannes Keller + !> @author Wolfgang Kurtz, Johannes Keller, Yorck Ewerdwalbesloh !> @date 20.11.2017 !> @brief Set local state vector STATE_L from global state vector STATE_P !> @details @@ -1418,6 +2483,8 @@ end subroutine init_dim_l_clm !> Source is STATE_P, the global (PE-local) state vector. subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) + use ColumnType , only : col + implicit none INTEGER, INTENT(in) :: domain_p ! Current local analysis domain @@ -1430,23 +2497,101 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) INTEGER :: n_domain INTEGER :: nshift_p + integer :: sub, g, j + ! call init_n_domains_clm(n_domain) ! DO i = 0, dim_l-1 ! nshift_p = domain_p + i * n_domain ! state_l(i+1) = state_p(nshift_p) ! ENDDO - + NOGRACE: if (clmupdate_tws/=1) then ! Column index inside gridcell index domain_p DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index: i state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) END DO + else NOGRACE + + if (clm_varsize_tws(5)/=0) then + sub=3 + else + sub=2 + end if + + select case (state_setup) + case(0) ! all compartmens in state vector + g = hactiveg_levels(domain_p,1) + do i = 1, dim_l-sub + do j = 1, num_layer(i) + if (g==hactiveg_levels(j,i)) then + if (i == 1) then + state_l(i) = state_p(j) + 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 + sum(clm_varsize_tws(1:2))) + state_l(dim_l-1) = state_p(j + sum(clm_varsize_tws(1:3))) + state_l(dim_l) = state_p(j + sum(clm_varsize_tws(1:4))) + else + state_l(dim_l-1) = state_p(j + sum(clm_varsize_tws(1:2))) + state_l(dim_l) = state_p(j + sum(clm_varsize_tws(1:3))) + end if + end if + end do + + case(1) ! 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(2) ! ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + + 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) ! surface SM + ! snow, same indexing as clm_varsize_tws(2:3) = 0 when only surface layers present + state_l(dim_l) = state_p(j + sum(clm_varsize_tws(1: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)) ! root zone SM + 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 + sum(clm_varsize_tws(1:2))) ! deep SM + end if + end do + end if + + end select + + end if NOGRACE end subroutine g2l_state_clm - !> @author Wolfgang Kurtz, Johannes Keller + !> @author Wolfgang Kurtz, Johannes Keller, Yorck Ewerdwalbesloh !> @date 20.11.2017 !> @brief Update global state vector STATE_P from local state vector STATE_L !> @details @@ -1455,6 +2600,8 @@ end subroutine g2l_state_clm !> Source is STATE_L, the local vector. subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) + use ColumnType , only : col + implicit none INTEGER, INTENT(in) :: domain_p ! Current local analysis domain @@ -1467,6 +2614,8 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) INTEGER :: n_domain INTEGER :: nshift_p + integer :: sub, j, g + ! ! beg and end gridcell for atm ! call init_n_domains_clm(n_domain) @@ -1474,13 +2623,93 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) ! nshift_p = domain_p + i * n_domain ! state_p(nshift_p) = state_l(i+1) ! ENDDO - + NOGRACE: if (clmupdate_tws==0) then ! Column index inside gridcell index domain_p DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index i state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) = state_l(i) END DO + else NOGRACE + + if (clm_varsize_tws(5)/=0) then + sub=3 + else + sub=2 + end if + + select case (state_setup) + case(0) ! all compartments in state vector + g = hactiveg_levels(domain_p,1) + do i = 1, dim_l-sub + do j = 1, num_layer(i) ! i is the layer that we are in right now + ! if the counter is the gridcell of the local domain, we know the position in the statevector + if (g==hactiveg_levels(j,i)) then + if (i == 1) then ! if first layer + state_p(j) = state_l(i) ! first liquid water as it is first in the statevector + else + state_p(j + sum(num_layer(1:i-1))) = state_l(i) + 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_p(j + sum(clm_varsize_tws(1:2))) = state_l(dim_l-2) + state_p(j + sum(clm_varsize_tws(1:3))) = state_l(dim_l-1) + state_p(j + sum(clm_varsize_tws(1:4))) = state_l(dim_l) + else + state_p(j + sum(clm_varsize_tws(1:2))) = state_l(dim_l-1) + state_p(j + sum(clm_varsize_tws(1:3))) = state_l(dim_l) + end if + + end if + end do + + case(1) ! TWS in statevector + + g = hactiveg_levels(domain_p,1) + do j = 1, num_layer(1) + if (g==hactiveg_levels(j,1)) then + state_p(j) = state_l(1) + end if + end do + + case(2) ! snow and soil moisture aggregated over surface, root zone and deep soil moisture in state vector + + g = hactiveg_levels(domain_p,1) + + do j = 1, num_layer(1) + if (g==hactiveg_levels(j,1)) then + + state_p(j) = state_l(1) + state_p(j + sum(clm_varsize_tws(1:3))) = state_l(dim_l) + + end if + end do + + do j = 1, num_layer(4) + if (g==hactiveg_levels(j,4)) then + + state_p(j + clm_varsize_tws(1)) = state_l(2) + + end if + end do + + do j = 1, num_layer(13) + if (g==hactiveg_levels(j,13)) then + + state_p(j + sum(clm_varsize_tws(1:2))) = state_l(3) + + end if + end do + + end select + + endif NOGRACE end subroutine l2g_state_clm #endif diff --git a/interface/model/eclm/print_update_clm_5.F90 b/interface/model/eclm/print_update_clm_5.F90 index 09dbc6810..5cbefccd0 100644 --- a/interface/model/eclm/print_update_clm_5.F90 +++ b/interface/model/eclm/print_update_clm_5.F90 @@ -219,6 +219,299 @@ subroutine print_update_clm(ts,ttot) bind(C,name="print_update_clm") ! deallocate(clmstate_tmp_local) end subroutine print_update_clm + +subroutine print_inc_clm() bind(C,name="print_inc_clm") + + ! use iso_c_binding + use shr_kind_mod , only : r8 => shr_kind_r8 + use domainMod , only : ldomain + use clm_varpar , only : nlevsoi + use clm_varcon , only : nameg, spval + use decompmod , only : get_proc_global, get_proc_bounds, ldecomp, get_proc_total + use spmdmod , only : masterproc, npes, mpicom, iam + use clm_time_manager , only : get_nstep + use clm_instMod, only : soilhydrology_inst, waterstate_inst, atm2lnd_inst + use netcdf, only : nf90_create + use netcdf, only : NF90_CLOBBER + use netcdf, only : nf90_def_dim + use netcdf, only : nf90_def_var + use netcdf, only : NF90_FLOAT + use netcdf, only : nf90_enddef + use netcdf, only : nf90_open + use netcdf, only : NF90_WRITE + use netcdf, only : nf90_inq_varid + use netcdf, only : nf90_put_var + use netcdf, only : nf90_close + ! use cime_comp_mod + use ColumnType , only : col + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use mpi, only: mpi_gatherv + use mpi, only: mpi_real8 + + implicit none + + ! local variables + + integer :: numg ! total number of gridcells across all processors + integer :: numl ! total number of landunits across all processors + integer :: numc ! total number of columns across all processors + integer :: nump ! total number of pfts across all processors + integer :: begg,endg ! local beg/end gridcells gdc + integer :: begl,endl ! local beg/end landunits + integer :: begc,endc ! local beg/end columns + integer :: begp,endp ! local beg/end pfts + integer :: ncells ! total number of gridcells on the processor + integer :: nlunits ! total number of landunits on the processor + integer :: ncols ! total number of columns on the processor + integer :: npfts ! total number of pfts on the processor + integer :: ncohorts + + integer :: isec, info, jn, jj, ji, g1, jx, c, l, j, g, index, p, count, count2 ! temporary integer + real(r8), pointer :: h2osoi_liq(:,:) + real(r8), pointer :: h2osoi_ice(:,:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: clmstate_tmp_local(:,:) + real(r8), pointer :: clmstate_tmp_global(:) + real(r8), allocatable :: clmstate_out(:,:,:) + integer ,dimension(3) :: dimids + integer ,dimension(2) :: dimids_1level + integer ,dimension(1) :: il_var_id + integer :: il_file_id, ncvarid(4), status + character(len = 300) :: inc_filename + integer :: nerror + integer :: ndlon,ndlat + + integer :: ier !return code + integer :: beg !temporary + integer :: numrecvv(0:npes-1) !vector of items to be received + integer :: displsv(0:npes-1) !displacement vector + integer :: numsend !number of items to be sent + integer :: pid ! processor id + integer :: count_columns + real(r8) :: sum_columns + real(r8), allocatable :: tws_inc(:) + + h2osoi_liq => waterstate_inst%h2osoi_liq_col_inc + h2osoi_ice => waterstate_inst%h2osoi_ice_col_inc + h2osno => waterstate_inst%h2osno_col_inc + + call get_proc_global(ng=numg,nl=numl,nc=numc,np=nump) + call get_proc_bounds(begg,endg,begl,endl,begc,endc,begp,endp) + allocate(clmstate_tmp_local(begg:endg,1:nlevsoi), stat=nerror) + allocate(tws_inc(begg:endg), stat=nerror) + tws_inc(begg:endg) = 0._r8 + + ndlon = ldomain%ni + ndlat = ldomain%nj + + if (masterproc) then + + allocate(clmstate_tmp_global(1:numg), stat=nerror) + allocate(clmstate_out(ndlon,ndlat,nlevsoi), stat=nerror) + clmstate_out(:,:,:) = nan + + end if + + call get_proc_total(iam, ncells, nlunits, ncols, npfts, ncohorts) + + numsend = ncells + + do pid = 0,npes-1 + call get_proc_total(pid, ncells, nlunits, ncols, npfts, ncohorts) + numrecvv(pid) = ncells + end do + beg = begg + displsv(0) = 0 + do pid = 1,npes-1 + displsv(pid) = displsv(pid-1) + numrecvv(pid-1) + end do + + if(masterproc) then + call get_inc_filename(inc_filename) + status = nf90_create(inc_filename, NF90_CLOBBER, il_file_id) + status = nf90_def_dim(il_file_id, "lon", ndlon, dimids(1)) + status = nf90_def_dim(il_file_id, "lat", ndlat, dimids(2)) + status = nf90_def_dim(il_file_id, "z", nlevsoi, dimids(3)) + + dimids_1level = [ dimids(1), dimids(2) ] + status = nf90_def_var(il_file_id, "SOILLIQ", NF90_FLOAT, dimids, ncvarid(1)) + status = nf90_def_var(il_file_id, "SOILICE", NF90_FLOAT, dimids, ncvarid(2)) + status = nf90_def_var(il_file_id, "H2OSNO", NF90_FLOAT, dimids_1level, ncvarid(3)) + status = nf90_def_var(il_file_id, "TWS", NF90_FLOAT, dimids_1level, ncvarid(3)) + status = nf90_enddef(il_file_id) + end if + + clmstate_tmp_local(begg:endg,:) = 0._r8 + + do j = 1, nlevsoi + do g = begg,endg + count_columns = 0 + sum_columns = 0 + do c=begc,endc + if (g==col%gridcell(c) .and. col%hydrologically_active(c).and.j<=col%nbedrock(c)) then + sum_columns = sum_columns+h2osoi_liq(c,j) + count_columns = count_columns+1 + end if + end do + if (count_columns == 0) then + !clmstate_tmp_local(g,j) = spval + else + clmstate_tmp_local(g,j) = sum_columns/count_columns + end if + if (j==1) then + tws_inc(g) = clmstate_tmp_local(g,j) + else + if (clmstate_tmp_local(g,j) /= spval) then + tws_inc(g) = tws_inc(g) + clmstate_tmp_local(g,j) + end if + end if + end do + end do + do jn = 1, nlevsoi + if (masterproc) then + call mpi_gatherv (clmstate_tmp_local(beg,jn), numsend, MPI_REAL8, & + clmstate_tmp_global, numrecvv, displsv, & + MPI_REAL8, 0, mpicom, ier) + else + call mpi_gatherv (clmstate_tmp_local(beg,jn), numsend , MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier) + end if + if(masterproc) then + do g1 = 1, numg + ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1 + jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1 + clmstate_out(ji,jj,jn) = clmstate_tmp_global(g1) + end do + end if + end do + if (masterproc) then + status = nf90_inq_varid(il_file_id, "SOILLIQ" , ncvarid(1)) + status = nf90_put_var( il_file_id, ncvarid(1), clmstate_out(:,:,:), & + start = [ 1, 1, 1 ], count = [ ndlon, ndlat, nlevsoi ] ) + end if + + + clmstate_tmp_local(begg:endg,:) = 0._r8 + + do j = 1, nlevsoi + do g=begg,endg + count_columns = 0 + sum_columns = 0 + do c=begc,endc + if (g==col%gridcell(c).and.col%hydrologically_active(c).and.j<=col%nbedrock(c)) then + sum_columns = sum_columns+h2osoi_ice(c,j) + count_columns = count_columns+1 + end if + end do + if (count_columns == 0) then + !clmstate_tmp_local(g,j) = spval + else + clmstate_tmp_local(g,j) = sum_columns/count_columns + end if + if (clmstate_tmp_local(g,j) /= spval) then + tws_inc(g) = tws_inc(g) + clmstate_tmp_local(g,j) + end if + end do + end do + do jn = 1, nlevsoi + if (masterproc) then + call mpi_gatherv (clmstate_tmp_local(beg,jn), numsend, MPI_REAL8, & + clmstate_tmp_global, numrecvv, displsv, & + MPI_REAL8, 0, mpicom, ier) + else + call mpi_gatherv (clmstate_tmp_local(beg,jn), numsend , MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier) + end if + if(masterproc) then + do g1 = 1, numg + ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1 + jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1 + clmstate_out(ji,jj,jn) = clmstate_tmp_global(g1) + end do + end if + end do + if (masterproc) then + status = nf90_inq_varid(il_file_id, "SOILICE" , ncvarid(2)) + status = nf90_put_var( il_file_id, ncvarid(2), clmstate_out(:,:,:), & + start = [ 1, 1, 1 ], count = [ ndlon, ndlat, nlevsoi ] ) + end if + + + + + clmstate_tmp_local(begg:endg,:) = 0._r8 + do g=begg,endg + count_columns = 0 + sum_columns = 0 + do c=begc,endc + if (g==col%gridcell(c).and.col%hydrologically_active(c)) then + sum_columns = sum_columns+h2osno(c) + count_columns = count_columns+1 + end if + end do + if (count_columns == 0) then + !clmstate_tmp_local(g,1) = spval + else + clmstate_tmp_local(g,1) = sum_columns/count_columns + tws_inc(g) = tws_inc(g) + clmstate_tmp_local(g,1) + end if + end do + + + if (masterproc) then + call mpi_gatherv (clmstate_tmp_local(beg,1), numsend, MPI_REAL8, & + clmstate_tmp_global, numrecvv, displsv, & + MPI_REAL8, 0, mpicom, ier) + else + call mpi_gatherv (clmstate_tmp_local(beg,1), numsend , MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier) + end if + if(masterproc) then + do g1 = 1, numg + ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1 + jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1 + clmstate_out(ji,jj,1) = clmstate_tmp_global(g1) + end do + + end if + if (masterproc) then + status = nf90_inq_varid(il_file_id, "H2OSNO" , ncvarid(3)) + status = nf90_put_var( il_file_id, ncvarid(3), clmstate_out(:,:,1), & + start = [ 1, 1 ], count = [ ndlon, ndlat ] ) + end if + + + + if (masterproc) then + call mpi_gatherv (tws_inc(beg), numsend , MPI_REAL8, clmstate_tmp_global, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier) + else + call mpi_gatherv (tws_inc(beg), numsend , MPI_REAL8, 0._r8, numrecvv, displsv, MPI_REAL8, 0, mpicom, ier) + end if + if(masterproc) then + do g1 = 1, numg + ji = mod(ldecomp%gdc2glo(g1)-1,ldomain%ni) + 1 + jj = (ldecomp%gdc2glo(g1) - 1)/ldomain%ni + 1 + clmstate_out(ji,jj,1) = clmstate_tmp_global(g1) + end do + end if + if (masterproc) then + status = nf90_inq_varid(il_file_id, "TWS" , ncvarid(4)) + status = nf90_put_var( il_file_id, ncvarid(4), clmstate_out(:,:,1), & + start = [ 1, 1 ], count = [ ndlon, ndlat ] ) + end if + + + + if(masterproc) then + status = nf90_close(il_file_id) + deallocate(clmstate_out) + deallocate(clmstate_tmp_global) + end if + deallocate(tws_inc) + deallocate(clmstate_tmp_local) + + + + + +end subroutine print_inc_clm #endif subroutine get_update_filename (iofile) @@ -243,3 +536,22 @@ subroutine get_update_filename (iofile) iofile = trim(caseid)//".update."//trim(cdate)//".nc" !iofile = trim(caseid)//".update.nc" end subroutine get_update_filename + +subroutine get_inc_filename (iofile) + use clm_varctl, only : caseid, inst_suffix + use clm_time_manager, only : get_curr_date, get_prev_date + ! !ARGUMENTS: + implicit none + character(len=300),intent(inout) :: iofile + ! LOCAL VARIABLES: + character(len=256) :: cdate !date char string + integer :: day !day (1 -> 31) + integer :: mon !month (1 -> 12) + integer :: yr !year (0 -> ...) + integer :: sec !seconds into current day + !----------------------------------------------------------------------- + call get_prev_date (yr, mon, day, sec) + write(cdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr,mon,day,sec + iofile = trim(caseid)//".inc"//trim(inst_suffix)//"."//trim(cdate)//".nc" + +end subroutine get_inc_filename diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 56687e032..f3f7edd67 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -198,11 +198,17 @@ void integrate_tsmp() { void update_tsmp(){ #if defined CLMSA - if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ + if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 0))){ update_clm(&tstartcycle, &mype_world); - if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ - print_update_clm(&tcycle, &total_steps); + /* if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ */ + /* print_update_clm(&tcycle, &total_steps); */ + /* } */ + /* TODO: enkfpf.par input switch "CLM:print_inc" */ +#ifdef CLMFIVE + if ((clmupdate_tws != 0) || (clmupdate_swc != 0)){ + print_inc_clm(); } +#endif } #endif diff --git a/src/Makefile b/src/Makefile index 76d481a52..d5ab4fa9b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -124,7 +124,8 @@ OBJ_PDAF_GEN = PDAF_analysis_utils.o \ PDAFlocalomi_put_state_nondiagR.o \ PDAFlocalomi_put_state_nondiagR_si.o \ PDAFlocalomi_put_state_si.o \ - PDAF_correlation_function.o + PDAF_correlation_function.o \ + PDAF_reset_dim_p.o # Specific PDAF-routines for SEIK diff --git a/src/PDAF_interfaces_module.F90 b/src/PDAF_interfaces_module.F90 index 68a685966..a652a7c99 100644 --- a/src/PDAF_interfaces_module.F90 +++ b/src/PDAF_interfaces_module.F90 @@ -923,6 +923,13 @@ SUBROUTINE PDAF_get_ensstats(skew_ptr, kurt_ptr, status) END SUBROUTINE PDAF_get_ensstats END INTERFACE + INTERFACE + SUBROUTINE PDAF_reset_dim_p(dim_p_in, outflag) + INTEGER, INTENT(in) :: dim_p_in ! Sub-type of filter + INTEGER, INTENT(inout):: outflag ! Status flag + END SUBROUTINE PDAF_reset_dim_p + END INTERFACE + INTERFACE SUBROUTINE PDAF_reset_forget(forget_in) REAL, INTENT(in) :: forget_in ! New value of forgetting factor diff --git a/src/PDAFomi_obs_f.F90 b/src/PDAFomi_obs_f.F90 index a86eba2ab..1c86c9fbc 100644 --- a/src/PDAFomi_obs_f.F90 +++ b/src/PDAFomi_obs_f.F90 @@ -110,6 +110,7 @@ MODULE PDAFomi_obs_f INTEGER :: off_obs_f !< Offset of this observation in overall full obs. vector INTEGER :: off_obs_g !< Offset of this observation in overall global obs. vector INTEGER :: obsid !< Index of observation over all assimilated observations + INTEGER :: infile !< Yorck addition for mutliple observations REAL, ALLOCATABLE :: obs_f(:) !< Full observed field REAL, ALLOCATABLE :: ocoord_f(:,:) !< Coordinates of full observation vector REAL, ALLOCATABLE :: ivar_obs_f(:) !< Inverse variance of full observations