From 437320093d68763644a4672525684fa736091ce4 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 12 Jul 2023 15:51:10 -0600 Subject: [PATCH 01/20] the intention --- src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml index cd3dbeca..210f321e 100644 --- a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -32,6 +32,7 @@ ! 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 = 11: Multivariate combination of 2 and 10; sfmc and TB obs update_type = 0 From a6e3341267d1438241363514de4aae2c04aaacbb Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 27 Jul 2023 11:13:04 -0600 Subject: [PATCH 02/20] new update_type --- .../clsm_ensupd_upd_routines.F90 | 180 +++++++++++++++++- 1 file changed, 179 insertions(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 603c5768..af132f0a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3426,7 +3426,7 @@ subroutine cat_enkf_increments( & integer :: n, n_e, kk, ii - 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_varnames real :: halo_minlon, halo_maxlon, halo_minlat, halo_maxlat real :: tmp_minlon, tmp_maxlon, tmp_minlat, tmp_maxlat @@ -3496,6 +3496,7 @@ subroutine cat_enkf_increments( & N_select_varnames = 0 N_select_species = 0 + N_varnames = 0 select_varnames = '' select_species = -8888 ! intentionally differs from init in get_select_species() @@ -4405,6 +4406,183 @@ subroutine cat_enkf_increments( & end do ! ---------------------------------- + + case (13) select_update_type ! All observation types from obs_param + + ! update each tile separately using all observations within customized halo around each tile + ! + ! amfox, 27 July 2023 + + if (logit) write (logunit,*) 'Get increments for all observation types in obs_param' + + N_select_varnames = 1 + + if (any(obs_param%varname == 'Tb')) then + N_varnames = N_varnames + 1 + select_varnames(N_varnames) = 'Tb' + end if + + if (any(obs_param%varname == 'sfds')) then + N_varnames = N_varnames + 1 + select_varnames(N_varnames) = 'sfds' + end if + + do ii = 1,N_varnames + + call get_select_species( & + N_select_varnames, select_varnames(ii), & + N_obs_param, obs_param, N_select_species, select_species ) + + if (select_varnames(ii)=='Tb') then + N_state_max = 7 + elseif (select_varnames(ii)=='sfds') then + N_state_max = 3 + N_state = 3 + else + err_msg = 'Dont know what state to use with this observation type.' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + 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 + 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 + + if (N_state_max==7 .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + + N_state = 7 ! srfexc, rzexc, catdef, tc1, tc2, tc4, ght1 + + elseif (N_state_max==7 .and. cat_param(kk)%poros= Date: Thu, 27 Jul 2023 13:32:33 -0400 Subject: [PATCH 03/20] fix if then --- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index af132f0a..6f0e4796 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4477,11 +4477,11 @@ subroutine cat_enkf_increments( & if (N_selected_obs>0) then - if (N_state_max==7 .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + if ((N_state_max==7 .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD)) then N_state = 7 ! srfexc, rzexc, catdef, tc1, tc2, tc4, ght1 - elseif (N_state_max==7 .and. cat_param(kk)%poros= Date: Thu, 27 Jul 2023 11:56:24 -0600 Subject: [PATCH 04/20] apply ASCAT _OR_ SMAP increment --- .../clsm_ensupd_enkf_update.F90 | 5 +- .../clsm_ensupd_upd_routines.F90 | 52 +++++++------------ 2 files changed, 23 insertions(+), 34 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 86d6449d..b49fc5e3 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1337,9 +1337,10 @@ 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 - ! for update_type 10, catdef increments may be zero by design + ! for update_type 10, catdef increments may be zero by design + ! for updagte_type 13, could be multiple zero increments if (logit) write (logunit,*) & 'apply_enkf_increments(): applying soil moisture and Tskin/ght1 increments' diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 6f0e4796..2f668050 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4416,6 +4416,11 @@ subroutine cat_enkf_increments( & if (logit) write (logunit,*) 'Get increments for all observation types in obs_param' N_select_varnames = 1 + N_state_max = 7 ! Needs to be constant size for applying increment, potential for lots of zeros + + allocate( State_incr(N_state_max,N_ens)) + allocate( State_lon( N_state_max )) + allocate( State_lat( N_state_max )) if (any(obs_param%varname == 'Tb')) then N_varnames = N_varnames + 1 @@ -4432,21 +4437,7 @@ subroutine cat_enkf_increments( & call get_select_species( & N_select_varnames, select_varnames(ii), & N_obs_param, obs_param, N_select_species, select_species ) - - if (select_varnames(ii)=='Tb') then - N_state_max = 7 - elseif (select_varnames(ii)=='sfds') then - N_state_max = 3 - N_state = 3 - else - err_msg = 'Dont know what state to use with this observation type.' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - - 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 @@ -4476,17 +4467,18 @@ subroutine cat_enkf_increments( & N_selected_obs, ind_obs ) if (N_selected_obs>0) then - - if ((N_state_max==7 .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD)) then - - N_state = 7 ! srfexc, rzexc, catdef, tc1, tc2, tc4, ght1 - - elseif ((N_state_max==7 .and. cat_param(kk)%poros=PEATCLSM_POROS_THRESHOLD) then + N_state = 7 + elseif (select_varnames(ii)=='Tb' .and. cat_param(kk)%poros Date: Fri, 11 Aug 2023 13:43:28 -0400 Subject: [PATCH 05/20] fix check_compact_support --- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 2f668050..d480a48b 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -5142,7 +5142,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 From 6e079f87060592c02749725681986a1543e90b52 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 17 Aug 2023 18:50:11 -0400 Subject: [PATCH 06/20] test read decode_ASCAT_ssom.x --- src/Applications/LDAS_App/CMakeLists.txt | 5 ++ .../LDAS_App/decode_ASCAT_ssom.F90 | 47 +++++++++++++++++++ 2 files changed, 52 insertions(+) create mode 100644 src/Applications/LDAS_App/decode_ASCAT_ssom.F90 diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index f2f5ed36..99d5028e 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -18,6 +18,11 @@ ecbuild_add_executable ( TARGET mwrtm_bin2nc4.x SOURCES util/inputs/mwRTM_params/mwrtm_bin2nc4.F90 LIBS GEOSlandassim_GridComp) + +ecbuild_add_executable ( + TARGET decode_ASCAT_ssom.x + SOURCES decode_ASCAT_ssom.F90 + LIBS NCEP_bufr_r4i4) set (scripts ldas_setup diff --git a/src/Applications/LDAS_App/decode_ASCAT_ssom.F90 b/src/Applications/LDAS_App/decode_ASCAT_ssom.F90 new file mode 100644 index 00000000..64d43ad5 --- /dev/null +++ b/src/Applications/LDAS_App/decode_ASCAT_ssom.F90 @@ -0,0 +1,47 @@ +program decode_ASCAT_ssom + + implicit none + + real*8, dimension(15) :: tmp_vdata + real*8, dimension(:,:), allocatable :: tmp_data + + integer, parameter :: lnbufr = 50 + integer, parameter :: max_obs = 250000 + integer :: idate,iret + integer :: ireadmg,ireadsb + integer :: N_obs + + character(8) :: subset + character(300) :: fname, mastertable_path + +! ------------------------------------------------------------------------- + + fname = '/home/amfox/smap/SMAP_Nature/ASCAT_EUMETSAT/Metop_C/Y2023/M03/' // & + 'M03-ASCA-ASCSMO02-NA-5.0-20230301090900.000000000Z-20230301105557-4839070.bfr' + mastertable_path = '/home/amfox/smap/SMAP_Nature/ASCAT_EUMETSAT' + +! Allocate the tmp_data array + allocate(tmp_data(max_obs, 15)) + + open(lnbufr, file=trim(fname), action='read',form='unformatted') + + call openbf(lnbufr,'SEC3', lnbufr) + call mtinfo( trim(mastertable_path) // '/BUFR_mastertable/', 51, 52) + call datelen(10) + + N_obs = 0 + msg_report: do while(ireadmg(lnbufr,subset,idate) ==0) + loop_report: do while(ireadsb(lnbufr) == 0) + call ufbint(lnbufr,tmp_vdata,15,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM DOMO SMPF SMCF ALFR TPCX IWFR CLATH CLONH') + N_obs = N_obs + 1 + tmp_data(N_obs,:) = tmp_vdata + end do loop_report + end do msg_report + + write(*,*) 'N_obs = ', N_obs + write(*,*) tmp_vdata + + call closbf(lnbufr) + close(lnbufr) + +end program decode_ASCAT_ssom From ef970abb705e008ef636f4296f155432d86dd8a7 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 13 Sep 2023 17:24:22 -0600 Subject: [PATCH 07/20] added State_incr_cum --- .../clsm_ensupd_upd_routines.F90 | 104 ++++++++++-------- 1 file changed, 59 insertions(+), 45 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index d480a48b..f3a27971 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3438,7 +3438,7 @@ subroutine cat_enkf_increments( & integer, dimension(N_obs) :: ind_obs - real, allocatable, dimension(:,:) :: State_incr + real, allocatable, dimension(:,:) :: State_incr, State_incr_cum real, allocatable, dimension(:,:) :: Obs_cov ! measurement error covariance real, allocatable, dimension(:) :: State_lon, State_lat @@ -3474,6 +3474,8 @@ subroutine cat_enkf_increments( & real, dimension( N_catd) :: tp1_ensavg type(obs_param_type) :: this_obs_param + + logical :: nonZeroFound ! ----------------------------------------------------------------------- @@ -4418,7 +4420,8 @@ subroutine cat_enkf_increments( & N_select_varnames = 1 N_state_max = 7 ! Needs to be constant size for applying increment, potential for lots of zeros - allocate( State_incr(N_state_max,N_ens)) + allocate( State_incr(N_state_max,N_ens)) + allocate( State_incr_cum(N_state_max,N_ens)) allocate( State_lon( N_state_max )) allocate( State_lat( N_state_max )) @@ -4432,14 +4435,16 @@ subroutine cat_enkf_increments( & select_varnames(N_varnames) = 'sfds' end if - do ii = 1,N_varnames + do kk=1,N_catd - call get_select_species( & - N_select_varnames, select_varnames(ii), & - N_obs_param, obs_param, N_select_species, select_species ) - - do kk=1,N_catd - + State_incr_cum = 0.0 + + do ii = 1,N_varnames + + call get_select_species( & + N_select_varnames, select_varnames(ii), & + N_obs_param, obs_param, N_select_species, select_species ) + ! compute increments only snow-free and non-frozen tiles if ( (SWE_ensavg(kk) < SWE_threshold) .and. & @@ -4460,7 +4465,7 @@ subroutine cat_enkf_increments( & halo_minlat = max(halo_minlat, -90.) halo_maxlat = min(halo_maxlat, 90.) - call get_ind_obs_lat_lon_box( & + 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), & @@ -4484,34 +4489,34 @@ subroutine cat_enkf_increments( & if (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 + State_incr(1,:) = (cat_progn( kk,:)%srfexc/scale_srfexc) + State_incr_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%catdef/scale_catdef) + State_incr_cum(3,:) elseif ( 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_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%catdef/scale_catdef) + State_incr_cum(3,:) ! catdef in State + State_incr(4,:) = (cat_progn( kk,:)%tc1 /scale_temp) + State_incr_cum(4,:) + State_incr(5,:) = (cat_progn( kk,:)%tc2 /scale_temp) + State_incr_cum(5,:) + State_incr(6,:) = (cat_progn( kk,:)%tc4 /scale_temp) + State_incr_cum(6,:) + State_incr(7,:) = (cat_progn( kk,:)%ght(1)/scale_ght1) + State_incr_cum(7,:) else ! N_state == 6 - State_incr(1,:) = cat_progn( kk,:)%srfexc/scale_srfexc - State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc - 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 + State_incr(1,:) = (cat_progn( kk,:)%srfexc/scale_srfexc) + State_incr_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%tc1 /scale_temp) + State_incr_cum(3,:) + State_incr(4,:) = (cat_progn( kk,:)%tc2 /scale_temp) + State_incr_cum(4,:) + State_incr(5,:) = (cat_progn( kk,:)%tc4 /scale_temp) + State_incr_cum(5,:) + State_incr(6,:) = (cat_progn( kk,:)%ght(1)/scale_ght1) + State_incr_cum(6,:) 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, & @@ -4532,33 +4537,41 @@ subroutine cat_enkf_increments( & fcsterr_inflation_fac ) deallocate(Obs_cov) + + State_incr_cum(1:N_state,:) = State_incr_cum(1:N_state,:) + State_incr(1:N_state,:) + + nonZeroFound = ANY(cat_progn_incr(kk,:)%srfexc /= 0.0) + if (nonZeroFound) then + write (logunit,*) "Non-zero values found. ii = ", ii, & + " kk = ", kk, "State_lon = ", State_lon(1), "State_lat = ", State_lat(1) + end if ! assemble cat_progn increments if (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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%catdef = State_incr_cum(3,:)*scale_catdef elseif (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 ! 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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%catdef = State_incr_cum(3,:)*scale_catdef ! catdef in State + cat_progn_incr(kk,:)%tc1 = State_incr_cum(4,:)*scale_temp + cat_progn_incr(kk,:)%tc2 = State_incr_cum(5,:)*scale_temp + cat_progn_incr(kk,:)%tc4 = State_incr_cum(6,:)*scale_temp + cat_progn_incr(kk,:)%ght(1) = State_incr_cum(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,:)%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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%tc1 = State_incr_cum(3,:)*scale_temp + cat_progn_incr(kk,:)%tc2 = State_incr_cum(4,:)*scale_temp + cat_progn_incr(kk,:)%tc4 = State_incr_cum(5,:)*scale_temp + cat_progn_incr(kk,:)%ght(1) = State_incr_cum(6,:)*scale_ght1 end if @@ -4566,9 +4579,9 @@ subroutine cat_enkf_increments( & end if ! thresholds - end do + end do ! varnames - end do !var_names + end do ! N_catd ! ---------------------------------- @@ -4581,6 +4594,7 @@ subroutine cat_enkf_increments( & ! clean up if (allocated( State_incr )) deallocate( State_incr ) + if (allocated( State_incr_cum )) deallocate( State_incr_cum ) if (allocated( State_lon )) deallocate( State_lon ) if (allocated( State_lat )) deallocate( State_lat ) if (allocated( select_tilenum )) deallocate( select_tilenum ) From 3d6615fab1ce5d349773a5d653b880a1854d1fba Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 4 Oct 2023 17:15:05 -0400 Subject: [PATCH 08/20] combining obs --- .../clsm_ensupd_upd_routines.F90 | 298 +++++++++++++++--- 1 file changed, 251 insertions(+), 47 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index d480a48b..f357a388 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3349,7 +3349,7 @@ end subroutine get_obs_pert ! ********************************************************************* - subroutine cat_enkf_increments( & + subroutine cat_enkf_increments( & N_ens, N_obs, N_catd, N_obs_param, & update_type, obs_param, & tile_coord, l2f, & @@ -3438,7 +3438,7 @@ subroutine cat_enkf_increments( & integer, dimension(N_obs) :: ind_obs - real, allocatable, dimension(:,:) :: State_incr + real, allocatable, dimension(:,:) :: State_incr, State_incr_cum real, allocatable, dimension(:,:) :: Obs_cov ! measurement error covariance real, allocatable, dimension(:) :: State_lon, State_lat @@ -3474,6 +3474,8 @@ subroutine cat_enkf_increments( & real, dimension( N_catd) :: tp1_ensavg type(obs_param_type) :: this_obs_param + + logical :: nonZeroFound ! ----------------------------------------------------------------------- @@ -4407,7 +4409,197 @@ subroutine cat_enkf_increments( & ! ---------------------------------- - case (13) select_update_type ! All observation types from obs_param + case (13) select_update_type ! 3d analysis + + ! update each tile separately using all observations within customized halo around each tile + ! + ! amfox, 14 September 2023 + + if (logit) write (logunit,*) 'Get increments for all observation types in obs_param' + + N_select_varnames = 0 + + if (any(obs_param%varname == 'Tb')) then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'Tb' + end if + + if (any(obs_param%varname == 'sfds')) then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'sfds' + end if + + ! 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 ) + + 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 = 0 + + ! 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 + 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 + + ! call get_select_species(1, 'sfds', N_obs_param, obs_param, N_select_species, select_species ) + + ! do ii = 1,N_select_species + ! if ( any(select_species(ii) == Observations(ind_obs(1:N_selected_obs))%species)) then + ! N_state = max(3, N_state) ! Keep 6 or 7 if we have Tb _and_ sfds obs) + ! if ( N_state > 3) then + ! write (logunit,*) "Two obs types in Catchment = ", kk + ! end if + ! end if + ! end do + + ! ! Determine which observation species are actually selected + ! call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species, select_species ) + + ! do ii = 1,N_select_species + ! if ( any(select_species(ii) == Observations(ind_obs(1:N_selected_obs))%species)) then + ! if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + ! N_state = 7 + ! else + ! N_state = 6 + ! end if + ! end if + ! end do + + N_state = 6 + + if ( N_state == 0 ) then + err_msg = 'Dont know what state to use with this observation type' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! assemble State_minus + ! (on input, cat_progn contains cat_progn_minus) + + if ( 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) + + elseif ( 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 + + 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,:)%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 ) + + ! EnKF update + write (logunit,*) "N_state = ", N_state + + 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==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 + + elseif ( 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 ! 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 + + 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(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 (33) select_update_type ! All observation types from obs_param ! update each tile separately using all observations within customized halo around each tile ! @@ -4418,7 +4610,8 @@ subroutine cat_enkf_increments( & N_select_varnames = 1 N_state_max = 7 ! Needs to be constant size for applying increment, potential for lots of zeros - allocate( State_incr(N_state_max,N_ens)) + allocate( State_incr(N_state_max,N_ens)) + allocate( State_incr_cum(N_state_max,N_ens)) allocate( State_lon( N_state_max )) allocate( State_lat( N_state_max )) @@ -4432,14 +4625,16 @@ subroutine cat_enkf_increments( & select_varnames(N_varnames) = 'sfds' end if - do ii = 1,N_varnames + do kk=1,N_catd - call get_select_species( & - N_select_varnames, select_varnames(ii), & - N_obs_param, obs_param, N_select_species, select_species ) - - do kk=1,N_catd - + State_incr_cum = 0.0 + + do ii = 1,N_varnames + + call get_select_species( & + N_select_varnames, select_varnames(ii), & + N_obs_param, obs_param, N_select_species, select_species ) + ! compute increments only snow-free and non-frozen tiles if ( (SWE_ensavg(kk) < SWE_threshold) .and. & @@ -4460,7 +4655,7 @@ subroutine cat_enkf_increments( & halo_minlat = max(halo_minlat, -90.) halo_maxlat = min(halo_maxlat, 90.) - call get_ind_obs_lat_lon_box( & + 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), & @@ -4484,34 +4679,34 @@ subroutine cat_enkf_increments( & if (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 + State_incr(1,:) = (cat_progn( kk,:)%srfexc/scale_srfexc) + State_incr_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%catdef/scale_catdef) + State_incr_cum(3,:) elseif ( 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_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%catdef/scale_catdef) + State_incr_cum(3,:) ! catdef in State + State_incr(4,:) = (cat_progn( kk,:)%tc1 /scale_temp) + State_incr_cum(4,:) + State_incr(5,:) = (cat_progn( kk,:)%tc2 /scale_temp) + State_incr_cum(5,:) + State_incr(6,:) = (cat_progn( kk,:)%tc4 /scale_temp) + State_incr_cum(6,:) + State_incr(7,:) = (cat_progn( kk,:)%ght(1)/scale_ght1) + State_incr_cum(7,:) else ! N_state == 6 - State_incr(1,:) = cat_progn( kk,:)%srfexc/scale_srfexc - State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc - 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 + State_incr(1,:) = (cat_progn( kk,:)%srfexc/scale_srfexc) + State_incr_cum(1,:) + State_incr(2,:) = (cat_progn( kk,:)%rzexc /scale_rzexc) + State_incr_cum(2,:) + State_incr(3,:) = (cat_progn( kk,:)%tc1 /scale_temp) + State_incr_cum(3,:) + State_incr(4,:) = (cat_progn( kk,:)%tc2 /scale_temp) + State_incr_cum(4,:) + State_incr(5,:) = (cat_progn( kk,:)%tc4 /scale_temp) + State_incr_cum(5,:) + State_incr(6,:) = (cat_progn( kk,:)%ght(1)/scale_ght1) + State_incr_cum(6,:) 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, & @@ -4532,33 +4727,41 @@ subroutine cat_enkf_increments( & fcsterr_inflation_fac ) deallocate(Obs_cov) + + State_incr_cum(1:N_state,:) = State_incr_cum(1:N_state,:) + State_incr(1:N_state,:) + + nonZeroFound = ANY(cat_progn_incr(kk,:)%srfexc /= 0.0) + if (nonZeroFound) then + write (logunit,*) "Non-zero values found. ii = ", ii, & + " kk = ", kk, "State_lon = ", State_lon(1), "State_lat = ", State_lat(1) + end if ! assemble cat_progn increments if (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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%catdef = State_incr_cum(3,:)*scale_catdef elseif (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 ! 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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%catdef = State_incr_cum(3,:)*scale_catdef ! catdef in State + cat_progn_incr(kk,:)%tc1 = State_incr_cum(4,:)*scale_temp + cat_progn_incr(kk,:)%tc2 = State_incr_cum(5,:)*scale_temp + cat_progn_incr(kk,:)%tc4 = State_incr_cum(6,:)*scale_temp + cat_progn_incr(kk,:)%ght(1) = State_incr_cum(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,:)%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 + cat_progn_incr(kk,:)%srfexc = State_incr_cum(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr_cum(2,:)*scale_rzexc + cat_progn_incr(kk,:)%tc1 = State_incr_cum(3,:)*scale_temp + cat_progn_incr(kk,:)%tc2 = State_incr_cum(4,:)*scale_temp + cat_progn_incr(kk,:)%tc4 = State_incr_cum(5,:)*scale_temp + cat_progn_incr(kk,:)%ght(1) = State_incr_cum(6,:)*scale_ght1 end if @@ -4566,9 +4769,9 @@ subroutine cat_enkf_increments( & end if ! thresholds - end do + end do ! varnames - end do !var_names + end do ! N_catd ! ---------------------------------- @@ -4581,6 +4784,7 @@ subroutine cat_enkf_increments( & ! clean up if (allocated( State_incr )) deallocate( State_incr ) + if (allocated( State_incr_cum )) deallocate( State_incr_cum ) if (allocated( State_lon )) deallocate( State_lon ) if (allocated( State_lat )) deallocate( State_lat ) if (allocated( select_tilenum )) deallocate( select_tilenum ) From 718a55af227d42d7f361566542bab8659d242507 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 6 Dec 2023 11:52:19 -0700 Subject: [PATCH 09/20] add select_species_Tb --- .../clsm_ensupd_upd_routines.F90 | 260 ++++-------------- 1 file changed, 46 insertions(+), 214 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 8188bdab..0640f998 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3422,9 +3422,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, N_varnames + integer :: N_state_max, N_state, N_selected_obs, N_select_varnames, N_select_species, N_varnames, N_select_species_Tb real :: halo_minlon, halo_maxlon, halo_minlat, halo_maxlat real :: tmp_minlon, tmp_maxlon, tmp_minlat, tmp_maxlat @@ -3441,7 +3441,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 @@ -3451,6 +3451,7 @@ subroutine cat_enkf_increments( & integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij => null() character(len=*), parameter :: Iam = 'cat_enkf_increments' + character(len=400) :: err_msg real, dimension( N_catd) :: r_x, tmp_dlon real :: r_y, tmp_dlat @@ -3472,7 +3473,7 @@ subroutine cat_enkf_increments( & type(obs_param_type) :: this_obs_param - logical :: nonZeroFound + logical :: nonZeroFound, found, smap_obs, ascat_obs ! ----------------------------------------------------------------------- @@ -4433,6 +4434,10 @@ subroutine cat_enkf_increments( & 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)) @@ -4471,37 +4476,30 @@ subroutine cat_enkf_increments( & if (N_selected_obs>0) then - ! call get_select_species(1, 'sfds', N_obs_param, obs_param, N_select_species, select_species ) - - ! do ii = 1,N_select_species - ! if ( any(select_species(ii) == Observations(ind_obs(1:N_selected_obs))%species)) then - ! N_state = max(3, N_state) ! Keep 6 or 7 if we have Tb _and_ sfds obs) - ! if ( N_state > 3) then - ! write (logunit,*) "Two obs types in Catchment = ", kk - ! end if - ! end if - ! end do - - ! ! Determine which observation species are actually selected - ! call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species, select_species ) - - ! do ii = 1,N_select_species - ! if ( any(select_species(ii) == Observations(ind_obs(1:N_selected_obs))%species)) then - ! if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then - ! N_state = 7 - ! else - ! N_state = 6 - ! end if - ! end if - ! end do - - N_state = 6 - - if ( N_state == 0 ) then - err_msg = 'Dont know what state to use with this observation type' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if - + ! Determine if Tb observations are present + + 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 = .true. + exit + end if + end do + if (found) then + exit + end if + end do + + if (found) then + if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + N_state = 7 + else + N_state = 6 + end if + else + N_state = 3 + end if + ! assemble State_minus ! (on input, cat_progn contains cat_progn_minus) @@ -4543,8 +4541,15 @@ subroutine cat_enkf_increments( & Observations(ind_obs(1:N_selected_obs)), Obs_cov ) ! EnKF update - write (logunit,*) "N_state = ", N_state + ! smap_obs = any(Observations(ind_obs(1:N_selected_obs))%species < 5) + ! ascat_obs = any(Observations(ind_obs(1:N_selected_obs))%species > 4) + + ! if (smap_obs .and. ascat_obs) then + ! write (logunit,*) "Found tile with both obs types. Tile = ", halo_minlon, halo_maxlon, halo_minlat, halo_maxlat + ! write (logunit,*) "Observations(ind_obs(1:N_selected_obs))%species = ", Observations(ind_obs(1:N_selected_obs))%species + ! end if + call enkf_increments( & N_state, N_selected_obs, N_ens, & Observations(ind_obs(1:N_selected_obs)), & @@ -4564,23 +4569,26 @@ subroutine cat_enkf_increments( & 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,:)%catdef = State_incr(3,:)*scale_catdef + write(logunit,*) 'N3 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc elseif ( 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 ! catdef in State + write(logunit,*) 'N7 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc 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 + write(logunit,*) 'N6 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc cat_progn_incr(kk,:)%tc1 = State_incr(3,:)*scale_temp cat_progn_incr(kk,:)%tc2 = State_incr(4,:)*scale_temp @@ -4597,182 +4605,6 @@ subroutine cat_enkf_increments( & ! ---------------------------------- - case (33) select_update_type ! All observation types from obs_param - - ! update each tile separately using all observations within customized halo around each tile - ! - ! amfox, 27 July 2023 - - if (logit) write (logunit,*) 'Get increments for all observation types in obs_param' - - N_select_varnames = 1 - N_state_max = 7 ! Needs to be constant size for applying increment, potential for lots of zeros - - allocate( State_incr(N_state_max,N_ens)) - allocate( State_incr_cum(N_state_max,N_ens)) - allocate( State_lon( N_state_max )) - allocate( State_lat( N_state_max )) - - if (any(obs_param%varname == 'Tb')) then - N_varnames = N_varnames + 1 - select_varnames(N_varnames) = 'Tb' - end if - - if (any(obs_param%varname == 'sfds')) then - N_varnames = N_varnames + 1 - select_varnames(N_varnames) = 'sfds' - end if - - do kk=1,N_catd - - State_incr_cum = 0.0 - - do ii = 1,N_varnames - - call get_select_species( & - N_select_varnames, select_varnames(ii), & - N_obs_param, obs_param, N_select_species, select_species ) - - ! 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 - 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 - - if (select_varnames(ii)=='Tb' .and. cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then - N_state = 7 - elseif (select_varnames(ii)=='Tb' .and. cat_param(kk)%poros 4) - - ! if (smap_obs .and. ascat_obs) then - ! write (logunit,*) "Found tile with both obs types. Tile = ", halo_minlon, halo_maxlon, halo_minlat, halo_maxlat - ! write (logunit,*) "Observations(ind_obs(1:N_selected_obs))%species = ", Observations(ind_obs(1:N_selected_obs))%species - ! end if call enkf_increments( & N_state, N_selected_obs, N_ens, & @@ -4614,7 +4604,6 @@ subroutine cat_enkf_increments( & ! clean up if (allocated( State_incr )) deallocate( State_incr ) - if (allocated( State_incr_cum )) deallocate( State_incr_cum ) if (allocated( State_lon )) deallocate( State_lon ) if (allocated( State_lat )) deallocate( State_lat ) if (allocated( select_tilenum )) deallocate( select_tilenum ) @@ -5176,7 +5165,7 @@ subroutine check_compact_support( & Iam // '(): reset for 1d update_type: ycompact = ', ycompact if (logit) write (logunit,*) - case (2,7,8,10, 13) ! "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 From 3f1842012c3518e9ca4c87647d16baf08503b0b1 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 6 Feb 2024 17:29:49 -0700 Subject: [PATCH 13/20] more tidying --- .../clsm_ensupd_upd_routines.F90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index cf3f39d0..17c91b9f 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3452,7 +3452,6 @@ subroutine cat_enkf_increments( & integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij => null() character(len=*), parameter :: Iam = 'cat_enkf_increments' - character(len=400) :: err_msg real, dimension( N_catd) :: r_x, tmp_dlon real :: r_y, tmp_dlat @@ -4559,16 +4558,14 @@ subroutine cat_enkf_increments( & 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 - write(logunit,*) 'N3 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc + cat_progn_incr(kk,:)%catdef = State_incr(3,:)*scale_catdef elseif ( 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 ! catdef in State - write(logunit,*) 'N7 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc - + 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 @@ -4578,7 +4575,6 @@ subroutine cat_enkf_increments( & cat_progn_incr(kk,:)%srfexc = State_incr(1,:)*scale_srfexc cat_progn_incr(kk,:)%rzexc = State_incr(2,:)*scale_rzexc - write(logunit,*) 'N6 After cat_progn_incr(kk,:)%srfexc = ', cat_progn_incr(kk,5)%srfexc cat_progn_incr(kk,:)%tc1 = State_incr(3,:)*scale_temp cat_progn_incr(kk,:)%tc2 = State_incr(4,:)*scale_temp From 96373047d2faf807ee17488683500bda7bb15a50 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 6 Feb 2024 17:32:44 -0700 Subject: [PATCH 14/20] remove white space --- .../GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 17c91b9f..ba42aafa 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4449,8 +4449,8 @@ subroutine cat_enkf_increments( & ! compute increments only snow-free and non-frozen tiles - if ( (SWE_ensavg(kk) < SWE_threshold) .and. & - (tp1_ensavg(kk) > tp1_threshold) ) then + if ( (SWE_ensavg(kk) < SWE_threshold) .and. & + (tp1_ensavg(kk) > tp1_threshold)) then ! find observations within halo around tile kk From bc0cef5f0ecf04ada191e77caebf55bee92d0871 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 9 Feb 2024 18:02:36 -0500 Subject: [PATCH 15/20] minor edits for new update_type=13 (clsm_ensupd_upd_routines.F90): - added 'sfmc' obs - added/edited comments - fixed indent - minor simplifications --- .../clsm_ensupd_upd_routines.F90 | 370 +++++++++--------- 1 file changed, 187 insertions(+), 183 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index ba42aafa..87db1a7f 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3363,7 +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 + ! amfox, 6 Feb 2024 - added update type 13 for combination of ASCAT SM and SMAP Tb ! ! -------------------------------------------------------------- @@ -3473,7 +3473,7 @@ subroutine cat_enkf_increments( & type(obs_param_type) :: this_obs_param - logical :: found + logical :: found_Tb_obs ! ----------------------------------------------------------------------- @@ -3666,7 +3666,7 @@ 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 @@ -4407,190 +4407,194 @@ subroutine cat_enkf_increments( & ! ---------------------------------- - case (13) select_update_type ! 3d analysis - - ! update each tile separately using all observations within customized halo around each tile - ! - ! amfox, 14 September 2023 - - if (logit) write (logunit,*) 'Get increments for all observation types in obs_param' - - N_select_varnames = 0 - - if (any(obs_param%varname == 'Tb')) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'Tb' - end if - - if (any(obs_param%varname == 'sfds')) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'sfds' - end if - - ! 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 = 0 - - ! 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 - 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 - - 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 = .true. - exit - end if - end do - if (found) then - exit - end if - end do - - if (found) then - if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then - N_state = 7 - else - N_state = 6 - end if - else - N_state = 3 - end if - - ! assemble State_minus - ! (on input, cat_progn contains cat_progn_minus) - - if ( 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) - - elseif ( 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 - - 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,:)%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==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 - - elseif ( 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 ! 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 + case (13) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs + + ! combines update_types 2 and 10 + ! + ! update each tile separately using all observations within customized halo around each tile + ! + ! amfox, 14 September 2023 + + if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' + + N_select_varnames = 0 + + if (any(obs_param%varname == 'Tb')) then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'Tb' + end if + + if (any(obs_param%varname == 'sfmc')) then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'sfmc' + end if - 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,:)%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 + if (any(obs_param%varname == 'sfds')) then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'sfds' + end if + + ! 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 = 3 ! for now, assume only sfmc and/or sfds obs + + ! 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 + 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 + + 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 (found_Tb_obs) then + if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + N_state = 7 + else + N_state = 6 + end if + end if + + ! assemble State_minus + ! (on input, cat_progn contains cat_progn_minus) + + if ( 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) + + elseif ( 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 + + 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,:)%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)) - end if ! thresholds - - end do + 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==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 + + elseif ( 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 ! 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 + + 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(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 call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown update_type') From 7e5b3c8e3b02352d4e4a698e39a898f93ebf10f0 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 19 Feb 2024 14:31:58 -0700 Subject: [PATCH 16/20] add trim when selecting by varname --- .../clsm_ensupd_upd_routines.F90 | 30 +++++++++++-------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 87db1a7f..c8710d6e 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4419,20 +4419,24 @@ subroutine cat_enkf_increments( & N_select_varnames = 0 - if (any(obs_param%varname == 'Tb')) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'Tb' - end if - - if (any(obs_param%varname == 'sfmc')) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'sfmc' - end if + 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' + end if + + if (trim(obs_param(ii)%varname) == 'sfmc') then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'sfmc' + end if + + if (trim(obs_param(ii)%varname) == 'sfds') then + N_select_varnames = N_select_varnames + 1 + select_varnames(N_select_varnames) = 'sfds' + end if - if (any(obs_param%varname == 'sfds')) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'sfds' - end if + end do ! Will get all species associated with Tb or sfds observations From f4780912f1912c26612f7504d0fbdb8fcff12c48 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 19 Feb 2024 14:38:37 -0700 Subject: [PATCH 17/20] Clarify comments in clsm_ensupd_enkf_update.F90 --- .../GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 30610618..16449a65 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1333,8 +1333,8 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & case (6,8,9,10,13) select_update_type ! soil moisture and temperature update ! for update_type 10, catdef increments may be zero by design - ! for update_type 13, could be multiple zero increments - + ! for update_type 13, could produce multiple zero increments (i.e. tc1, tc2, tc4 will be unchanged if soil moisture observations only) + if (logit) write (logunit,*) & 'apply_enkf_increments(): applying soil moisture and Tskin/ght1 increments' From e78cc6bedb1c0a6834c3e7ef6bde4cef5a5a1702 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 19 Feb 2024 20:59:55 -0700 Subject: [PATCH 18/20] Fix issue with select_varnames --- .../clsm_ensupd_upd_routines.F90 | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index c8710d6e..8ff56681 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4420,22 +4420,27 @@ subroutine cat_enkf_increments( & 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' + 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' + 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' - end if - + select_varnames(N_select_varnames) = 'sfds' + exit + end if end do ! Will get all species associated with Tb or sfds observations From 47d4bb98203011c98702c15f4a9282bf280b339f Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 26 Feb 2024 12:46:18 -0700 Subject: [PATCH 19/20] remove catdef from state for non-peatland tiles --- .../clsm_ensupd_upd_routines.F90 | 30 ++++++++++++++----- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 8ff56681..05474126 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4461,7 +4461,7 @@ subroutine cat_enkf_increments( & do kk=1,N_catd - N_state = 3 ! for now, assume only sfmc and/or sfds obs + N_state = 2 ! for now, assume only sfmc and/or sfds observations an no peatland catchments ! compute increments only snow-free and non-frozen tiles @@ -4490,6 +4490,12 @@ subroutine cat_enkf_increments( & N_selected_obs, ind_obs ) if (N_selected_obs>0) then + + ! Determine if peatland catchment and set N_state accordingly for SM observations + + if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then + N_state = 3 + end if ! Determine if Tb observations are present @@ -4514,11 +4520,16 @@ subroutine cat_enkf_increments( & ! assemble State_minus ! (on input, cat_progn contains cat_progn_minus) - if ( N_state==3 ) then + 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) + 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==7 ) then @@ -4567,11 +4578,16 @@ subroutine cat_enkf_increments( & ! assemble cat_progn increments - if ( N_state==3 ) then + 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 + cat_progn_incr(kk,:)%catdef = State_incr(3,:)*scale_catdef ! catdef in State elseif ( N_state==7 ) then From 9f02392ce6dc5766477a5fc5ff82619e8e2aa032 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 27 Feb 2024 07:51:06 -0500 Subject: [PATCH 20/20] disabled updated_type=[1,2]; cleanup of new update_type=13 (clsm_ensupd_upd_routines.F90, clsm_ensupd_enkf_update.F90, LDASsa_DEFAULT_inputs_ensupd.nml) --- .../LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml | 30 +++++----- .../clsm_ensupd_enkf_update.F90 | 5 +- .../clsm_ensupd_upd_routines.F90 | 59 +++++++++++-------- 3 files changed, 54 insertions(+), 40 deletions(-) diff --git a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml index b849d51d..afacd0f3 100644 --- a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -19,20 +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 -! update_type = 13: Multivariate combination of 2 and 10; sfmc and 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 diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 16449a65..d90ae7cc 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1332,8 +1332,9 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & case (6,8,9,10,13) select_update_type ! soil moisture and temperature update - ! for update_type 10, catdef increments may be zero by design - ! for update_type 13, could produce multiple zero increments (i.e. tc1, tc2, tc4 will be unchanged if soil moisture observations only) + ! 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) if (logit) write (logunit,*) & 'apply_enkf_increments(): applying soil moisture and Tskin/ght1 increments' diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 05474126..c402cbf8 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3604,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' @@ -3673,6 +3677,10 @@ subroutine cat_enkf_increments( & 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' @@ -4409,11 +4417,18 @@ subroutine cat_enkf_increments( & case (13) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs - ! combines update_types 2 and 10 - ! ! update each tile separately using all observations within customized halo around each tile ! - ! amfox, 14 September 2023 + ! 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' @@ -4461,12 +4476,12 @@ subroutine cat_enkf_increments( & do kk=1,N_catd - N_state = 2 ! for now, assume only sfmc and/or sfds observations an no peatland catchments + N_state = 2 ! initialize (always have srfexc and rzexc in state vector) - ! compute increments only snow-free and non-frozen tiles + ! compute increments only for snow-free and non-frozen tiles - if ( (SWE_ensavg(kk) < SWE_threshold) .and. & - (tp1_ensavg(kk) > tp1_threshold)) then + if ( (SWE_ensavg(kk) < SWE_threshold) .and. & + (tp1_ensavg(kk) > tp1_threshold) ) then ! find observations within halo around tile kk @@ -4491,14 +4506,10 @@ subroutine cat_enkf_increments( & if (N_selected_obs>0) then - ! Determine if peatland catchment and set N_state accordingly for SM observations - - if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then - N_state = 3 - end if - ! 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 @@ -4509,13 +4520,13 @@ subroutine cat_enkf_increments( & if (found_Tb_obs) exit end do - if (found_Tb_obs) then - if (cat_param(kk)%poros>=PEATCLSM_POROS_THRESHOLD) then - N_state = 7 - else - N_state = 6 - end if - end if + ! 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) @@ -4531,11 +4542,10 @@ subroutine cat_enkf_increments( & State_incr(2,:) = cat_progn( kk,:)%rzexc /scale_rzexc State_incr(3,:) = cat_progn( kk,:)%catdef/scale_catdef ! catdef in State - elseif ( N_state==7 ) then + 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(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 @@ -4546,6 +4556,7 @@ subroutine cat_enkf_increments( & 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 @@ -4589,11 +4600,10 @@ subroutine cat_enkf_increments( & 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==7 ) then + 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,:)%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 @@ -4604,6 +4614,7 @@ subroutine cat_enkf_increments( & 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