diff --git a/interface/framework/obs_GRACE_pdafomi.F90 b/interface/framework/obs_GRACE_pdafomi.F90 index d7930f5b..149a4581 100644 --- a/interface/framework/obs_GRACE_pdafomi.F90 +++ b/interface/framework/obs_GRACE_pdafomi.F90 @@ -217,7 +217,8 @@ SUBROUTINE init_dim_obs_GRACE(step, dim_obs) 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 + 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 @@ -565,7 +566,8 @@ SUBROUTINE init_dim_obs_GRACE(step, dim_obs) 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 (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)) @@ -680,6 +682,8 @@ SUBROUTINE obs_op_GRACE(dim_p, dim_obs, state_p, ostate) use mod_parallel_pdaf, & only: comm_filter + use mod_parallel_pdaf, only: abort_parallel + use mod_parallel_pdaf, only: mype_world use clm_varpar , only : nlevsoi @@ -805,6 +809,12 @@ SUBROUTINE obs_op_GRACE(dim_p, dim_obs, state_p, ostate) tws_from_statevector(g) = tws_from_statevector(g) + state_p(count + clm_varsize_tws(1) + clm_varsize_tws(2)) end do + case default + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR", & + " invalid `state_setup` in obs_GRACE_pdafomi.F90." + call abort_parallel() + end select @@ -986,6 +996,8 @@ subroutine add_obs_err_GRACE(step, dim_obs, C) use mod_read_obs, only: multierr USE mod_parallel_pdaf, & ONLY: npes_filter + USE mod_parallel_pdaf, ONLY: abort_parallel + USE mod_parallel_pdaf, ONLY: mype_world use PDAFomi, only: obsdims @@ -1026,6 +1038,12 @@ subroutine add_obs_err_GRACE(step, dim_obs, C) C(i,j) = C(i,j) + obscov(obs_pdaf2nc(i),obs_pdaf2nc(j)) end do end do + case default + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR", & + " invalid `multierr` in obs_GRACE_pdafomi.F90." + call abort_parallel() + end select @@ -1040,6 +1058,8 @@ subroutine init_obscovar_GRACE(step, dim_obs, dim_obs_p, covar, m_state_p, isdia USE mod_parallel_pdaf, & ONLY: npes_filter + USE mod_parallel_pdaf, ONLY: abort_parallel + USE mod_parallel_pdaf, ONLY: mype_world use PDAFomi, only: obsdims, map_obs_id @@ -1105,6 +1125,12 @@ subroutine init_obscovar_GRACE(step, dim_obs, dim_obs_p, covar, m_state_p, isdia isdiag = .FALSE. + case default + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR", & + " invalid `multierr` in obs_GRACE_pdafomi.F90." + call abort_parallel() + end select @@ -1119,6 +1145,8 @@ 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 + USE mod_parallel_pdaf, ONLY: abort_parallel + USE mod_parallel_pdaf, ONLY: mype_world INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector @@ -1149,6 +1177,12 @@ subroutine prodRinvA_GRACE(step, dim_obs_p, rank, obs_p, A_p, C_p) 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,:)) + case default + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR", & + " invalid `multierr` in obs_GRACE_pdafomi.F90." + call abort_parallel() + end select @@ -1161,12 +1195,15 @@ subroutine prodRinvA_l_GRACE(domain_p, step, dim_obs, rank, obs_l, A_l, C_l) use mod_assimilation, only: obscov, obs_pdaf2nc, cradius_GRACE, locweight use mod_read_obs, only: multierr use PDAFomi, only: PDAFomi_observation_localization_weights + USE mod_parallel_pdaf, ONLY: abort_parallel + USE mod_parallel_pdaf, ONLY: mype_world 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) :: 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 @@ -1284,6 +1321,12 @@ subroutine prodRinvA_l_GRACE(domain_p, step, dim_obs, rank, obs_l, A_l, C_l) C_l(off+1:off+thisobs_l%dim_obs_l,:) = matmul(obscov_inv_l,A_l(off+1:off+thisobs_l%dim_obs_l,:)) + case default + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR", & + " invalid `multierr` in obs_GRACE_pdafomi.F90." + call abort_parallel() + end select deallocate(weight) diff --git a/interface/framework/obs_SM_pdafomi.F90 b/interface/framework/obs_SM_pdafomi.F90 index 8c1e25a7..be0d8abe 100644 --- a/interface/framework/obs_SM_pdafomi.F90 +++ b/interface/framework/obs_SM_pdafomi.F90 @@ -442,15 +442,42 @@ SUBROUTINE init_dim_obs_SM(step, dim_obs) ! Use index array for setting the correct state vector index in `obs_id_p` - ! id_obs_p has to be allocated but it is not really used when a self implemented observation operator is used. Else, the structure has to point for each observation - ! on the state vector element that is used to predict that observation. The dimension is then not the number of gridcells but dim_obs - ! Normally, this does not yield any error. But when the number of observations is larger than the process-local state vector size, the model will crash. - ! This is caused by an internal check in the PDAFomi_obs_f script: IF (MAXVAL(thisobs%id_obs_p) > dim_p .AND. dim_obs_p>0) - ! Here, I will just set the id_obs_p to 1, not to i to not cause this error for large observation dimensions - ! Note that id_obs_p could be used correctly if an internal observation operator is used but in this case, it has to be allocated with dim_obs_p, so after this loop - ! This could then be filled with the state vector element index used to predict each observation, for now I set this everywhere to 1. - ! Note that I keep the initial implementation for GRACE as I use the obs_id_p there for the observation operator. But the number of observations is there usually pretty - ! low, so this should not make a problem + ! id_obs_p has to be allocated but it is not really + ! used when a self implemented observation operator + ! is used. Else, the structure has to point for each + ! observation on the state vector element that is + ! used to predict that observation. + ! + ! The dimension is then not the number of gridcells + ! but dim_obs + ! + ! Normally, this does not yield any error. But when + ! the number of observations is larger than the + ! process-local state vector size, the model will + ! crash. + ! + ! This is caused by an internal check in the + ! PDAFomi_obs_f script: IF (MAXVAL(thisobs%id_obs_p) + ! > dim_p .AND. dim_obs_p>0) + ! + ! Here, I will just set the id_obs_p to 1, not to i + ! to not cause this error for large observation + ! dimensions + ! + ! Note that id_obs_p could be used correctly if an + ! internal observation operator is used but in this + ! case, it has to be allocated with dim_obs_p, so + ! after this loop + ! + ! This could then be filled with the state vector + ! element index used to predict each observation, + ! for now I set this everywhere to 1. + ! + ! Note that I keep the initial implementation for + ! GRACE as I use the obs_id_p there for the + ! observation operator. But the number of + ! observations is there usually pretty low, so this + ! should not make a problem thisobs%id_obs_p(1,state_clm2pdaf_p(c,layer_obs(i))) = 1 if (obs_snapped) then @@ -1134,11 +1161,20 @@ subroutine add_obs_err_SM(step, dim_obs, C) 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! + ! Initialize indices --> we only have information about local + ! obs. dims per PE, so we get the global indices, more + ! generalizable than using the afrrays 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 @@ -1185,7 +1221,9 @@ subroutine init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) 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 + ! 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)) @@ -1196,7 +1234,9 @@ subroutine init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) 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! + ! 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 @@ -1209,8 +1249,12 @@ subroutine init_obscovar_SM(step, dim_obs, dim_obs_p, covar, m_state_p, isdiag) 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 + ! the inverse of the observation variance is saved for + ! each observation, so we do not need any other + covar(i, i) = covar(i, i) + 1.0/thisobs%ivar_obs_f(cnt) + ! 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 @@ -1260,7 +1304,8 @@ subroutine prodRinvA_l_SM(domain_p, step, dim_obs, rank, obs_l, A_l, C_l) 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) :: 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 diff --git a/interface/framework/prepoststep_ens_pdaf.F90 b/interface/framework/prepoststep_ens_pdaf.F90 index a1188fc8..7ae7fac0 100644 --- a/interface/framework/prepoststep_ens_pdaf.F90 +++ b/interface/framework/prepoststep_ens_pdaf.F90 @@ -110,7 +110,9 @@ SUBROUTINE prepoststep_ens_pdaf(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, & INTEGER, SAVE :: allocflag = 0 ! Flag for memory counting LOGICAL, SAVE :: firstio = .TRUE. ! File output is peformed for first time? LOGICAL, SAVE :: firsttime = .TRUE. ! Routine is called for first time? - LOGICAL, SAVE :: firsttime_omi = .TRUE. ! Routine is called for first time? --> for PDAF OMI, for backward compatibility if the statement in 1==2 should be used at one point + LOGICAL, SAVE :: firsttime_omi = .TRUE. ! Routine is called for first time? + ! --> for PDAF OMI, for backward compatibility + ! if the statement in 1==2 should be used at one point REAL :: invdim_ens ! Inverse ensemble size REAL :: invdim_ensm1 ! Inverse of ensemble size minus 1 REAL :: rmserror_est ! estimated RMS error diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index deb745f5..fb7a796a 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -212,9 +212,17 @@ subroutine define_clm_statevec(mype) 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). + ! 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 @@ -671,6 +679,10 @@ subroutine define_clm_statevec_tws(mype) clm_varsize_tws(4) = num_layer(1) clm_statevecsize = clm_statevecsize + num_layer(1) + case default + + error stop "Unsupported state_setup" + end select end subroutine define_clm_statevec_tws @@ -900,6 +912,10 @@ subroutine set_clm_statevec_tws() sfc => waterstate_inst%h2osfc_col_mean can => waterstate_inst%h2ocan_patch_mean TWS => waterstate_inst%tws_hactive_mean + case default + + error stop "Unsupported TWS_smoother" + end select select case (state_setup) @@ -1113,6 +1129,10 @@ subroutine set_clm_statevec_tws() end do + case default + + error stop "Unsupported state_setup" + end select @@ -1503,6 +1523,8 @@ subroutine clm_update_tws() 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() + case default + error stop "Unsupported state_setup" end select call finalize_increments() @@ -2469,6 +2491,10 @@ subroutine init_dim_l_clm(domain_p, dim_l) dim_l = dim_l+1 end if end do + case default + + error stop "Unsupported state_setup" + end select endif @@ -2585,6 +2611,10 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) end do end if + case default + + error stop "Unsupported state_setup" + end select end if NOGRACE @@ -2707,6 +2737,10 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) end if end do + case default + + error stop "Unsupported state_setup" + end select endif NOGRACE