Skip to content

Commit

Permalink
Add "update_type" for joint assimilation of ASCAT SM and SMAP Tb (#703)
Browse files Browse the repository at this point in the history
  • Loading branch information
@gmao-rreichle
gmao-rreichle committed Feb 27, 2024
2 parents 14ca604 + 175ccb5 commit b175515
Show file tree
Hide file tree
Showing 3 changed files with 253 additions and 21 deletions.
29 changes: 16 additions & 13 deletions src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,19 +19,22 @@
! update type - for details see subroutine cat_enkf_update()
! (note: all 3d updates use compact support)
!
! local = "1d", regional = "3d"
!
! update_type = 0: NO assimilation, NO bias correction
! update_type = 1: 1d soil moisture analysis; sfmc obs
! update_type = 2: 3d soil moisture analysis; sfmc obs
! update_type = 3: 1d Tskin (assim incr NOT applied, use w/ bias corr) analysis; Tskin obs
! update_type = 4: 1d Tskin/ght1 (assim incr applied, use w/ or w/o bias corr) analysis; Tskin obs
! update_type = 5: 1d Tskin/ght1 (assim incr NOT applied, use w/ bias corr) analysis; Tskin obs
! update_type = 6: 1d soil moisture/Tskin/ght(1); TB obs
! 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/Tskin/ght(1) excl. catdef unless PEATCLSM tile; TB obs
! local = "1d", regional = "3d"
!
! # = no longer supported
!
! update_type = 0: NO assimilation, NO bias correction
! # update_type = 1: 1d soil moisture analysis; sfmc obs
! # update_type = 2: 3d soil moisture analysis; sfmc obs
! update_type = 3: 1d Tskin (assim incr NOT applied, use w/ bias corr) analysis; Tskin obs
! update_type = 4: 1d Tskin/ght1 (assim incr applied, use w/ or w/o bias corr) analysis; Tskin obs
! update_type = 5: 1d Tskin/ght1 (assim incr NOT applied, use w/ bias corr) analysis; Tskin obs
! update_type = 6: 1d soil moisture/Tskin/ght(1); TB obs
! 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/Tskin/ght(1) excl. catdef unless PEATCLSM tile; TB obs
! update_type = 13: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile; sfmc and TB obs

update_type = 0

Expand Down
8 changes: 5 additions & 3 deletions src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1330,10 +1330,12 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, &

cat_progn_has_changed = .false.

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

! some of the increments fields below may be zero by design
! (e.g., tc[X]=ght(1)=0 in update_type=13 when only sfmc or sfds obs are assimilated;
! or catdef=0 in update_type 10 or 13 when tile has mineral soil)

! 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'

Expand Down
237 changes: 232 additions & 5 deletions src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3363,6 +3363,7 @@ subroutine cat_enkf_increments( &
! 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
! amfox, 6 Feb 2024 - added update type 13 for combination of ASCAT SM and SMAP Tb
!
! --------------------------------------------------------------

Expand Down Expand Up @@ -3422,9 +3423,9 @@ subroutine cat_enkf_increments( &

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

integer :: n, n_e, kk, ii
integer :: n, n_e, kk, ii, jj

integer :: N_state_max, 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, N_select_species_Tb

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

real, allocatable, dimension(:) :: State_lon, State_lat

integer, dimension(N_obs_param) :: select_species ! alloc max possible length
integer, dimension(N_obs_param) :: select_species, select_species_Tb ! alloc max possible length

character(40), dimension(N_obs_param) :: select_varnames ! alloc max possible length

Expand Down Expand Up @@ -3471,6 +3472,8 @@ subroutine cat_enkf_increments( &
real, dimension( N_catd) :: tp1_ensavg

type(obs_param_type) :: this_obs_param

logical :: found_Tb_obs

! -----------------------------------------------------------------------

Expand Down Expand Up @@ -3601,6 +3604,10 @@ subroutine cat_enkf_increments( &

if (logit) write (logunit,*) 'get 1d soil moisture increments; sfmc obs'

! disable update_type=1 (b/c it includes catdef in state vector for mineral soil)

call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'update_type=1 no longer supported; use update_type=13 instead')

N_select_varnames = 2

select_varnames(1) = 'sfmc'
Expand Down Expand Up @@ -3663,13 +3670,17 @@ subroutine cat_enkf_increments( &

end do

case (2) select_update_type ! 3d soil moisture analysis; sfmc and/or sfds obs
case (2) select_update_type ! 3d soil moisture analysis; sfmc+sfds obs

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

if (logit) write (logunit,*) 'get 3d soil moisture increments; sfmc obs'

! disable update_type=2 (b/c it includes catdef in state vector for mineral soil)

call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'update_type=2 no longer supported; use update_type=13 instead')

N_select_varnames = 2

select_varnames(1) = 'sfmc'
Expand Down Expand Up @@ -4403,6 +4414,222 @@ subroutine cat_enkf_increments( &
end do

! ----------------------------------

case (13) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs

! update each tile separately using all observations within customized halo around each tile
!
! state vector differs for each tile depending on assimilated obs and soil type
!
! obs | soil | N_state | state vector
! ----------------------------------------------------------------------
! sfcm/sfds only | mineral | 2 | srfexc, rzexc
! sfcm/sfds only | peat | 3 | srfexc, rzexc, catdef,
! sfcm/sfds & Tb | mineral | 6 | srfexc, rzexc, tc[x], ght(1)
! sfcm/sfds & Tb | peat | 7 | srfexc, rzexc, catdef, tc[x], ght(1)
!
! amfox+rreichle, 26 Feb 2024

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

N_select_varnames = 0

do ii = 1,N_obs_param
if (trim(obs_param(ii)%varname) == 'Tb') then
N_select_varnames = N_select_varnames + 1
select_varnames(N_select_varnames) = 'Tb'
exit
end if
end do

do ii = 1,N_obs_param
if (trim(obs_param(ii)%varname) == 'sfmc') then
N_select_varnames = N_select_varnames + 1
select_varnames(N_select_varnames) = 'sfmc'
exit
end if
end do

do ii = 1,N_obs_param
if (trim(obs_param(ii)%varname) == 'sfds') then
N_select_varnames = N_select_varnames + 1
select_varnames(N_select_varnames) = 'sfds'
exit
end if
end do

! Will get all species associated with Tb or sfds observations

call get_select_species( &
N_select_varnames, select_varnames(1:N_select_varnames), &
N_obs_param, obs_param, N_select_species, select_species )

! Determine which species are Tb

call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb )

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

N_state = 2 ! initialize (always have srfexc and rzexc in state vector)

! compute increments only for 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
halo_maxlon = tile_coord(kk)%com_lon + xcompact
halo_minlat = tile_coord(kk)%com_lat - ycompact
halo_maxlat = tile_coord(kk)%com_lat + ycompact

! simple approach to dateline issue (cut halo back to at most -180:180, -90:90)
! - reichle, 28 May 2013

halo_minlon = max(halo_minlon,-180.)
halo_maxlon = min(halo_maxlon, 180.)
halo_minlat = max(halo_minlat, -90.)
halo_maxlat = min(halo_maxlat, 90.)

call get_ind_obs_lat_lon_box( &
N_obs, Observations, &
halo_minlon, halo_maxlon, halo_minlat, halo_maxlat, &
N_select_species, select_species(1:N_select_species), &
N_selected_obs, ind_obs )

if (N_selected_obs>0) then

! Determine if Tb observations are present

found_Tb_obs = .false.

do ii = 1,N_select_species_Tb
do jj = 1,N_selected_obs
if (select_species_Tb(ii) == Observations(ind_obs(jj))%species) then
found_Tb_obs = .true.
exit
end if
end do
if (found_Tb_obs) exit
end do

! if Tb_obs are present, add tc[X] and ght(1) to state vector

if (found_Tb_obs) N_state = N_state + 4

! for peatland tile, add catdef to state vector

if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) N_state = N_state + 1

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

if ( N_state==2 ) then

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

elseif ( N_state==3 ) 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

elseif ( N_state==6 ) then

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

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
State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc
State_incr(3,:) = cat_progn( kk,:)%catdef/scale_catdef ! catdef in State

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

end if

State_lon( :) = tile_coord(kk )%com_lon
State_lat( :) = tile_coord(kk )%com_lat

allocate(Obs_cov(N_selected_obs,N_selected_obs))

call assemble_obs_cov( N_selected_obs, N_obs_param, obs_param, &
Observations(ind_obs(1:N_selected_obs)), Obs_cov )

call enkf_increments( &
N_state, N_selected_obs, N_ens, &
Observations(ind_obs(1:N_selected_obs)), &
Obs_pred(ind_obs(1:N_selected_obs),:), &
Obs_pert(ind_obs(1:N_selected_obs),:), &
Obs_cov, &
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 ( N_state==2 ) then

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

elseif ( N_state==3 ) 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 ! catdef in State

elseif ( N_state==6 ) 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,:)%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

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(3,:)*scale_temp
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

end do

! ----------------------------------

case default

Expand Down Expand Up @@ -4974,7 +5201,7 @@ subroutine check_compact_support( &
Iam // '(): reset for 1d update_type: ycompact = ', ycompact
if (logit) write (logunit,*)

case (2,7,8,10) ! "3d" updates, check consistency of xcompact, ycompact
case (2,7,8,10,13) ! "3d" updates, check consistency of xcompact, ycompact

! check xcompact/ycompact against corr scales of model error

Expand Down

0 comments on commit b175515

Please sign in to comment.