From 5cec2190f59624e36079602f6cc24e5b7fca4582 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 31 Oct 2024 10:39:30 -0600 Subject: [PATCH 01/21] add snow to update type 13 --- .../clsm_ensupd_upd_routines.F90 | 245 ++++++++++++++++-- 1 file changed, 222 insertions(+), 23 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index c8398b6..e7dd848 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3524,6 +3524,7 @@ subroutine cat_enkf_increments( & real, dimension( N_catd) :: SWE_ensavg real, dimension( N_catd) :: tp1_ensavg real, dimension( N_catd) :: asnow_ensavg + real, dimension( N_catd) :: swe_incr_ensavg type(obs_param_type) :: this_obs_param @@ -3537,7 +3538,7 @@ subroutine cat_enkf_increments( & real, dimension(N_snow) :: tpsn, fice_snow_vec ! for snow model relayer real, dimension(N_snow,N_constit) :: rconstit - logical :: found_Tb_obs + logical :: found_Tb_obs, unique, has_asnow ! ----------------------------------------------------------------------- @@ -4690,31 +4691,228 @@ subroutine cat_enkf_increments( & 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 = 0 + + do ii = 1, N_obs_param + unique = .true. + do jj = 1, N_select_varnames + if (trim(obs_param(ii)%varname) == trim(select_varnames(jj))) then + unique = .false. + exit + end if + end do + if (unique) then N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'Tb' - exit + select_varnames(N_select_varnames) = trim(obs_param(ii)%varname) 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' + end do + + ! Check if we potentially have snow fraction observations by looking for 'asnow' in select_varnames + + do ii = 1, N_select_varnames + if (trim(select_varnames(ii)) == 'asnow') then + has_asnow = .true. exit end if - end do + 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 + ! If we have snow fraction observations, we need to calculate the increments for the snow layers + + if (has_asnow) then + + if (logit) write (logunit, *) 'get 1d snow increments (Toure et al. 2018 empirical gain); snow cover fraction obs' + + ! ensure that max SWE increment parameter is less than WEMIN; larger increments make no sense because + ! at SWE=WEMIN, the tile is fully snow covered (asnow=1) + + if (SCF_ANA_MAXINCRSWE>WEMIN) call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'must use SCF_ANA_MAXINCRSWE<=WEMIN') + + ! identify the obs species of interest + + N_select_varnames = 1 + + select_varnames(1) = 'asnow' + + call get_select_species( & + N_select_varnames, select_varnames(1:N_select_varnames), & + N_obs_param, obs_param, N_select_species, select_species ) + + allocate(select_tilenum(1)) + + swe_incr = 0. ! total SWE increment; initialize to NO CHANGE + swe_incr_ensavg = 0. ! total SWE ensemble average increment; initialize to NO CHANGE + + ! loop through tiles and compute increments + + do kk=1,N_catd + + ! find observations for tile kk + + select_tilenum(1) = l2f(kk) + + call get_ind_obs( & + N_obs, Observations, & + 1, select_tilenum, & + N_select_species, select_species(1:N_select_species), & + N_selected_obs, ind_obs ) + + if (N_selected_obs > 0) then + + ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) + + tmp_obs = sum(Observations(ind_obs(1:N_selected_obs))%obs) + + if (N_selected_obs > 1) tmp_obs = tmp_obs/real(N_selected_obs) + + do n_e=1,N_ens ! compute analysis separately for each ensemble member + + ! 1. Diagnose model forecast snow cover area fraction and total SWE + + asnow_fcst = asnow(kk,n_e) + swe_fcst = sum(cat_progn(kk,n_e)%wesn(1:N_snow)) + + ! 2. Calculate SWE increment based on modified eq 1 of Toure et al (2018) + + if (asnow_fcst .lt. tmp_obs * SCF_ANA_ALPHA) then + + ! ADD SNOW: Forecast SCF is less than observed SCF (after "bias" adjustment with alpha) + + swe_incr(kk,n_e) = SCF_ANA_MAXINCRSWE * (tmp_obs - asnow_fcst/SCF_ANA_ALPHA) + + elseif (tmp_obs .lt. SCF_ANA_BETA) then + + ! REMOVE SNOW: Simulated SCF is greater than observed SCF (after "bias" adjustment) + ! and observed SCF is less than beta threshold + + swe_incr(kk,n_e) = (-1.) * SCF_ANA_MAXINCRSWE * asnow_fcst * (1. - tmp_obs/SCF_ANA_BETA) + + else + + cycle ! NO CHANGE, skip rest of increment calcs and go straight to next ens member + + endif ! (Toure et al. 2018 Equation 1) + + ! 3. Derive SWE, snow heat content, and snow depth increments for each layer from total SWE increment + + swe_ana = max(swe_fcst + swe_incr(kk,n_e), 0.0) ! total SWE after analysis + + call StieglitzSnow_calc_asnow( swe_ana, asnow_ana ) ! asnow after analysis + + if (swe_fcst>=StieglitzSnow_MINSWE) then + swe_ratio = swe_ana / swe_fcst + else + swe_ratio = MAPL_UNDEF ! swe_ratio unreliable; set to MAPL_UNDEF to expose inadvertent use + end if + + ! loop through snow layers and compute SWE, snow heat content, and snow depth analysis for each layer + + do isnow=1,N_snow + + if (asnow_ana == 0.0) then + + ! no snow in analysis, remove all snow + + tmp_wesn(kk,n_e,isnow) = 0.0 + tmp_htsn(kk,n_e,isnow) = 0.0 + tmp_sndz(kk,n_e,isnow) = 0.0 + + elseif (swe_fcst < StieglitzSnow_MINSWE) then + + ! too little snow in forecast, use generic properties for added snow + + tmp_wesn(kk,n_e,isnow) = swe_ana / N_snow ! distribute SWE evenly across layers + + ! assign heat content for snow at 0 deg C and without liquid water content (100% frozen) + ! (based on StieglitzSnow: htsn = (CPW*tsnow - fice*MAPL_ALHF)*swe ) + + tmp_htsn(kk,n_e,isnow) = (0.0 - MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) + + ! assign snow depth consistent with density of freshly fallen snow (must have SCF_ANA_MAXINCRSWE<=WEMIN) + + tmp_sndz(kk,n_e,isnow) = (WEMIN / RHOFS) / N_snow + + else + + ! snow in forecast and analysis, derive properties of analysis snow from properties of forecast snow + + ! update SWE: + + tmp_wesn(kk,n_e,isnow) = cat_progn(kk,n_e)%wesn(isnow) * swe_ratio + + ! update snow heat content (keep snow temperature constant): + + call StieglitzSnow_calc_tpsnow( cat_progn(kk,n_e)%htsn(isnow), cat_progn(kk,n_e)%wesn(isnow), & + snow_temp, fice_snow, log_dum, log_dum2, .false. ) + + tmp_htsn(kk,n_e,isnow) = (StieglitzSnow_CPW*snow_temp - fice_snow*MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) + + ! update snow depth: + + if (asnow_ana < 1. .and. asnow_fcst < 1.) then + + ! keep snow depth constant when less than full snow cover in fcst and ana + + tmp_sndz(kk,n_e,isnow) = cat_progn(kk,n_e)%sndz(isnow) + + else + + ! compute analysis snow depth by keeping snow density constant + ! + ! in this case, it is possible that either asnow_fcst<1 or asnow_ana<1; + ! when computing density or depth, make sure that SWE value (which is per unit area) is + ! adjusted to reflect SWE value (per unit area) in the snow-covered fraction of the tile + + ! i) diagnose (layer-specific) forecast snow density + + snow_dens = ( cat_progn(kk,n_e)%wesn(isnow)/asnow_fcst ) / cat_progn(kk,n_e)%sndz(isnow) + + ! ii) diagnose analysis snow depth using forecast density + + tmp_sndz(kk,n_e,isnow) = ( tmp_wesn(kk,n_e,isnow)/asnow_ana ) / snow_dens + + end if + + end if + + end do ! isnow=1,N_snow (compute SWE, snow heat content, and snow depth analysis for each layer) + + ! 4. Relayer to balance the snow column (call with optional args for adjustment of htsnn) + + call StieglitzSnow_relayer( N_snow, N_constit, & + MAPL_LAND, CATCH_SNOW_DZPARAM, & + tmp_htsn(kk,n_e,1:N_snow), & + tmp_wesn(kk,n_e,1:N_snow), & + tmp_sndz(kk,n_e,1:N_snow), & + rconstit, tpsn, fice_snow_vec ) + + ! print the old and new swe, heat content and snow density + + !if (logit) write (logunit, *) & + ! 'fcst_wesn = ', cat_progn(kk, n_e)%wesn(1:N_snow), & + ! 'tmp_wesn = ', tmp_wesn( kk,n_e, 1:N_snow), & + ! 'fcst_htsn = ', cat_progn(kk, n_e)%htsn(1:N_snow), & + ! 'tmp_htsn = ', tmp_htsn( kk, n_e, 1:N_snow), & + ! 'fcst_sndz = ', cat_progn(kk, n_e)%sndz(1:N_snow), & + ! 'tmp_sndz = ', tmp_sndz( kk ,n_e, 1:N_snow), & + ! '--------------------------------------' + + ! 5. Diagnose increments + + cat_progn_incr(kk,n_e)%wesn(1:N_snow) = tmp_wesn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%wesn(1:N_snow) + cat_progn_incr(kk,n_e)%htsn(1:N_snow) = tmp_htsn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%htsn(1:N_snow) + cat_progn_incr(kk,n_e)%sndz(1:N_snow) = tmp_sndz(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%sndz(1:N_snow) + + end do ! n_e=1,N_ens + + end if ! if (N_selected_obs > 0) + + end do ! kk=1,N_catd + + end if ! if (has_asnow) + + ! Calculate swe_incr_ensavg as the average of the ensemble members to add to check if tile is snow-free + + swe_incr_ensavg = sum(swe_incr, dim=2) / real(N_ens) ! Will get all species associated with Tb or sfds observations @@ -4738,8 +4936,9 @@ subroutine cat_enkf_increments( & ! 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. & + (swe_incr_ensavg(kk) < SWE_threshold) .and. & + (tp1_ensavg(kk) > tp1_threshold) ) then ! find observations within halo around tile kk From b4a29c85ae4d495ece49a26321131f6753d0f938 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Sun, 3 Nov 2024 19:22:01 -0500 Subject: [PATCH 02/21] working code --- .../clsm_ensupd_enkf_update.F90 | 16 +++++++- .../clsm_ensupd_upd_routines.F90 | 41 +++++++++++++++++-- 2 files changed, 53 insertions(+), 4 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index ea0be94..c7fce5a 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1396,7 +1396,21 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn(n,n_e)%ght(1) = & cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1) - + + do ii=1,N_snow ! for each snow layer + + cat_progn(n,n_e)%wesn(ii) = & + cat_progn(n,n_e)%wesn(ii) + cat_progn_incr(n,n_e)%wesn(ii) + + cat_progn(n,n_e)%sndz(ii) = & + cat_progn(n,n_e)%sndz(ii) + cat_progn_incr(n,n_e)%sndz(ii) + + cat_progn(n,n_e)%htsn(ii) = & + cat_progn(n,n_e)%htsn(ii) + cat_progn_incr(n,n_e)%htsn(ii) + + end do + + end do end do diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index e52aac6..5e9e8c5 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4689,7 +4689,6 @@ subroutine cat_enkf_increments( & ! ! amfox+rreichle, 26 Feb 2024 - if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' N_select_varnames = 0 @@ -4757,7 +4756,9 @@ subroutine cat_enkf_increments( & N_selected_obs, ind_obs ) if (N_selected_obs > 0) then - + + write(*,*) 'Doing snow DA' + ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) tmp_obs = sum(Observations(ind_obs(1:N_selected_obs))%obs) @@ -4795,7 +4796,9 @@ subroutine cat_enkf_increments( & ! 3. Derive SWE, snow heat content, and snow depth increments for each layer from total SWE increment swe_ana = max(swe_fcst + swe_incr(kk,n_e), 0.0) ! total SWE after analysis - + + write (*,*) 'swe_incr(kk,n_e): ',swe_incr(kk,n_e) + call StieglitzSnow_calc_asnow( swe_ana, asnow_ana ) ! asnow after analysis if (swe_fcst>=StieglitzSnow_MINSWE) then @@ -4915,12 +4918,40 @@ subroutine cat_enkf_increments( & swe_incr_ensavg = sum(swe_incr, dim=2) / real(N_ens) ! Will get all species associated with Tb or sfds observations + + 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 + 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 + + if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb ) @@ -4963,6 +4994,8 @@ subroutine cat_enkf_increments( & if (N_selected_obs>0) then + write(*,*) 'Doing SM DA' + ! Determine if Tb observations are present found_Tb_obs = .false. @@ -5079,6 +5112,8 @@ subroutine cat_enkf_increments( & cat_progn_incr(kk,:)%ght(1) = State_incr(7,:)*scale_ght1 end if + + write (*,*) 'State_incr: ', State_incr(1,1) end if From 995a414d2e487c3eaf5cb6385d3944d491d58413 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 6 Nov 2024 17:48:36 -0700 Subject: [PATCH 03/21] single do loop --- .../clsm_ensupd_upd_routines.F90 | 495 +++++++++--------- 1 file changed, 240 insertions(+), 255 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 5e9e8c5..dbb144e 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3478,7 +3478,7 @@ subroutine cat_enkf_increments( & integer :: n, n_e, kk, ii, jj - integer :: N_state_max, N_state, N_selected_obs, N_select_varnames, N_select_species, N_select_species_Tb + integer :: N_state_max, N_state, N_selected_obs, N_select_varnames, N_select_species, N_select_species_Tb, N_select_species_asnow real :: halo_minlon, halo_maxlon, halo_minlat, halo_maxlat real :: tmp_minlon, tmp_maxlon, tmp_minlat, tmp_maxlat @@ -3495,7 +3495,7 @@ subroutine cat_enkf_increments( & real, allocatable, dimension(:) :: State_lon, State_lat - integer, dimension(N_obs_param) :: select_species, select_species_Tb ! alloc max possible length + integer, dimension(N_obs_param) :: select_species, select_species_Tb, select_species_asnow ! alloc max possible length character(40), dimension(N_obs_param) :: select_varnames ! alloc max possible length @@ -3538,7 +3538,7 @@ subroutine cat_enkf_increments( & real, dimension(N_snow) :: tpsn, fice_snow_vec ! for snow model relayer real, dimension(N_snow,N_constit) :: rconstit - logical :: found_Tb_obs, unique, has_asnow + logical :: found_Tb_obs, unique, has_asnow, has_Tb, has_sfmc, has_sfds ! ----------------------------------------------------------------------- @@ -3564,6 +3564,8 @@ subroutine cat_enkf_increments( & select_varnames = '' select_species = -8888 ! intentionally differs from init in get_select_species() + select_species_Tb = -8888 ! intentionally differs from init in get_select_species() + select_species_asnow = -8888 ! intentionally differs from init in get_select_species() ! ---------------------------------------------------------------------- ! @@ -4688,34 +4690,67 @@ subroutine cat_enkf_increments( & ! sfcm/sfds & Tb | peat | 7 | srfexc, rzexc, catdef, tc[x], ght(1) ! ! amfox+rreichle, 26 Feb 2024 + + ! Figure out what types of observtion we have - - N_select_varnames = 0 - - do ii = 1, N_obs_param - unique = .true. - do jj = 1, N_select_varnames - if (trim(obs_param(ii)%varname) == trim(select_varnames(jj))) then - unique = .false. - exit - end if - end do - if (unique) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = trim(obs_param(ii)%varname) - 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' + has_Tb = .true. + 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' + has_sfmc = .true. + exit + end if end do - ! Check if we potentially have snow fraction observations by looking for 'asnow' in select_varnames + 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' + has_sfds = .true. + exit + end if + end do - do ii = 1, N_select_varnames - if (trim(select_varnames(ii)) == 'asnow') then + do ii = 1,N_obs_param + if (trim(obs_param(ii)%varname) == 'asnow') then has_asnow = .true. - exit - end if - end do + exit + end if + end do - ! If we have snow fraction observations, we need to calculate the increments for the snow layers + ! If we have any of the soil moisture observtions we need to get the species ID + + if (N_select_varnames > 0) then ! We have soil moisture observations + + if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' + + call get_select_species( & + N_select_varnames, select_varnames(1:N_select_varnames), & + N_obs_param, obs_param, N_select_species, select_species ) + + ! If we have Tb observations we need to get the species ID + if (has_Tb) then + call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb ) + end if + + N_state_max = 7 + + allocate( State_incr(N_state_max,N_ens)) + allocate( State_lon( N_state_max )) + allocate( State_lat( N_state_max )) + + end if + + ! If we have snow observations we need to get the species ID if (has_asnow) then @@ -4725,40 +4760,34 @@ subroutine cat_enkf_increments( & ! at SWE=WEMIN, the tile is fully snow covered (asnow=1) if (SCF_ANA_MAXINCRSWE>WEMIN) call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'must use SCF_ANA_MAXINCRSWE<=WEMIN') - - ! identify the obs species of interest - - N_select_varnames = 1 - - select_varnames(1) = 'asnow' - - call get_select_species( & - N_select_varnames, select_varnames(1:N_select_varnames), & - N_obs_param, obs_param, N_select_species, select_species ) - + allocate(select_tilenum(1)) swe_incr = 0. ! total SWE increment; initialize to NO CHANGE swe_incr_ensavg = 0. ! total SWE ensemble average increment; initialize to NO CHANGE - + + call get_select_species(1, 'asnow', N_obs_param, obs_param, N_select_species_asnow, select_species_asnow ) + end if ! loop through tiles and compute increments - do kk=1,N_catd - + do kk=1,N_catd + + if (has_asnow) then + ! find observations for tile kk select_tilenum(1) = l2f(kk) call get_ind_obs( & - N_obs, Observations, & - 1, select_tilenum, & - N_select_species, select_species(1:N_select_species), & - N_selected_obs, ind_obs ) + N_obs, Observations, & + 1, select_tilenum, & + N_select_species_asnow, select_species_asnow(1:N_select_species_asnow), & + N_selected_obs, ind_obs ) if (N_selected_obs > 0) then - + write(*,*) 'Doing snow DA' - + ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) tmp_obs = sum(Observations(ind_obs(1:N_selected_obs))%obs) @@ -4824,14 +4853,14 @@ subroutine cat_enkf_increments( & ! too little snow in forecast, use generic properties for added snow tmp_wesn(kk,n_e,isnow) = swe_ana / N_snow ! distribute SWE evenly across layers - + ! assign heat content for snow at 0 deg C and without liquid water content (100% frozen) ! (based on StieglitzSnow: htsn = (CPW*tsnow - fice*MAPL_ALHF)*swe ) - + tmp_htsn(kk,n_e,isnow) = (0.0 - MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) - + ! assign snow depth consistent with density of freshly fallen snow (must have SCF_ANA_MAXINCRSWE<=WEMIN) - + tmp_sndz(kk,n_e,isnow) = (WEMIN / RHOFS) / N_snow else @@ -4841,16 +4870,16 @@ subroutine cat_enkf_increments( & ! update SWE: tmp_wesn(kk,n_e,isnow) = cat_progn(kk,n_e)%wesn(isnow) * swe_ratio - + ! update snow heat content (keep snow temperature constant): - + call StieglitzSnow_calc_tpsnow( cat_progn(kk,n_e)%htsn(isnow), cat_progn(kk,n_e)%wesn(isnow), & - snow_temp, fice_snow, log_dum, log_dum2, .false. ) + snow_temp, fice_snow, log_dum, log_dum2, .false. ) tmp_htsn(kk,n_e,isnow) = (StieglitzSnow_CPW*snow_temp - fice_snow*MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) ! update snow depth: - + if (asnow_ana < 1. .and. asnow_fcst < 1.) then ! keep snow depth constant when less than full snow cover in fcst and ana @@ -4870,7 +4899,7 @@ subroutine cat_enkf_increments( & snow_dens = ( cat_progn(kk,n_e)%wesn(isnow)/asnow_fcst ) / cat_progn(kk,n_e)%sndz(isnow) ! ii) diagnose analysis snow depth using forecast density - + tmp_sndz(kk,n_e,isnow) = ( tmp_wesn(kk,n_e,isnow)/asnow_ana ) / snow_dens end if @@ -4882,11 +4911,11 @@ subroutine cat_enkf_increments( & ! 4. Relayer to balance the snow column (call with optional args for adjustment of htsnn) call StieglitzSnow_relayer( N_snow, N_constit, & - MAPL_LAND, CATCH_SNOW_DZPARAM, & - tmp_htsn(kk,n_e,1:N_snow), & - tmp_wesn(kk,n_e,1:N_snow), & - tmp_sndz(kk,n_e,1:N_snow), & - rconstit, tpsn, fice_snow_vec ) + MAPL_LAND, CATCH_SNOW_DZPARAM, & + tmp_htsn(kk,n_e,1:N_snow), & + tmp_wesn(kk,n_e,1:N_snow), & + tmp_sndz(kk,n_e,1:N_snow), & + rconstit, tpsn, fice_snow_vec ) ! print the old and new swe, heat content and snow density @@ -4909,217 +4938,173 @@ subroutine cat_enkf_increments( & end if ! if (N_selected_obs > 0) - end do ! kk=1,N_catd + end if ! if (has_asnow) - end if ! if (has_asnow) + ! Calculate swe_incr_ensavg as the average of the ensemble members to add to check if tile is snow-free - ! Calculate swe_incr_ensavg as the average of the ensemble members to add to check if tile is snow-free + swe_incr_ensavg(kk) = sum((swe_incr(kk,:))) / real(N_ens) - swe_incr_ensavg = sum(swe_incr, dim=2) / real(N_ens) + if (N_select_species > 0) then ! We have soil moisture observations - ! Will get all species associated with Tb or sfds observations - - 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 + 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. & + (swe_incr_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 - 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 + write(*,*) 'Doing SM DA' - 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 + ! Determine if Tb observations are present + + found_Tb_obs = .false. - if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' - - 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. & - (swe_incr_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 + 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 - write(*,*) 'Doing SM DA' + 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 - ! Determine if Tb observations are present - - found_Tb_obs = .false. + 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 - 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 ( N_state==2 ) then + + cat_progn_incr(kk,:)%srfexc = State_incr(1,:)*scale_srfexc + cat_progn_incr(kk,:)%rzexc = State_incr(2,:)*scale_rzexc - 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 + + 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 - 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)) + write (*,*) 'State_incr: ', State_incr(1,1) + + end if ! if (N_selected_obs > 0) - 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 + end if ! thresholds - 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 - - write (*,*) 'State_incr: ', State_incr(1,1) - - end if - - end if ! thresholds + end if ! if (N_select_species > 0) - end do + end do ! kk=1,N_catd ! ---------------------------------- From d3546c8e6666b5663fac71917aae45818fface69 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 8 Nov 2024 15:39:33 -0700 Subject: [PATCH 04/21] rename to update type 14 --- .../clsm_ensupd_enkf_update.F90 | 66 ++++-- .../clsm_ensupd_upd_routines.F90 | 218 +++++++++++++++++- 2 files changed, 267 insertions(+), 17 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index c7fce5a..17b598a 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1395,22 +1395,7 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn(n,n_e)%tc4 + cat_progn_incr(n,n_e)%tc4 cat_progn(n,n_e)%ght(1) = & - cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1) - - do ii=1,N_snow ! for each snow layer - - cat_progn(n,n_e)%wesn(ii) = & - cat_progn(n,n_e)%wesn(ii) + cat_progn_incr(n,n_e)%wesn(ii) - - cat_progn(n,n_e)%sndz(ii) = & - cat_progn(n,n_e)%sndz(ii) + cat_progn_incr(n,n_e)%sndz(ii) - - cat_progn(n,n_e)%htsn(ii) = & - cat_progn(n,n_e)%htsn(ii) + cat_progn_incr(n,n_e)%htsn(ii) - - end do - - + cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1) end do end do @@ -1441,6 +1426,55 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .true. + case (14) select_update_type ! soil moisture, temperature and snow cover 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) + + if (logit) write (logunit,*) & + 'apply_enkf_increments(): applying soil moisture and Tskin/ght1 increments' + + do n=1,N_catd + do n_e=1,N_ens + + cat_progn(n,n_e)%srfexc = & + cat_progn(n,n_e)%srfexc + cat_progn_incr(n,n_e)%srfexc + cat_progn(n,n_e)%rzexc = & + cat_progn(n,n_e)%rzexc + cat_progn_incr(n,n_e)%rzexc + cat_progn(n,n_e)%catdef = & + cat_progn(n,n_e)%catdef + cat_progn_incr(n,n_e)%catdef + + cat_progn(n,n_e)%tc1 = & + cat_progn(n,n_e)%tc1 + cat_progn_incr(n,n_e)%tc1 + cat_progn(n,n_e)%tc2 = & + cat_progn(n,n_e)%tc2 + cat_progn_incr(n,n_e)%tc2 + cat_progn(n,n_e)%tc4 = & + cat_progn(n,n_e)%tc4 + cat_progn_incr(n,n_e)%tc4 + + cat_progn(n,n_e)%ght(1) = & + cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1) + + do ii=1,N_snow ! for each snow layer + + cat_progn(n,n_e)%wesn(ii) = & + cat_progn(n,n_e)%wesn(ii) + cat_progn_incr(n,n_e)%wesn(ii) + + cat_progn(n,n_e)%sndz(ii) = & + cat_progn(n,n_e)%sndz(ii) + cat_progn_incr(n,n_e)%sndz(ii) + + cat_progn(n,n_e)%htsn(ii) = & + cat_progn(n,n_e)%htsn(ii) + cat_progn_incr(n,n_e)%htsn(ii) + + end do + + end do + end do + + cat_progn_has_changed = .true. + + check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case + case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown update_type') diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index dbb144e..ed0293b 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4676,7 +4676,223 @@ subroutine cat_enkf_increments( & ! ---------------------------------- - case (13) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs + 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(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 (14) 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 ! From 7d8c26b82127dc9e213c4a191fb8a246e16f741e Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 8 Nov 2024 19:38:15 -0500 Subject: [PATCH 05/21] add update type to check_compact_support --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index ed0293b..331fd83 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -5894,7 +5894,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,14) ! "3d" updates, check consistency of xcompact, ycompact ! check xcompact/ycompact against corr scales of model error From 1a678a386f2abaa9fec7d534d46c031b9ffaa6b1 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Nov 2024 15:00:17 -0700 Subject: [PATCH 06/21] minor edits --- .../clsm_ensupd_upd_routines.F90 | 60 ++++++++++--------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index ed0293b..8a40786 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4893,9 +4893,20 @@ subroutine cat_enkf_increments( & ! ---------------------------------- case (14) 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 + ! this case is a combination of cases 11 and 13 + ! first, obs_param%varname is checked to determine if snow and/or soil moisture obs are potentially in the analysis + ! then, if in the analysis the species ID is determined for each type of soil moisture observation + ! if snow observations are in the analysis, the species ID is determined and the snow analysis is performed + ! finally, if soil moisture observations are in the analysis, the soil moisture analysis is performed + + ! in the case of the snow analysis: + ! update each tile separately using all observations associated with that tile + ! loops through ensemble members and compute analysis separately for each ensemble member + + ! in the case of the 3d soil moisture/Tskin/ght(1) analysis: ! 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 @@ -4941,13 +4952,11 @@ subroutine cat_enkf_increments( & has_asnow = .true. exit end if - end do + end do ! If we have any of the soil moisture observtions we need to get the species ID - if (N_select_varnames > 0) then ! We have soil moisture observations - - if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' + if (has_Tb .or. has_sfmc .or. has_sfds) then ! We have soil moisture observations call get_select_species( & N_select_varnames, select_varnames(1:N_select_varnames), & @@ -4982,27 +4991,28 @@ subroutine cat_enkf_increments( & swe_incr = 0. ! total SWE increment; initialize to NO CHANGE swe_incr_ensavg = 0. ! total SWE ensemble average increment; initialize to NO CHANGE - call get_select_species(1, 'asnow', N_obs_param, obs_param, N_select_species_asnow, select_species_asnow ) + call get_select_species(1, 'asnow', N_obs_param, obs_param, N_select_species_asnow, select_species_asnow ) + end if - ! loop through tiles and compute increments + + ! loop through tiles and compute increments do kk=1,N_catd if (has_asnow) then - ! find observations for tile kk + ! find snow observations for tile kk select_tilenum(1) = l2f(kk) - call get_ind_obs( & - N_obs, Observations, & - 1, select_tilenum, & - N_select_species_asnow, select_species_asnow(1:N_select_species_asnow), & + call get_ind_obs( & + N_obs, Observations, & + 1, select_tilenum, & + N_select_species_asnow, & + select_species_asnow(1:N_select_species_asnow), & N_selected_obs, ind_obs ) if (N_selected_obs > 0) then - - write(*,*) 'Doing snow DA' ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) @@ -5041,8 +5051,6 @@ subroutine cat_enkf_increments( & ! 3. Derive SWE, snow heat content, and snow depth increments for each layer from total SWE increment swe_ana = max(swe_fcst + swe_incr(kk,n_e), 0.0) ! total SWE after analysis - - write (*,*) 'swe_incr(kk,n_e): ',swe_incr(kk,n_e) call StieglitzSnow_calc_asnow( swe_ana, asnow_ana ) ! asnow after analysis @@ -5160,13 +5168,15 @@ subroutine cat_enkf_increments( & swe_incr_ensavg(kk) = sum((swe_incr(kk,:))) / real(N_ens) - if (N_select_species > 0) then ! We have soil moisture observations + if (has_Tb .or. has_sfmc .or. has_sfds) then ! Have soil moisture observations + + if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' 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. & + if ( (SWE_ensavg(kk) < SWE_threshold) .and. & (swe_incr_ensavg(kk) < SWE_threshold) .and. & (tp1_ensavg(kk) > tp1_threshold) ) then @@ -5185,16 +5195,14 @@ subroutine cat_enkf_increments( & halo_minlat = max(halo_minlat, -90.) halo_maxlat = min(halo_maxlat, 90.) - call get_ind_obs_lat_lon_box( & - N_obs, Observations, & + 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 - write(*,*) 'Doing SM DA' - ! Determine if Tb observations are present found_Tb_obs = .false. @@ -5262,7 +5270,7 @@ subroutine cat_enkf_increments( & call assemble_obs_cov( N_selected_obs, N_obs_param, obs_param, & Observations(ind_obs(1:N_selected_obs)), Obs_cov ) - call enkf_increments( & + 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),:), & @@ -5311,14 +5319,12 @@ subroutine cat_enkf_increments( & cat_progn_incr(kk,:)%ght(1) = State_incr(7,:)*scale_ght1 end if - - write (*,*) 'State_incr: ', State_incr(1,1) - end if ! if (N_selected_obs > 0) + end if ! if (N_selected_obs > 0) end if ! thresholds - end if ! if (N_select_species > 0) + end if ! if (have soil moisture observations) end do ! kk=1,N_catd From d94d405030fbf51a4aebf374899516f8ee16d8d7 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Nov 2024 16:24:54 -0700 Subject: [PATCH 07/21] edit comments --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 3418ffa..012305d 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4896,10 +4896,12 @@ subroutine cat_enkf_increments( & ! 1d snow analysis (Toure et al. 2018 empirical gain); snow cover fraction obs ! this case is a combination of cases 11 and 13 - ! first, obs_param%varname is checked to determine if snow and/or soil moisture obs are potentially in the analysis - ! then, if in the analysis the species ID is determined for each type of soil moisture observation - ! if snow observations are in the analysis, the species ID is determined and the snow analysis is performed - ! finally, if soil moisture observations are in the analysis, the soil moisture analysis is performed + ! 1. obs_param%varname is checked to determine if snow and/or soil moisture obs are potentially in the analysis + ! 2. if sm obs are in the analysis the species ID is determined for each type + ! 3. if snow observations are in the analysis, their species ID is determined + ! 4. start looping through all the tiles + ! 5. if snow obs are in the analysis, check for snow obs with each tile, and if found perform the snow analysis + ! 6. if sm obs are in the analysis, check for SM obs around each tile, and if found perform sm analysis ! in the case of the snow analysis: ! update each tile separately using all observations associated with that tile From aa4c6013ebf5ba1a954996f76f73f10082a9f427 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Nov 2024 17:19:47 -0700 Subject: [PATCH 08/21] update log message when adding increments --- GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 17b598a..1929ccd 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1433,7 +1433,7 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & ! 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' + 'apply_enkf_increments(): applying soil moisture, Tskin/ght1 and snow increments' do n=1,N_catd do n_e=1,N_ens From fdf417fa1e7c95938cecf4df6138be40fe910735 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Nov 2024 17:25:23 -0700 Subject: [PATCH 09/21] try to fix white space change --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 012305d..10a2a32 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4675,7 +4675,7 @@ subroutine cat_enkf_increments( & 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 From 321f5e9e3353627748d0fbb0982268bf8ab1e661 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Nov 2024 18:07:55 -0700 Subject: [PATCH 10/21] try to fix to enable easy Github diff --- .../clsm_ensupd_upd_routines.F90 | 436 +++++++++--------- 1 file changed, 218 insertions(+), 218 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 10a2a32..12b6a8a 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3524,7 +3524,7 @@ subroutine cat_enkf_increments( & real, dimension( N_catd) :: SWE_ensavg real, dimension( N_catd) :: tp1_ensavg real, dimension( N_catd) :: asnow_ensavg - real, dimension( N_catd) :: swe_incr_ensavg + real, dimension( N_catd) :: swe_incr_ensavg type(obs_param_type) :: this_obs_param @@ -3538,7 +3538,7 @@ subroutine cat_enkf_increments( & real, dimension(N_snow) :: tpsn, fice_snow_vec ! for snow model relayer real, dimension(N_snow,N_constit) :: rconstit - logical :: found_Tb_obs, unique, has_asnow, has_Tb, has_sfmc, has_sfds + logical :: found_Tb_obs, has_asnow, has_Tb, has_sfmc, has_sfds ! ----------------------------------------------------------------------- @@ -3564,8 +3564,8 @@ subroutine cat_enkf_increments( & select_varnames = '' select_species = -8888 ! intentionally differs from init in get_select_species() - select_species_Tb = -8888 ! intentionally differs from init in get_select_species() - select_species_asnow = -8888 ! intentionally differs from init in get_select_species() + select_species_Tb = -8888 + select_species_asnow = -8888 ! ---------------------------------------------------------------------- ! @@ -4675,221 +4675,221 @@ subroutine cat_enkf_increments( & 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 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 )) + + 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 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 + 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(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)) - 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==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 (14) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs @@ -5331,7 +5331,7 @@ subroutine cat_enkf_increments( & end do ! kk=1,N_catd ! ---------------------------------- - + case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown update_type') From 541508919d37107bd2c749b9dec0a4fe13713515 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 26 Nov 2024 16:51:50 -0500 Subject: [PATCH 11/21] add obs_param%assim logical --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 12b6a8a..5104e6f 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4918,12 +4918,12 @@ subroutine cat_enkf_increments( & ! 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 + ! amfox+rreichle, Dec 2024 ! Figure out what types of observtion we have do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'Tb') then + if (trim(obs_param(ii)%varname) == 'Tb' .and. obs_param(ii)%assim) then N_select_varnames = N_select_varnames + 1 select_varnames(N_select_varnames) = 'Tb' has_Tb = .true. @@ -4932,7 +4932,7 @@ subroutine cat_enkf_increments( & end do do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'sfmc') then + if (trim(obs_param(ii)%varname) == 'sfmc' .and. obs_param(ii)%assim) then N_select_varnames = N_select_varnames + 1 select_varnames(N_select_varnames) = 'sfmc' has_sfmc = .true. @@ -4941,7 +4941,7 @@ subroutine cat_enkf_increments( & end do do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'sfds') then + if (trim(obs_param(ii)%varname) == 'sfds' .and. obs_param(ii)%assim) then N_select_varnames = N_select_varnames + 1 select_varnames(N_select_varnames) = 'sfds' has_sfds = .true. @@ -4950,7 +4950,7 @@ subroutine cat_enkf_increments( & end do do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'asnow') then + if (trim(obs_param(ii)%varname) == 'asnow' .and. obs_param(ii)%assim) then has_asnow = .true. exit end if From 1f7f6a6bdb0419f40ba26f4fc57d0be9252eb1c1 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 26 Nov 2024 16:01:36 -0700 Subject: [PATCH 12/21] add update types 11 and 14 description to default nml --- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index ece206c..53e428f 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -34,7 +34,10 @@ ! 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: 1d snow analysis (Toure et al. 2018 empirical gain); SCF obs ! update_type = 13: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile; sfmc and TB obs +! update_type = 14: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile & +! 1d snow analysis (Toure et al. 2018 empirical gain); sfmc, TB and SCF obs update_type = 0 From 13e5f63062aa72fd38117551f1c9dda47627a0b6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 29 Nov 2024 14:37:08 -0500 Subject: [PATCH 13/21] edited new update_type for soil moisture, soil temperature, and snow analysis: - changed number associated with update_type 14 --> 12 - modified "check_snow" (maintains 0-diff for update_type=10 only!) - simplified how species of interest are obtained - added check of consistency between update_type and obs_param%assim - additional clean-up (log messages, indentation) --- .../clsm_ensupd_enkf_update.F90 | 54 +- .../clsm_ensupd_upd_routines.F90 | 890 +++++++++--------- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 35 +- 3 files changed, 489 insertions(+), 490 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 1929ccd..65e762a 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1339,6 +1339,8 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .false. + ! ------------------------------------------------------------------ + case (4,7) select_update_type ! Tskin/ght1 update if (logit) write (logunit,*) & @@ -1361,6 +1363,8 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .true. + ! ------------------------------------------------------------------ + case (5) select_update_type ! Tskin/ght1 update if (logit) write (logunit,*) & @@ -1368,6 +1372,8 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .false. + ! ------------------------------------------------------------------ + 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 @@ -1401,13 +1407,21 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .true. - check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case + if (select_update_type==10) then + + check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case - case(11) select_update_type ! empirical MODIS SCF update + end if + + ! ------------------------------------------------------------------ + + case(11) select_update_type ! snow update - do n=1,N_catd ! for each tile + if (logit) write (logunit,*) 'apply_enkf_increments(): applying snow increments' + + do n=1,N_catd ! for each tile - do n_e=1,N_ens ! for each ensemble member + do n_e=1,N_ens ! for each ensemble member do ii=1,N_snow ! for each snow layer @@ -1426,7 +1440,9 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .true. - case (14) select_update_type ! soil moisture, temperature and snow cover update + ! ------------------------------------------------------------------ + + case (12) select_update_type ! soil moisture, temperature, and snow 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; @@ -1455,26 +1471,26 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn(n,n_e)%ght(1) = & cat_progn(n,n_e)%ght(1) + cat_progn_incr(n,n_e)%ght(1) - do ii=1,N_snow ! for each snow layer - - cat_progn(n,n_e)%wesn(ii) = & - cat_progn(n,n_e)%wesn(ii) + cat_progn_incr(n,n_e)%wesn(ii) - - cat_progn(n,n_e)%sndz(ii) = & - cat_progn(n,n_e)%sndz(ii) + cat_progn_incr(n,n_e)%sndz(ii) - - cat_progn(n,n_e)%htsn(ii) = & - cat_progn(n,n_e)%htsn(ii) + cat_progn_incr(n,n_e)%htsn(ii) - - end do + do ii=1,N_snow ! for each snow layer + + cat_progn(n,n_e)%wesn(ii) = & + cat_progn(n,n_e)%wesn(ii) + cat_progn_incr(n,n_e)%wesn(ii) + + cat_progn(n,n_e)%sndz(ii) = & + cat_progn(n,n_e)%sndz(ii) + cat_progn_incr(n,n_e)%sndz(ii) + cat_progn(n,n_e)%htsn(ii) = & + cat_progn(n,n_e)%htsn(ii) + cat_progn_incr(n,n_e)%htsn(ii) + + end do + end do end do cat_progn_has_changed = .true. - check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case - + ! ------------------------------------------------------------------ + case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown update_type') diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 5104e6f..95bad86 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3478,7 +3478,8 @@ subroutine cat_enkf_increments( & integer :: n, n_e, kk, ii, jj - integer :: N_state_max, N_state, N_selected_obs, N_select_varnames, N_select_species, N_select_species_Tb, N_select_species_asnow + integer :: N_state_max, N_state, N_selected_obs, N_select_varnames + integer :: N_select_species, N_select_species_smTb, N_select_species_Tb, N_select_species_asnow real :: halo_minlon, halo_maxlon, halo_minlat, halo_maxlat real :: tmp_minlon, tmp_maxlon, tmp_minlat, tmp_maxlat @@ -3495,7 +3496,8 @@ subroutine cat_enkf_increments( & real, allocatable, dimension(:) :: State_lon, State_lat - integer, dimension(N_obs_param) :: select_species, select_species_Tb, select_species_asnow ! alloc max possible length + integer, dimension(N_obs_param) :: select_species ! alloc max possible length + integer, dimension(N_obs_param) :: select_species_smTb, select_species_Tb, select_species_asnow ! alloc max possible length character(40), dimension(N_obs_param) :: select_varnames ! alloc max possible length @@ -3538,7 +3540,7 @@ subroutine cat_enkf_increments( & real, dimension(N_snow) :: tpsn, fice_snow_vec ! for snow model relayer real, dimension(N_snow,N_constit) :: rconstit - logical :: found_Tb_obs, has_asnow, has_Tb, has_sfmc, has_sfds + logical :: found_Tb_obs ! ----------------------------------------------------------------------- @@ -3559,13 +3561,19 @@ subroutine cat_enkf_increments( & ! more initializations - N_select_varnames = 0 - N_select_species = 0 + N_select_varnames = 0 - select_varnames = '' - select_species = -8888 ! intentionally differs from init in get_select_species() - select_species_Tb = -8888 - select_species_asnow = -8888 + N_select_species = 0 + N_select_species_smTb = 0 + N_select_species_Tb = 0 + N_select_species_asnow = 0 + + select_varnames = '' + + select_species = -8888 ! intentionally differs from init in get_select_species() + select_species_smTb = -8888 + select_species_Tb = -8888 + select_species_asnow = -8888 ! ---------------------------------------------------------------------- ! @@ -3740,6 +3748,8 @@ subroutine cat_enkf_increments( & end do + ! ---------------------------------------------------------------------------------------------------------------------- + case (2) select_update_type ! 3d soil moisture analysis; sfmc+sfds obs ! update each tile separately using all observations within @@ -3829,7 +3839,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (3) select_update_type ! 1d Tskin analysis; tskin obs @@ -3900,7 +3910,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (4,5) select_update_type ! 1d Tskin/ght(1) analysis; tskin obs @@ -3974,7 +3984,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (6) select_update_type ! 1d soil moisture/Tskin/ght(1) analysis; Tb obs @@ -4053,7 +4063,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (7) select_update_type ! 3d Tskin/ght(1) analysis; tskin obs @@ -4144,7 +4154,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (8,10) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb obs @@ -4302,7 +4312,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (9) select_update_type ! 1d Tskin/ght(1) analysis; FT obs @@ -4483,7 +4493,7 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (11) select_update_type ! 1d snow analysis (Toure et al. 2018 empirical gain); snow cover fraction obs @@ -4674,7 +4684,7 @@ subroutine cat_enkf_increments( & end do ! kk=1,N_catd - ! ---------------------------------- + ! ---------------------------------------------------------------------------------------------------------------------- case (13) select_update_type ! 3d soil moisture/Tskin/ght(1) analysis; Tb+sfmc+sfds obs @@ -4691,35 +4701,15 @@ subroutine cat_enkf_increments( & ! ! amfox+rreichle, 26 Feb 2024 - if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' + if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc+sfds 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 + ! Get all species associated with *assimilated* Tb, sfmc, and sfds observations - 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 + N_select_varnames = 3 - ! Will get all species associated with Tb or sfds observations + 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), & @@ -4890,389 +4880,374 @@ subroutine cat_enkf_increments( & end do - ! ---------------------------------- - - case (14) 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 + ! ---------------------------------------------------------------------------------------------------------------------- - ! this case is a combination of cases 11 and 13 - ! 1. obs_param%varname is checked to determine if snow and/or soil moisture obs are potentially in the analysis - ! 2. if sm obs are in the analysis the species ID is determined for each type - ! 3. if snow observations are in the analysis, their species ID is determined - ! 4. start looping through all the tiles - ! 5. if snow obs are in the analysis, check for snow obs with each tile, and if found perform the snow analysis - ! 6. if sm obs are in the analysis, check for SM obs around each tile, and if found perform sm analysis + 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 - ! in the case of the snow analysis: - ! update each tile separately using all observations associated with that tile - ! loops through ensemble members and compute analysis separately for each ensemble member - - ! in the case of the 3d soil moisture/Tskin/ght(1) analysis: - ! 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 + ! update_type 12 is a combination of update_type 11 and update_type 13: ! - ! 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) + ! Loop through tiles, first apply 1d snow update, then apply 3d soil moisture update. + ! + ! 1d snow analysis: + ! Update each tile separately using all observations associated with the tile. + ! Loop through ensemble members and compute analysis separately for each ensemble member. + ! + ! 3d soil moisture/Tskin/ght(1) analysis: + ! Update each tile separately using all observations within customized halo around the tile. + ! State vector of each tile depends 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, Dec 2024 + ! + ! ------------------------------------------------------------------------------------------------------- - ! Figure out what types of observtion we have + if (logit) write (logunit, *) 'get increments: update_type = ', select_update_type - do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'Tb' .and. obs_param(ii)%assim) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'Tb' - has_Tb = .true. - exit - end if - end do - - do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'sfmc' .and. obs_param(ii)%assim) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'sfmc' - has_sfmc = .true. - exit - end if - end do - - do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'sfds' .and. obs_param(ii)%assim) then - N_select_varnames = N_select_varnames + 1 - select_varnames(N_select_varnames) = 'sfds' - has_sfds = .true. - exit - end if - end do - - do ii = 1,N_obs_param - if (trim(obs_param(ii)%varname) == 'asnow' .and. obs_param(ii)%assim) then - has_asnow = .true. - exit - end if - end do - - ! If we have any of the soil moisture observtions we need to get the species ID - - if (has_Tb .or. has_sfmc .or. has_sfds) then ! We have soil moisture observations - - call get_select_species( & - N_select_varnames, select_varnames(1:N_select_varnames), & - N_obs_param, obs_param, N_select_species, select_species ) - - ! If we have Tb observations we need to get the species ID - if (has_Tb) then - call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb ) - end if - - N_state_max = 7 - - allocate( State_incr(N_state_max,N_ens)) - allocate( State_lon( N_state_max )) - allocate( State_lat( N_state_max )) - - end if - - ! If we have snow observations we need to get the species ID - - if (has_asnow) then - - if (logit) write (logunit, *) 'get 1d snow increments (Toure et al. 2018 empirical gain); snow cover fraction obs' + ! determine species of assimilated obs associated with snow analysis - ! ensure that max SWE increment parameter is less than WEMIN; larger increments make no sense because - ! at SWE=WEMIN, the tile is fully snow covered (asnow=1) - - if (SCF_ANA_MAXINCRSWE>WEMIN) call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'must use SCF_ANA_MAXINCRSWE<=WEMIN') - - allocate(select_tilenum(1)) + call get_select_species(1, 'asnow', N_obs_param, obs_param, N_select_species_asnow, select_species_asnow ) + + if (N_select_species_asnow>0) then + + if (logit) write (logunit, *) '- get 1d snow increments (Toure et al. 2018 empirical gain); snow cover fraction obs' + + ! ensure that max SWE increment parameter is less than WEMIN; larger increments make no sense because + ! at SWE=WEMIN, the tile is fully snow covered (asnow=1) + + if (SCF_ANA_MAXINCRSWE>WEMIN) call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'must use SCF_ANA_MAXINCRSWE<=WEMIN') + + allocate(select_tilenum(1)) + + swe_incr = 0. ! total SWE increment; initialize to NO CHANGE + swe_incr_ensavg = 0. ! total SWE ensemble average increment; initialize to NO CHANGE - swe_incr = 0. ! total SWE increment; initialize to NO CHANGE - swe_incr_ensavg = 0. ! total SWE ensemble average increment; initialize to NO CHANGE - - call get_select_species(1, 'asnow', N_obs_param, obs_param, N_select_species_asnow, select_species_asnow ) - - end if + + end if + + ! determine species of assimilated obs associated with soil moisture (+Tskin/ght(1)) analysis: + + 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_smTb, select_species_smTb ) + + ! determine which species are Tb, if any + + call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb ) + + if (N_select_species_smTb>0) then + + if (logit) write (logunit, *) '- get 1d snow increments (Toure et al. 2018 empirical gain); snow cover fraction obs' + + end if - ! loop through tiles and compute increments - - do kk=1,N_catd + ! check consistency of update_type and obs_param%assim config + + if ( N_select_species_asnow==0 .and. N_select_species_smTb==0 ) then + + call ldas_abort(LDAS_GENERIC_ERROR, Iam, & + 'update_type inconsistent with obs_param%assim==.false. for all asnow/Tb/sfmc/sfds species') + + end if + + ! allocate State_* - if (has_asnow) then - - ! find snow observations for tile kk - - select_tilenum(1) = l2f(kk) - - call get_ind_obs( & + N_state_max = 7 + + allocate( State_incr(N_state_max,N_ens)) + allocate( State_lon( N_state_max )) + allocate( State_lat( N_state_max )) + + ! ---------------------------------------------------- + ! + ! loop through tiles and compute increments + + do kk=1,N_catd + + if (N_select_species_asnow>0) then ! assimilate snow observations + + ! find snow observations for tile kk + + select_tilenum(1) = l2f(kk) + + call get_ind_obs( & N_obs, Observations, & 1, select_tilenum, & N_select_species_asnow, & select_species_asnow(1:N_select_species_asnow), & N_selected_obs, ind_obs ) - - if (N_selected_obs > 0) then - - ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) - - tmp_obs = sum(Observations(ind_obs(1:N_selected_obs))%obs) - - if (N_selected_obs > 1) tmp_obs = tmp_obs/real(N_selected_obs) - - do n_e=1,N_ens ! compute analysis separately for each ensemble member - - ! 1. Diagnose model forecast snow cover area fraction and total SWE - - asnow_fcst = asnow(kk,n_e) - swe_fcst = sum(cat_progn(kk,n_e)%wesn(1:N_snow)) - - ! 2. Calculate SWE increment based on modified eq 1 of Toure et al (2018) - - if (asnow_fcst .lt. tmp_obs * SCF_ANA_ALPHA) then - - ! ADD SNOW: Forecast SCF is less than observed SCF (after "bias" adjustment with alpha) - - swe_incr(kk,n_e) = SCF_ANA_MAXINCRSWE * (tmp_obs - asnow_fcst/SCF_ANA_ALPHA) - - elseif (tmp_obs .lt. SCF_ANA_BETA) then - - ! REMOVE SNOW: Simulated SCF is greater than observed SCF (after "bias" adjustment) - ! and observed SCF is less than beta threshold - - swe_incr(kk,n_e) = (-1.) * SCF_ANA_MAXINCRSWE * asnow_fcst * (1. - tmp_obs/SCF_ANA_BETA) - - else - - cycle ! NO CHANGE, skip rest of increment calcs and go straight to next ens member - - endif ! (Toure et al. 2018 Equation 1) - - ! 3. Derive SWE, snow heat content, and snow depth increments for each layer from total SWE increment - - swe_ana = max(swe_fcst + swe_incr(kk,n_e), 0.0) ! total SWE after analysis - - call StieglitzSnow_calc_asnow( swe_ana, asnow_ana ) ! asnow after analysis - - if (swe_fcst>=StieglitzSnow_MINSWE) then - swe_ratio = swe_ana / swe_fcst - else - swe_ratio = MAPL_UNDEF ! swe_ratio unreliable; set to MAPL_UNDEF to expose inadvertent use - end if - - ! loop through snow layers and compute SWE, snow heat content, and snow depth analysis for each layer - - do isnow=1,N_snow - - if (asnow_ana == 0.0) then - - ! no snow in analysis, remove all snow - - tmp_wesn(kk,n_e,isnow) = 0.0 - tmp_htsn(kk,n_e,isnow) = 0.0 - tmp_sndz(kk,n_e,isnow) = 0.0 - - elseif (swe_fcst < StieglitzSnow_MINSWE) then - - ! too little snow in forecast, use generic properties for added snow - - tmp_wesn(kk,n_e,isnow) = swe_ana / N_snow ! distribute SWE evenly across layers - - ! assign heat content for snow at 0 deg C and without liquid water content (100% frozen) - ! (based on StieglitzSnow: htsn = (CPW*tsnow - fice*MAPL_ALHF)*swe ) - - tmp_htsn(kk,n_e,isnow) = (0.0 - MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) - - ! assign snow depth consistent with density of freshly fallen snow (must have SCF_ANA_MAXINCRSWE<=WEMIN) - - tmp_sndz(kk,n_e,isnow) = (WEMIN / RHOFS) / N_snow - - else - - ! snow in forecast and analysis, derive properties of analysis snow from properties of forecast snow - - ! update SWE: - - tmp_wesn(kk,n_e,isnow) = cat_progn(kk,n_e)%wesn(isnow) * swe_ratio - - ! update snow heat content (keep snow temperature constant): - - call StieglitzSnow_calc_tpsnow( cat_progn(kk,n_e)%htsn(isnow), cat_progn(kk,n_e)%wesn(isnow), & + + if (N_selected_obs > 0) then + + ! average in case there are multiple "asnow" obs (e.g., from MODIS and VIIRS) + + tmp_obs = sum(Observations(ind_obs(1:N_selected_obs))%obs) + + if (N_selected_obs > 1) tmp_obs = tmp_obs/real(N_selected_obs) + + do n_e=1,N_ens ! compute analysis separately for each ensemble member + + ! 1. Diagnose model forecast snow cover area fraction and total SWE + + asnow_fcst = asnow(kk,n_e) + swe_fcst = sum(cat_progn(kk,n_e)%wesn(1:N_snow)) + + ! 2. Calculate SWE increment based on modified eq 1 of Toure et al (2018) + + if (asnow_fcst .lt. tmp_obs * SCF_ANA_ALPHA) then + + ! ADD SNOW: Forecast SCF is less than observed SCF (after "bias" adjustment with alpha) + + swe_incr(kk,n_e) = SCF_ANA_MAXINCRSWE * (tmp_obs - asnow_fcst/SCF_ANA_ALPHA) + + elseif (tmp_obs .lt. SCF_ANA_BETA) then + + ! REMOVE SNOW: Simulated SCF is greater than observed SCF (after "bias" adjustment) + ! and observed SCF is less than beta threshold + + swe_incr(kk,n_e) = (-1.) * SCF_ANA_MAXINCRSWE * asnow_fcst * (1. - tmp_obs/SCF_ANA_BETA) + + else + + cycle ! NO CHANGE, skip rest of increment calcs and go straight to next ens member + + endif ! (Toure et al. 2018 Equation 1) + + ! 3. Derive SWE, snow heat content, and snow depth increments for each layer from total SWE increment + + swe_ana = max(swe_fcst + swe_incr(kk,n_e), 0.0) ! total SWE after analysis + + call StieglitzSnow_calc_asnow( swe_ana, asnow_ana ) ! asnow after analysis + + if (swe_fcst>=StieglitzSnow_MINSWE) then + swe_ratio = swe_ana / swe_fcst + else + swe_ratio = MAPL_UNDEF ! swe_ratio unreliable; set to MAPL_UNDEF to expose inadvertent use + end if + + ! loop through snow layers and compute SWE, snow heat content, and snow depth analysis for each layer + + do isnow=1,N_snow + + if (asnow_ana == 0.0) then + + ! no snow in analysis, remove all snow + + tmp_wesn(kk,n_e,isnow) = 0.0 + tmp_htsn(kk,n_e,isnow) = 0.0 + tmp_sndz(kk,n_e,isnow) = 0.0 + + elseif (swe_fcst < StieglitzSnow_MINSWE) then + + ! too little snow in forecast, use generic properties for added snow + + tmp_wesn(kk,n_e,isnow) = swe_ana / N_snow ! distribute SWE evenly across layers + + ! assign heat content for snow at 0 deg C and without liquid water content (100% frozen) + ! (based on StieglitzSnow: htsn = (CPW*tsnow - fice*MAPL_ALHF)*swe ) + + tmp_htsn(kk,n_e,isnow) = (0.0 - MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) + + ! assign snow depth consistent with density of freshly fallen snow (must have SCF_ANA_MAXINCRSWE<=WEMIN) + + tmp_sndz(kk,n_e,isnow) = (WEMIN / RHOFS) / N_snow + + else + + ! snow in forecast and analysis, derive properties of analysis snow from properties of forecast snow + + ! update SWE: + + tmp_wesn(kk,n_e,isnow) = cat_progn(kk,n_e)%wesn(isnow) * swe_ratio + + ! update snow heat content (keep snow temperature constant): + + call StieglitzSnow_calc_tpsnow( cat_progn(kk,n_e)%htsn(isnow), cat_progn(kk,n_e)%wesn(isnow), & snow_temp, fice_snow, log_dum, log_dum2, .false. ) - - tmp_htsn(kk,n_e,isnow) = (StieglitzSnow_CPW*snow_temp - fice_snow*MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) - - ! update snow depth: - - if (asnow_ana < 1. .and. asnow_fcst < 1.) then - - ! keep snow depth constant when less than full snow cover in fcst and ana - - tmp_sndz(kk,n_e,isnow) = cat_progn(kk,n_e)%sndz(isnow) - - else - - ! compute analysis snow depth by keeping snow density constant - ! - ! in this case, it is possible that either asnow_fcst<1 or asnow_ana<1; - ! when computing density or depth, make sure that SWE value (which is per unit area) is - ! adjusted to reflect SWE value (per unit area) in the snow-covered fraction of the tile - - ! i) diagnose (layer-specific) forecast snow density - - snow_dens = ( cat_progn(kk,n_e)%wesn(isnow)/asnow_fcst ) / cat_progn(kk,n_e)%sndz(isnow) - - ! ii) diagnose analysis snow depth using forecast density - - tmp_sndz(kk,n_e,isnow) = ( tmp_wesn(kk,n_e,isnow)/asnow_ana ) / snow_dens - - end if - - end if - - end do ! isnow=1,N_snow (compute SWE, snow heat content, and snow depth analysis for each layer) - - ! 4. Relayer to balance the snow column (call with optional args for adjustment of htsnn) - - call StieglitzSnow_relayer( N_snow, N_constit, & + + tmp_htsn(kk,n_e,isnow) = (StieglitzSnow_CPW*snow_temp - fice_snow*MAPL_ALHF)*tmp_wesn(kk,n_e,isnow) + + ! update snow depth: + + if (asnow_ana < 1. .and. asnow_fcst < 1.) then + + ! keep snow depth constant when less than full snow cover in fcst and ana + + tmp_sndz(kk,n_e,isnow) = cat_progn(kk,n_e)%sndz(isnow) + + else + + ! compute analysis snow depth by keeping snow density constant + ! + ! in this case, it is possible that either asnow_fcst<1 or asnow_ana<1; + ! when computing density or depth, make sure that SWE value (which is per unit area) is + ! adjusted to reflect SWE value (per unit area) in the snow-covered fraction of the tile + + ! i) diagnose (layer-specific) forecast snow density + + snow_dens = ( cat_progn(kk,n_e)%wesn(isnow)/asnow_fcst ) / cat_progn(kk,n_e)%sndz(isnow) + + ! ii) diagnose analysis snow depth using forecast density + + tmp_sndz(kk,n_e,isnow) = ( tmp_wesn(kk,n_e,isnow)/asnow_ana ) / snow_dens + + end if + + end if + + end do ! isnow=1,N_snow (compute SWE, snow heat content, and snow depth analysis for each layer) + + ! 4. Relayer to balance the snow column (call with optional args for adjustment of htsnn) + + call StieglitzSnow_relayer( N_snow, N_constit, & MAPL_LAND, CATCH_SNOW_DZPARAM, & tmp_htsn(kk,n_e,1:N_snow), & tmp_wesn(kk,n_e,1:N_snow), & tmp_sndz(kk,n_e,1:N_snow), & rconstit, tpsn, fice_snow_vec ) - - ! print the old and new swe, heat content and snow density - - !if (logit) write (logunit, *) & - ! 'fcst_wesn = ', cat_progn(kk, n_e)%wesn(1:N_snow), & - ! 'tmp_wesn = ', tmp_wesn( kk,n_e, 1:N_snow), & - ! 'fcst_htsn = ', cat_progn(kk, n_e)%htsn(1:N_snow), & - ! 'tmp_htsn = ', tmp_htsn( kk, n_e, 1:N_snow), & - ! 'fcst_sndz = ', cat_progn(kk, n_e)%sndz(1:N_snow), & - ! 'tmp_sndz = ', tmp_sndz( kk ,n_e, 1:N_snow), & - ! '--------------------------------------' - - ! 5. Diagnose increments - - cat_progn_incr(kk,n_e)%wesn(1:N_snow) = tmp_wesn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%wesn(1:N_snow) - cat_progn_incr(kk,n_e)%htsn(1:N_snow) = tmp_htsn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%htsn(1:N_snow) - cat_progn_incr(kk,n_e)%sndz(1:N_snow) = tmp_sndz(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%sndz(1:N_snow) - - end do ! n_e=1,N_ens - - end if ! if (N_selected_obs > 0) - - end if ! if (has_asnow) - - ! Calculate swe_incr_ensavg as the average of the ensemble members to add to check if tile is snow-free - + + ! print the old and new swe, heat content and snow density + + !if (logit) write (logunit, *) & + ! 'fcst_wesn = ', cat_progn(kk, n_e)%wesn(1:N_snow), & + ! 'tmp_wesn = ', tmp_wesn( kk,n_e, 1:N_snow), & + ! 'fcst_htsn = ', cat_progn(kk, n_e)%htsn(1:N_snow), & + ! 'tmp_htsn = ', tmp_htsn( kk, n_e, 1:N_snow), & + ! 'fcst_sndz = ', cat_progn(kk, n_e)%sndz(1:N_snow), & + ! 'tmp_sndz = ', tmp_sndz( kk ,n_e, 1:N_snow), & + ! '--------------------------------------' + + ! 5. Diagnose increments + + cat_progn_incr(kk,n_e)%wesn(1:N_snow) = tmp_wesn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%wesn(1:N_snow) + cat_progn_incr(kk,n_e)%htsn(1:N_snow) = tmp_htsn(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%htsn(1:N_snow) + cat_progn_incr(kk,n_e)%sndz(1:N_snow) = tmp_sndz(kk,n_e,1:N_snow) - cat_progn(kk,n_e)%sndz(1:N_snow) + + end do ! n_e=1,N_ens + + end if ! if (N_selected_obs > 0) + + end if ! if (N_select_species_asnow>0) [assimilate snow observations for tile kk] + + ! Calculate swe_incr_ensavg as the average of the ensemble members to add to check if tile is snow-free + swe_incr_ensavg(kk) = sum((swe_incr(kk,:))) / real(N_ens) + + ! ------------------------------------------------------------------------------------------ - if (has_Tb .or. has_sfmc .or. has_sfds) then ! Have soil moisture observations - - if (logit) write (logunit,*) 'get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc obs' - - 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. & + if (N_select_species_smTb>0) then ! assimilate soil moisture observations + + 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. & (swe_incr_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), & + + ! 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_smTb, select_species_smTb(1:N_select_species_smTb), & 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, & + + 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( & + + 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),:), & @@ -5283,55 +5258,55 @@ subroutine cat_enkf_increments( & 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 ! if (N_selected_obs > 0) - - end if ! thresholds - - end if ! if (have soil moisture observations) + + 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 ! if (N_selected_obs > 0) + + end if ! thresholds + + 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 default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown update_type') @@ -5362,9 +5337,18 @@ subroutine get_select_species( & N_select_varnames, select_varnames, N_obs_param, obs_param, & N_select_species, select_species ) - ! find out obs species ID numbers ("obs_param%species") that correspond - ! to a set of obs species variables names ("obs_param%varname") - ! - reichle, 16 Oct 2014 + ! find out obs species ID numbers ("obs_param%species") that correspond to a set of + ! obs variables names ("obs_param%varname") with obs_param(species)%assim==.true. + + ! original version -- reichle, 16 Oct 2014 + ! + ! added constraint obs_param%assim==.true. + ! --> excludes obs species that are categorically not assimilated (per obs_param); + ! before this change, this constraint was applied at a later stage of + ! selecting the assimilated obs; + ! with this change, N_select_species==0 can be used to check if + ! any obs with select_varnames are assimilated. + ! -reichle, 29 Nov 2024 implicit none @@ -5394,7 +5378,7 @@ subroutine get_select_species( & do ii=1,N_obs_param - if (any(trim(obs_param(ii)%varname)==select_varnames)) then + if ( any(trim(obs_param(ii)%varname)==select_varnames) .and. obs_param(ii)%assim ) then kk = kk+1 @@ -5455,9 +5439,7 @@ subroutine get_ind_obs( & ! locals integer :: i, k - - - + ! -------------------------------------------------------------- if (N_select_species==0 .and. N_select_tilenum==0) then @@ -5528,7 +5510,7 @@ end subroutine get_ind_obs ! ********************************************************************** - + subroutine get_ind_obs_assim( N_obs, assim_flag, N_obs_assim, ind_obs_assim ) ! loop through Observations%assim and construct an index vector that maps to @@ -5902,7 +5884,7 @@ subroutine check_compact_support( & Iam // '(): reset for 1d update_type: ycompact = ', ycompact if (logit) write (logunit,*) - case (2,7,8,10,13,14) ! "3d" updates, check consistency of xcompact, ycompact + case (2,7,8,10,12,13) ! "3d" updates, check consistency of xcompact, ycompact ! check xcompact/ycompact against corr scales of model error diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index 53e428f..6f4951a 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -9,11 +9,10 @@ ! 7 May 2021 - removed "dtstep_assim" and "centered_update"; replaced with MAPL ! resource parameters "LANDASSIM_DT" and "LANDASSIM_T0" (in LDAS.rc) ! -! -------------------------------------------------------------------- +! ---------------------------------------------------------------------- &ens_upd_inputs - ! ---------------------------------------------------------------------- ! ! update type - for details see subroutine cat_enkf_update() @@ -23,21 +22,23 @@ ! ! # = 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 = 11: 1d snow analysis (Toure et al. 2018 empirical gain); SCF obs -! update_type = 13: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile; sfmc and TB obs -! update_type = 14: 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile & -! 1d snow analysis (Toure et al. 2018 empirical gain); sfmc, TB and SCF obs +! update_type | analysis state vector | assimilated obs +! ------------------------------------------------------------------------------------------------- +! 0 | NO assimilation, NO bias correction | n/a +! # 1 | 1d soil moisture | sfmc +! # 2 | 3d soil moisture | sfmc +! 3 | 1d Tskin (assim incr NOT applied, use w/ bias corr) | Tskin +! 4 | 1d Tskin/ght1 (assim incr applied, use w/ or w/o bias corr) | Tskin +! 5 | 1d Tskin/ght1 (assim incr NOT applied, use w/ bias corr) | Tskin +! 6 | 1d soil moisture/Tskin/ght(1) | Tb +! 7 | 3d Tskin/ght1 update | Tskin +! 8 | 3d soil moisture/Tskin/ght(1) | Tb +! 9 | 1d Tskin/ght1 update | FT +! 10 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | Tb +! 11 | 1d snow analysis (Toure et al. 2018 empirical gain) | SCF +! 12 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | sfmc/sfds, Tb, SCF +! | & 1d snow analysis (Toure et al. 2018 empirical gain) | +! 13 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | sfmc/sfds, Tb update_type = 0 From a1b9ca2240c61312f137663433cbd3d25943067a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 29 Nov 2024 14:51:45 -0500 Subject: [PATCH 14/21] changed error message to warning (inconsistent update_type and obs_param%assim) --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 95bad86..af66955 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -3507,6 +3507,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 @@ -3542,6 +3543,8 @@ subroutine cat_enkf_increments( & logical :: found_Tb_obs + + ! ----------------------------------------------------------------------- if (logit) write (logunit,*) & @@ -4957,8 +4960,12 @@ subroutine cat_enkf_increments( & if ( N_select_species_asnow==0 .and. N_select_species_smTb==0 ) then - call ldas_abort(LDAS_GENERIC_ERROR, Iam, & - 'update_type inconsistent with obs_param%assim==.false. for all asnow/Tb/sfmc/sfds species') + err_msg = & + 'update_type not consistent with obs_param%assim==.false. for all asnow/Tb/sfmc/sfds species ' & + // '(may be intentional for "innovations" run, i.e., with obs_param%innov==.true.)' + + call ldas_warn(LDAS_GENERIC_ERROR, Iam, & + end if From c88b2a2c2da1d7f7bc0f609368effc93fbe616d2 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 29 Nov 2024 17:08:35 -0500 Subject: [PATCH 15/21] updated CHANGELOG.md --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d07a99f..55fce11 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- New update_type for joint 3d soil moisture and 1d snow analysis (Tb+sfmc+sfds+SCF obs) + ### Changed ### Fixed From f4e9626f545fecc7e913a4d450f91352c625362a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 29 Nov 2024 18:59:21 -0500 Subject: [PATCH 16/21] fixed syntax error in previous commit --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index af66955..13f2ace 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4964,7 +4964,7 @@ subroutine cat_enkf_increments( & 'update_type not consistent with obs_param%assim==.false. for all asnow/Tb/sfmc/sfds species ' & // '(may be intentional for "innovations" run, i.e., with obs_param%innov==.true.)' - call ldas_warn(LDAS_GENERIC_ERROR, Iam, & + call ldas_warn(LDAS_GENERIC_ERROR, Iam, err_msg) end if From e7c91e11718b466a3f0c1db661aa5bb9dca4aa48 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 30 Nov 2024 08:55:33 -0500 Subject: [PATCH 17/21] fixed build error in previous commit --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 13f2ace..52958e0 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4911,7 +4911,7 @@ subroutine cat_enkf_increments( & ! ! ------------------------------------------------------------------------------------------------------- - if (logit) write (logunit, *) 'get increments: update_type = ', select_update_type + if (logit) write (logunit, *) 'get increments: update_type = ', update_type ! determine species of assimilated obs associated with snow analysis From bced10ad21585e06f05d36689c8ac66e6b88960b Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 30 Nov 2024 09:27:51 -0500 Subject: [PATCH 18/21] fixed build error in previous commit --- GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 index 65e762a..eedd28c 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_enkf_update.F90 @@ -1407,11 +1407,7 @@ subroutine apply_enkf_increments( N_catd, N_ens, update_type, & cat_progn_has_changed = .true. - if (select_update_type==10) then - - check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case - - end if + if (update_type==10) check_snow = .false. ! turn off for now to maintain 0-diff w/ SMAP Tb DA test case ! ------------------------------------------------------------------ From c4a1826e68fc16cb7893455e1f9d60fadd39b718 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 30 Nov 2024 10:25:22 -0500 Subject: [PATCH 19/21] fixing linking error in previous commit --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 52958e0..3a9e182 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -135,6 +135,7 @@ module clsm_ensupd_upd_routines use LDAS_ExceptionsMod, ONLY: & ldas_abort, & + ldas_warn, & LDAS_GENERIC_ERROR use enkf_general, ONLY: & From 067305fd0ebd16508e27a26a32d8e3efc32f4779 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 2 Dec 2024 14:25:39 -0500 Subject: [PATCH 20/21] fixed log statement for soil moisture analysis in new update_type=12 (clsm_ensupd_upd_routines.F90) --- GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 3a9e182..85e7c1c 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -4951,12 +4951,9 @@ subroutine cat_enkf_increments( & call get_select_species(1, 'Tb', N_obs_param, obs_param, N_select_species_Tb, select_species_Tb ) - if (N_select_species_smTb>0) then - - if (logit) write (logunit, *) '- get 1d snow increments (Toure et al. 2018 empirical gain); snow cover fraction obs' + if ((N_select_species_smTb>0) .and. logit) & + write (logunit, *) '- get 3d soil moisture/Tskin/ght(1) increments; Tb+sfmc+sfds obs' - end if - ! check consistency of update_type and obs_param%assim config if ( N_select_species_asnow==0 .and. N_select_species_smTb==0 ) then From 829f229baa6b7bd1613c7120465d286059787708 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 2 Dec 2024 14:29:55 -0500 Subject: [PATCH 21/21] 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 ! ----------------------------------------------------------------------------------------------------------------------