Skip to content

fortitude fixes for OMI code #46

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jjokella wants to merge 1 commit into
base: tsmp-pdaf-patched-omi
Choose a base branch
from
Open

fortitude fixes for OMI code #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 46 additions & 3 deletions interface/framework/obs_GRACE_pdafomi.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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


Expand All @@ -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

Expand Down Expand Up @@ -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


Expand All @@ -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
Expand Down Expand Up @@ -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


Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down
83 changes: 64 additions & 19 deletions interface/framework/obs_SM_pdafomi.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion interface/framework/prepoststep_ens_pdaf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
40 changes: 37 additions & 3 deletions interface/model/eclm/enkf_clm_mod_5.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1113,6 +1129,10 @@ subroutine set_clm_statevec_tws()

end do

case default

error stop "Unsupported state_setup"

end select


Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down