Skip to content
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.

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

modified EnKF update_type 10 for PEATCLSM (include catdef in EnKF State) #529

Merged
merged 3 commits into from
Feb 26, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
! update_type = 7: 3d Tskin/ght1 update; Tskin obs
! update_type = 8: 3d soil moisture/Tskin/ght(1); TB obs
! update_type = 9: 1d Tskin/ght1 update; FT obs
! update_type = 10: 3d soil moisture excl catdef/Tskin/ght(1); TB obs
! update_type = 10: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile; TB obs

update_type = 0

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -435,6 +435,11 @@ subroutine check_cat_bias_inputs( update_type, cat_bias_param )
end if

case (10)

! For PEATCLSM tiles, update_type 10 may have non-zero catdef increments,
! but catdef increments still vanish for non-PEATCLSM tiles.
! Leaving code here unchanged for now, that is, cannot request bias
! corr for catdef. - reichle, 20 Feb 2022

if ( &
!(cat_bias_param%Nparam%tc1>0) .or. &
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1340,25 +1340,27 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, &

case (6,8,9,10) select_update_type ! soil moisture and temperature update

! for update_type 10, catdef increments may be zero by design

if (logit) write (logunit,*) &
'apply_enkf_increments(): applying soil moisture and Tskin/ght1 increments'

do n=1,N_catd
do n_e=1,N_ens

cat_progn(n,n_e)%srfexc = &
cat_progn(n,n_e)%srfexc + cat_progn_incr(n,n_e)%srfexc
cat_progn(n,n_e)%rzexc = &
cat_progn(n,n_e)%rzexc + cat_progn_incr(n,n_e)%rzexc
cat_progn(n,n_e)%catdef = &
cat_progn(n,n_e)%catdef + cat_progn_incr(n,n_e)%catdef

cat_progn(n,n_e)%tc1 = &
cat_progn(n,n_e)%tc1 + cat_progn_incr(n,n_e)%tc1
cat_progn(n,n_e)%tc1 + cat_progn_incr(n,n_e)%tc1
cat_progn(n,n_e)%tc2 = &
cat_progn(n,n_e)%tc2 + cat_progn_incr(n,n_e)%tc2
cat_progn(n,n_e)%tc2 + cat_progn_incr(n,n_e)%tc2
cat_progn(n,n_e)%tc4 = &
cat_progn(n,n_e)%tc4 + cat_progn_incr(n,n_e)%tc4
cat_progn(n,n_e)%tc4 + cat_progn_incr(n,n_e)%tc4

cat_progn(n,n_e)%ght(1) = &
cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,8 @@ module clsm_ensupd_upd_routines

use catch_constants, ONLY: &
N_snow => CATCH_N_SNOW, &
N_gt => CATCH_N_GT
N_gt => CATCH_N_GT, &
PEATCLSM_POROS_THRESHOLD

use StieglitzSnow, ONLY: &
StieglitzSnow_calc_asnow
Expand Down Expand Up @@ -3389,6 +3390,7 @@ subroutine cat_enkf_increments( &
! reichle, 27 Jul 2005
! reichle, 18 Oct 2005 - return increments (instead of updated cat_progn)
! reichle, 17 Oct 2011 - added "l2f" for revised (MPI) analysis
! reichle, 20 Feb 2022 - modified update_type 10 for PEATCLSM
!
! --------------------------------------------------------------

Expand Down Expand Up @@ -3449,10 +3451,10 @@ subroutine cat_enkf_increments( &
real, parameter :: SWE_threshold = +HUGE(1.) ! = 1.e-4 ! [kg/m2]

real, parameter :: tp1_threshold = -HUGE(1.) ! = 0.2 ! [CELSIUS]

integer :: n, n_e, kk, ii

integer :: N_state, N_selected_obs, N_select_varnames, N_select_species
integer :: N_state_max, N_state, N_selected_obs, N_select_varnames, N_select_species

real :: halo_minlon, halo_maxlon, halo_minlat, halo_maxlat
real :: tmp_minlon, tmp_maxlon, tmp_minlat, tmp_maxlat
Expand All @@ -3461,7 +3463,7 @@ subroutine cat_enkf_increments( &

real :: fice_minus, tp1_minus, ght1_minus
real :: fice_plus, tp1_plus, ght1_plus

integer, dimension(N_obs) :: ind_obs

real, allocatable, dimension(:,:) :: State_incr
Expand Down Expand Up @@ -4094,16 +4096,21 @@ subroutine cat_enkf_increments( &

case (8,10) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb obs

! update each tile separately using all observations within
! the customized halo around each tile

! update each tile separately using all observations within customized halo around each tile
!
! state vector includes different subsets of Catchment model soil moisture prognostics:
!
! update_type = 8 -- state vector: srfexc, rzexc, catdef, tc1, tc2, tc4, ght1
! update_type = 10 -- state vector: srfexc, rzexc, tc1, tc2, tc4, ght1
! update_type | subset of tiles | state vector
! ===================================================================================================
! 8 | all | srfexc, rzexc, catdef, tc1, tc2, tc4, ght1
! ---------------------------------------------------------------------------------------------------
! 10 | PEATCLSM tiles | srfexc, rzexc, catdef, tc1, tc2, tc4, ght1
! | otherwise | srfexc, rzexc, tc1, tc2, tc4, ght1 (incl. NLv4 peat tiles)
! ---------------------------------------------------------------------------------------------------
!
! reichle, 27 Nov 2017

! reichle, 20 Feb 2022 - modified for PEATCLSM

if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb obs'

N_select_varnames = 1
Expand All @@ -4113,28 +4120,20 @@ subroutine cat_enkf_increments( &
call get_select_species( &
N_select_varnames, select_varnames(1:N_select_varnames), &
N_obs_param, obs_param, N_select_species, select_species )

if (update_type==8) then

N_state = 7

else

N_state = 6

end if

allocate( State_incr(N_state,N_ens))
allocate( State_lon( N_state ))
allocate( State_lat( N_state ))
N_state_max = 7

allocate( State_incr(N_state_max,N_ens))
allocate( State_lon( N_state_max ))
allocate( State_lat( N_state_max ))

do kk=1,N_catd

! compute increments only snow-free and non-frozen tiles

if ( (SWE_ensavg(kk) < SWE_threshold) .and. &
(tp1_ensavg(kk) > tp1_threshold) ) then

! find observations within halo around tile kk

halo_minlon = tile_coord(kk)%com_lon - xcompact
Expand All @@ -4158,20 +4157,32 @@ subroutine cat_enkf_increments( &

if (N_selected_obs>0) then

if ( (update_type== 8 ) .or. &
(update_type==10 .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) &
) then

N_state = 7 ! srfexc, rzexc, catdef, tc1, tc2, tc4, ght1

else

N_state = 6 ! srfexc, rzexc, tc1, tc2, tc4, ght1

end if

! assemble State_minus
! (on input, cat_progn contains cat_progn_minus)

if (update_type==8) then
if ( N_state==7 ) then

State_incr(1,:) = cat_progn( kk,:)%srfexc/scale_srfexc
State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc
State_incr(3,:) = cat_progn( kk,:)%catdef/scale_catdef ! catdef in State

State_incr(4,:) = cat_progn( kk,:)%tc1 /scale_temp
State_incr(5,:) = cat_progn( kk,:)%tc2 /scale_temp
State_incr(6,:) = cat_progn( kk,:)%tc4 /scale_temp
State_incr(7,:) = cat_progn( kk,:)%ght(1)/scale_ght1

State_incr(1,:) = cat_progn( kk,:)%srfexc/scale_srfexc
State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc
State_incr(3,:) = cat_progn( kk,:)%catdef/scale_catdef

State_incr(4,:) = cat_progn( kk,:)%tc1 /scale_temp
State_incr(5,:) = cat_progn( kk,:)%tc2 /scale_temp
State_incr(6,:) = cat_progn( kk,:)%tc4 /scale_temp
State_incr(7,:) = cat_progn( kk,:)%ght(1)/scale_ght1

else

State_incr(1,:) = cat_progn( kk,:)%srfexc/scale_srfexc
Expand Down Expand Up @@ -4200,23 +4211,26 @@ subroutine cat_enkf_increments( &
Obs_pred(ind_obs(1:N_selected_obs),:), &
Obs_pert(ind_obs(1:N_selected_obs),:), &
Obs_cov, &
State_incr, State_lon, State_lat, xcompact, ycompact, &
State_incr(1:N_state,:), &
State_lon( 1:N_state ), &
State_lat( 1:N_state ), &
xcompact, ycompact, &
fcsterr_inflation_fac )

deallocate(Obs_cov)

! assemble cat_progn increments

if (update_type==8) then
if (N_state==7) then

cat_progn_incr(kk,:)%srfexc = State_incr(1,:)*scale_srfexc
cat_progn_incr(kk,:)%rzexc = State_incr(2,:)*scale_rzexc
cat_progn_incr(kk,:)%catdef = State_incr(3,:)*scale_catdef

cat_progn_incr(kk,:)%tc1 = State_incr(4,:)*scale_temp
cat_progn_incr(kk,:)%tc2 = State_incr(5,:)*scale_temp
cat_progn_incr(kk,:)%tc4 = State_incr(6,:)*scale_temp
cat_progn_incr(kk,:)%ght(1) = State_incr(7,:)*scale_ght1
cat_progn_incr(kk,:)%srfexc = State_incr(1,:)*scale_srfexc
cat_progn_incr(kk,:)%rzexc = State_incr(2,:)*scale_rzexc
cat_progn_incr(kk,:)%catdef = State_incr(3,:)*scale_catdef ! catdef in State
cat_progn_incr(kk,:)%tc1 = State_incr(4,:)*scale_temp
cat_progn_incr(kk,:)%tc2 = State_incr(5,:)*scale_temp
cat_progn_incr(kk,:)%tc4 = State_incr(6,:)*scale_temp
cat_progn_incr(kk,:)%ght(1) = State_incr(7,:)*scale_ght1

else

Expand All @@ -4227,9 +4241,9 @@ subroutine cat_enkf_increments( &
cat_progn_incr(kk,:)%tc2 = State_incr(4,:)*scale_temp
cat_progn_incr(kk,:)%tc4 = State_incr(5,:)*scale_temp
cat_progn_incr(kk,:)%ght(1) = State_incr(6,:)*scale_ght1

end if

end if

end if ! thresholds
Expand Down