From 829f229baa6b7bd1613c7120465d286059787708 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 2 Dec 2024 14:29:55 -0500 Subject: [PATCH] switched order of code blocks for update_type=13 and update_type=12 in "select case" (clsm_ensupd_upd_routines.F90) --- .../clsm_ensupd_upd_routines.F90 | 394 +++++++++--------- 1 file changed, 197 insertions(+), 197 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 85e7c1c..0ab23b8 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4690,202 +4690,6 @@ subroutine cat_enkf_increments( & ! ---------------------------------------------------------------------------------------------------------------------- - 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+sfds obs' - - ! Get all species associated with *assimilated* Tb, sfmc, and sfds observations - - N_select_varnames = 3 - - select_varnames(1) = 'Tb' - select_varnames(2) = 'sfmc' - select_varnames(3) = 'sfds' - - 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(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 - - 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(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 - - 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(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 - - 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(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 - - end if - - end if - - end if ! thresholds - - end do - - ! ---------------------------------------------------------------------------------------------------------------------- - case (12) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs ! & 1d snow analysis (Toure et al. 2018 empirical gain); snow cover fraction obs @@ -5308,7 +5112,203 @@ subroutine cat_enkf_increments( & end if ! if (N_select_species_smTb>0) [assimilate soil moisture observations for tile kk] - end do ! kk=1,N_catd + end do ! kk=1,N_catd + + ! ---------------------------------------------------------------------------------------------------------------------- + + 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+sfds obs' + + ! Get all species associated with *assimilated* Tb, sfmc, and sfds observations + + N_select_varnames = 3 + + select_varnames(1) = 'Tb' + select_varnames(2) = 'sfmc' + select_varnames(3) = 'sfds' + + 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(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 + + 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(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 + + 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(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 + + 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(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 + + end if + + end if + + end if ! thresholds + + end do ! kk=1,N_catd ! ----------------------------------------------------------------------------------------------------------------------