From 5cec2190f59624e36079602f6cc24e5b7fca4582 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 31 Oct 2024 10:39:30 -0600 Subject: [PATCH 01/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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')