From ff870eb45fe9ab8abd2f09727296222d8adc993e Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 8 Nov 2024 15:12:42 +0100 Subject: [PATCH 01/32] Snow-DA: initial implementation from CLM5-PDAF --- interface/framework/obs_op_pdaf.F90 | 11 +- interface/model/clm5_0/enkf_clm_mod_5.F90 | 355 +++++++++++++++++++++- interface/model/common/enkf.h | 2 + interface/model/common/read_enkfpar.c | 2 + interface/model/eclm/enkf_clm_mod_5.F90 | 355 +++++++++++++++++++++- interface/model/wrapper_tsmp.c | 2 +- 6 files changed, 721 insertions(+), 6 deletions(-) diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 23a755d4f..29bf1e8cb 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -211,9 +211,14 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) if(lpointobs) then - DO i = 1, dim_obs_p - m_state_p(i) = state_p(obs_index_p(i)) - END DO + ! WORKAROUND FOR SUDDEN SEG 11 error + ! DO NOT LEAVE IN FINAL VERSION + ! JUST DURING DEBUG PROCESS + ! JUST VALID FOR SINGLE GRID CELLS WITH SINGLE VALUE OBS + m_state_p(1) = state_p(1) +! DO i = 1, dim_obs_p +! m_state_p(i) = state_p(obs_index_p(i)) +! END DO end if diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index c0b89c299..1e7cea324 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -50,6 +50,8 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture + integer(c_int),bind(C,name="clmupdate_snow") :: clmupdate_snow + integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et @@ -270,6 +272,19 @@ subroutine define_clm_statevec(mype) clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) endif + ! Snow assimilation + ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 + ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell + if(clmupdate_snow.eq.1) then + clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA + clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority + endif + ! Case 2: Assimilation of snow water equivalent same as above + if(clmupdate_snow.eq.2) then + clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA + clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority + endif + !hcp LST DA if(clmupdate_T.eq.1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -322,6 +337,8 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -331,6 +348,9 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -384,6 +404,60 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! Snow assimilation + ! Case 1: Snow depth + if(clmupdate_snow.eq.1) then + cc = 1 + do j=clm_begg,clm_endg + ! Only get the snow_depth from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = snow_depth(jj) + endif + endif + end do + cc = cc + 1 + end do + endif + ! Case 2: SWE + if(clmupdate_snow.eq.2) then + cc = 1 + + if(clmstatevec_allcol.eq.0) then + + do j=clm_begg,clm_endg + ! Only get the SWE from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = h2osno(jj) + endif + endif + end do + cc = cc + 1 + end do + + else + + do jj=clm_begc,clm_endc + ! SWC from all columns of each gridcell + clm_statevec(cc+offset) = h2osno(jj) + cc = cc + 1 + end do + + end if + + endif + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -419,10 +493,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice + + real(r8) :: rliq,rice,incr_h2osno + real(r8) :: rsnow(clm_statevecsize) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -435,6 +514,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu + integer, pointer :: snlsno(:) logical :: swc_zero_before_update = .false. #ifdef PDAF_DEBUG @@ -455,9 +535,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + dz => col%dz h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + snlsno => col%snl ! number of snow layers (negative) #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -634,8 +718,277 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif + ! Snow assimilation: + ! Case 1: Snow depth + ! Write updated snow depth back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.eq.1) then + cc = 1 + do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do jj=clm_begc,clm_endc + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", j, jj, clm_statevec(cc+offset) + else + rsnow(jj) = snow_depth(jj) - clm_statevec(cc+offset) + if ( ABS(SUM(rsnow(:))) .gt.0.000001) then + snow_depth(jj) = clm_statevec(cc+offset) + endif + endif + end do + cc = cc + 1 + end do + if ( clmupdate_snow_repartitioning.ne.0) then + if ( ABS(SUM(rsnow(:))).gt.0.000001) then + call clm_repartition_snow() + end if + end if + endif + ! Case 2: Snow water equivalent + ! Write updated snow depth back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.eq.2) then + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if + + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif + + end if + endif + end do + ! cc = cc + 1 + ! end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + endif + end subroutine update_clm + subroutine clm_repartition_snow(h2osno_in) + use ColumnType, only : col + use clm_instMod, only : waterstate_inst + use clm_varpar, only : nlevsno, nlevsoi + use clm_varcon, only : bdsno, denice + use shr_kind_mod, only : r8 => shr_kind_r8 + + implicit none + + real(r8), intent(in), optional :: h2osno_in(clm_begc:clm_endc) + + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: h2oliq(:,:) + real(r8), pointer :: h2oice(:,:) + real(r8), pointer :: frac_sno(:) + real(r8), pointer :: snowdp(:) + integer, pointer :: snlsno(:) + + real(r8) :: dzsno(clm_begc:clm_endc,-nlevsno+1:0) + real(r8) :: h2osno_po(clm_begc:clm_endc) + real(r8) :: h2osno_pr(clm_begc:clm_endc) + real(r8) :: snowden, frac_swe, frac_liq, frac_ice + real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice + real(r8) :: gain_dzsno(-nlevsno+1:0) + real(r8) :: rho_avg, z_avg + integer :: i,ii,j,jj,g,cc=1,offset=0 + + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + snowdp => waterstate_inst%snowdp_col ! area-averaged snow height (m) + h2osno => waterstate_inst%h2osno_col ! col snow water (mm H2O) + h2oliq => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) + h2oice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) + + snlsno => col%snl ! number of snow layers (negative) + + frac_sno => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow + ! dz for snow layers is defined like in the history output as col%dz for the snow layers + dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0) + ! Iterate through all columns + do jj=clm_begc,clm_endc + if (h2osno(jj).lt.0.0) then ! No snow in column + print *, "WARNING: negative snow in col: ", jj, h2osno +! ! Set existing layers to near zero and let CLM do the layer aggregation + do i=0,snlsno(jj)+1,-1 + h2oliq(jj,i) = 0.0_r8 + h2oice(jj,i) = 0.00000001_r8 + dzsno(jj,i) = 0.00000001_r8 + end do + snow_depth(jj) = sum(dzsno(jj,:)) + h2osno(jj) = sum(h2oice(jj,:)) + + if (isnan(h2osno(jj))) then + print *, "WARNING: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING: snow_depth at jj is nan: ", jj + endif + + else ! snow (h2osno) present + if (snlsno(jj).lt.0) then ! snow layers in the column + if (clmupdate_snow .eq. 1) then + ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90 + ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE + ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth + ! Therefore need to have a transform to get h2osno_po + ! v1 init + ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD + ! v2 SoilTemperatureMod + if (snowdp(jj).gt.0.0_r8) then + rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj)) + else + rho_avg = 200.0_r8 + end if + if (frac_sno(jj).gt.0.0_r8 .and. snlsno(jj).lt.0.0_r8) then + h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj)) + else + h2osno_po(jj) = snow_depth(jj) * denice + end if + h2osno_pr(jj) = h2osno(jj) + else if (clmupdate_snow .eq. 2) then + ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already + h2osno_po(jj) = h2osno(jj) + h2osno_pr(jj) = h2osno_in(jj) + end if + + do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers + ! DART VERSION: ii = nlevsoi - i + 1 + ! snow density prior for each layer + if (dzsno(jj,ii).gt.0.0_r8) then + snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii)) + else + snowden = 0.0_r8 + endif + ! fraction of SWE in each active layers + if(h2osno_pr(jj).gt.0.0_r8) then + ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else) + if(clmupdate_snow_repartitioning.eq.1) then + if (ii .eq. 0) then ! DART version indexing check against nlevsno but for us 0 + frac_swe = 1.0_r8 + ! JK: Let CLM repartitioning do the job + ! afterwards. Provide all the snow in a single layer + else + frac_swe = 0.0_r8 + end if + ! repartition Option 2: Adjust all active layers + else if (clmupdate_snow_repartitioning.eq.2) then + frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj) + end if + else + frac_swe = 0.0_r8 ! no fraction SWE if no snow is present in column + end if ! end SWE fraction if + + ! fraction of liquid and ice + if ((h2oliq(jj,ii) + h2oice(jj,ii)).gt.0.0_r8) then + frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii)) + frac_ice = 1.0_r8 - frac_liq + else + frac_liq = 0.0_r8 + frac_ice = 0.0_r8 + end if + + ! SWE adjustment per layer + ! assumes identical layer distribution of liq and ice than before DA (frac_*) + if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe + gain_h2oliq = gain_h2osno * frac_liq + gain_h2oice = gain_h2osno * frac_ice + else + gain_h2osno = 0.0_r8 + gain_h2oliq = 0.0_r8 + gain_h2oice = 0.0_r8 + end if + ! layer level adjustments + if (snowden.gt.0.0_r8) then + gain_dzsno(ii) = gain_h2osno / snowden + else + gain_dzsno(ii) = 0.0_r8 + end if + h2oliq(jj,ii) = h2oliq(jj,ii) + gain_h2oliq + h2oice(jj,ii) = h2oice(jj,ii) + gain_h2oice + + ! Adjust snow layer dimensions so that CLM5 can calculate compaction / aggregation + ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic + ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod + ! therefore we adjust dz(:, snow_layer) here + if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii) + ! mid point and interface adjustments + ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers)) + ! DART version the sum goes from ilevel:nlevsno to fit with our indexing: + col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8 + ! In DART the check is ilevel == nlevsno but here + if (ii.eq.0) then + col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8 + else + col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8 + end if + end if + + end do ! End iteration of snow layers + endif ! End of snow layers present check + + ! Column level variables + if (clmupdate_snow .eq. 1) then + ! Finally adjust SWE (h2osno) since the prior value is no longer needed + ! column level variables so we can adjust it outside the layer loop + h2osno(jj) = h2osno_po(jj) + else if (clmupdate_snow .eq. 2) then + ! Update the total snow depth to match udpates to layers for active snow layers + if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then + snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) + end if + end if + + if (isnan(h2osno(jj))) then + print *, "WARNING2: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING2: snow_depth at jj is nan: ", jj + endif + + end if ! End of snow present check + end do ! End of column iteration + + end subroutine clm_repartition_snow + subroutine clm_correct_texture() use clm_varpar , only : nlevsoi diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index c3826684c..4c4067d85 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -87,6 +87,8 @@ GLOBAL int nx_local,ny_local,nz_local; GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; +GLOBAL int clmupdate_snow; +GLOBAL int clmupdate_snow_repartitioning; GLOBAL int clmprint_swc; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index 2adcfa156..a8a7f4761 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -78,6 +78,8 @@ void read_enkfpar(char *parname) clmupdate_swc = iniparser_getint(pardict,"CLM:update_swc",1); clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); + clmupdate_snow = iniparser_getint(pardict,"CLM:update_snow",0); + clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",1); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index c0b89c299..1e7cea324 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -50,6 +50,8 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture + integer(c_int),bind(C,name="clmupdate_snow") :: clmupdate_snow + integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et @@ -270,6 +272,19 @@ subroutine define_clm_statevec(mype) clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) endif + ! Snow assimilation + ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 + ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell + if(clmupdate_snow.eq.1) then + clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA + clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority + endif + ! Case 2: Assimilation of snow water equivalent same as above + if(clmupdate_snow.eq.2) then + clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA + clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority + endif + !hcp LST DA if(clmupdate_T.eq.1) then error stop "Not implemented: clmupdate_T.eq.1" @@ -322,6 +337,8 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -331,6 +348,9 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -384,6 +404,60 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! Snow assimilation + ! Case 1: Snow depth + if(clmupdate_snow.eq.1) then + cc = 1 + do j=clm_begg,clm_endg + ! Only get the snow_depth from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = snow_depth(jj) + endif + endif + end do + cc = cc + 1 + end do + endif + ! Case 2: SWE + if(clmupdate_snow.eq.2) then + cc = 1 + + if(clmstatevec_allcol.eq.0) then + + do j=clm_begg,clm_endg + ! Only get the SWE from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = h2osno(jj) + endif + endif + end do + cc = cc + 1 + end do + + else + + do jj=clm_begc,clm_endc + ! SWC from all columns of each gridcell + clm_statevec(cc+offset) = h2osno(jj) + cc = cc + 1 + end do + + end if + + endif + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -419,10 +493,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice + + real(r8) :: rliq,rice,incr_h2osno + real(r8) :: rsnow(clm_statevecsize) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -435,6 +514,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu + integer, pointer :: snlsno(:) logical :: swc_zero_before_update = .false. #ifdef PDAF_DEBUG @@ -455,9 +535,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + dz => col%dz h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + snlsno => col%snl ! number of snow layers (negative) #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -634,8 +718,277 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif + ! Snow assimilation: + ! Case 1: Snow depth + ! Write updated snow depth back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.eq.1) then + cc = 1 + do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do jj=clm_begc,clm_endc + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", j, jj, clm_statevec(cc+offset) + else + rsnow(jj) = snow_depth(jj) - clm_statevec(cc+offset) + if ( ABS(SUM(rsnow(:))) .gt.0.000001) then + snow_depth(jj) = clm_statevec(cc+offset) + endif + endif + end do + cc = cc + 1 + end do + if ( clmupdate_snow_repartitioning.ne.0) then + if ( ABS(SUM(rsnow(:))).gt.0.000001) then + call clm_repartition_snow() + end if + end if + endif + ! Case 2: Snow water equivalent + ! Write updated snow depth back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.eq.2) then + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if + + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif + + end if + endif + end do + ! cc = cc + 1 + ! end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + endif + end subroutine update_clm + subroutine clm_repartition_snow(h2osno_in) + use ColumnType, only : col + use clm_instMod, only : waterstate_inst + use clm_varpar, only : nlevsno, nlevsoi + use clm_varcon, only : bdsno, denice + use shr_kind_mod, only : r8 => shr_kind_r8 + + implicit none + + real(r8), intent(in), optional :: h2osno_in(clm_begc:clm_endc) + + real(r8), pointer :: snow_depth(:) + real(r8), pointer :: h2osno(:) + real(r8), pointer :: h2oliq(:,:) + real(r8), pointer :: h2oice(:,:) + real(r8), pointer :: frac_sno(:) + real(r8), pointer :: snowdp(:) + integer, pointer :: snlsno(:) + + real(r8) :: dzsno(clm_begc:clm_endc,-nlevsno+1:0) + real(r8) :: h2osno_po(clm_begc:clm_endc) + real(r8) :: h2osno_pr(clm_begc:clm_endc) + real(r8) :: snowden, frac_swe, frac_liq, frac_ice + real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice + real(r8) :: gain_dzsno(-nlevsno+1:0) + real(r8) :: rho_avg, z_avg + integer :: i,ii,j,jj,g,cc=1,offset=0 + + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + snowdp => waterstate_inst%snowdp_col ! area-averaged snow height (m) + h2osno => waterstate_inst%h2osno_col ! col snow water (mm H2O) + h2oliq => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) + h2oice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) + + snlsno => col%snl ! number of snow layers (negative) + + frac_sno => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow + ! dz for snow layers is defined like in the history output as col%dz for the snow layers + dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0) + ! Iterate through all columns + do jj=clm_begc,clm_endc + if (h2osno(jj).lt.0.0) then ! No snow in column + print *, "WARNING: negative snow in col: ", jj, h2osno +! ! Set existing layers to near zero and let CLM do the layer aggregation + do i=0,snlsno(jj)+1,-1 + h2oliq(jj,i) = 0.0_r8 + h2oice(jj,i) = 0.00000001_r8 + dzsno(jj,i) = 0.00000001_r8 + end do + snow_depth(jj) = sum(dzsno(jj,:)) + h2osno(jj) = sum(h2oice(jj,:)) + + if (isnan(h2osno(jj))) then + print *, "WARNING: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING: snow_depth at jj is nan: ", jj + endif + + else ! snow (h2osno) present + if (snlsno(jj).lt.0) then ! snow layers in the column + if (clmupdate_snow .eq. 1) then + ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90 + ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE + ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth + ! Therefore need to have a transform to get h2osno_po + ! v1 init + ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD + ! v2 SoilTemperatureMod + if (snowdp(jj).gt.0.0_r8) then + rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj)) + else + rho_avg = 200.0_r8 + end if + if (frac_sno(jj).gt.0.0_r8 .and. snlsno(jj).lt.0.0_r8) then + h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj)) + else + h2osno_po(jj) = snow_depth(jj) * denice + end if + h2osno_pr(jj) = h2osno(jj) + else if (clmupdate_snow .eq. 2) then + ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already + h2osno_po(jj) = h2osno(jj) + h2osno_pr(jj) = h2osno_in(jj) + end if + + do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers + ! DART VERSION: ii = nlevsoi - i + 1 + ! snow density prior for each layer + if (dzsno(jj,ii).gt.0.0_r8) then + snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii)) + else + snowden = 0.0_r8 + endif + ! fraction of SWE in each active layers + if(h2osno_pr(jj).gt.0.0_r8) then + ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else) + if(clmupdate_snow_repartitioning.eq.1) then + if (ii .eq. 0) then ! DART version indexing check against nlevsno but for us 0 + frac_swe = 1.0_r8 + ! JK: Let CLM repartitioning do the job + ! afterwards. Provide all the snow in a single layer + else + frac_swe = 0.0_r8 + end if + ! repartition Option 2: Adjust all active layers + else if (clmupdate_snow_repartitioning.eq.2) then + frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj) + end if + else + frac_swe = 0.0_r8 ! no fraction SWE if no snow is present in column + end if ! end SWE fraction if + + ! fraction of liquid and ice + if ((h2oliq(jj,ii) + h2oice(jj,ii)).gt.0.0_r8) then + frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii)) + frac_ice = 1.0_r8 - frac_liq + else + frac_liq = 0.0_r8 + frac_ice = 0.0_r8 + end if + + ! SWE adjustment per layer + ! assumes identical layer distribution of liq and ice than before DA (frac_*) + if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe + gain_h2oliq = gain_h2osno * frac_liq + gain_h2oice = gain_h2osno * frac_ice + else + gain_h2osno = 0.0_r8 + gain_h2oliq = 0.0_r8 + gain_h2oice = 0.0_r8 + end if + ! layer level adjustments + if (snowden.gt.0.0_r8) then + gain_dzsno(ii) = gain_h2osno / snowden + else + gain_dzsno(ii) = 0.0_r8 + end if + h2oliq(jj,ii) = h2oliq(jj,ii) + gain_h2oliq + h2oice(jj,ii) = h2oice(jj,ii) + gain_h2oice + + ! Adjust snow layer dimensions so that CLM5 can calculate compaction / aggregation + ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic + ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod + ! therefore we adjust dz(:, snow_layer) here + if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii) + ! mid point and interface adjustments + ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers)) + ! DART version the sum goes from ilevel:nlevsno to fit with our indexing: + col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8 + ! In DART the check is ilevel == nlevsno but here + if (ii.eq.0) then + col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8 + else + col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8 + end if + end if + + end do ! End iteration of snow layers + endif ! End of snow layers present check + + ! Column level variables + if (clmupdate_snow .eq. 1) then + ! Finally adjust SWE (h2osno) since the prior value is no longer needed + ! column level variables so we can adjust it outside the layer loop + h2osno(jj) = h2osno_po(jj) + else if (clmupdate_snow .eq. 2) then + ! Update the total snow depth to match udpates to layers for active snow layers + if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then + snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) + end if + end if + + if (isnan(h2osno(jj))) then + print *, "WARNING2: h2osno at jj is nan: ", jj + endif + if (isnan(snow_depth(jj))) then + print *, "WARNING2: snow_depth at jj is nan: ", jj + endif + + end if ! End of snow present check + end do ! End of column iteration + + end subroutine clm_repartition_snow + subroutine clm_correct_texture() use clm_varpar , only : nlevsoi diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 1af2bb6ae..d3e42b9de 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -198,7 +198,7 @@ void integrate_tsmp() { void update_tsmp(){ #if defined CLMSA - if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ + if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_snow!=0))){ update_clm(&tstartcycle, &mype_world); print_update_clm(&tcycle, &total_steps); } From 867d7155a3c33982ee2d6fe4ab4bcf123d524de6 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 7 Feb 2025 14:56:54 +0100 Subject: [PATCH 02/32] Snow-DA: Add MPI-barrier after `integrate_tsmp` --- interface/framework/pdaf_terrsysmp.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/interface/framework/pdaf_terrsysmp.F90 b/interface/framework/pdaf_terrsysmp.F90 index 4ed09922b..c28606fa0 100644 --- a/interface/framework/pdaf_terrsysmp.F90 +++ b/interface/framework/pdaf_terrsysmp.F90 @@ -88,6 +88,9 @@ PROGRAM pdaf_terrsysmp ! forward simulation of component models CALL integrate_tsmp() + ! barrier before model integration starts + CALL MPI_BARRIER(MPI_COMM_WORLD, MPIerr) + ! assimilation step CALL assimilate_pdaf() From c9fb4f66028919761ecd75ddf564492fa62f10b4 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 4 Mar 2025 15:56:44 +0100 Subject: [PATCH 03/32] Update enkf_clm_mod_5.F90 (#265) * Update enkf_clm_mod_5.F90 Updated snow depth (clmupdate_snow.eq.1) to include the clmupdate_snow_repartitioning.eq.3 option and updated the column iteriation counter to be identical to clmupdate_snow.eq.2 . * copy updates to eCLM-path --------- Co-authored-by: Johannes Keller --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 49 +++++++++++++++-------- interface/model/eclm/enkf_clm_mod_5.F90 | 49 +++++++++++++++-------- 2 files changed, 64 insertions(+), 34 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 1e7cea324..9a9131f76 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -722,28 +722,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Case 1: Snow depth ! Write updated snow depth back to CLM and then repartition snow and adjust related variables if(clmupdate_snow.eq.1) then - cc = 1 - do j=clm_begg,clm_endg + do j=clm_begc,clm_endc ! iterate through the columns and copy from the same gridcell ! i.e. statevec position (cc) for each column - do jj=clm_begc,clm_endc - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", j, jj, clm_statevec(cc+offset) - else - rsnow(jj) = snow_depth(jj) - clm_statevec(cc+offset) - if ( ABS(SUM(rsnow(:))) .gt.0.000001) then - snow_depth(jj) = clm_statevec(cc+offset) - endif - endif - end do - cc = cc + 1 - end do - if ( clmupdate_snow_repartitioning.ne.0) then + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) + else + rsnow(j) = snow_depth(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))) .gt.0.000001) then + snow_depth(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + endif + endif + endif + end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then if ( ABS(SUM(rsnow(:))).gt.0.000001) then call clm_repartition_snow() end if - end if + end if endif ! Case 2: Snow water equivalent ! Write updated snow depth back to CLM and then repartition snow and adjust related variables diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 1e7cea324..9a9131f76 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -722,28 +722,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Case 1: Snow depth ! Write updated snow depth back to CLM and then repartition snow and adjust related variables if(clmupdate_snow.eq.1) then - cc = 1 - do j=clm_begg,clm_endg + do j=clm_begc,clm_endc ! iterate through the columns and copy from the same gridcell ! i.e. statevec position (cc) for each column - do jj=clm_begc,clm_endc - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", j, jj, clm_statevec(cc+offset) - else - rsnow(jj) = snow_depth(jj) - clm_statevec(cc+offset) - if ( ABS(SUM(rsnow(:))) .gt.0.000001) then - snow_depth(jj) = clm_statevec(cc+offset) - endif - endif - end do - cc = cc + 1 - end do - if ( clmupdate_snow_repartitioning.ne.0) then + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) + else + rsnow(j) = snow_depth(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))) .gt.0.000001) then + snow_depth(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + endif + endif + endif + end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then if ( ABS(SUM(rsnow(:))).gt.0.000001) then call clm_repartition_snow() end if - end if + end if endif ! Case 2: Snow water equivalent ! Write updated snow depth back to CLM and then repartition snow and adjust related variables From 6dd9f7f328e1afb0f25e245759589d8a528c7776 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Mar 2025 14:47:09 +0100 Subject: [PATCH 04/32] clmupdate_snow.eq.3: snow depth observations, h2osoi_ice in state --- interface/model/eclm/enkf_clm_mod_5.F90 | 378 +++++++++++++++--------- 1 file changed, 233 insertions(+), 145 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9a9131f76..3a5b20c03 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -86,7 +86,7 @@ module enkf_clm_mod subroutine define_clm_statevec(mype) use shr_kind_mod, only: r8 => shr_kind_r8 use decompMod , only : get_proc_bounds - use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevsoi, nlevsno, nlevgrnd use clm_varcon , only : ispval use ColumnType , only : col @@ -284,6 +284,11 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif + ! Case 3: Assimilation of snow depth: adding h2osoi_ice in the state vector + if(clmupdate_snow.eq.3) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*((clm_endg-clm_begg+1) * (nlevsno + nlevgrnd)) !2D h2osoi_ice + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -325,7 +330,7 @@ end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst - use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevsoi, nlevsno, nlevgrnd ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col @@ -339,6 +344,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: porgm(:,:) real(r8), pointer :: snow_depth(:) real(r8), pointer :: h2osno(:) + real(r8), pointer :: h2oice(:,:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -350,6 +356,7 @@ subroutine set_clm_statevec(tstartcycle, mype) snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + h2osoi_ice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -457,6 +464,49 @@ subroutine set_clm_statevec(tstartcycle, mype) end if endif + ! Case 3: Snow Depth with h2osoi_ice in state vector + if(clmupdate_snow.eq.3) then + + ! snow_depth into state vector + cc = 1 + do j=clm_begg,clm_endg + ! Only get the snow_depth from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = snow_depth(jj) + endif + endif + end do + cc = cc + 1 + end do + + ! h2osoi_ice into state vector + cc = 1 + ! Additional loop over layers + do i=-nlevsno+1,nlevgrnd + do j=clm_begg,clm_endg + ! Only get the h2osoi_ice from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + ! Shift by clm_varsize (size of snow-depth in state vector) + clm_statevec(cc+clm_varsize+offset) = h2osoi_ice(jj,i) + endif + endif + end do + cc = cc + 1 + end do + end do + endif #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -576,113 +626,113 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! write updated swc back to CLM if(clmupdate_swc.ne.0) then - ! Set minimum soil moisture for checking the state vector and - ! for setting minimum swc for CLM - if(clmwatmin_switch.eq.3) then - ! CLM3.5 type watmin - watmin_check = 0.00 - watmin_set = 0.05 - else if(clmwatmin_switch.eq.5) then - ! CLM5.0 type watmin - watmin_check = watmin - watmin_set = watmin - else - ! Default - watmin_check = 0.0 - watmin_set = 0.0 - end if + ! Set minimum soil moisture for checking the state vector and + ! for setting minimum swc for CLM + if(clmwatmin_switch.eq.3) then + ! CLM3.5 type watmin + watmin_check = 0.00 + watmin_set = 0.05 + else if(clmwatmin_switch.eq.5) then + ! CLM5.0 type watmin + watmin_check = watmin + watmin_set = watmin + else + ! Default + watmin_check = 0.0 + watmin_set = 0.0 + end if - ! cc = 0 - do i=1,nlevsoi - ! CLM3.5: iterate over grid cells - ! CLM5.0: iterate over columns - ! do j=clm_begg,clm_endg - do j=clm_begc,clm_endc - - ! Update only those SWCs that are not excluded by ispval - if(state_clm2pdaf_p(j,i) .ne. ispval) then - - if(swc(j,i).eq.0.0) then - swc_zero_before_update = .true. - - ! Zero-SWC leads to zero denominator in computation of - ! rliq/rice, therefore setting rliq/rice to special - ! value - rliq = spval - rice = spval - else - swc_zero_before_update = .false. + ! cc = 0 + do i=1,nlevsoi + ! CLM3.5: iterate over grid cells + ! CLM5.0: iterate over columns + ! do j=clm_begg,clm_endg + do j=clm_begc,clm_endc - rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) - rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) - !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) - end if + ! Update only those SWCs that are not excluded by ispval + if(state_clm2pdaf_p(j,i) .ne. ispval) then - swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + if(swc(j,i).eq.0.0) then + swc_zero_before_update = .true. - if(swc_update.le.watmin_check) then - swc(j,i) = watmin_set - else if(swc_update.ge.watsat(j,i)) then - swc(j,i) = watsat(j,i) - else - swc(j,i) = swc_update - endif + ! Zero-SWC leads to zero denominator in computation of + ! rliq/rice, therefore setting rliq/rice to special + ! value + rliq = spval + rice = spval + else + swc_zero_before_update = .false. - if (isnan(swc(j,i))) then - swc(j,i) = watmin_set - print *, "WARNING: swc at j,i is nan: ", j, i - endif + rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) + rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) + !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) + end if + + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + + if(swc_update.le.watmin_check) then + swc(j,i) = watmin_set + else if(swc_update.ge.watsat(j,i)) then + swc(j,i) = watsat(j,i) + else + swc(j,i) = swc_update + endif + + if (isnan(swc(j,i))) then + swc(j,i) = watmin_set + print *, "WARNING: swc at j,i is nan: ", j, i + endif - if(swc_zero_before_update) then - ! This case should not appear for hydrologically - ! active columns/layers, where always: swc > watmin - ! - ! If you want to make sure that no zero SWCs appear in - ! the code, comment out the error stop + if(swc_zero_before_update) then + ! This case should not appear for hydrologically + ! active columns/layers, where always: swc > watmin + ! + ! If you want to make sure that no zero SWCs appear in + ! the code, comment out the error stop #ifdef PDAF_DEBUG - ! error stop "ERROR: Update of zero-swc" - print *, "WARNING: Update of zero-swc" - print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i + ! error stop "ERROR: Update of zero-swc" + print *, "WARNING: Update of zero-swc" + print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i #endif - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o - h2osoi_ice(j,i) = 0.0 - else - ! update liquid water content - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq - ! update ice content - h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice - end if + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o + h2osoi_ice(j,i) = 0.0 + else + ! update liquid water content + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq + ! update ice content + h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice + end if - end if - ! cc = cc + 1 - end do + end if + ! cc = cc + 1 end do + end do #ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN - - IF(clmupdate_swc.NE.0) THEN - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn3, action="write") - WRITE (71,"(es22.15)") h2osoi_liq(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn4, action="write") - WRITE (71,"(es22.15)") h2osoi_ice(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") swc(:,:) - CLOSE(71) - END IF - + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + + IF(clmupdate_swc.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn3, action="write") + WRITE (71,"(es22.15)") h2osoi_liq(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") swc(:,:) + CLOSE(71) END IF + + END IF #endif endif @@ -746,73 +796,111 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 if ( clmupdate_snow_repartitioning.eq.3) then incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do endif endif endif end do if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:))).gt.0.000001) then - call clm_repartition_snow() - end if + if ( ABS(SUM(rsnow(:))).gt.0.000001) then + call clm_repartition_snow() + end if end if endif ! Case 2: Snow water equivalent ! Write updated snow depth back to CLM and then repartition snow and adjust related variables if(clmupdate_snow.eq.2) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if - - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if - end if - endif - end do - ! cc = cc + 1 - ! end do + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif + + end if + endif + end do + ! cc = cc + 1 + ! end do - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + endif + ! Case 3: Snow depth with h2osoi_ice in state vector + ! Write updated h2osoi_ice back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.eq.3) then + do i=-nlevsno+1,nlevsoi + do j=clm_begc,clm_endc + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + (i + nlevsno -1)*(clm_endg - clm_begg + 1) + else + cc = (j - clm_begc + 1) + (i + nlevsno -1)*(clm_endc - clm_begc + 1) + end if + ! Catch negative or 0 values from DA + if (i .eq. -nlevsno+1) then + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow_depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) + end if end if + if (clm_statevec(cc+clm_varsize+offset).lt.0.0) then + print *, "WARNING: h2osoi_ice at g,c is negative: ", cc, j, i, clm_statevec(cc+clm_varsize+offset) + else + ! UPDATE + h2osoi_ice(j,i) = clm_statevec(cc+clm_varsize+offset) + endif + end do + end do + + ! How to interact with repartitioning? + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:))).gt.0.000001) then + call clm_repartition_snow() end if + end if endif end subroutine update_clm From ddad8394bab5b50744f22a1830deed7137408a68 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Mar 2025 16:44:39 +0100 Subject: [PATCH 05/32] snow: update eclm --- interface/model/eclm/enkf_clm_mod_5.F90 | 186 ++++++++++++------------ 1 file changed, 93 insertions(+), 93 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 543f4494b..a661072b8 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -627,113 +627,113 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! write updated swc back to CLM if(clmupdate_swc.ne.0) then - ! Set minimum soil moisture for checking the state vector and - ! for setting minimum swc for CLM - if(clmwatmin_switch.eq.3) then - ! CLM3.5 type watmin - watmin_check = 0.00 - watmin_set = 0.05 - else if(clmwatmin_switch.eq.5) then - ! CLM5.0 type watmin - watmin_check = watmin - watmin_set = watmin - else - ! Default - watmin_check = 0.0 - watmin_set = 0.0 - end if - - ! cc = 0 - do i=1,nlevsoi - ! CLM3.5: iterate over grid cells - ! CLM5.0: iterate over columns - ! do j=clm_begg,clm_endg - do j=clm_begc,clm_endc - - ! Update only those SWCs that are not excluded by ispval - if(state_clm2pdaf_p(j,i) .ne. ispval) then - - if(swc(j,i).eq.0.0) then - swc_zero_before_update = .true. + ! Set minimum soil moisture for checking the state vector and + ! for setting minimum swc for CLM + if(clmwatmin_switch.eq.3) then + ! CLM3.5 type watmin + watmin_check = 0.00 + watmin_set = 0.05 + else if(clmwatmin_switch.eq.5) then + ! CLM5.0 type watmin + watmin_check = watmin + watmin_set = watmin + else + ! Default + watmin_check = 0.0 + watmin_set = 0.0 + end if - ! Zero-SWC leads to zero denominator in computation of - ! rliq/rice, therefore setting rliq/rice to special - ! value - rliq = spval - rice = spval - else - swc_zero_before_update = .false. + ! cc = 0 + do i=1,nlevsoi + ! CLM3.5: iterate over grid cells + ! CLM5.0: iterate over columns + ! do j=clm_begg,clm_endg + do j=clm_begc,clm_endc + + ! Update only those SWCs that are not excluded by ispval + if(state_clm2pdaf_p(j,i) .ne. ispval) then + + if(swc(j,i).eq.0.0) then + swc_zero_before_update = .true. + + ! Zero-SWC leads to zero denominator in computation of + ! rliq/rice, therefore setting rliq/rice to special + ! value + rliq = spval + rice = spval + else + swc_zero_before_update = .false. - rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) - rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) - !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) - end if + rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) + rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) + !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) + end if - swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) - if(swc_update.le.watmin_check) then - swc(j,i) = watmin_set - else if(swc_update.ge.watsat(j,i)) then - swc(j,i) = watsat(j,i) - else - swc(j,i) = swc_update - endif + if(swc_update.le.watmin_check) then + swc(j,i) = watmin_set + else if(swc_update.ge.watsat(j,i)) then + swc(j,i) = watsat(j,i) + else + swc(j,i) = swc_update + endif - if (isnan(swc(j,i))) then - swc(j,i) = watmin_set - print *, "WARNING: swc at j,i is nan: ", j, i - endif + if (isnan(swc(j,i))) then + swc(j,i) = watmin_set + print *, "WARNING: swc at j,i is nan: ", j, i + endif - if(swc_zero_before_update) then - ! This case should not appear for hydrologically - ! active columns/layers, where always: swc > watmin - ! - ! If you want to make sure that no zero SWCs appear in - ! the code, comment out the error stop + if(swc_zero_before_update) then + ! This case should not appear for hydrologically + ! active columns/layers, where always: swc > watmin + ! + ! If you want to make sure that no zero SWCs appear in + ! the code, comment out the error stop #ifdef PDAF_DEBUG - ! error stop "ERROR: Update of zero-swc" - print *, "WARNING: Update of zero-swc" - print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i + ! error stop "ERROR: Update of zero-swc" + print *, "WARNING: Update of zero-swc" + print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i #endif - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o - h2osoi_ice(j,i) = 0.0 - else - ! update liquid water content - h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq - ! update ice content - h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice - end if + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o + h2osoi_ice(j,i) = 0.0 + else + ! update liquid water content + h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq + ! update ice content + h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice + end if - end if - ! cc = cc + 1 + end if + ! cc = cc + 1 + end do end do - end do #ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN - - IF(clmupdate_swc.NE.0) THEN - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn3, action="write") - WRITE (71,"(es22.15)") h2osoi_liq(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn4, action="write") - WRITE (71,"(es22.15)") h2osoi_ice(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") swc(:,:) - CLOSE(71) - END IF + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + + IF(clmupdate_swc.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn3, action="write") + WRITE (71,"(es22.15)") h2osoi_liq(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") swc(:,:) + CLOSE(71) + END IF - END IF + END IF #endif endif From 8654093d1984bb08b1c51898e06956550651b7ee Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Mar 2025 16:44:53 +0100 Subject: [PATCH 06/32] clmupdate_snow.eq.3: snow_depth and h2osno in state vector Currently: - snow_depth observations - h2osno updated in state vector through correlations - updated h2osno used for updating h2osoi_ice in CLM - repartitioning: use same as h2osno-state-vector --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 89 ++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 229 +++++++++++----------- 2 files changed, 200 insertions(+), 118 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 172c86297..cf7fff0e9 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -285,6 +285,11 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif + ! Case 3: Assimilation of snow depth: adding swe in the state vector + if(clmupdate_snow.eq.3) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -458,6 +463,29 @@ subroutine set_clm_statevec(tstartcycle, mype) end if endif + ! Case 3: Snow Depth with swe in state vector + if(clmupdate_snow.eq.3) then + + ! snow_depth and swe into state vector + cc = 1 + do j=clm_begg,clm_endg + ! Only get the snow_depth/swe from the first column of each gridcell + ! and add it to the clm_statevec at the position of the gridcell (cc) + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + endif + endif + end do + cc = cc + 1 + end do + + endif #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -815,6 +843,63 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end if end if endif + ! Case 3: Snow depth with swe in state vector + ! Use updated swe (from snow_depth observations) to update h2osoi_ice by increment + if(clmupdate_snow.eq.3) then + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc + + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) + end if + if (clm_statevec(cc+clm_varsize+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+clm_varsize+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+clm_varsize+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+clm_varsize+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if + + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif + + end if + endif + end do + ! cc = cc + 1 + ! end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + endif end subroutine update_clm @@ -898,7 +983,7 @@ subroutine clm_repartition_snow(h2osno_in) h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2) then + else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already h2osno_po(jj) = h2osno(jj) h2osno_pr(jj) = h2osno_in(jj) @@ -986,7 +1071,7 @@ subroutine clm_repartition_snow(h2osno_in) ! Finally adjust SWE (h2osno) since the prior value is no longer needed ! column level variables so we can adjust it outside the layer loop h2osno(jj) = h2osno_po(jj) - else if (clmupdate_snow .eq. 2) then + else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then ! Update the total snow depth to match udpates to layers for active snow layers if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index a661072b8..cf7fff0e9 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -87,7 +87,7 @@ module enkf_clm_mod subroutine define_clm_statevec(mype) use shr_kind_mod, only: r8 => shr_kind_r8 use decompMod , only : get_proc_bounds - use clm_varpar , only : nlevsoi, nlevsno, nlevgrnd + use clm_varpar , only : nlevsoi use clm_varcon , only : ispval use ColumnType , only : col @@ -285,10 +285,10 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif - ! Case 3: Assimilation of snow depth: adding h2osoi_ice in the state vector + ! Case 3: Assimilation of snow depth: adding swe in the state vector if(clmupdate_snow.eq.3) then clm_varsize = (clm_endg-clm_begg+1) - clm_statevecsize = 2*((clm_endg-clm_begg+1) * (nlevsno + nlevgrnd)) !2D h2osoi_ice + clm_statevecsize = 2*(clm_endg-clm_begg+1) endif !hcp LST DA @@ -331,7 +331,7 @@ end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst - use clm_varpar , only : nlevsoi, nlevsno, nlevgrnd + use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col @@ -345,7 +345,6 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: porgm(:,:) real(r8), pointer :: snow_depth(:) real(r8), pointer :: h2osno(:) - real(r8), pointer :: h2oice(:,:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -357,7 +356,6 @@ subroutine set_clm_statevec(tstartcycle, mype) snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) - h2osoi_ice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -465,13 +463,13 @@ subroutine set_clm_statevec(tstartcycle, mype) end if endif - ! Case 3: Snow Depth with h2osoi_ice in state vector + ! Case 3: Snow Depth with swe in state vector if(clmupdate_snow.eq.3) then - ! snow_depth into state vector + ! snow_depth and swe into state vector cc = 1 do j=clm_begg,clm_endg - ! Only get the snow_depth from the first column of each gridcell + ! Only get the snow_depth/swe from the first column of each gridcell ! and add it to the clm_statevec at the position of the gridcell (cc) newgridcell = .true. do jj=clm_begc,clm_endc @@ -480,33 +478,13 @@ subroutine set_clm_statevec(tstartcycle, mype) if (newgridcell) then newgridcell = .false. clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) endif endif end do cc = cc + 1 end do - ! h2osoi_ice into state vector - cc = 1 - ! Additional loop over layers - do i=-nlevsno+1,nlevgrnd - do j=clm_begg,clm_endg - ! Only get the h2osoi_ice from the first column of each gridcell - ! and add it to the clm_statevec at the position of the gridcell (cc) - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - ! Shift by clm_varsize (size of snow-depth in state vector) - clm_statevec(cc+clm_varsize+offset) = h2osoi_ice(jj,i) - endif - endif - end do - cc = cc + 1 - end do - end do endif #ifdef PDAF_DEBUG @@ -797,111 +775,130 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 if ( clmupdate_snow_repartitioning.eq.3) then incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do endif endif endif end do if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:))).gt.0.000001) then - call clm_repartition_snow() - end if + if ( ABS(SUM(rsnow(:))).gt.0.000001) then + call clm_repartition_snow() + end if end if endif ! Case 2: Snow water equivalent ! Write updated snow depth back to CLM and then repartition snow and adjust related variables if(clmupdate_snow.eq.2) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if + + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + end if + endif + end do + ! cc = cc + 1 + ! end do + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) end if - endif - end do - ! cc = cc + 1 - ! end do - - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) end if - end if endif - ! Case 3: Snow depth with h2osoi_ice in state vector - ! Write updated h2osoi_ice back to CLM and then repartition snow and adjust related variables + ! Case 3: Snow depth with swe in state vector + ! Use updated swe (from snow_depth observations) to update h2osoi_ice by increment if(clmupdate_snow.eq.3) then - do i=-nlevsno+1,nlevsoi - do j=clm_begc,clm_endc - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) + (i + nlevsno -1)*(clm_endg - clm_begg + 1) - else - cc = (j - clm_begc + 1) + (i + nlevsno -1)*(clm_endc - clm_begc + 1) - end if - ! Catch negative or 0 values from DA - if (i .eq. -nlevsno+1) then - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow_depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) - end if - end if - if (clm_statevec(cc+clm_varsize+offset).lt.0.0) then - print *, "WARNING: h2osoi_ice at g,c is negative: ", cc, j, i, clm_statevec(cc+clm_varsize+offset) - else - ! UPDATE - h2osoi_ice(j,i) = clm_statevec(cc+clm_varsize+offset) - endif - end do - end do + ! cc = 1 + ! do j=clm_begg,clm_endg + ! iterate through the columns and copy from the same gridcell + ! i.e. statevec position (cc) for each column + do j=clm_begc,clm_endc - ! How to interact with repartitioning? - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:))).gt.0.000001) then - call clm_repartition_snow() + ! Set cc (the state vector index) from the + ! CLM5-grid-index and the `CLM5-layer-index times + ! num_gridcells` + if(clmstatevec_allcol.eq.0) then + cc = (col%gridcell(j) - clm_begg + 1) + else + cc = (j - clm_begc + 1) + end if + + ! Catch negative or 0 values from DA + if (clm_statevec(cc+offset).lt.0.0) then + print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) + end if + if (clm_statevec(cc+clm_varsize+offset).le.0.0) then + print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+clm_varsize+offset) + else + rsnow(j) = h2osno(j) + if ( ABS(SUM(rsnow(:) - clm_statevec(cc+clm_varsize+offset))).gt.0.000001) then + h2osno(j) = clm_statevec(cc+clm_varsize+offset) + ! JK: clmupdate_snow_repartitioning.eq.3 is experimental + ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) + ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 + if ( clmupdate_snow_repartitioning.eq.3) then + incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if + + if (isnan(rsnow(j))) then + print *, "WARNING: rsnow at j is nan: ", j + endif + if (isnan(h2osno(j))) then + print *, "WARNING: h2osno at j is nan: ", j + endif + + end if + endif + end do + ! cc = cc + 1 + ! end do + + if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then + if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if end if - end if endif end subroutine update_clm @@ -986,7 +983,7 @@ subroutine clm_repartition_snow(h2osno_in) h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2) then + else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already h2osno_po(jj) = h2osno(jj) h2osno_pr(jj) = h2osno_in(jj) @@ -1074,7 +1071,7 @@ subroutine clm_repartition_snow(h2osno_in) ! Finally adjust SWE (h2osno) since the prior value is no longer needed ! column level variables so we can adjust it outside the layer loop h2osno(jj) = h2osno_po(jj) - else if (clmupdate_snow .eq. 2) then + else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then ! Update the total snow depth to match udpates to layers for active snow layers if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) From 7c1a2c0b87ccf134ff796ef4862da90e04caa829 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 17 Mar 2025 07:09:41 +0100 Subject: [PATCH 07/32] Snow-DA: Re-introduce standard "point-obs" observation operator --- interface/framework/obs_op_pdaf.F90 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 172ec7d6a..97becfc66 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -320,14 +320,9 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) if(lpointobs) then - ! WORKAROUND FOR SUDDEN SEG 11 error - ! DO NOT LEAVE IN FINAL VERSION - ! JUST DURING DEBUG PROCESS - ! JUST VALID FOR SINGLE GRID CELLS WITH SINGLE VALUE OBS - m_state_p(1) = state_p(1) -! DO i = 1, dim_obs_p -! m_state_p(i) = state_p(obs_index_p(i)) -! END DO + DO i = 1, dim_obs_p + m_state_p(i) = state_p(obs_index_p(i)) + END DO end if From 12ce4fff41cba2909ba98a7869a91075f7cb0922 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 17 Mar 2025 11:28:01 +0100 Subject: [PATCH 08/32] Snow-DA: Refactoring to prepare for PR into `master` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 275 ++++++---------------- interface/model/eclm/enkf_clm_mod_5.F90 | 275 ++++++---------------- 2 files changed, 150 insertions(+), 400 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index cf7fff0e9..9490f3561 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -276,16 +276,13 @@ subroutine define_clm_statevec(mype) ! Snow assimilation ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell - if(clmupdate_snow.eq.1) then - clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA - clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority - endif ! Case 2: Assimilation of snow water equivalent same as above - if(clmupdate_snow.eq.2) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif - ! Case 3: Assimilation of snow depth: adding swe in the state vector + ! Case 3: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector if(clmupdate_snow.eq.3) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) @@ -410,63 +407,11 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif - ! Snow assimilation + ! Snow assimilation state vector ! Case 1: Snow depth - if(clmupdate_snow.eq.1) then - cc = 1 - do j=clm_begg,clm_endg - ! Only get the snow_depth from the first column of each gridcell - ! and add it to the clm_statevec at the position of the gridcell (cc) - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = snow_depth(jj) - endif - endif - end do - cc = cc + 1 - end do - endif ! Case 2: SWE - if(clmupdate_snow.eq.2) then - cc = 1 - - if(clmstatevec_allcol.eq.0) then - - do j=clm_begg,clm_endg - ! Only get the SWE from the first column of each gridcell - ! and add it to the clm_statevec at the position of the gridcell (cc) - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = h2osno(jj) - endif - endif - end do - cc = cc + 1 - end do - - else - - do jj=clm_begc,clm_endc - ! SWC from all columns of each gridcell - clm_statevec(cc+offset) = h2osno(jj) - cc = cc + 1 - end do - - end if - - endif - ! Case 3: Snow Depth with swe in state vector - if(clmupdate_snow.eq.3) then - - ! snow_depth and swe into state vector + ! Case 3: Snow depth + SWE + if(clmupdate_snow.ne.0) then cc = 1 do j=clm_begg,clm_endg ! Only get the snow_depth/swe from the first column of each gridcell @@ -477,8 +422,18 @@ subroutine set_clm_statevec(tstartcycle, mype) if (g .eq. j) then if (newgridcell) then newgridcell = .false. - clm_statevec(cc+offset) = snow_depth(jj) - clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + + if(clmupdate_snow.eq.1) then + clm_statevec(cc+offset) = snow_depth(jj) + else if(clmupdate_snow.eq.2) then + clm_statevec(cc+offset) = h2osno(jj) + else if(clmupdate_snow.eq.3) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else + error stop "Wrong input for clmupdate_snow" + end if + endif endif end do @@ -531,6 +486,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: rliq,rice,incr_h2osno real(r8) :: rsnow(clm_statevecsize) + real(r8) :: nsnow(clm_statevecsize) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -749,156 +705,75 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Snow assimilation: ! Case 1: Snow depth - ! Write updated snow depth back to CLM and then repartition snow and adjust related variables - if(clmupdate_snow.eq.1) then + ! Case 2: Snow water equivalent + ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector + ! Write updated snow variable back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.ne.0) then do j=clm_begc,clm_endc - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column + ! Iterate through the columns + + ! For each column index, copy the gridcell-linked state vector + ! value at state vector index `cc` ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + ! CLM5-grid-index + cc = (col%gridcell(j) - clm_begg + 1) + ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) - else - rsnow(j) = snow_depth(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))) .gt.0.000001) then - snow_depth(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - endif - endif - endif - end do + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) + end if - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:))).gt.0.000001) then - call clm_repartition_snow() - end if - end if - endif - ! Case 2: Snow water equivalent - ! Write updated snow depth back to CLM and then repartition snow and adjust related variables - if(clmupdate_snow.eq.2) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc - - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + if (clmupdate_snow.eq.1) then + rsnow(j) = snow_depth(j) + else if (clmupdate_snow.eq.2) then + rsnow(j) = h2osno(j) + else if (clmupdate_snow.eq.3) then + rsnow(j) = h2osno(j) + end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if - - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + nsnow(j) = clm_statevec(cc+offset) + else if (clmupdate_snow.eq.3) then + nsnow(j) = clm_statevec(cc+clm_varsize+offset) + end if + + ! Update state variable to CLM + ! Not needed for repartioning-case 3? + if(clmupdate_snow.eq.1) then + snow_depth(j) = nsnow(j) + end if + if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + h2osno(j) = nsnow(j) + end if - end if - endif - end do - ! cc = cc + 1 - ! end do + end do - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) - end if + ! Repartitioning + if ( clmupdate_snow_repartitioning.eq.1 .and. clmupdate_snow_repartitioning.eq.2) then + + if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) end if - endif - ! Case 3: Snow depth with swe in state vector - ! Use updated swe (from snow_depth observations) to update h2osoi_ice by increment - if(clmupdate_snow.eq.3) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc - - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) - end if - if (clm_statevec(cc+clm_varsize+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+clm_varsize+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+clm_varsize+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+clm_varsize+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if - - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + end if - end if - endif - end do - ! cc = cc + 1 - ! end do + if ( clmupdate_snow_repartitioning.eq.3) then - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + do j=clm_begc,clm_endc + if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then + if ( ABS(rsnow(j)).gt.0.000001) then + ! Update h2osoi_ice with increment + incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if end if - end if + end do + + end if + endif end subroutine update_clm diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index cf7fff0e9..9490f3561 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -276,16 +276,13 @@ subroutine define_clm_statevec(mype) ! Snow assimilation ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell - if(clmupdate_snow.eq.1) then - clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA - clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority - endif ! Case 2: Assimilation of snow water equivalent same as above - if(clmupdate_snow.eq.2) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif - ! Case 3: Assimilation of snow depth: adding swe in the state vector + ! Case 3: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector if(clmupdate_snow.eq.3) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) @@ -410,63 +407,11 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif - ! Snow assimilation + ! Snow assimilation state vector ! Case 1: Snow depth - if(clmupdate_snow.eq.1) then - cc = 1 - do j=clm_begg,clm_endg - ! Only get the snow_depth from the first column of each gridcell - ! and add it to the clm_statevec at the position of the gridcell (cc) - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = snow_depth(jj) - endif - endif - end do - cc = cc + 1 - end do - endif ! Case 2: SWE - if(clmupdate_snow.eq.2) then - cc = 1 - - if(clmstatevec_allcol.eq.0) then - - do j=clm_begg,clm_endg - ! Only get the SWE from the first column of each gridcell - ! and add it to the clm_statevec at the position of the gridcell (cc) - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = h2osno(jj) - endif - endif - end do - cc = cc + 1 - end do - - else - - do jj=clm_begc,clm_endc - ! SWC from all columns of each gridcell - clm_statevec(cc+offset) = h2osno(jj) - cc = cc + 1 - end do - - end if - - endif - ! Case 3: Snow Depth with swe in state vector - if(clmupdate_snow.eq.3) then - - ! snow_depth and swe into state vector + ! Case 3: Snow depth + SWE + if(clmupdate_snow.ne.0) then cc = 1 do j=clm_begg,clm_endg ! Only get the snow_depth/swe from the first column of each gridcell @@ -477,8 +422,18 @@ subroutine set_clm_statevec(tstartcycle, mype) if (g .eq. j) then if (newgridcell) then newgridcell = .false. - clm_statevec(cc+offset) = snow_depth(jj) - clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + + if(clmupdate_snow.eq.1) then + clm_statevec(cc+offset) = snow_depth(jj) + else if(clmupdate_snow.eq.2) then + clm_statevec(cc+offset) = h2osno(jj) + else if(clmupdate_snow.eq.3) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else + error stop "Wrong input for clmupdate_snow" + end if + endif endif end do @@ -531,6 +486,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: rliq,rice,incr_h2osno real(r8) :: rsnow(clm_statevecsize) + real(r8) :: nsnow(clm_statevecsize) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -749,156 +705,75 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Snow assimilation: ! Case 1: Snow depth - ! Write updated snow depth back to CLM and then repartition snow and adjust related variables - if(clmupdate_snow.eq.1) then + ! Case 2: Snow water equivalent + ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector + ! Write updated snow variable back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.ne.0) then do j=clm_begc,clm_endc - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column + ! Iterate through the columns + + ! For each column index, copy the gridcell-linked state vector + ! value at state vector index `cc` ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + ! CLM5-grid-index + cc = (col%gridcell(j) - clm_begg + 1) + ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) - else - rsnow(j) = snow_depth(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))) .gt.0.000001) then - snow_depth(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = snow_depth(j) / rsnow(j) ! INC = New SD / OLD SD - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - endif - endif - endif - end do + if (clm_statevec(cc+offset).le.0.0) then + print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) + end if - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:))).gt.0.000001) then - call clm_repartition_snow() - end if - end if - endif - ! Case 2: Snow water equivalent - ! Write updated snow depth back to CLM and then repartition snow and adjust related variables - if(clmupdate_snow.eq.2) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc - - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if + if (clmupdate_snow.eq.1) then + rsnow(j) = snow_depth(j) + else if (clmupdate_snow.eq.2) then + rsnow(j) = h2osno(j) + else if (clmupdate_snow.eq.3) then + rsnow(j) = h2osno(j) + end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if - - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + nsnow(j) = clm_statevec(cc+offset) + else if (clmupdate_snow.eq.3) then + nsnow(j) = clm_statevec(cc+clm_varsize+offset) + end if + + ! Update state variable to CLM + ! Not needed for repartioning-case 3? + if(clmupdate_snow.eq.1) then + snow_depth(j) = nsnow(j) + end if + if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + h2osno(j) = nsnow(j) + end if - end if - endif - end do - ! cc = cc + 1 - ! end do + end do - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) - end if + ! Repartitioning + if ( clmupdate_snow_repartitioning.eq.1 .and. clmupdate_snow_repartitioning.eq.2) then + + if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) end if - endif - ! Case 3: Snow depth with swe in state vector - ! Use updated swe (from snow_depth observations) to update h2osoi_ice by increment - if(clmupdate_snow.eq.3) then - ! cc = 1 - ! do j=clm_begg,clm_endg - ! iterate through the columns and copy from the same gridcell - ! i.e. statevec position (cc) for each column - do j=clm_begc,clm_endc - - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.0) then - cc = (col%gridcell(j) - clm_begg + 1) - else - cc = (j - clm_begc + 1) - end if - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).lt.0.0) then - print *, "WARNING: snow depth at g,c is negative: ", cc, j, clm_statevec(cc+offset) - end if - if (clm_statevec(cc+clm_varsize+offset).le.0.0) then - print *, "WARNING: SWE at g,c is negative or zero: ", j, clm_statevec(cc+clm_varsize+offset) - else - rsnow(j) = h2osno(j) - if ( ABS(SUM(rsnow(:) - clm_statevec(cc+clm_varsize+offset))).gt.0.000001) then - h2osno(j) = clm_statevec(cc+clm_varsize+offset) - ! JK: clmupdate_snow_repartitioning.eq.3 is experimental - ! JK: clmupdate_snow_repartitioning.eq.3 from NASA-Code (based on older CLM3.5 version) - ! https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90 - if ( clmupdate_snow_repartitioning.eq.3) then - incr_h2osno = h2osno(j) / rsnow(j) ! INC = New SWE / OLD SWE - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do - end if - - if (isnan(rsnow(j))) then - print *, "WARNING: rsnow at j is nan: ", j - endif - if (isnan(h2osno(j))) then - print *, "WARNING: h2osno at j is nan: ", j - endif + end if - end if - endif - end do - ! cc = cc + 1 - ! end do + if ( clmupdate_snow_repartitioning.eq.3) then - if ( clmupdate_snow_repartitioning.ne.0 .and. clmupdate_snow_repartitioning.ne.3) then - if ( ABS(SUM(rsnow(:) - h2osno(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + do j=clm_begc,clm_endc + if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then + if ( ABS(rsnow(j)).gt.0.000001) then + ! Update h2osoi_ice with increment + incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if end if - end if + end do + + end if + endif end subroutine update_clm From dc4f8c17af4dd22fb0fc101d35c63e8a3df7ea8d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 17 Mar 2025 18:37:13 +0100 Subject: [PATCH 09/32] Snow-DA: adapt dimension of `rsnow`, `nsnow` --- interface/model/eclm/enkf_clm_mod_5.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9490f3561..13bd051b4 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -485,8 +485,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_ice(:,:) real(r8) :: rliq,rice,incr_h2osno - real(r8) :: rsnow(clm_statevecsize) - real(r8) :: nsnow(clm_statevecsize) + real(r8) :: rsnow(clm_begc:clm_endc) + real(r8) :: nsnow(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop From 33252adcbe5f603b6880944dc309d9a56f7e77af Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 17 Mar 2025 18:39:39 +0100 Subject: [PATCH 10/32] Snow-DA: try to fix update conditions --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 46 ++++++++++++----------- interface/model/eclm/enkf_clm_mod_5.F90 | 42 +++++++++++---------- 2 files changed, 46 insertions(+), 42 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 9490f3561..39711a7e3 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -485,8 +485,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_ice(:,:) real(r8) :: rliq,rice,incr_h2osno - real(r8) :: rsnow(clm_statevecsize) - real(r8) :: nsnow(clm_statevecsize) + real(r8) :: rsnow(clm_begc:clm_endc) + real(r8) :: nsnow(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -719,11 +719,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! CLM5-grid-index cc = (col%gridcell(j) - clm_begg + 1) - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) - end if - if (clmupdate_snow.eq.1) then rsnow(j) = snow_depth(j) else if (clmupdate_snow.eq.2) then @@ -737,14 +732,19 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.3) then nsnow(j) = clm_statevec(cc+clm_varsize+offset) end if - - ! Update state variable to CLM - ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) - end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) + + if (nsnow(j).gt.0.0) then + ! Update state variable to CLM + ! Not needed for repartioning-case 3? + if(clmupdate_snow.eq.1) then + snow_depth(j) = nsnow(j) + end if + if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + h2osno(j) = nsnow(j) + end if + else + ! Catch negative or 0 values from DA + print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) end if end do @@ -761,13 +761,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then do j=clm_begc,clm_endc - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.000001) then - ! Update h2osoi_ice with increment - incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do + if (nsnow(j).gt.0.0) then + if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then + if ( rsnow(j).ne.0.0) then + ! Update h2osoi_ice with increment + incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if end if end if end do diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 13bd051b4..39711a7e3 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -719,11 +719,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! CLM5-grid-index cc = (col%gridcell(j) - clm_begg + 1) - ! Catch negative or 0 values from DA - if (clm_statevec(cc+offset).le.0.0) then - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) - end if - if (clmupdate_snow.eq.1) then rsnow(j) = snow_depth(j) else if (clmupdate_snow.eq.2) then @@ -737,14 +732,19 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.3) then nsnow(j) = clm_statevec(cc+clm_varsize+offset) end if - - ! Update state variable to CLM - ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) - end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) + + if (nsnow(j).gt.0.0) then + ! Update state variable to CLM + ! Not needed for repartioning-case 3? + if(clmupdate_snow.eq.1) then + snow_depth(j) = nsnow(j) + end if + if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + h2osno(j) = nsnow(j) + end if + else + ! Catch negative or 0 values from DA + print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) end if end do @@ -761,13 +761,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then do j=clm_begc,clm_endc - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.000001) then - ! Update h2osoi_ice with increment - incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno - end do + if (nsnow(j).gt.0.0) then + if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then + if ( rsnow(j).ne.0.0) then + ! Update h2osoi_ice with increment + incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + end do + end if end if end if end do From f1c1145dc578a206e92ffe6f544e5bae15e0c9eb Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 17 Mar 2025 18:48:28 +0100 Subject: [PATCH 11/32] Snow-DA: debug output of `h2osoi_ice` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 13 +++++++++++++ interface/model/eclm/enkf_clm_mod_5.F90 | 13 +++++++++++++ 2 files changed, 26 insertions(+) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 39711a7e3..78ff66690 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -537,7 +537,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") OPEN(unit=71, file=fn5, action="write") WRITE (71,"(es22.15)") h2osoi_liq(:,:) CLOSE(71) + END IF + IF(clmupdate_swc.NE.0 .OR. clmupdate_snow.NE.0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn6, action="write") @@ -776,6 +778,17 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end if +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + END IF +#endif endif end subroutine update_clm diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 39711a7e3..78ff66690 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -537,7 +537,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") OPEN(unit=71, file=fn5, action="write") WRITE (71,"(es22.15)") h2osoi_liq(:,:) CLOSE(71) + END IF + IF(clmupdate_swc.NE.0 .OR. clmupdate_snow.NE.0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn6, action="write") @@ -776,6 +778,17 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end if +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + END IF +#endif endif end subroutine update_clm From b0385e4a84295103e94ed466460cedb1f4ae434c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Mar 2025 17:06:32 +0100 Subject: [PATCH 12/32] Snow-DA: remove experimental barrier --- interface/framework/pdaf_terrsysmp.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/interface/framework/pdaf_terrsysmp.F90 b/interface/framework/pdaf_terrsysmp.F90 index c28606fa0..4ed09922b 100644 --- a/interface/framework/pdaf_terrsysmp.F90 +++ b/interface/framework/pdaf_terrsysmp.F90 @@ -88,9 +88,6 @@ PROGRAM pdaf_terrsysmp ! forward simulation of component models CALL integrate_tsmp() - ! barrier before model integration starts - CALL MPI_BARRIER(MPI_COMM_WORLD, MPIerr) - ! assimilation step CALL assimilate_pdaf() From 62e0381017b9e744612c5b888ac83ae28a1abf53 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Mar 2025 17:06:45 +0100 Subject: [PATCH 13/32] Snow-DA: No layer variable used for setting `obs_index_p` --- interface/framework/init_dim_obs_f_pdaf.F90 | 8 +++++++- interface/framework/init_dim_obs_pdaf.F90 | 8 +++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index d3f1484d3..a34125f73 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -135,6 +135,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) USE enkf_clm_mod, only: domain_def_clm USE enkf_clm_mod, only: get_interp_idx use enkf_clm_mod, only: clmstatevec_allcol + use enkf_clm_mod, only: clmupdate_snow !hcp end #endif #endif @@ -966,7 +967,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end if #endif else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + if(clmupdate_snow.ne.0) then + ! Snow-DA: no layer in state vector variables + obs_index_p(cnt) = g-begg+1 + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index e1a6ee4e0..beece1221 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -131,6 +131,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) USE enkf_clm_mod, only: domain_def_clm USE enkf_clm_mod, only: get_interp_idx use enkf_clm_mod, only: clmstatevec_allcol + use enkf_clm_mod, only: clmupdate_snow !hcp end #endif #endif @@ -959,7 +960,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if #endif else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + if(clmupdate_snow.ne.0) then + ! Snow-DA: no layer in state vector variables + obs_index_p(cnt) = g-begg+1 + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) From e9113fe75ab499cc3e47b88b26bf619af089e2b1 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Mar 2025 17:15:33 +0100 Subject: [PATCH 14/32] Snow-DA: Default for `CLM:update_snow_repartitioning` to `3` --- interface/model/common/read_enkfpar.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index fcadb26a7..22958091f 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -79,7 +79,7 @@ void read_enkfpar(char *parname) clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); clmupdate_snow = iniparser_getint(pardict,"CLM:update_snow",0); - clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",1); + clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",3); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); From 96a895fe22b5a66d0961dcff1c1f759e7fc278b9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Mar 2025 17:22:17 +0100 Subject: [PATCH 15/32] Snow-DA: CLM3.5 variable declarations and error messages --- interface/model/clm3_5/enkf_clm_mod.F90 | 32 +++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 5a8114de0..7918ce69d 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -65,6 +65,8 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture + integer(c_int),bind(C,name="clmupdate_snow") :: clmupdate_snow + integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et @@ -147,6 +149,19 @@ subroutine define_clm_statevec(mype) error stop "Not implemented: clmupdate_texture.eq.2" endif + ! Snow assimilation + ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 + ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell + ! Case 2: Assimilation of snow water equivalent same as above + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + error stop "Not implemented: clmupdate_snow.eq.1 or clmupdate_snow.eq.1" + endif + ! Case 3: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector + if(clmupdate_snow.eq.3) then + error stop "Not implemented: clmupdate_snow.eq.3" + endif + !hcp LST DA if(clmupdate_T.eq.1) then clm_varsize = endg-begg+1 @@ -287,6 +302,14 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! Snow assimilation state vector + ! Case 1: Snow depth + ! Case 2: SWE + ! Case 3: Snow depth + SWE + if(clmupdate_snow.ne.0) then + error stop "Not implemented: clmupdate_snow.ne.0" + endif + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files @@ -534,6 +557,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif + ! Snow assimilation: + ! Case 1: Snow depth + ! Case 2: Snow water equivalent + ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector + ! Write updated snow variable back to CLM and then repartition snow and adjust related variables + if(clmupdate_snow.ne.0) then + error stop "Not implemented: clmupdate_snow.ne.0" + endif + end subroutine update_clm subroutine clm_correct_texture() From b3f23370836fefac5b0d01d7caf9143caa4596d9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 09:31:23 +0100 Subject: [PATCH 16/32] Snow-DA: bugfix repartitioning condition --- interface/model/eclm/enkf_clm_mod_5.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 78ff66690..2f2d2470d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -752,7 +752,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do ! Repartitioning - if ( clmupdate_snow_repartitioning.eq.1 .and. clmupdate_snow_repartitioning.eq.2) then + if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then call clm_repartition_snow(rsnow(:)) From a75d8418303ae8d2ac4e7259901988de0a526b2c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 09:34:19 +0100 Subject: [PATCH 17/32] Snow-DA: bugfix repartitioning condition II --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 78ff66690..2f2d2470d 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -752,7 +752,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do ! Repartitioning - if ( clmupdate_snow_repartitioning.eq.1 .and. clmupdate_snow_repartitioning.eq.2) then + if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then call clm_repartition_snow(rsnow(:)) From 1d16b75add06bec5a1d514408cf3dbd98bd7a25e Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 09:36:18 +0100 Subject: [PATCH 18/32] Snow-DA: numerically better check for `rsnow.ne.0.0` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 2f2d2470d..ca765b79f 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -765,7 +765,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (nsnow(j).gt.0.0) then if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( rsnow(j).ne.0.0) then + if ( ABS(rsnow(j)).gt.0.0) then ! Update h2osoi_ice with increment incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var do i=snlsno(j)+1,0 diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 2f2d2470d..ca765b79f 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -765,7 +765,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (nsnow(j).gt.0.0) then if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( rsnow(j).ne.0.0) then + if ( ABS(rsnow(j)).gt.0.0) then ! Update h2osoi_ice with increment incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var do i=snlsno(j)+1,0 From 2dde41f962b26eeb27562615d22c0925269e2c77 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 10:06:20 +0100 Subject: [PATCH 19/32] Snow-DA: repartitioning function only for `clmupdate_snow.eq.1,2` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 10 ++++++++-- interface/model/eclm/enkf_clm_mod_5.F90 | 12 +++++++++--- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index ca765b79f..7cb5ccc9d 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -754,8 +754,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + + if (clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then + error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.eq.3/4" end if end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index ca765b79f..67d37c34d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -754,8 +754,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then + call clm_repartition_snow(rsnow(:)) + end if + end if + + if (clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then + error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.eq.3/4" end if end if @@ -873,7 +879,7 @@ subroutine clm_repartition_snow(h2osno_in) h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then + else if (clmupdate_snow .eq. 2) then ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already h2osno_po(jj) = h2osno(jj) h2osno_pr(jj) = h2osno_in(jj) From d2ae20d4dc9baa10f17642479a2575c98f8c7028 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 10:19:43 +0100 Subject: [PATCH 20/32] Snow-DA: rename `incr_h2osno` to `incr_sno` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 6 +++--- interface/model/eclm/enkf_clm_mod_5.F90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 7cb5ccc9d..a6eff6c69 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -484,7 +484,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_h2osno + real(r8) :: rliq,rice,incr_sno real(r8) :: rsnow(clm_begc:clm_endc) real(r8) :: nsnow(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) @@ -773,9 +773,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then if ( ABS(rsnow(j)).gt.0.0) then ! Update h2osoi_ice with increment - incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 67d37c34d..f8c3c94e6 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -484,7 +484,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_h2osno + real(r8) :: rliq,rice,incr_sno real(r8) :: rsnow(clm_begc:clm_endc) real(r8) :: nsnow(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) @@ -773,9 +773,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then if ( ABS(rsnow(j)).gt.0.0) then ! Update h2osoi_ice with increment - incr_h2osno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var + incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_h2osno + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if end if From 37a10d097318998106bae93cae6bbcc16d53189f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 27 Mar 2025 16:09:18 +0100 Subject: [PATCH 21/32] Snow-DA: Refactor Snow-DA preparing for introduction of `CLM:update_snow=4` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 101 +++++++++++++++------- interface/model/eclm/enkf_clm_mod_5.F90 | 99 ++++++++++++++------- 2 files changed, 135 insertions(+), 65 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index a6eff6c69..3592696ab 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -282,7 +282,7 @@ subroutine define_clm_statevec(mype) clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif ! Case 3: Assimilation of snow depth: Snow depth and snow water - ! equivalent in the state vector + ! equivalent in the state vector. Update of h2osoi_ice if(clmupdate_snow.eq.3) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) @@ -484,9 +484,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_sno - real(r8) :: rsnow(clm_begc:clm_endc) - real(r8) :: nsnow(clm_begc:clm_endc) + real(r8) :: rliq,rice + real(r8) :: incr_sno + real(r8) :: incr_sd + real(r8) :: incr_swe + real(r8) :: h2osno_in(clm_begc:clm_endc) + real(r8) :: snow_depth_in(clm_begc:clm_endc) + real(r8) :: h2osno_out(clm_begc:clm_endc) + real(r8) :: snow_depth_out(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -722,31 +727,37 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") cc = (col%gridcell(j) - clm_begg + 1) if (clmupdate_snow.eq.1) then - rsnow(j) = snow_depth(j) + snow_depth_in(j) = snow_depth(j) else if (clmupdate_snow.eq.2) then - rsnow(j) = h2osno(j) + h2osno_in(j) = h2osno(j) else if (clmupdate_snow.eq.3) then - rsnow(j) = h2osno(j) + h2osno_in(j) = h2osno(j) end if - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - nsnow(j) = clm_statevec(cc+offset) + if (clmupdate_snow.eq.1) then + snow_depth_out(j) = clm_statevec(cc+offset) + else if (clmupdate_snow.eq.2) then + h2osno_out(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then - nsnow(j) = clm_statevec(cc+clm_varsize+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if (nsnow(j).gt.0.0) then + if(clmupdate_snow.eq.1) then ! Update state variable to CLM ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) + if (snow_depth_out(j).gt.0.0) then + snow_depth(j) = snow_depth_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + if (h2osno_out(j).gt.0.0) then + h2osno(j) = h2osno_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: SWE is negative/zero at cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j) end if - else - ! Catch negative or 0 values from DA - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) end if end do @@ -754,9 +765,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:))).gt.0.000001) then + call clm_repartition_snow(snow_depth_in(:)) + end if + end if + + if (clmupdate_snow.eq.2) then + if ( SUM(ABS(h2osno_in(:) - h2osno_out(:))).gt.0.000001) then + call clm_repartition_snow(h2osno_in(:)) end if end if @@ -768,19 +785,37 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then - do j=clm_begc,clm_endc - if (nsnow(j).gt.0.0) then - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.0) then - ! Update h2osoi_ice with increment - incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno - end do + if (clmupdate_snow.eq.1) then + do j=clm_begc,clm_endc + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if end if end if - end if - end do + end do + end if + + if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if + end if + end if + end do + end if end if @@ -879,7 +914,7 @@ subroutine clm_repartition_snow(h2osno_in) h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then + else if (clmupdate_snow .eq. 2) then ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already h2osno_po(jj) = h2osno(jj) h2osno_pr(jj) = h2osno_in(jj) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index f8c3c94e6..3592696ab 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -282,7 +282,7 @@ subroutine define_clm_statevec(mype) clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif ! Case 3: Assimilation of snow depth: Snow depth and snow water - ! equivalent in the state vector + ! equivalent in the state vector. Update of h2osoi_ice if(clmupdate_snow.eq.3) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) @@ -484,9 +484,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_sno - real(r8) :: rsnow(clm_begc:clm_endc) - real(r8) :: nsnow(clm_begc:clm_endc) + real(r8) :: rliq,rice + real(r8) :: incr_sno + real(r8) :: incr_sd + real(r8) :: incr_swe + real(r8) :: h2osno_in(clm_begc:clm_endc) + real(r8) :: snow_depth_in(clm_begc:clm_endc) + real(r8) :: h2osno_out(clm_begc:clm_endc) + real(r8) :: snow_depth_out(clm_begc:clm_endc) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop @@ -722,31 +727,37 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") cc = (col%gridcell(j) - clm_begg + 1) if (clmupdate_snow.eq.1) then - rsnow(j) = snow_depth(j) + snow_depth_in(j) = snow_depth(j) else if (clmupdate_snow.eq.2) then - rsnow(j) = h2osno(j) + h2osno_in(j) = h2osno(j) else if (clmupdate_snow.eq.3) then - rsnow(j) = h2osno(j) + h2osno_in(j) = h2osno(j) end if - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - nsnow(j) = clm_statevec(cc+offset) + if (clmupdate_snow.eq.1) then + snow_depth_out(j) = clm_statevec(cc+offset) + else if (clmupdate_snow.eq.2) then + h2osno_out(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then - nsnow(j) = clm_statevec(cc+clm_varsize+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if (nsnow(j).gt.0.0) then + if(clmupdate_snow.eq.1) then ! Update state variable to CLM ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) + if (snow_depth_out(j).gt.0.0) then + snow_depth(j) = snow_depth_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + if (h2osno_out(j).gt.0.0) then + h2osno(j) = h2osno_out(j) + else + ! Catch negative or 0 values from DA + print *, "WARNING: SWE is negative/zero at cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j) end if - else - ! Catch negative or 0 values from DA - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) end if end do @@ -754,9 +765,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:))).gt.0.000001) then + call clm_repartition_snow(snow_depth_in(:)) + end if + end if + + if (clmupdate_snow.eq.2) then + if ( SUM(ABS(h2osno_in(:) - h2osno_out(:))).gt.0.000001) then + call clm_repartition_snow(h2osno_in(:)) end if end if @@ -768,19 +785,37 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then - do j=clm_begc,clm_endc - if (nsnow(j).gt.0.0) then - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.0) then - ! Update h2osoi_ice with increment - incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno - end do + if (clmupdate_snow.eq.1) then + do j=clm_begc,clm_endc + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if end if end if - end if - end do + end do + end if + + if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno + end do + end if + end if + end if + end do + end if end if From e4101a32e5868d507889ad2c715654ee81c407e9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 27 Mar 2025 16:20:03 +0100 Subject: [PATCH 22/32] Snow-DA: Introduce `CLM:update_snow.eq.4` --- interface/model/clm3_5/enkf_clm_mod.F90 | 8 +++- interface/model/clm5_0/enkf_clm_mod_5.F90 | 57 ++++++++++++++++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 57 ++++++++++++++++++++++- 3 files changed, 117 insertions(+), 5 deletions(-) diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7918ce69d..c056f183f 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -157,10 +157,16 @@ subroutine define_clm_statevec(mype) error stop "Not implemented: clmupdate_snow.eq.1 or clmupdate_snow.eq.1" endif ! Case 3: Assimilation of snow depth: Snow depth and snow water - ! equivalent in the state vector + ! equivalent in the state vector. Update of h2osoi_ice if(clmupdate_snow.eq.3) then error stop "Not implemented: clmupdate_snow.eq.3" endif + ! Case 4: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + ! and dz + if(clmupdate_snow.eq.4) then + error stop "Not implemented: clmupdate_snow.eq.4" + endif !hcp LST DA if(clmupdate_T.eq.1) then diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 3592696ab..da15ed1b1 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -287,6 +287,13 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif + ! Case 4: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + ! and dz + if(clmupdate_snow.eq.4) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -430,6 +437,9 @@ subroutine set_clm_statevec(tstartcycle, mype) else if(clmupdate_snow.eq.3) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.4) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) else error stop "Wrong input for clmupdate_snow" end if @@ -732,6 +742,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_in(j) = h2osno(j) else if (clmupdate_snow.eq.3) then h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.4) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) end if if (clmupdate_snow.eq.1) then @@ -740,9 +753,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_out(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.4) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then ! Update state variable to CLM ! Not needed for repartioning-case 3? if (snow_depth_out(j).gt.0.0) then @@ -751,7 +767,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else @@ -817,6 +833,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if + if (clmupdate_snow.eq.4) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end do + end if + end if + end if + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end do + end if + end if + end if + end do + end if + end if #ifdef PDAF_DEBUG diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 3592696ab..da15ed1b1 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -287,6 +287,13 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif + ! Case 4: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + ! and dz + if(clmupdate_snow.eq.4) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -430,6 +437,9 @@ subroutine set_clm_statevec(tstartcycle, mype) else if(clmupdate_snow.eq.3) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.4) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) else error stop "Wrong input for clmupdate_snow" end if @@ -732,6 +742,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_in(j) = h2osno(j) else if (clmupdate_snow.eq.3) then h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.4) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) end if if (clmupdate_snow.eq.1) then @@ -740,9 +753,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_out(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.4) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then ! Update state variable to CLM ! Not needed for repartioning-case 3? if (snow_depth_out(j).gt.0.0) then @@ -751,7 +767,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else @@ -817,6 +833,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if + if (clmupdate_snow.eq.4) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end do + end if + end if + end if + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end do + end if + end if + end if + end do + end if + end if #ifdef PDAF_DEBUG From 664eed33e7bcb3bc1fbcf35adfdf9551d561c207 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 12:51:08 +0100 Subject: [PATCH 23/32] Snow-DA: `CLM:update_snow=5,6,7` intermediate updates between `3` and `4`. --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 162 +++++++++++++++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 162 +++++++++++++++++++++- 2 files changed, 320 insertions(+), 4 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index da15ed1b1..0c034d3c0 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -294,6 +294,26 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif + ! Case 5: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! and dz + if(clmupdate_snow.eq.5) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 6: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + if(clmupdate_snow.eq.6) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 7: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! Should reproduce case 3 + if(clmupdate_snow.eq.6) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -440,6 +460,15 @@ subroutine set_clm_statevec(tstartcycle, mype) else if(clmupdate_snow.eq.4) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.5) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.6) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.7) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) else error stop "Wrong input for clmupdate_snow" end if @@ -745,6 +774,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.4) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.5) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.6) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.7) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) end if if (clmupdate_snow.eq.1) then @@ -756,9 +794,18 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.4) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.5) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.6) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.7) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then ! Update state variable to CLM ! Not needed for repartioning-case 3? if (snow_depth_out(j).gt.0.0) then @@ -767,7 +814,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else @@ -870,6 +917,117 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if + if (clmupdate_snow.eq.5) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + ! if (isnan(h2osoi_liq(j,i))) then + ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + ! endif + end do + end if + end if + end if + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end do + end if + end if + end if + end do + end if + + if (clmupdate_snow.eq.6) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end do + end if + end if + end if + ! if (snow_depth_out(j).gt.0.0) then + ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + ! if (snow_depth_in(j).gt.0.0) then + ! ! Update snow_depth with increment + ! incr_sd = snow_depth_out(j) / snow_depth_in(j) + ! do i=snlsno(j)+1,0 + ! dz(j,i) = dz(j,i) * incr_sd + ! if (isnan(dz(j,i))) then + ! print *, "WARNING: dz at j,i is nan: ", j, i + ! endif + ! end do + ! end if + ! end if + ! end if + end do + end if + + if (clmupdate_snow.eq.7) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + ! if (isnan(h2osoi_liq(j,i))) then + ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + ! endif + end do + end if + end if + end if + ! if (snow_depth_out(j).gt.0.0) then + ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + ! if (snow_depth_in(j).gt.0.0) then + ! ! Update snow_depth with increment + ! incr_sd = snow_depth_out(j) / snow_depth_in(j) + ! do i=snlsno(j)+1,0 + ! dz(j,i) = dz(j,i) * incr_sd + ! if (isnan(dz(j,i))) then + ! print *, "WARNING: dz at j,i is nan: ", j, i + ! endif + ! end do + ! end if + ! end if + ! end if + end do + end if + end if #ifdef PDAF_DEBUG diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index da15ed1b1..0c034d3c0 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -294,6 +294,26 @@ subroutine define_clm_statevec(mype) clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif + ! Case 5: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! and dz + if(clmupdate_snow.eq.5) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 6: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq + if(clmupdate_snow.eq.6) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif + ! Case 7: Assimilation of snow depth: Snow depth and snow water + ! equivalent in the state vector. Update of h2osoi_ice + ! Should reproduce case 3 + if(clmupdate_snow.eq.6) then + clm_varsize = (clm_endg-clm_begg+1) + clm_statevecsize = 2*(clm_endg-clm_begg+1) + endif !hcp LST DA if(clmupdate_T.eq.1) then @@ -440,6 +460,15 @@ subroutine set_clm_statevec(tstartcycle, mype) else if(clmupdate_snow.eq.4) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.5) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.6) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) + else if(clmupdate_snow.eq.7) then + clm_statevec(cc+offset) = snow_depth(jj) + clm_statevec(cc+clm_varsize+offset) = h2osno(jj) else error stop "Wrong input for clmupdate_snow" end if @@ -745,6 +774,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.4) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.5) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.6) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + else if (clmupdate_snow.eq.7) then + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) end if if (clmupdate_snow.eq.1) then @@ -756,9 +794,18 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") else if (clmupdate_snow.eq.4) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.5) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.6) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.7) then + snow_depth_out(j) = clm_statevec(cc+offset) + h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then + if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then ! Update state variable to CLM ! Not needed for repartioning-case 3? if (snow_depth_out(j).gt.0.0) then @@ -767,7 +814,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then + else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else @@ -870,6 +917,117 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if + if (clmupdate_snow.eq.5) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + ! if (isnan(h2osoi_liq(j,i))) then + ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + ! endif + end do + end if + end if + end if + if (snow_depth_out(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + if (snow_depth_in(j).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end do + end if + end if + end if + end do + end if + + if (clmupdate_snow.eq.6) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end do + end if + end if + end if + ! if (snow_depth_out(j).gt.0.0) then + ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + ! if (snow_depth_in(j).gt.0.0) then + ! ! Update snow_depth with increment + ! incr_sd = snow_depth_out(j) / snow_depth_in(j) + ! do i=snlsno(j)+1,0 + ! dz(j,i) = dz(j,i) * incr_sd + ! if (isnan(dz(j,i))) then + ! print *, "WARNING: dz at j,i is nan: ", j, i + ! endif + ! end do + ! end if + ! end if + ! end if + end do + end if + + if (clmupdate_snow.eq.7) then + do j=clm_begc,clm_endc + if (h2osno_out(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then + if (h2osno_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + do i=snlsno(j)+1,0 + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + ! if (isnan(h2osoi_liq(j,i))) then + ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + ! endif + end do + end if + end if + end if + ! if (snow_depth_out(j).gt.0.0) then + ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then + ! if (snow_depth_in(j).gt.0.0) then + ! ! Update snow_depth with increment + ! incr_sd = snow_depth_out(j) / snow_depth_in(j) + ! do i=snlsno(j)+1,0 + ! dz(j,i) = dz(j,i) * incr_sd + ! if (isnan(dz(j,i))) then + ! print *, "WARNING: dz at j,i is nan: ", j, i + ! endif + ! end do + ! end if + ! end if + ! end if + end do + end if + end if #ifdef PDAF_DEBUG From 3a4fda7ba92df10a0b302af7b5a5b6a3c7cb5553 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 14:40:18 +0100 Subject: [PATCH 24/32] Snow-DA: Fix error stop condition --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 4 ++-- interface/model/eclm/enkf_clm_mod_5.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 0c034d3c0..487adf8c5 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -840,8 +840,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end if end if - if (clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then - error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.eq.3/4" + if (clmupdate_snow.ge.3) then + error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3" end if end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 0c034d3c0..487adf8c5 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -840,8 +840,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end if end if - if (clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then - error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.eq.3/4" + if (clmupdate_snow.ge.3) then + error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3" end if end if From 6a87e079c4c8b79fdd456fd92d1698db665c877b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 14:42:59 +0100 Subject: [PATCH 25/32] Snow-DA: No SWE/SD update for `CLM:update_snow >= 3` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 6 +++--- interface/model/eclm/enkf_clm_mod_5.F90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 487adf8c5..c3ac4c84d 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -805,16 +805,16 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + if(clmupdate_snow.eq.1) then ! Update state variable to CLM - ! Not needed for repartioning-case 3? + ! Needed for Case 1/2 if they use repartioning function if (snow_depth_out(j).gt.0.0) then snow_depth(j) = snow_depth_out(j) else ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + else if(clmupdate_snow.eq.2) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 487adf8c5..c3ac4c84d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -805,16 +805,16 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + if(clmupdate_snow.eq.1) then ! Update state variable to CLM - ! Not needed for repartioning-case 3? + ! Needed for Case 1/2 if they use repartioning function if (snow_depth_out(j).gt.0.0) then snow_depth(j) = snow_depth_out(j) else ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + else if(clmupdate_snow.eq.2) then if (h2osno_out(j).gt.0.0) then h2osno(j) = h2osno_out(j) else From 81b1a6a7caf2d70d0d47e043c9f7ee08f27ee459 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 14:47:10 +0100 Subject: [PATCH 26/32] Snow-DA: Remove difference-check for increment-updates --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 20 -------------------- interface/model/eclm/enkf_clm_mod_5.F90 | 20 -------------------- 2 files changed, 40 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index c3ac4c84d..751683906 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -851,7 +851,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.1) then do j=clm_begc,clm_endc if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = snow_depth_out(j) / snow_depth_in(j) @@ -859,7 +858,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if - end if end if end do end if @@ -867,7 +865,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = h2osno_out(j) / h2osno_in(j) @@ -875,7 +872,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if - end if end if end do end if @@ -883,7 +879,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.4) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -898,10 +893,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update snow_depth with increment incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -912,7 +905,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if end do end if @@ -920,7 +912,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.5) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -935,10 +926,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif end do end if - end if end if if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update snow_depth with increment incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -949,7 +938,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if end do end if @@ -957,7 +945,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.6) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -972,10 +959,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if ! if (snow_depth_out(j).gt.0.0) then - ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then ! if (snow_depth_in(j).gt.0.0) then ! ! Update snow_depth with increment ! incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -986,7 +971,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif ! end do ! end if - ! end if ! end if end do end if @@ -994,7 +978,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.7) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -1009,10 +992,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif end do end if - end if end if ! if (snow_depth_out(j).gt.0.0) then - ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then ! if (snow_depth_in(j).gt.0.0) then ! ! Update snow_depth with increment ! incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -1023,7 +1004,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif ! end do ! end if - ! end if ! end if end do end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index c3ac4c84d..751683906 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -851,7 +851,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.1) then do j=clm_begc,clm_endc if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = snow_depth_out(j) / snow_depth_in(j) @@ -859,7 +858,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if - end if end if end do end if @@ -867,7 +865,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = h2osno_out(j) / h2osno_in(j) @@ -875,7 +872,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno end do end if - end if end if end do end if @@ -883,7 +879,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.4) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -898,10 +893,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update snow_depth with increment incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -912,7 +905,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if end do end if @@ -920,7 +912,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.5) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -935,10 +926,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif end do end if - end if end if if (snow_depth_out(j).gt.0.0) then - if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then if (snow_depth_in(j).gt.0.0) then ! Update snow_depth with increment incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -949,7 +938,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if end do end if @@ -957,7 +945,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.6) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -972,10 +959,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif end do end if - end if end if ! if (snow_depth_out(j).gt.0.0) then - ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then ! if (snow_depth_in(j).gt.0.0) then ! ! Update snow_depth with increment ! incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -986,7 +971,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif ! end do ! end if - ! end if ! end if end do end if @@ -994,7 +978,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.7) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then if (h2osno_in(j).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) @@ -1009,10 +992,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif end do end if - end if end if ! if (snow_depth_out(j).gt.0.0) then - ! if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then ! if (snow_depth_in(j).gt.0.0) then ! ! Update snow_depth with increment ! incr_sd = snow_depth_out(j) / snow_depth_in(j) @@ -1023,7 +1004,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! endif ! end do ! end if - ! end if ! end if end do end if From 306688be85d1a67897859e39755c38d1e6d2f90a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 14:58:18 +0100 Subject: [PATCH 27/32] Snow-DA: `CLM:update_snow=4..7`, check both SD and SWE update with both increments or don't update at all --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 159 +++++----------------- interface/model/eclm/enkf_clm_mod_5.F90 | 159 +++++----------------- 2 files changed, 68 insertions(+), 250 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 751683906..7f1574691 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -876,135 +876,44 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if - if (clmupdate_snow.eq.4) then + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - if (isnan(h2osoi_liq(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - end do - end if - end if - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then - ! Update snow_depth with increment - incr_sd = snow_depth_out(j) / snow_depth_in(j) - do i=snlsno(j)+1,0 - dz(j,i) = dz(j,i) * incr_sd - if (isnan(dz(j,i))) then - print *, "WARNING: dz at j,i is nan: ", j, i - endif - end do - end if - end if - end do - end if - - if (clmupdate_snow.eq.5) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - ! if (isnan(h2osoi_liq(j,i))) then - ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - ! endif - end do - end if - end if - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then - ! Update snow_depth with increment - incr_sd = snow_depth_out(j) / snow_depth_in(j) - do i=snlsno(j)+1,0 - dz(j,i) = dz(j,i) * incr_sd - if (isnan(dz(j,i))) then - print *, "WARNING: dz at j,i is nan: ", j, i - endif - end do - end if - end if - end do - end if - - if (clmupdate_snow.eq.6) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - if (isnan(h2osoi_liq(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - end do - end if - end if - ! if (snow_depth_out(j).gt.0.0) then - ! if (snow_depth_in(j).gt.0.0) then - ! ! Update snow_depth with increment - ! incr_sd = snow_depth_out(j) / snow_depth_in(j) - ! do i=snlsno(j)+1,0 - ! dz(j,i) = dz(j,i) * incr_sd - ! if (isnan(dz(j,i))) then - ! print *, "WARNING: dz at j,i is nan: ", j, i - ! endif - ! end do - ! end if - ! end if - end do - end if - - if (clmupdate_snow.eq.7) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - ! if (isnan(h2osoi_liq(j,i))) then - ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - ! endif - end do + if (h2osno_in(j).gt.0.0) then + if (snow_depth_out(j).gt.0.0) then + if (snow_depth_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.6) then + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5) then + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end if + + end do + end if end if + end if end if - ! if (snow_depth_out(j).gt.0.0) then - ! if (snow_depth_in(j).gt.0.0) then - ! ! Update snow_depth with increment - ! incr_sd = snow_depth_out(j) / snow_depth_in(j) - ! do i=snlsno(j)+1,0 - ! dz(j,i) = dz(j,i) * incr_sd - ! if (isnan(dz(j,i))) then - ! print *, "WARNING: dz at j,i is nan: ", j, i - ! endif - ! end do - ! end if - ! end if end do end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 751683906..7f1574691 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -876,135 +876,44 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if - if (clmupdate_snow.eq.4) then + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then do j=clm_begc,clm_endc if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - if (isnan(h2osoi_liq(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - end do - end if - end if - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then - ! Update snow_depth with increment - incr_sd = snow_depth_out(j) / snow_depth_in(j) - do i=snlsno(j)+1,0 - dz(j,i) = dz(j,i) * incr_sd - if (isnan(dz(j,i))) then - print *, "WARNING: dz at j,i is nan: ", j, i - endif - end do - end if - end if - end do - end if - - if (clmupdate_snow.eq.5) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - ! if (isnan(h2osoi_liq(j,i))) then - ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - ! endif - end do - end if - end if - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then - ! Update snow_depth with increment - incr_sd = snow_depth_out(j) / snow_depth_in(j) - do i=snlsno(j)+1,0 - dz(j,i) = dz(j,i) * incr_sd - if (isnan(dz(j,i))) then - print *, "WARNING: dz at j,i is nan: ", j, i - endif - end do - end if - end if - end do - end if - - if (clmupdate_snow.eq.6) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - if (isnan(h2osoi_liq(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - end do - end if - end if - ! if (snow_depth_out(j).gt.0.0) then - ! if (snow_depth_in(j).gt.0.0) then - ! ! Update snow_depth with increment - ! incr_sd = snow_depth_out(j) / snow_depth_in(j) - ! do i=snlsno(j)+1,0 - ! dz(j,i) = dz(j,i) * incr_sd - ! if (isnan(dz(j,i))) then - ! print *, "WARNING: dz at j,i is nan: ", j, i - ! endif - ! end do - ! end if - ! end if - end do - end if - - if (clmupdate_snow.eq.7) then - do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - ! Update h2osoi_ice/h2osoi_liq with increment - incr_swe = h2osno_out(j) / h2osno_in(j) - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe - ! h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe - if (isnan(h2osoi_ice(j,i))) then - print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - endif - ! if (isnan(h2osoi_liq(j,i))) then - ! print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i - ! endif - end do + if (h2osno_in(j).gt.0.0) then + if (snow_depth_out(j).gt.0.0) then + if (snow_depth_in(j).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno_out(j) / h2osno_in(j) + ! Update snow_depth with increment + incr_sd = snow_depth_out(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe + if (isnan(h2osoi_ice(j,i))) then + print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.6) then + h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe + if (isnan(h2osoi_liq(j,i))) then + print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i + endif + end if + + if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5) then + dz(j,i) = dz(j,i) * incr_sd + if (isnan(dz(j,i))) then + print *, "WARNING: dz at j,i is nan: ", j, i + endif + end if + + end do + end if end if + end if end if - ! if (snow_depth_out(j).gt.0.0) then - ! if (snow_depth_in(j).gt.0.0) then - ! ! Update snow_depth with increment - ! incr_sd = snow_depth_out(j) / snow_depth_in(j) - ! do i=snlsno(j)+1,0 - ! dz(j,i) = dz(j,i) * incr_sd - ! if (isnan(dz(j,i))) then - ! print *, "WARNING: dz at j,i is nan: ", j, i - ! endif - ! end do - ! end if - ! end if end do end if From f435467a1a171c8bb10ef936ab10ee73daace059 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Mar 2025 14:59:57 +0100 Subject: [PATCH 28/32] Snow-DA: Rule out very small SD/SWE for updates --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 16 ++++++++-------- interface/model/eclm/enkf_clm_mod_5.F90 | 16 ++++++++-------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 7f1574691..9ef496877 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -850,8 +850,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.1) then do j=clm_begc,clm_endc - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then + if (snow_depth_out(j).gt.1.0d-6) then + if (snow_depth_in(j).gt.1.0d-6) then ! Update h2osoi_ice with increment incr_sno = snow_depth_out(j) / snow_depth_in(j) do i=snlsno(j)+1,0 @@ -864,8 +864,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then + if (h2osno_out(j).gt.1.0d-6) then + if (h2osno_in(j).gt.1.0d-6) then ! Update h2osoi_ice with increment incr_sno = h2osno_out(j) / h2osno_in(j) do i=snlsno(j)+1,0 @@ -878,10 +878,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then + if (h2osno_out(j).gt.1.0d-6) then + if (h2osno_in(j).gt.1.0d-6) then + if (snow_depth_out(j).gt.1.0d-6) then + if (snow_depth_in(j).gt.1.0d-6) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) ! Update snow_depth with increment diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 7f1574691..9ef496877 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -850,8 +850,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.1) then do j=clm_begc,clm_endc - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then + if (snow_depth_out(j).gt.1.0d-6) then + if (snow_depth_in(j).gt.1.0d-6) then ! Update h2osoi_ice with increment incr_sno = snow_depth_out(j) / snow_depth_in(j) do i=snlsno(j)+1,0 @@ -864,8 +864,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then + if (h2osno_out(j).gt.1.0d-6) then + if (h2osno_in(j).gt.1.0d-6) then ! Update h2osoi_ice with increment incr_sno = h2osno_out(j) / h2osno_in(j) do i=snlsno(j)+1,0 @@ -878,10 +878,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.0.0) then - if (h2osno_in(j).gt.0.0) then - if (snow_depth_out(j).gt.0.0) then - if (snow_depth_in(j).gt.0.0) then + if (h2osno_out(j).gt.1.0d-6) then + if (h2osno_in(j).gt.1.0d-6) then + if (snow_depth_out(j).gt.1.0d-6) then + if (snow_depth_in(j).gt.1.0d-6) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) ! Update snow_depth with increment From 83588932905eaec40259da9f22ae68db6e4fc2c1 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Sat, 29 Mar 2025 07:05:14 +0100 Subject: [PATCH 29/32] Snow-DA: Bugfix for `CLM:update_snow=7` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 9ef496877..58611cd7b 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -310,7 +310,7 @@ subroutine define_clm_statevec(mype) ! Case 7: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice ! Should reproduce case 3 - if(clmupdate_snow.eq.6) then + if(clmupdate_snow.eq.7) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9ef496877..58611cd7b 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -310,7 +310,7 @@ subroutine define_clm_statevec(mype) ! Case 7: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice ! Should reproduce case 3 - if(clmupdate_snow.eq.6) then + if(clmupdate_snow.eq.7) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif From 4b0c24f37e24ed798672f6d17d97d46acd31001d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 31 Oct 2025 13:15:22 +0100 Subject: [PATCH 30/32] fortitude style updates --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/framework/init_dim_obs_pdaf.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 171 ++++++++++---------- 3 files changed, 89 insertions(+), 86 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index d3e865c4e..bc6eddad2 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -987,7 +987,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) if(clmstatevec_only_active==1) then obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) else - if(clmupdate_snow.ne.0) then + if(clmupdate_snow/=0) then ! Snow-DA: no layer in state vector variables obs_index_p(cnt) = g-begg+1 else diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index f9736df65..aea74e227 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -980,7 +980,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) if(clmstatevec_only_active==1) then obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) else - if(clmupdate_snow.ne.0) then + if(clmupdate_snow/=0) then ! Snow-DA: no layer in state vector variables obs_index_p(cnt) = g-begg+1 else diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 3e4273758..fd4da863f 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -143,40 +143,40 @@ subroutine define_clm_statevec(mype) ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5 ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell ! Case 2: Assimilation of snow water equivalent same as above - if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then + if(clmupdate_snow==1 .or. clmupdate_snow==2) then clm_varsize = (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA clm_statevecsize = (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority endif ! Case 3: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice - if(clmupdate_snow.eq.3) then + if(clmupdate_snow==3) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif ! Case 4: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq ! and dz - if(clmupdate_snow.eq.4) then + if(clmupdate_snow==4) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif ! Case 5: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice ! and dz - if(clmupdate_snow.eq.5) then + if(clmupdate_snow==5) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif ! Case 6: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq - if(clmupdate_snow.eq.6) then + if(clmupdate_snow==6) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif ! Case 7: Assimilation of snow depth: Snow depth and snow water ! equivalent in the state vector. Update of h2osoi_ice ! Should reproduce case 3 - if(clmupdate_snow.eq.7) then + if(clmupdate_snow==7) then clm_varsize = (clm_endg-clm_begg+1) clm_statevecsize = 2*(clm_endg-clm_begg+1) endif @@ -452,7 +452,7 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col - snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) + snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) #ifdef PDAF_DEBUG @@ -509,7 +509,7 @@ subroutine set_clm_statevec(tstartcycle, mype) ! Case 1: Snow depth ! Case 2: SWE ! Case 3: Snow depth + SWE - if(clmupdate_snow.ne.0) then + if(clmupdate_snow/=0) then cc = 1 do j=clm_begg,clm_endg ! Only get the snow_depth/swe from the first column of each gridcell @@ -517,27 +517,27 @@ subroutine set_clm_statevec(tstartcycle, mype) newgridcell = .true. do jj=clm_begc,clm_endc g = col%gridcell(jj) - if (g .eq. j) then + if (g == j) then if (newgridcell) then newgridcell = .false. - if(clmupdate_snow.eq.1) then + if(clmupdate_snow==1) then clm_statevec(cc+offset) = snow_depth(jj) - else if(clmupdate_snow.eq.2) then + else if(clmupdate_snow==2) then clm_statevec(cc+offset) = h2osno(jj) - else if(clmupdate_snow.eq.3) then + else if(clmupdate_snow==3) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) - else if(clmupdate_snow.eq.4) then + else if(clmupdate_snow==4) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) - else if(clmupdate_snow.eq.5) then + else if(clmupdate_snow==5) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) - else if(clmupdate_snow.eq.6) then + else if(clmupdate_snow==6) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) - else if(clmupdate_snow.eq.7) then + else if(clmupdate_snow==7) then clm_statevec(cc+offset) = snow_depth(jj) clm_statevec(cc+clm_varsize+offset) = h2osno(jj) else @@ -545,7 +545,7 @@ subroutine set_clm_statevec(tstartcycle, mype) end if endif - endif + endif end do cc = cc + 1 end do @@ -673,7 +673,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") #endif snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) - h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) + h2osno => waterstate_inst%h2osno_col ! snow water equivalent (mm) h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col @@ -690,7 +690,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") CLOSE(71) END IF - IF(clmupdate_swc.NE.0 .OR. clmupdate_snow.NE.0) THEN + IF(clmupdate_swc/=0 .OR. clmupdate_snow/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn6, action="write") @@ -732,7 +732,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Case 2: Snow water equivalent ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector ! Write updated snow variable back to CLM and then repartition snow and adjust related variables - if(clmupdate_snow.ne.0) then + if(clmupdate_snow/=0) then do j=clm_begc,clm_endc ! Iterate through the columns @@ -743,57 +743,57 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! CLM5-grid-index cc = (col%gridcell(j) - clm_begg + 1) - if (clmupdate_snow.eq.1) then + if (clmupdate_snow==1) then snow_depth_in(j) = snow_depth(j) - else if (clmupdate_snow.eq.2) then + else if (clmupdate_snow==2) then h2osno_in(j) = h2osno(j) - else if (clmupdate_snow.eq.3) then + else if (clmupdate_snow==3) then h2osno_in(j) = h2osno(j) - else if (clmupdate_snow.eq.4) then + else if (clmupdate_snow==4) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) - else if (clmupdate_snow.eq.5) then + else if (clmupdate_snow==5) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) - else if (clmupdate_snow.eq.6) then + else if (clmupdate_snow==6) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) - else if (clmupdate_snow.eq.7) then + else if (clmupdate_snow==7) then snow_depth_in(j) = snow_depth(j) h2osno_in(j) = h2osno(j) end if - if (clmupdate_snow.eq.1) then + if (clmupdate_snow==1) then snow_depth_out(j) = clm_statevec(cc+offset) - else if (clmupdate_snow.eq.2) then + else if (clmupdate_snow==2) then h2osno_out(j) = clm_statevec(cc+offset) - else if (clmupdate_snow.eq.3) then + else if (clmupdate_snow==3) then h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) - else if (clmupdate_snow.eq.4) then + else if (clmupdate_snow==4) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) - else if (clmupdate_snow.eq.5) then + else if (clmupdate_snow==5) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) - else if (clmupdate_snow.eq.6) then + else if (clmupdate_snow==6) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) - else if (clmupdate_snow.eq.7) then + else if (clmupdate_snow==7) then snow_depth_out(j) = clm_statevec(cc+offset) h2osno_out(j) = clm_statevec(cc+clm_varsize+offset) end if - if(clmupdate_snow.eq.1) then + if(clmupdate_snow==1) then ! Update state variable to CLM ! Needed for Case 1/2 if they use repartioning function - if (snow_depth_out(j).gt.0.0) then + if (snow_depth_out(j)>0.0) then snow_depth(j) = snow_depth_out(j) else ! Catch negative or 0 values from DA print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) end if - else if(clmupdate_snow.eq.2) then - if (h2osno_out(j).gt.0.0) then + else if(clmupdate_snow==2) then + if (h2osno_out(j)>0.0) then h2osno(j) = h2osno_out(j) else ! Catch negative or 0 values from DA @@ -804,32 +804,32 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do ! Repartitioning - if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then + if ( clmupdate_snow_repartitioning==1 .or. clmupdate_snow_repartitioning==2) then - if (clmupdate_snow.eq.1) then - if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:))).gt.0.000001) then + if (clmupdate_snow==1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:)))>0.000001) then call clm_repartition_snow(snow_depth_in(:)) end if end if - if (clmupdate_snow.eq.2) then - if ( SUM(ABS(h2osno_in(:) - h2osno_out(:))).gt.0.000001) then + if (clmupdate_snow==2) then + if ( SUM(ABS(h2osno_in(:) - h2osno_out(:)))>0.000001) then call clm_repartition_snow(h2osno_in(:)) end if end if - if (clmupdate_snow.ge.3) then + if (clmupdate_snow>=3) then error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3" end if end if - if ( clmupdate_snow_repartitioning.eq.3) then + if ( clmupdate_snow_repartitioning==3) then - if (clmupdate_snow.eq.1) then + if (clmupdate_snow==1) then do j=clm_begc,clm_endc - if (snow_depth_out(j).gt.1.0d-6) then - if (snow_depth_in(j).gt.1.0d-6) then + if (snow_depth_out(j)>1.0d-6) then + if (snow_depth_in(j)>1.0d-6) then ! Update h2osoi_ice with increment incr_sno = snow_depth_out(j) / snow_depth_in(j) do i=snlsno(j)+1,0 @@ -840,10 +840,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if - if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + if (clmupdate_snow==2 .or. clmupdate_snow==3) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.1.0d-6) then - if (h2osno_in(j).gt.1.0d-6) then + if (h2osno_out(j)>1.0d-6) then + if (h2osno_in(j)>1.0d-6) then ! Update h2osoi_ice with increment incr_sno = h2osno_out(j) / h2osno_in(j) do i=snlsno(j)+1,0 @@ -854,33 +854,33 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end if - if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + if (clmupdate_snow==4 .or. clmupdate_snow==5 .or. clmupdate_snow==6 .or. clmupdate_snow==7) then do j=clm_begc,clm_endc - if (h2osno_out(j).gt.1.0d-6) then - if (h2osno_in(j).gt.1.0d-6) then - if (snow_depth_out(j).gt.1.0d-6) then - if (snow_depth_in(j).gt.1.0d-6) then + if (h2osno_out(j)>1.0d-6) then + if (h2osno_in(j)>1.0d-6) then + if (snow_depth_out(j)>1.0d-6) then + if (snow_depth_in(j)>1.0d-6) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno_out(j) / h2osno_in(j) ! Update snow_depth with increment incr_sd = snow_depth_out(j) / snow_depth_in(j) do i=snlsno(j)+1,0 - if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then + if (clmupdate_snow==4 .or. clmupdate_snow==5 .or. clmupdate_snow==6 .or. clmupdate_snow==7) then h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe if (isnan(h2osoi_ice(j,i))) then print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i endif end if - if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.6) then + if (clmupdate_snow==4 .or. clmupdate_snow==6) then h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe if (isnan(h2osoi_liq(j,i))) then print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i endif end if - if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5) then + if (clmupdate_snow==4 .or. clmupdate_snow==5) then dz(j,i) = dz(j,i) * incr_sd if (isnan(dz(j,i))) then print *, "WARNING: dz at j,i is nan: ", j, i @@ -1157,22 +1157,25 @@ subroutine clm_repartition_snow(h2osno_in) real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice real(r8) :: gain_dzsno(-nlevsno+1:0) real(r8) :: rho_avg, z_avg - integer :: i,ii,j,jj,g,cc=1,offset=0 + integer :: i,ii,j,jj,g,cc,offset + + cc = 1 + offset = 0 snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) - snowdp => waterstate_inst%snowdp_col ! area-averaged snow height (m) + snowdp => waterstate_inst%snowdp_col ! area-averaged snow height (m) h2osno => waterstate_inst%h2osno_col ! col snow water (mm H2O) - h2oliq => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) + h2oliq => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) h2oice => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd) snlsno => col%snl ! number of snow layers (negative) - frac_sno => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow + frac_sno => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow ! dz for snow layers is defined like in the history output as col%dz for the snow layers dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0) ! Iterate through all columns do jj=clm_begc,clm_endc - if (h2osno(jj).lt.0.0) then ! No snow in column + if (h2osno(jj)<0.0) then ! No snow in column print *, "WARNING: negative snow in col: ", jj, h2osno ! ! Set existing layers to near zero and let CLM do the layer aggregation do i=0,snlsno(jj)+1,-1 @@ -1191,8 +1194,8 @@ subroutine clm_repartition_snow(h2osno_in) endif else ! snow (h2osno) present - if (snlsno(jj).lt.0) then ! snow layers in the column - if (clmupdate_snow .eq. 1) then + if (snlsno(jj)<0) then ! snow layers in the column + if (clmupdate_snow == 1) then ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90 ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth @@ -1200,18 +1203,18 @@ subroutine clm_repartition_snow(h2osno_in) ! v1 init ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD ! v2 SoilTemperatureMod - if (snowdp(jj).gt.0.0_r8) then + if (snowdp(jj)>0.0_r8) then rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj)) else rho_avg = 200.0_r8 end if - if (frac_sno(jj).gt.0.0_r8 .and. snlsno(jj).lt.0.0_r8) then + if (frac_sno(jj)>0.0_r8 .and. snlsno(jj)<0.0_r8) then h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj)) else h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2) then + else if (clmupdate_snow == 2) then ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already h2osno_po(jj) = h2osno(jj) h2osno_pr(jj) = h2osno_in(jj) @@ -1220,16 +1223,16 @@ subroutine clm_repartition_snow(h2osno_in) do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers ! DART VERSION: ii = nlevsoi - i + 1 ! snow density prior for each layer - if (dzsno(jj,ii).gt.0.0_r8) then + if (dzsno(jj,ii)>0.0_r8) then snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii)) else snowden = 0.0_r8 endif ! fraction of SWE in each active layers - if(h2osno_pr(jj).gt.0.0_r8) then + if(h2osno_pr(jj)>0.0_r8) then ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else) - if(clmupdate_snow_repartitioning.eq.1) then - if (ii .eq. 0) then ! DART version indexing check against nlevsno but for us 0 + if(clmupdate_snow_repartitioning==1) then + if (ii == 0) then ! DART version indexing check against nlevsno but for us 0 frac_swe = 1.0_r8 ! JK: Let CLM repartitioning do the job ! afterwards. Provide all the snow in a single layer @@ -1237,7 +1240,7 @@ subroutine clm_repartition_snow(h2osno_in) frac_swe = 0.0_r8 end if ! repartition Option 2: Adjust all active layers - else if (clmupdate_snow_repartitioning.eq.2) then + else if (clmupdate_snow_repartitioning==2) then frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj) end if else @@ -1245,7 +1248,7 @@ subroutine clm_repartition_snow(h2osno_in) end if ! end SWE fraction if ! fraction of liquid and ice - if ((h2oliq(jj,ii) + h2oice(jj,ii)).gt.0.0_r8) then + if ((h2oliq(jj,ii) + h2oice(jj,ii))>0.0_r8) then frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii)) frac_ice = 1.0_r8 - frac_liq else @@ -1253,9 +1256,9 @@ subroutine clm_repartition_snow(h2osno_in) frac_ice = 0.0_r8 end if - ! SWE adjustment per layer + ! SWE adjustment per layer ! assumes identical layer distribution of liq and ice than before DA (frac_*) - if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + if (abs(h2osno_po(jj) - h2osno_pr(jj))>0.0_r8) then gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe gain_h2oliq = gain_h2osno * frac_liq gain_h2oice = gain_h2osno * frac_ice @@ -1265,7 +1268,7 @@ subroutine clm_repartition_snow(h2osno_in) gain_h2oice = 0.0_r8 end if ! layer level adjustments - if (snowden.gt.0.0_r8) then + if (snowden>0.0_r8) then gain_dzsno(ii) = gain_h2osno / snowden else gain_dzsno(ii) = 0.0_r8 @@ -1277,14 +1280,14 @@ subroutine clm_repartition_snow(h2osno_in) ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod ! therefore we adjust dz(:, snow_layer) here - if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then + if (abs(h2osno_po(jj) - h2osno_pr(jj))>0.0_r8) then col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii) ! mid point and interface adjustments ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers)) ! DART version the sum goes from ilevel:nlevsno to fit with our indexing: col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8 ! In DART the check is ilevel == nlevsno but here - if (ii.eq.0) then + if (ii==0) then col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8 else col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8 @@ -1295,13 +1298,13 @@ subroutine clm_repartition_snow(h2osno_in) endif ! End of snow layers present check ! Column level variables - if (clmupdate_snow .eq. 1) then + if (clmupdate_snow == 1) then ! Finally adjust SWE (h2osno) since the prior value is no longer needed ! column level variables so we can adjust it outside the layer loop - h2osno(jj) = h2osno_po(jj) - else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then + h2osno(jj) = h2osno_po(jj) + else if (clmupdate_snow == 2 .or. clmupdate_snow == 3) then ! Update the total snow depth to match udpates to layers for active snow layers - if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then + if (abs(h2osno_po(jj) - h2osno_pr(jj)) > 0.0_r8 .and. snlsno(jj) < 0.0_r8) then snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0)) end if end if From ddbc17da39973fba277f478e866bc596ef084d17 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 31 Oct 2025 13:44:10 +0100 Subject: [PATCH 31/32] fortitude style updates II --- interface/model/eclm/enkf_clm_mod_5.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index fd4da863f..b12733cae 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -790,7 +790,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") snow_depth(j) = snow_depth_out(j) else ! Catch negative or 0 values from DA - print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j) + print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", & + cc, j, offset, snow_depth_out(j) end if else if(clmupdate_snow==2) then if (h2osno_out(j)>0.0) then From 25c1b588a2c68475a0fd22194fb125fc029d73ac Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 31 Oct 2025 14:43:50 +0100 Subject: [PATCH 32/32] compilation errors: re-introduce necessary declarations --- interface/model/eclm/enkf_clm_mod_5.F90 | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index b12733cae..bc098c11d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -635,6 +635,7 @@ end subroutine set_clm_statevec_swc subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use ColumnType , only : col use clm_instMod, only : waterstate_inst implicit none @@ -650,14 +651,28 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) + real(r8) :: incr_sno + real(r8) :: incr_sd + real(r8) :: incr_swe + real(r8) :: h2osno_in(clm_begc:clm_endc) + real(r8) :: snow_depth_in(clm_begc:clm_endc) + real(r8) :: h2osno_out(clm_begc:clm_endc) + real(r8) :: snow_depth_out(clm_begc:clm_endc) + integer :: i + integer :: j + integer :: cc + integer :: offset character (len = 31) :: fn !TSMP-PDAF: function name for state vector output + character (len = 32) :: fn4 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu integer, pointer :: snlsno(:) logical :: swc_zero_before_update + cc = 0 + offset = 0 swc_zero_before_update = .false. #ifdef PDAF_DEBUG