From c49edbe9f95e458fde5220e7e664a72cbc2ef459 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Dec 2024 09:37:21 +0100 Subject: [PATCH 01/36] first implementation of LST-DA from TSMP repo, branch `dev-lstda`, commit 143ea800 --- interface/framework/init_dim_obs_f_pdaf.F90 | 159 +++++++++++++++----- interface/framework/init_dim_obs_pdaf.F90 | 159 +++++++++++++++----- interface/framework/obs_op_pdaf.F90 | 14 +- interface/model/clm3_5/enkf_clm_mod.F90 | 2 +- interface/model/clm5_0/enkf_clm_mod_5.F90 | 86 ++++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 86 ++++++++++- 6 files changed, 405 insertions(+), 101 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index d3f1484d3..60c2e75c3 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -119,6 +119,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) #ifdef CLMFIVE use GridcellType, only: grc use ColumnType, only : col + use PatchType, only : patch ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg use enkf_clm_mod, only: state_clm2pdaf_p @@ -135,6 +136,8 @@ 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_swc + use enkf_clm_mod, only: clmupdate_T !hcp end #endif #endif @@ -156,9 +159,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) integer :: ierror INTEGER :: max_var_id ! Multi-scale DA INTEGER :: sum_dim_obs_p + INTEGER :: p ! CLM Patch index INTEGER :: c ! CLM Column index INTEGER :: g ! CLM Gridcell index INTEGER :: cg + INTEGER :: pg + INTEGER :: pc INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -921,68 +927,141 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) cnt = 1 do i = 1, dim_obs - obs(i) = clm_obs(i) - do g = begg,endg - newgridcell = .true. + obs(i) = clm_obs(i) - do c = begc,endc + if(clmupdate_swc.eq.1) then - cg = mycgridcell(c) + do g = begg,endg + newgridcell = .true. - if(cg .eq. g) then + do c = begc,endc - if(newgridcell) then + cg = mycgridcell(c) - if(is_use_dr) then - deltax = abs(lon(g)-clmobs_lon(i)) - deltay = abs(lat(g)-clmobs_lat(i)) - end if + if(cg .eq. g) then + + if(newgridcell) then + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then - ! Different settings of observation-location-index in - ! state vector depending on the method of state - ! vector assembling. - if(clmstatevec_allcol.eq.1) then + ! Different settings of observation-location-index in + ! state vector depending on the method of state + ! vector assembling. + if(clmstatevec_allcol.eq.1) then #ifdef CLMFIVE - if(clmstatevec_only_active.eq.1) then - - ! Error if observation deeper than clmstatevec_max_layer - if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." - print *, "i=", i - print *, "c=", c - print *, "clmobs_layer(i)=", clmobs_layer(i) - print *, "col%nbedrock(c)=", col%nbedrock(c) - print *, "clmstatevec_max_layer=", clmstatevec_max_layer - call abort_parallel() - end if - obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) - else + if(clmstatevec_only_active.eq.1) then + + ! Error if observation deeper than clmstatevec_max_layer + if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "i=", i + print *, "c=", c + print *, "clmobs_layer(i)=", clmobs_layer(i) + print *, "col%nbedrock(c)=", col%nbedrock(c) + print *, "clmstatevec_max_layer=", clmstatevec_max_layer + call abort_parallel() + end if + obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) + else #endif - obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) #ifdef CLMFIVE - end if + end if #endif - else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 end if - !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) - obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) - cnt = cnt + 1 + newgridcell = .false. + end if - newgridcell = .false. + end if + + end do + end do + + else if(clmupdate_T.eq.1) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + pg = patch%gridcell(p) + pc = patch%column(p) + + if(pg .eq. g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! Set index in state vector, LST will be computed + ! for first patch appearing here + obs_index_p(cnt) = state_clm2pdaf_p(p,1) + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end if + + end if end if + newgridcell = .false. + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 end if + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + end do - end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR unsupported update in setting obs_index_p." + call abort_parallel() + + end if + end do if(obs_interp_switch.eq.1) then diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index e1a6ee4e0..4a5635a23 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -115,6 +115,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) #ifdef CLMFIVE use GridcellType, only: grc use ColumnType, only : col + use PatchType, only : patch ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg use enkf_clm_mod, only: state_clm2pdaf_p @@ -131,6 +132,8 @@ 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_swc + use enkf_clm_mod, only: clmupdate_T !hcp end #endif #endif @@ -152,9 +155,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) integer :: ierror INTEGER :: max_var_id ! Multi-scale DA INTEGER :: sum_dim_obs_p + INTEGER :: p ! CLM Patch index INTEGER :: c ! CLM Column index INTEGER :: g ! CLM Gridcell index INTEGER :: cg + INTEGER :: pg + INTEGER :: pc INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -914,68 +920,141 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) cnt = 1 do i = 1, dim_obs - obs(i) = clm_obs(i) - do g = begg,endg - newgridcell = .true. + obs(i) = clm_obs(i) - do c = begc,endc + if(clmupdate_swc.eq.1) then - cg = mycgridcell(c) + do g = begg,endg + newgridcell = .true. - if(cg .eq. g) then + do c = begc,endc - if(newgridcell) then + cg = mycgridcell(c) - if(is_use_dr) then - deltax = abs(lon(g)-clmobs_lon(i)) - deltay = abs(lat(g)-clmobs_lat(i)) - end if + if(cg .eq. g) then + + if(newgridcell) then + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then - ! Different settings of observation-location-index in - ! state vector depending on the method of state - ! vector assembling. - if(clmstatevec_allcol.eq.1) then + ! Different settings of observation-location-index in + ! state vector depending on the method of state + ! vector assembling. + if(clmstatevec_allcol.eq.1) then #ifdef CLMFIVE - if(clmstatevec_only_active.eq.1) then - - ! Error if observation deeper than clmstatevec_max_layer - if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." - print *, "i=", i - print *, "c=", c - print *, "clmobs_layer(i)=", clmobs_layer(i) - print *, "col%nbedrock(c)=", col%nbedrock(c) - print *, "clmstatevec_max_layer=", clmstatevec_max_layer - call abort_parallel() - end if - obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) - else + if(clmstatevec_only_active.eq.1) then + + ! Error if observation deeper than clmstatevec_max_layer + if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "i=", i + print *, "c=", c + print *, "clmobs_layer(i)=", clmobs_layer(i) + print *, "col%nbedrock(c)=", col%nbedrock(c) + print *, "clmstatevec_max_layer=", clmstatevec_max_layer + call abort_parallel() + end if + obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) + else #endif - obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) #ifdef CLMFIVE - end if + end if #endif - else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + else + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 end if - !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) - obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) - cnt = cnt + 1 + newgridcell = .false. + end if - newgridcell = .false. + end if + + end do + end do + + else if(clmupdate_T.eq.1) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + pg = patch%gridcell(p) + pc = patch%column(p) + + if(pg .eq. g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! Set index in state vector, LST will be computed + ! for first patch appearing here + obs_index_p(cnt) = state_clm2pdaf_p(p,1) + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end if + + end if end if + newgridcell = .false. + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 end if + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + end do - end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR unsupported update in setting obs_index_p." + call abort_parallel() + + end if + end do if(obs_interp_switch.eq.1) then diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 23a755d4f..265a1e472 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -114,13 +114,14 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) lpointobs = .false. + ! Equation for LST computation from Kustas2009, Eq(7) + ! http://dx.doi.org/10.1016/j.agrformet.2009.05.016 + ! + ! Comment: Fractional vegetation cover (Eq(8) from Kustas2009) + ! currently implemented with simplified settings: Vegetation + ! clumping parameter `Omega=1`; radiometer view angle `phi=0` + DO i = 1, dim_obs_p - ! Equation for LST computation from Kustas2009, Eq(7) - ! http://dx.doi.org/10.1016/j.agrformet.2009.05.016 - ! - ! Comment: Fractional vegetation cover (Eq(8) from Kustas2009) - ! currently implemented with simplified settings: Vegetation - ! clumping parameter `Omega=1`; radiometer view angle `phi=0` m_state_p(i) & = (exp(-0.5*clm_paramarr(obs_index_p(i))) & *state_p(obs_index_p(i))**4 & @@ -132,6 +133,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! write(*,*) 'model LST', m_state_p(:) ! write(*,*) 'TG', state_p(obs_index_p(:)) ! write(*,*) 'TV', state_p(clm_varsize+obs_index_p(:)) + endif #endif diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 0e9d832cc..db5719b47 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -70,7 +70,7 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol integer(c_int),bind(C,name="clmt_printensemble") :: clmt_printensemble - integer(c_int),bind(C,name="clmwatmin_switch") :: clmwatmin_switch + integer(c_int),bind(C,name="clmwatmin_switch") :: clmwatmin_switch integer :: nstep ! time step index real(r8) :: dtime ! time step increment (sec) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index c0b89c299..3e5d3b2ec 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -42,6 +42,7 @@ module enkf_clm_mod integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) integer,allocatable :: state_pdaf2clm_c_p(:) + integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) @@ -87,6 +88,7 @@ subroutine define_clm_statevec(mype) use clm_varpar , only : nlevsoi use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -95,6 +97,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: p integer :: c integer :: g integer :: cc @@ -272,7 +275,37 @@ subroutine define_clm_statevec(mype) !hcp LST DA if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = (endp-begp+1)*2 !TG, then TV + + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + endif !end hcp @@ -291,7 +324,7 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp if ((clmupdate_T.ne.0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec @@ -309,11 +342,15 @@ subroutine cleanup_clm_statevec() end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer,intent(in) :: tstartcycle @@ -322,6 +359,9 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: tlai(:) 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 +371,12 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + ! LST variables (t_veg is patch variable...) + t_grnd => temperature_inst%t_grnd_col + t_veg => temperature_inst%t_veg_patch + tlai => canopystate_inst%tlai_patch + + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -359,7 +405,18 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + do cc = 1, clm_varsize + ! t_grnd iterated over cols + ! t_veg iterated over patches + clm_statevec(cc) = t_grnd(state_pdaf2clm_c_p(cc)) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) + end do + + do cc = 1, clm_paramsize + ! Works only if clm_paramsize corresponds to clm_varsize (also + ! the order) + clm_paramarr(cc) = tlai(state_pdaf2clm_c_p(cc)) + end do endif !end hcp LAI @@ -402,8 +459,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use PatchType , only : patch use ColumnType , only : col - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst use clm_varcon , only : denh2o, denice, watmin use clm_varcon , only : ispval use clm_varcon , only : spval @@ -419,6 +479,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) @@ -427,7 +490,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,p,cc=0,offset=0 character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -459,6 +522,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + ! LST + t_grnd => temperature_inst%t_grnd_col + t_veg => temperature_inst%t_veg_patch + ! tlai => canopystate_inst%tlai_patch + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -605,7 +673,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !hcp: TG, TV if(clmupdate_T.EQ.1) then - error stop "Not implemented: clmupdate_T.eq.1" + do p = clm_begp, clm_endp + c = patch%column(p) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + end do endif ! end hcp TG, TV diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index c0b89c299..3e5d3b2ec 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -42,6 +42,7 @@ module enkf_clm_mod integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) integer,allocatable :: state_pdaf2clm_c_p(:) + integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) @@ -87,6 +88,7 @@ subroutine define_clm_statevec(mype) use clm_varpar , only : nlevsoi use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -95,6 +97,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: p integer :: c integer :: g integer :: cc @@ -272,7 +275,37 @@ subroutine define_clm_statevec(mype) !hcp LST DA if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = (endp-begp+1)*2 !TG, then TV + + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + endif !end hcp @@ -291,7 +324,7 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp if ((clmupdate_T.ne.0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec @@ -309,11 +342,15 @@ subroutine cleanup_clm_statevec() end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer,intent(in) :: tstartcycle @@ -322,6 +359,9 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: tlai(:) 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 +371,12 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + ! LST variables (t_veg is patch variable...) + t_grnd => temperature_inst%t_grnd_col + t_veg => temperature_inst%t_veg_patch + tlai => canopystate_inst%tlai_patch + + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -359,7 +405,18 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + do cc = 1, clm_varsize + ! t_grnd iterated over cols + ! t_veg iterated over patches + clm_statevec(cc) = t_grnd(state_pdaf2clm_c_p(cc)) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) + end do + + do cc = 1, clm_paramsize + ! Works only if clm_paramsize corresponds to clm_varsize (also + ! the order) + clm_paramarr(cc) = tlai(state_pdaf2clm_c_p(cc)) + end do endif !end hcp LAI @@ -402,8 +459,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use PatchType , only : patch use ColumnType , only : col - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst use clm_varcon , only : denh2o, denice, watmin use clm_varcon , only : ispval use clm_varcon , only : spval @@ -419,6 +479,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) @@ -427,7 +490,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,p,cc=0,offset=0 character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -459,6 +522,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + ! LST + t_grnd => temperature_inst%t_grnd_col + t_veg => temperature_inst%t_veg_patch + ! tlai => canopystate_inst%tlai_patch + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -605,7 +673,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !hcp: TG, TV if(clmupdate_T.EQ.1) then - error stop "Not implemented: clmupdate_T.eq.1" + do p = clm_begp, clm_endp + c = patch%column(p) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + end do endif ! end hcp TG, TV From 0fb734a4d7c9eb8823ba5b13366d580e5fd58d8b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 13:51:49 +0100 Subject: [PATCH 02/36] bugfix: re-introduce default for general update --- interface/framework/init_dim_obs_f_pdaf.F90 | 24 +++++++++++++++++++-- interface/framework/init_dim_obs_pdaf.F90 | 24 +++++++++++++++++++-- 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 60c2e75c3..589aa4085 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1057,8 +1057,28 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) #endif else - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR unsupported update in setting obs_index_p." - call abort_parallel() + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do end if diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 4a5635a23..90af36521 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1050,8 +1050,28 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) #endif else - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR unsupported update in setting obs_index_p." - call abort_parallel() + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do end if From b186564e3ab4d6ff5c5a92399bc55c7266f5b87f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 13:52:06 +0100 Subject: [PATCH 03/36] bugfix: lai is a patch-array --- 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 3e5d3b2ec..bc6502874 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -415,7 +415,7 @@ subroutine set_clm_statevec(tstartcycle, mype) do cc = 1, clm_paramsize ! Works only if clm_paramsize corresponds to clm_varsize (also ! the order) - clm_paramarr(cc) = tlai(state_pdaf2clm_c_p(cc)) + clm_paramarr(cc) = tlai(state_pdaf2clm_p_p(cc)) end do endif !end hcp LAI diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 3e5d3b2ec..bc6502874 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -415,7 +415,7 @@ subroutine set_clm_statevec(tstartcycle, mype) do cc = 1, clm_paramsize ! Works only if clm_paramsize corresponds to clm_varsize (also ! the order) - clm_paramarr(cc) = tlai(state_pdaf2clm_c_p(cc)) + clm_paramarr(cc) = tlai(state_pdaf2clm_p_p(cc)) end do endif !end hcp LAI From 35c9c9341ec6c6cfbb8519065fca76212c8cd446 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 13:58:24 +0100 Subject: [PATCH 04/36] bugfix: typo --- 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 bc6502874..76ee9b480 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -287,7 +287,7 @@ subroutine define_clm_statevec(mype) clm_paramsize = endp-begp+1 !LAI clm_statevecsize = (endp-begp+1)*2 !TG, then TV - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) allocate(state_pdaf2clm_c_p(clm_statevecsize)) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index bc6502874..76ee9b480 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -287,7 +287,7 @@ subroutine define_clm_statevec(mype) clm_paramsize = endp-begp+1 !LAI clm_statevecsize = (endp-begp+1)*2 !TG, then TV - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) allocate(state_pdaf2clm_c_p(clm_statevecsize)) From 929471f1a31d3bf71d1e43b4770fac395f9ab225 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 14:19:30 +0100 Subject: [PATCH 05/36] dev: skin temperature state vector, first implementation --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 56 ++++++++++++++++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 56 ++++++++++++++++++++++- 2 files changed, 110 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 76ee9b480..ed07cb36e 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -309,6 +309,42 @@ subroutine define_clm_statevec(mype) endif !end hcp + ! skin temperature state vector + if(clmupdate_T.eq.2) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = (endp-begp+1) !TSKIN + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + + endif + + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, "define_clm_statevec: clm_statevecsize=", clm_statevecsize @@ -323,7 +359,7 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T.ne.0)) then !hcp + if ((clmupdate_T.eq.1)) then !hcp allocate(clm_paramarr(clm_paramsize)) end if @@ -374,6 +410,7 @@ subroutine set_clm_statevec(tstartcycle, mype) ! LST variables (t_veg is patch variable...) t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch tlai => canopystate_inst%tlai_patch @@ -420,6 +457,14 @@ subroutine set_clm_statevec(tstartcycle, mype) endif !end hcp LAI + ! skin temperature state vector + if(clmupdate_T.eq.2) then + do cc = 1, clm_statevecsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -525,6 +570,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! LST t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch ! tlai => canopystate_inst%tlai_patch #ifdef PDAF_DEBUG @@ -681,6 +727,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif ! end hcp TG, TV + ! skin temperature state vector + if(clmupdate_T.EQ.2) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 76ee9b480..ed07cb36e 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -309,6 +309,42 @@ subroutine define_clm_statevec(mype) endif !end hcp + ! skin temperature state vector + if(clmupdate_T.eq.2) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = (endp-begp+1) !TSKIN + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + + endif + + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, "define_clm_statevec: clm_statevecsize=", clm_statevecsize @@ -323,7 +359,7 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T.ne.0)) then !hcp + if ((clmupdate_T.eq.1)) then !hcp allocate(clm_paramarr(clm_paramsize)) end if @@ -374,6 +410,7 @@ subroutine set_clm_statevec(tstartcycle, mype) ! LST variables (t_veg is patch variable...) t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch tlai => canopystate_inst%tlai_patch @@ -420,6 +457,14 @@ subroutine set_clm_statevec(tstartcycle, mype) endif !end hcp LAI + ! skin temperature state vector + if(clmupdate_T.eq.2) then + do cc = 1, clm_statevecsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -525,6 +570,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! LST t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch ! tlai => canopystate_inst%tlai_patch #ifdef PDAF_DEBUG @@ -681,6 +727,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif ! end hcp TG, TV + ! skin temperature state vector + if(clmupdate_T.EQ.2) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi From 0f3738b7a5196889a978a52e8c3e53ec43e836fa Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 14:23:41 +0100 Subject: [PATCH 06/36] dev: skin temperature state vector, first implementation II --- interface/framework/obs_op_pdaf.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 265a1e472..0137e30c9 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -134,6 +134,17 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! write(*,*) 'TG', state_p(obs_index_p(:)) ! write(*,*) 'TV', state_p(clm_varsize+obs_index_p(:)) +endif + +if (clmupdate_T.EQ.2) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + endif #endif From 53998fb638694ef9848e50f331eaa5617294621f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 13 Dec 2024 09:56:05 +0100 Subject: [PATCH 07/36] bugfix: t_skin declarations ``` eclm/enkf_clm_mod_5.F90(413): error #6404: This name does not have a type, and must have an explicit type. [T_SKIN] ``` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 2 ++ interface/model/eclm/enkf_clm_mod_5.F90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index ed07cb36e..f9a886878 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -397,6 +397,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output @@ -526,6 +527,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: t_grnd(:) real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index ed07cb36e..f9a886878 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -397,6 +397,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) integer :: i,j,jj,g,cc=0,offset=0 character (len = 34) :: fn !TSMP-PDAF: function name for state vector output @@ -526,6 +527,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: t_grnd(:) real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) From 75ac19fb1ba4a5e3d46dac46206bbd9ce5e11c6f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 08:23:29 +0100 Subject: [PATCH 08/36] bugfix: location of `newgridcell = .false` this was wrongly placed outside the condition checking `newgridcell` leading to `obs_index_p` not being set at all, if not at the very first patch. --- interface/framework/init_dim_obs_f_pdaf.F90 | 4 ++-- interface/framework/init_dim_obs_pdaf.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 589aa4085..9581edfd8 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1028,11 +1028,11 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end if + newgridcell = .false. + end if end if - newgridcell = .false. - end do end do #else diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 90af36521..6ab7fc2cc 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1021,11 +1021,11 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if + newgridcell = .false. + end if end if - newgridcell = .false. - end do end do #else From 34a47ca95f42b3457387fb1bf986a3dc1b439442 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 08:49:07 +0100 Subject: [PATCH 09/36] LST-DA with TSKIN: use same obs_index_p setting as for TG/TV --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/framework/init_dim_obs_pdaf.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 9581edfd8..98d5abb3b 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -993,7 +993,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end do end do - else if(clmupdate_T.eq.1) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 6ab7fc2cc..9220d106c 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -986,7 +986,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T.eq.1) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2) then #ifdef CLMFIVE ! patch loop do g = begg,endg From 43d1ba30eb47e754bed0345f6cf083af3cfa3096 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 09:14:43 +0100 Subject: [PATCH 10/36] LST-DA: state-vector with TSKIN, plus TG and TV --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 21 ++++++++++++++------- interface/model/eclm/enkf_clm_mod_5.F90 | 21 ++++++++++++++------- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index f9a886878..8df584cfa 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -321,7 +321,7 @@ subroutine define_clm_statevec(mype) clm_varsize = endp-begp+1 ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = (endp-begp+1) !TSKIN + clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TG and TV IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) @@ -334,12 +334,15 @@ subroutine define_clm_statevec(mype) do p=clm_begp,clm_endp cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TG - state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_p_p(cc+clm_varsize) = p !TG + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TG state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 end do endif @@ -460,9 +463,11 @@ subroutine set_clm_statevec(tstartcycle, mype) ! skin temperature state vector if(clmupdate_T.eq.2) then - do cc = 1, clm_statevecsize + do cc = 1, clm_varsize ! t_skin iterated over patches - clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_grnd(state_pdaf2clm_c_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) end do endif @@ -734,6 +739,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) end do endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index f9a886878..8df584cfa 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -321,7 +321,7 @@ subroutine define_clm_statevec(mype) clm_varsize = endp-begp+1 ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = (endp-begp+1) !TSKIN + clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TG and TV IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) @@ -334,12 +334,15 @@ subroutine define_clm_statevec(mype) do p=clm_begp,clm_endp cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TG - state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_p_p(cc+clm_varsize) = p !TG + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TG state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 end do endif @@ -460,9 +463,11 @@ subroutine set_clm_statevec(tstartcycle, mype) ! skin temperature state vector if(clmupdate_T.eq.2) then - do cc = 1, clm_statevecsize + do cc = 1, clm_varsize ! t_skin iterated over patches - clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_grnd(state_pdaf2clm_c_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) end do endif @@ -734,6 +739,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) end do endif From 3df6fd0530d7edda2e9dfe1c52e574212db73487 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 09:39:14 +0100 Subject: [PATCH 11/36] LST-DA: state-vector with TSKIN, plus TSOIL and TV `clmupdate_T.eq.3` --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/framework/init_dim_obs_pdaf.F90 | 2 +- interface/framework/obs_op_pdaf.F90 | 2 +- interface/model/clm5_0/enkf_clm_mod_5.F90 | 63 ++++++++++++++++++++- interface/model/eclm/enkf_clm_mod_5.F90 | 63 ++++++++++++++++++++- 5 files changed, 127 insertions(+), 5 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 98d5abb3b..070482313 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -993,7 +993,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end do end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 9220d106c..0fe391a9e 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -986,7 +986,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 0137e30c9..cb0b2ec43 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -136,7 +136,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) endif -if (clmupdate_T.EQ.2) then +if (clmupdate_T.EQ.2 .OR. clmupdate_T.EQ.3) then lpointobs = .false. diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 8df584cfa..e5bc3ed68 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -347,6 +347,43 @@ subroutine define_clm_statevec(mype) endif + if(clmupdate_T.eq.3) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TSOIL and TV + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 + end do + + endif + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize @@ -399,6 +436,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) @@ -411,10 +449,11 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col - ! LST variables (t_veg is patch variable...) + ! LST variables t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch t_skin => temperature_inst%t_skin_patch + t_soisno => temperature_inst%t_soisno_col tlai => canopystate_inst%tlai_patch @@ -471,6 +510,16 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! skin temperature state vector updating soil temperature + if(clmupdate_T.eq.3) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+clm_varsize), state_pdaf2clm_j_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -531,6 +580,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) @@ -576,6 +626,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! LST t_grnd => temperature_inst%t_grnd_col + t_soisno => temperature_inst%t_soisno_col t_veg => temperature_inst%t_veg_patch t_skin => temperature_inst%t_skin_patch ! tlai => canopystate_inst%tlai_patch @@ -744,6 +795,16 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do endif + ! skin temperature state vector updating soil temperature + if(clmupdate_T.EQ.3) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 8df584cfa..e5bc3ed68 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -347,6 +347,43 @@ subroutine define_clm_statevec(mype) endif + if(clmupdate_T.eq.3) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TSOIL and TV + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 + end do + + endif + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize @@ -399,6 +436,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) @@ -411,10 +449,11 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col - ! LST variables (t_veg is patch variable...) + ! LST variables t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch t_skin => temperature_inst%t_skin_patch + t_soisno => temperature_inst%t_soisno_col tlai => canopystate_inst%tlai_patch @@ -471,6 +510,16 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! skin temperature state vector updating soil temperature + if(clmupdate_T.eq.3) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+clm_varsize), state_pdaf2clm_j_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -531,6 +580,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: porgm(:,:) real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) @@ -576,6 +626,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! LST t_grnd => temperature_inst%t_grnd_col + t_soisno => temperature_inst%t_soisno_col t_veg => temperature_inst%t_veg_patch t_skin => temperature_inst%t_skin_patch ! tlai => canopystate_inst%tlai_patch @@ -744,6 +795,16 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do endif + ! skin temperature state vector updating soil temperature + if(clmupdate_T.EQ.3) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi From 9e88c2a31174dd24cc6bc053aa9b4094f457dfe1 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 24 Jan 2025 14:51:28 +0100 Subject: [PATCH 12/36] LST-DA: add debug output for first layer `t_soisno` --- interface/model/clm3_5/enkf_clm_mod.F90 | 67 +++++++++++++--------- interface/model/clm5_0/enkf_clm_mod_5.F90 | 68 ++++++++++++++--------- interface/model/eclm/enkf_clm_mod_5.F90 | 68 ++++++++++++++--------- 3 files changed, 125 insertions(+), 78 deletions(-) diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index db5719b47..b4bc2ba8a 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -213,6 +213,14 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + END IF #endif @@ -472,32 +480,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") 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 - - END IF -#endif - endif !hcp: TG, TV @@ -533,6 +515,39 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif +#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(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index e5bc3ed68..2f22db725 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -468,6 +468,14 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + END IF #endif @@ -747,32 +755,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") 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 - - END IF -#endif - endif !hcp: TG, TV @@ -830,6 +812,40 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif +#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(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index e5bc3ed68..2f22db725 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -468,6 +468,14 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + END IF #endif @@ -747,32 +755,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") 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 - - END IF -#endif - endif !hcp: TG, TV @@ -830,6 +812,40 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif +#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(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() From ff0373a4cdc6180a8b40d3b6138eaf5446e637eb Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 28 Jan 2025 20:47:32 +0100 Subject: [PATCH 13/36] LST-DA: all temperature layers updated --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 36 +++++++++++++++-------- interface/model/eclm/enkf_clm_mod_5.F90 | 36 +++++++++++++++-------- 2 files changed, 48 insertions(+), 24 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 2f22db725..40712623a 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -86,6 +86,7 @@ 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 : nlevgrnd use clm_varcon , only : ispval use ColumnType , only : col use PatchType , only : patch @@ -97,6 +98,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev integer :: p integer :: c integer :: g @@ -358,7 +360,7 @@ subroutine define_clm_statevec(mype) clm_varsize = endp-begp+1 ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TSOIL and TV + clm_statevecsize = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and TV IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) @@ -374,13 +376,15 @@ subroutine define_clm_statevec(mype) state_pdaf2clm_p_p(cc) = p !TSKIN state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TSOIL - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TSOIL - state_pdaf2clm_j_p(cc+clm_varsize) = 1 - state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV - state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 - end do + do lev=1,nlevgrnd + ! ivar = 2-26: TSOIL + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) + state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + end do + state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 endif @@ -423,6 +427,7 @@ subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : temperature_inst use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col @@ -441,6 +446,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) integer :: i,j,jj,g,cc=0,offset=0 + integer :: lev character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -523,8 +529,10 @@ subroutine set_clm_statevec(tstartcycle, mype) do cc = 1, clm_varsize ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) - clm_statevec(cc+clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+clm_varsize), state_pdaf2clm_j_p(cc+clm_varsize)) - clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + do lev=1,nlevgrnd + clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), state_pdaf2clm_j_p(cc+lev*clm_varsize)) + end do + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) end do endif @@ -565,6 +573,7 @@ end subroutine set_clm_statevec subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 use PatchType , only : patch @@ -601,6 +610,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: swc_update ! updated SWC in loop integer :: i,j,jj,g,c,p,cc=0,offset=0 + integer :: lev character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -782,8 +792,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) - t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + do lev=1,nlevgrnd + t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1) + (lev)*clm_varsize) + end do + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + (1+nlevgrnd)*clm_varsize) end do endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 2f22db725..40712623a 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -86,6 +86,7 @@ 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 : nlevgrnd use clm_varcon , only : ispval use ColumnType , only : col use PatchType , only : patch @@ -97,6 +98,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev integer :: p integer :: c integer :: g @@ -358,7 +360,7 @@ subroutine define_clm_statevec(mype) clm_varsize = endp-begp+1 ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TSOIL and TV + clm_statevecsize = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and TV IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) @@ -374,13 +376,15 @@ subroutine define_clm_statevec(mype) state_pdaf2clm_p_p(cc) = p !TSKIN state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TSOIL - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TSOIL - state_pdaf2clm_j_p(cc+clm_varsize) = 1 - state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV - state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 - end do + do lev=1,nlevgrnd + ! ivar = 2-26: TSOIL + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) + state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + end do + state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 endif @@ -423,6 +427,7 @@ subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : temperature_inst use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col @@ -441,6 +446,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) integer :: i,j,jj,g,cc=0,offset=0 + integer :: lev character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -523,8 +529,10 @@ subroutine set_clm_statevec(tstartcycle, mype) do cc = 1, clm_varsize ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) - clm_statevec(cc+clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+clm_varsize), state_pdaf2clm_j_p(cc+clm_varsize)) - clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + do lev=1,nlevgrnd + clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), state_pdaf2clm_j_p(cc+lev*clm_varsize)) + end do + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) end do endif @@ -565,6 +573,7 @@ end subroutine set_clm_statevec subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 use PatchType , only : patch @@ -601,6 +610,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: swc_update ! updated SWC in loop integer :: i,j,jj,g,c,p,cc=0,offset=0 + integer :: lev character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -782,8 +792,10 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) - t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + do lev=1,nlevgrnd + t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1) + (lev)*clm_varsize) + end do + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + (1+nlevgrnd)*clm_varsize) end do endif From 4bd276d60ed5b172e7f2d775d690369b977b30ed Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 29 Jan 2025 09:25:27 +0100 Subject: [PATCH 14/36] bugfix: LST-DA missing `end do` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 1 + interface/model/eclm/enkf_clm_mod_5.F90 | 1 + 2 files changed, 2 insertions(+) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index 40712623a..3febe498e 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -385,6 +385,7 @@ subroutine define_clm_statevec(mype) state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 + end do endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 40712623a..3febe498e 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -385,6 +385,7 @@ subroutine define_clm_statevec(mype) state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 + end do endif From 0778801bd13b879580cdab74da81f519a532f80c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 25 Apr 2025 09:55:42 +0200 Subject: [PATCH 15/36] introduce `clmupdate_T.eq.4` --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/framework/init_dim_obs_pdaf.F90 | 2 +- interface/framework/obs_op_pdaf.F90 | 2 +- interface/model/clm5_0/enkf_clm_mod_5.F90 | 46 +++++++++++++++++++++ interface/model/eclm/enkf_clm_mod_5.F90 | 46 +++++++++++++++++++++ 5 files changed, 95 insertions(+), 3 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 070482313..98b44597c 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -993,7 +993,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end do end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.4) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 0fe391a9e..c17dc28f8 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -986,7 +986,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3) then + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.4) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index abc9f0d57..62ef12e5a 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -155,7 +155,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) endif -if (clmupdate_T.EQ.2 .OR. clmupdate_T.EQ.3) then +if (clmupdate_T.EQ.2 .OR. clmupdate_T.EQ.3 .OR. clmupdate_T.EQ.4) then lpointobs = .false. diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index a58ecce74..bdb7df336 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -389,6 +389,37 @@ subroutine define_clm_statevec(mype) endif + if(clmupdate_T.eq.4) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = 1* (endp-begp+1) !TSOIL + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSOIL + state_pdaf2clm_c_p(cc) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc) = 1 + end do + + endif + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize @@ -537,6 +568,13 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! soil temperature state vector updating soil temperature + if(clmupdate_T.eq.4) then + do cc = 1, clm_varsize + clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -800,6 +838,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do endif + ! soil temperature state vector updating soil temperature + if(clmupdate_T.EQ.4) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index a58ecce74..bdb7df336 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -389,6 +389,37 @@ subroutine define_clm_statevec(mype) endif + if(clmupdate_T.eq.4) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = 1* (endp-begp+1) !TSOIL + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSOIL + state_pdaf2clm_c_p(cc) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc) = 1 + end do + + endif + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize @@ -537,6 +568,13 @@ subroutine set_clm_statevec(tstartcycle, mype) end do endif + ! soil temperature state vector updating soil temperature + if(clmupdate_T.eq.4) then + do cc = 1, clm_varsize + clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -800,6 +838,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do endif + ! soil temperature state vector updating soil temperature + if(clmupdate_T.EQ.4) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi From fb67c10cfd63480c276dff51a30fb1abf060c1ba Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 9 Oct 2025 22:53:43 +0200 Subject: [PATCH 16/36] fortitude fixes --- interface/framework/init_dim_obs_pdaf.F90 | 20 +++++++------- interface/framework/obs_op_pdaf.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 32 +++++++++++------------ 3 files changed, 27 insertions(+), 27 deletions(-) diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index a87f91820..026dbb775 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -948,7 +948,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) obs(i) = clm_obs(i) - if(clmupdate_swc.eq.1) then + if(clmupdate_swc==1) then do g = begg,endg newgridcell = .true. @@ -957,7 +957,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) cg = mycgridcell(c) - if(cg .eq. g) then + if(cg == g) then if(newgridcell) then @@ -1010,7 +1010,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.4) then + else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then #ifdef CLMFIVE ! patch loop do g = begg,endg @@ -1021,7 +1021,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) pg = patch%gridcell(p) pc = patch%column(p) - if(pg .eq. g) then + if(pg == g) then if(newgridcell) then ! Sets first patch/column in a gridcell. TODO: Make ! patch / column information part of the observation @@ -1032,7 +1032,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then ! Set index in state vector, LST will be computed ! for first patch appearing here @@ -1040,7 +1040,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end if @@ -1061,13 +1061,13 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end do @@ -1086,13 +1086,13 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end do diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index d76ce581e..9d74ff996 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -155,7 +155,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) endif -if (clmupdate_T.EQ.2 .OR. clmupdate_T.EQ.3 .OR. clmupdate_T.EQ.4) then +if (clmupdate_T==2 .OR. clmupdate_T==3 .OR. clmupdate_T==4) then lpointobs = .false. diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 961b84267..5a8e5cb35 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -303,7 +303,7 @@ subroutine define_clm_statevec(mype) endif !hcp LST DA - if(clmupdate_T.eq.1) then + if(clmupdate_T==1) then IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) allocate(state_clm2pdaf_p(begp:endp,1)) @@ -339,7 +339,7 @@ subroutine define_clm_statevec(mype) !end hcp ! skin temperature state vector - if(clmupdate_T.eq.2) then + if(clmupdate_T==2) then IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) allocate(state_clm2pdaf_p(begp:endp,1)) @@ -376,7 +376,7 @@ subroutine define_clm_statevec(mype) endif - if(clmupdate_T.eq.3) then + if(clmupdate_T==3) then IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) allocate(state_clm2pdaf_p(begp:endp,1)) @@ -416,7 +416,7 @@ subroutine define_clm_statevec(mype) endif - if(clmupdate_T.eq.4) then + if(clmupdate_T==4) then IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) allocate(state_clm2pdaf_p(begp:endp,1)) @@ -470,7 +470,7 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T.eq.1)) then !hcp + if ((clmupdate_T==1)) then !hcp allocate(clm_paramarr(clm_paramsize)) end if @@ -545,7 +545,7 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF - IF(clmupdate_T.NE.0) THEN + IF(clmupdate_T/=0) THEN ! TSMP-PDAF: Debug output of CLM t_soisno, first layer WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" OPEN(unit=71, file=fn2, action="write") @@ -611,7 +611,7 @@ subroutine set_clm_statevec(tstartcycle, mype) endif !hcp LAI - if(clmupdate_T.eq.1) then + if(clmupdate_T==1) then do cc = 1, clm_varsize ! t_grnd iterated over cols ! t_veg iterated over patches @@ -628,7 +628,7 @@ subroutine set_clm_statevec(tstartcycle, mype) !end hcp LAI ! skin temperature state vector - if(clmupdate_T.eq.2) then + if(clmupdate_T==2) then do cc = 1, clm_varsize ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) @@ -638,7 +638,7 @@ subroutine set_clm_statevec(tstartcycle, mype) endif ! skin temperature state vector updating soil temperature - if(clmupdate_T.eq.3) then + if(clmupdate_T==3) then do cc = 1, clm_varsize ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) @@ -650,7 +650,7 @@ subroutine set_clm_statevec(tstartcycle, mype) endif ! soil temperature state vector updating soil temperature - if(clmupdate_T.eq.4) then + if(clmupdate_T==4) then do cc = 1, clm_varsize clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) end do @@ -944,7 +944,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif !hcp: TG, TV - if(clmupdate_T.EQ.1) then + if(clmupdate_T==1) then do p = clm_begp, clm_endp c = patch%column(p) t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1)) @@ -954,7 +954,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! end hcp TG, TV ! skin temperature state vector - if(clmupdate_T.EQ.2) then + if(clmupdate_T==2) then do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) @@ -964,7 +964,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif ! skin temperature state vector updating soil temperature - if(clmupdate_T.EQ.3) then + if(clmupdate_T==3) then do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) @@ -976,7 +976,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif ! soil temperature state vector updating soil temperature - if(clmupdate_T.EQ.4) then + if(clmupdate_T==4) then do p = clm_begp, clm_endp c = patch%column(p) t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) @@ -1011,7 +1011,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN - IF(clmupdate_swc.NE.0) THEN + IF(clmupdate_swc/=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") @@ -1031,7 +1031,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") CLOSE(71) END IF - IF(clmupdate_T.NE.0) THEN + IF(clmupdate_T/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" OPEN(unit=71, file=fn2, action="write") From 75f2b8d4b26da8528318b970a5371bd1bd63e93a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 9 Oct 2025 23:06:35 +0200 Subject: [PATCH 17/36] keep indentation of SM-case --- interface/framework/init_dim_obs_f_pdaf.F90 | 36 +++++++++++---------- interface/framework/init_dim_obs_pdaf.F90 | 28 ++++++++++------ 2 files changed, 38 insertions(+), 26 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index fbb500562..d9222e96b 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -952,17 +952,20 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) cnt = 1 do i = 1, dim_obs + obs(i) = clm_obs(i) - obs(i) = clm_obs(i) + if(clmupdate_swc==1) then - if(clmupdate_swc.eq.1) then + do g = begg,endg + newgridcell = .true. - do g = begg,endg - newgridcell = .true. + do c = begc,endc - do c = begc,endc + cg = mycgridcell(c) - cg = mycgridcell(c) + if(cg == g) then + + if(newgridcell) then if(is_use_dr) then deltax = abs(lon(g)-clmobs_lon(i)) @@ -1016,16 +1019,15 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end if - newgridcell = .false. - end if + newgridcell = .false. end if - end do end do + end do - else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.4) then + else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then #ifdef CLMFIVE ! patch loop do g = begg,endg @@ -1036,7 +1038,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) pg = patch%gridcell(p) pc = patch%column(p) - if(pg .eq. g) then + if(pg == g) then if(newgridcell) then ! Sets first patch/column in a gridcell. TODO: Make ! patch / column information part of the observation @@ -1047,7 +1049,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then ! Set index in state vector, LST will be computed ! for first patch appearing here @@ -1055,7 +1057,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end if @@ -1076,13 +1078,13 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end do @@ -1101,13 +1103,13 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) end if !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) cnt = cnt + 1 end do diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 026dbb775..8f3f80f97 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -945,21 +945,31 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) cnt = 1 do i = 1, dim_obs - - obs(i) = clm_obs(i) + obs(i) = clm_obs(i) if(clmupdate_swc==1) then - do g = begg,endg - newgridcell = .true. + do g = begg,endg + newgridcell = .true. - do c = begc,endc + do c = begc,endc - cg = mycgridcell(c) + cg = mycgridcell(c) - if(cg == g) then + if(cg == g) then - if(newgridcell) then + if(newgridcell) then + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + if (deltax > 180.0) then + deltax = 360.0 - deltax + end if + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then ! Different settings of observation-location-index in ! state vector depending on the method of state @@ -1007,8 +1017,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if - end do end do + end do else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then #ifdef CLMFIVE From 2b7f9570829947e33111cba6fa804ca194403550 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 9 Oct 2025 23:16:53 +0200 Subject: [PATCH 18/36] syntax fix --- interface/framework/init_dim_obs_f_pdaf.F90 | 1 + interface/framework/init_dim_obs_pdaf.F90 | 1 + 2 files changed, 2 insertions(+) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index d9222e96b..7abaf0cac 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1024,6 +1024,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end if + end if end do end do diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 8f3f80f97..d509bd2d8 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1017,6 +1017,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end if + end if end do end do From 6a6706a23a7c63d944cc799393c11dc982efd158 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 9 Oct 2025 23:16:57 +0200 Subject: [PATCH 19/36] line length fixes --- interface/framework/init_dim_obs_f_pdaf.F90 | 9 ++++++--- interface/framework/init_dim_obs_pdaf.F90 | 9 ++++++--- interface/model/eclm/enkf_clm_mod_5.F90 | 3 ++- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 7abaf0cac..ac25b42fe 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1050,7 +1050,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then ! Set index in state vector, LST will be computed ! for first patch appearing here @@ -1079,7 +1080,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 end if @@ -1104,7 +1106,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) end if diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index d509bd2d8..ffddc4832 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1043,7 +1043,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then ! Set index in state vector, LST will be computed ! for first patch appearing here @@ -1072,7 +1073,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 end if @@ -1097,7 +1099,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) deltay = abs(lat(g)-clmobs_lat(i)) end if - if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or.((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) end if diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 5a8e5cb35..ef199dd49 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -643,7 +643,8 @@ subroutine set_clm_statevec(tstartcycle, mype) ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) do lev=1,nlevgrnd - clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), state_pdaf2clm_j_p(cc+lev*clm_varsize)) + clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), & + state_pdaf2clm_j_p(cc+lev*clm_varsize)) end do clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) end do From 993641496dfd904376835e79767fd338fe2e3a7b Mon Sep 17 00:00:00 2001 From: Johannes Keller <16795031+jjokella@users.noreply.github.com> Date: Wed, 29 Oct 2025 10:36:07 +0100 Subject: [PATCH 20/36] CI-fix: handle clmswc_mask_snow as integer (#29) --- 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 ef199dd49..dba9f4044 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -834,7 +834,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc ! If snow is masked, update only, when snow depth is less than 1mm - if( (.not. clmswc_mask_snow) .or. snow_depth(j) < 0.001 ) then + if( (clmswc_mask_snow == 0) .or. snow_depth(j) < 0.001 ) then ! Update only those SWCs that are not excluded by ispval if(state_clm2pdaf_p(j,i) /= ispval) then From 319bd8c6ce0e3599fa050d33ad064c67e19523ea Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Nov 2025 08:17:34 +0100 Subject: [PATCH 21/36] style changes --- interface/model/eclm/enkf_clm_mod_5.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index f701f6386..9c3786a03 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -142,7 +142,7 @@ subroutine define_clm_statevec(mype) if(clmupdate_T/=0) then call define_clm_statevec_T(mype) end if - + !end hcp #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize @@ -170,7 +170,6 @@ subroutine define_clm_statevec(mype) allocate(clm_paramarr(clm_paramsize)) end if - end subroutine define_clm_statevec @@ -561,8 +560,7 @@ subroutine cleanup_clm_statevec() end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) - use clm_instMod, only : soilstate_inst - use clm_instMod, only : waterstate_inst + use clm_instMod, only : soilstate_inst, waterstate_inst use clm_varpar , only : nlevsoi ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite From 35b1dc0c28cb26673ab2139399a30028ea315eff Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Nov 2025 08:18:41 +0100 Subject: [PATCH 22/36] set_clm_statevec_T: debug output to dedicated subroutine --- interface/model/eclm/enkf_clm_mod_5.F90 | 39 +++++++++++++------------ 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9c3786a03..4140227d2 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -573,13 +573,7 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) - real(r8), pointer :: t_grnd(:) - real(r8), pointer :: t_soisno(:,:) - real(r8), pointer :: t_veg(:) - real(r8), pointer :: t_skin(:) - real(r8), pointer :: tlai(:) integer :: i,j,jj,g,c,cc,offset - integer :: lev integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -592,8 +586,6 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col - t_soisno => temperature_inst%t_soisno_col - #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN @@ -605,14 +597,6 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF - IF(clmupdate_T/=0) THEN - ! TSMP-PDAF: Debug output of CLM t_soisno, first layer - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") t_soisno(:,1) - CLOSE(71) - END IF - END IF #endif @@ -627,7 +611,7 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T/=0) then - call set_clm_statevec_T + call set_clm_statevec_T(tstartcycle,mype) endif !end hcp LAI @@ -731,7 +715,7 @@ subroutine set_clm_statevec_swc() end subroutine set_clm_statevec_swc - subroutine set_clm_statevec_T + subroutine set_clm_statevec_T(tstartcycle,mype) use clm_instMod, only : temperature_inst use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevgrnd @@ -739,6 +723,10 @@ subroutine set_clm_statevec_T use shr_kind_mod, only: r8 => shr_kind_r8 implicit none + + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + real(r8), pointer :: t_grnd(:) real(r8), pointer :: t_soisno(:,:) real(r8), pointer :: t_veg(:) @@ -746,6 +734,7 @@ subroutine set_clm_statevec_T real(r8), pointer :: tlai(:) integer :: j,g,cc,c integer :: n_c + integer :: lev ! LST variables t_grnd => temperature_inst%t_grnd_col @@ -755,6 +744,20 @@ subroutine set_clm_statevec_T tlai => canopystate_inst%tlai_patch +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN + + IF(clmupdate_T/=0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + + END IF +#endif + !hcp LAI if(clmupdate_T==1) then do cc = 1, clm_varsize From 58501ab2997e67c79a288128826af8779d8208c7 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Nov 2025 08:19:11 +0100 Subject: [PATCH 23/36] update_clm_T: all declarations to dedicated subroutine --- interface/model/eclm/enkf_clm_mod_5.F90 | 55 +++++++++++++------------ 1 file changed, 28 insertions(+), 27 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 4140227d2..8bc1862f1 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -808,36 +808,19 @@ subroutine set_clm_statevec_T(tstartcycle,mype) end subroutine set_clm_statevec_T subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") - use clm_varpar , only : nlevsoi - use clm_varpar , only : nlevgrnd use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 - use PatchType , only : patch - use ColumnType , only : col - use clm_instMod, only : soilstate_inst use clm_instMod, only : waterstate_inst - use clm_instMod, only : temperature_inst - use clm_varcon , only : denh2o, denice, watmin - use clm_varcon , only : ispval - use clm_varcon , only : spval implicit none integer,intent(in) :: tstartcycle integer,intent(in) :: mype - real(r8), pointer :: t_grnd(:) - real(r8), pointer :: t_soisno(:,:) - real(r8), pointer :: t_veg(:) - real(r8), pointer :: t_skin(:) - - real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - ! integer :: i,j,jj,g,c,p,cc,offset integer :: i - integer :: lev character (len = 31) :: fn !TSMP-PDAF: function name for state vector output character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu @@ -861,13 +844,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col - ! LST - t_grnd => temperature_inst%t_grnd_col - t_soisno => temperature_inst%t_soisno_col - t_veg => temperature_inst%t_veg_patch - t_skin => temperature_inst%t_skin_patch - ! tlai => canopystate_inst%tlai_patch - #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -904,8 +880,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif !hcp: TG, TV - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call update_clm_T(tstartcycle, mype) endif ! end hcp TG, TV @@ -1093,10 +1069,35 @@ subroutine update_clm_swc(tstartcycle, mype) end subroutine update_clm_swc - subroutine update_clm_T + subroutine update_clm_T(tstartcycle,mype) + + use PatchType , only : patch + use clm_varpar , only : nlevgrnd + use clm_instMod, only : temperature_inst implicit none + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + + integer :: i + integer :: j + integer :: c + integer :: p + integer :: lev + + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + + ! LST + t_grnd => temperature_inst%t_grnd_col + t_soisno => temperature_inst%t_soisno_col + t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch + ! tlai => canopystate_inst%tlai_patch + !hcp: TG, TV if(clmupdate_T==1) then do p = clm_begp, clm_endp From 3847db7e48d22a4bdae5bdf5e38807fe3e0c4c8a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Nov 2025 09:42:33 +0100 Subject: [PATCH 24/36] compilation fixes --- interface/model/eclm/enkf_clm_mod_5.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 8bc1862f1..03956ac15 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -367,17 +367,20 @@ subroutine define_clm_statevec_swc(mype) end subroutine define_clm_statevec_swc - subroutine define_clm_statevec_T + subroutine define_clm_statevec_T(mype) use decompMod , only : get_proc_bounds use clm_varpar , only : nlevgrnd use PatchType , only : patch implicit none + integer,intent(in) :: mype + integer :: j integer :: jj integer :: lev integer :: p + integer :: cc integer :: begp, endp ! per-proc beginning and ending pft indices integer :: begc, endc ! per-proc beginning and ending column indices From ee6cc3e308e687707030a033c8d51b6b796c5a8c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 14 Nov 2025 11:27:11 +0100 Subject: [PATCH 25/36] compilation fixes for DEBUG --- interface/model/eclm/enkf_clm_mod_5.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 03956ac15..a2d9eed75 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -739,6 +739,8 @@ subroutine set_clm_statevec_T(tstartcycle,mype) integer :: n_c integer :: lev + character (len = 34) :: fn !TSMP-PDAF: function name for swc output + ! LST variables t_grnd => temperature_inst%t_grnd_col t_veg => temperature_inst%t_veg_patch @@ -752,8 +754,8 @@ subroutine set_clm_statevec_T(tstartcycle,mype) IF(clmupdate_T/=0) THEN ! TSMP-PDAF: Debug output of CLM t_soisno, first layer - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" - OPEN(unit=71, file=fn2, action="write") + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn, action="write") WRITE (71,"(es22.15)") t_soisno(:,1) CLOSE(71) END IF @@ -1094,6 +1096,8 @@ subroutine update_clm_T(tstartcycle,mype) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) + character (len = 31) :: fn !TSMP-PDAF: function name for swc output + ! LST t_grnd => temperature_inst%t_grnd_col t_soisno => temperature_inst%t_soisno_col @@ -1144,8 +1148,8 @@ subroutine update_clm_T(tstartcycle,mype) #ifdef PDAF_DEBUG IF(clmupdate_T/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn, action="write") WRITE (71,"(es22.15)") t_soisno(:,1) CLOSE(71) END IF From 7174f1ba9d61c7e83d4c65af6799f510a9e4faa5 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 20 Nov 2025 16:12:33 +0100 Subject: [PATCH 26/36] add LSTDA of TSKIN/TVEG (clmupdate_T==5) --- interface/framework/init_dim_obs_pdaf.F90 | 2 +- interface/framework/obs_op_pdaf.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 66 +++++++++++++++++++++-- 3 files changed, 64 insertions(+), 6 deletions(-) diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 259bfde80..4855fe470 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1008,7 +1008,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then + else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 9d74ff996..c58595151 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -155,7 +155,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) endif -if (clmupdate_T==2 .OR. clmupdate_T==3 .OR. clmupdate_T==4) then +if (clmupdate_T==2 .OR. clmupdate_T==3 .OR. clmupdate_T==4 .OR. clmupdate_T==5) then lpointobs = .false. diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index a2d9eed75..89ffaa1bc 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -547,6 +547,47 @@ subroutine define_clm_statevec_T(mype) endif + if(clmupdate_T==5) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + ! clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = 2 * (endp-begp+1) !TSKIN and TV + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + + ! Check that cc is state_clm2pdaf_p + if (cc /= state_clm2pdaf_p(p,1)) then + print *, "ERROR: Index problem in clmupdate_T==5" + print *, "ERROR: Differences cc,state_clm2pdaf_p(p,1) = ", cc, state_clm2pdaf_p(p,1) + end if + + end do + + endif + end subroutine define_clm_statevec_T subroutine cleanup_clm_statevec() @@ -790,7 +831,7 @@ subroutine set_clm_statevec_T(tstartcycle,mype) end do endif - ! skin temperature state vector updating soil temperature + ! skin temperature updating state vector with skin, soil and vegetation temperature if(clmupdate_T==3) then do cc = 1, clm_varsize ! t_skin iterated over patches @@ -803,13 +844,22 @@ subroutine set_clm_statevec_T(tstartcycle,mype) end do endif - ! soil temperature state vector updating soil temperature + ! soil temperature updating state vector with soil temperature if(clmupdate_T==4) then do cc = 1, clm_varsize clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) end do endif + ! skin temperature updating state vector with skin and vegetation temperature + if(clmupdate_T==5) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+clm_varsize)) + end do + endif + end subroutine set_clm_statevec_T subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") @@ -1125,7 +1175,7 @@ subroutine update_clm_T(tstartcycle,mype) end do endif - ! skin temperature state vector updating soil temperature + ! skin temperature updating skin, soil and vegetation temperature if(clmupdate_T==3) then do p = clm_begp, clm_endp c = patch%column(p) @@ -1137,7 +1187,7 @@ subroutine update_clm_T(tstartcycle,mype) end do endif - ! soil temperature state vector updating soil temperature + ! soil temperature updating soil temperature if(clmupdate_T==4) then do p = clm_begp, clm_endp c = patch%column(p) @@ -1145,6 +1195,14 @@ subroutine update_clm_T(tstartcycle,mype) end do endif + ! skin temperature updating skin and vegetation temperature + if(clmupdate_T==5) then + do p = clm_begp, clm_endp + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + end do + endif + #ifdef PDAF_DEBUG IF(clmupdate_T/=0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files From 5c307e59b35666ee8b13a7952399d520060ea31c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 11 Dec 2025 12:42:18 +0100 Subject: [PATCH 27/36] LSTDA for local filters: first implementation --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/model/clm3_5/enkf_clm_mod.F90 | 1 + interface/model/eclm/enkf_clm_mod_5.F90 | 123 +++++++++++++++++++- 3 files changed, 120 insertions(+), 6 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 252bcc359..0f17eb18e 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1015,7 +1015,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end do end do - else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then + else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index c05f5a3b0..ad5badf3b 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -60,6 +60,7 @@ module enkf_clm_mod integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) integer,allocatable :: state_loc2clm_c_p(:) + integer,allocatable :: state_loc2clm_p_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 89ffaa1bc..ae3d35e77 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -47,7 +47,9 @@ module enkf_clm_mod integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) + integer,allocatable :: state_pdaf2clm_v_p(:) integer,allocatable :: state_loc2clm_c_p(:) + integer,allocatable :: state_loc2clm_p_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) @@ -478,11 +480,32 @@ subroutine define_clm_statevec_T(mype) if(clmupdate_T==3) then + ! Allocate with full dimension for all variables/layers + ! + ! Here, the second index of STATE_CLM2PDAF_P is NOT simply the + ! layer index of CLM, but rather a variable index of the state + ! vector, in order tomake the index mapping 1:1. IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begp:endp,1)) + allocate(state_clm2pdaf_p(begp:endp,1:(2+nlevgrnd))) + ! ^ + ! dimension layout: + ! 1: TSKIN + ! 2:(1+nlevgrnd): TSOIL layers + ! (2+nlevgrnd): TVEG do p=clm_begp,clm_endp + cc = (p - clm_begp + 1) + + ! TSKIN (variable index 1) state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + + ! TSOIL layers (variable indices 2 to 1+nlevgrnd) + do lev=1,nlevgrnd + state_clm2pdaf_p(p, 1+lev) = cc + lev*clm_varsize + end do + + ! TVEG (variable index 2+nlevgrnd) + state_clm2pdaf_p(p, 2+nlevgrnd) = cc + (1+nlevgrnd)*clm_varsize end do clm_varsize = endp-begp+1 @@ -495,6 +518,8 @@ subroutine define_clm_statevec_T(mype) allocate(state_pdaf2clm_c_p(clm_statevecsize)) IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) allocate(state_pdaf2clm_j_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_v_p)) deallocate(state_pdaf2clm_v_p) + allocate(state_pdaf2clm_v_p(clm_statevecsize)) cc = 0 @@ -503,15 +528,18 @@ subroutine define_clm_statevec_T(mype) state_pdaf2clm_p_p(cc) = p !TSKIN state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_v_p(cc) = 1 do lev=1,nlevgrnd ! ivar = 2-26: TSOIL state_pdaf2clm_p_p(cc + lev*clm_varsize) = p state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + state_pdaf2clm_v_p(cc + lev*clm_varsize) = 1+lev end do state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 + state_pdaf2clm_v_p(cc+(1+nlevgrnd)*clm_varsize) = 2+nlevgrnd end do endif @@ -1181,9 +1209,9 @@ subroutine update_clm_T(tstartcycle,mype) c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) do lev=1,nlevgrnd - t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1) + (lev)*clm_varsize) + t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1+lev)) end do - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + (1+nlevgrnd)*clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,2+levgrnd)) end do endif @@ -1664,6 +1692,7 @@ subroutine init_n_domains_clm(n_domains_p) use decompMod, only : get_proc_bounds use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -1671,6 +1700,7 @@ subroutine init_n_domains_clm(n_domains_p) integer :: domain_p integer :: begg, endg ! per-proc gridcell ending gridcell indices integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begp, endp ! per-proc beginning and ending pft indices integer :: g integer :: c @@ -1690,6 +1720,11 @@ subroutine init_n_domains_clm(n_domains_p) ! -> DIM_L: number of layers in gridcell n_domains_p = endg - begg + 1 end if + else if (clmupdate_T/=0) then + ! For LSTDA: patches are domains + ! Each patch is a local domain + ! -> DIM_L: number of layers in patch + n_domains_p = endp - begp + 1 else ! Process-local number of gridcells Default, possibly not tested ! for other updates except SWC @@ -1704,6 +1739,7 @@ subroutine init_n_domains_clm(n_domains_p) ! local domain domain_p (from the column, the gridcell can be ! derived) + if(clmupdate_swc==1) then ! Allocate state_loc2clm_c_p with preliminary n_domains_p IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) allocate(state_loc2clm_c_p(n_domains_p)) @@ -1720,7 +1756,7 @@ subroutine init_n_domains_clm(n_domains_p) if(clmstatevec_allcol == 1) then ! COLUMNS - ! Each hydrologically active layer is a local domain + ! Each hydrologically active column is a local domain ! -> DIM_L: number of layers in hydrologically active column do c=clm_begc,clm_endc ! Skip state vector loop and directly check if column is @@ -1769,15 +1805,46 @@ subroutine init_n_domains_clm(n_domains_p) else ! GRIDCELLS do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begg + domain_p - 1 + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) end do end if end if + end if ! Possibly: Warning when final n_domains_p actually excludes ! hydrologically inactive gridcells + if(clmupdate_T/=0) then + ! For LSTDA: Patches are domains + + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do + + ! Columns corresponding to patches + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = patch%column(clm_begp + domain_p - 1) + end do + + ! Allocate state_loc2clm_p_p + IF (allocated(state_loc2clm_p_p)) deallocate(state_loc2clm_p_p) + allocate(state_loc2clm_p_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_p_p(domain_p) = ispval + end do + + ! Patches + do domain_p=1,n_domains_p + state_loc2clm_p_p(domain_p) = clm_begp + domain_p - 1 + end do + + end if + + end subroutine init_n_domains_clm @@ -1788,6 +1855,7 @@ end subroutine init_n_domains_clm !> This routine sets DIM_L, the local state vector dimension. subroutine init_dim_l_clm(domain_p, dim_l) use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use ColumnType , only : col implicit none @@ -1822,6 +1890,31 @@ subroutine init_dim_l_clm(domain_p, dim_l) dim_l = 3*nlevsoi + nshift endif + if(clmupdate_T==1) then + ! TG + TV: 2 temperatures per patch + dim_l = 2 + endif + + if(clmupdate_T==2) then + ! TSKIN + TG + TV: 3 temperatures per patch + dim_l = 3 + endif + + if(clmupdate_T==3) then + ! TSKIN + TSOIL(nlevgrnd layers) + TV + dim_l = 1 + nlevgrnd + 1 + endif + + if(clmupdate_T==4) then + ! TSOIL (first layer only) + dim_l = 1 + endif + + if(clmupdate_T==5) then + ! TSKIN + TV (first layer only) + dim_l = 2 + endif + end subroutine init_dim_l_clm !> @author Wolfgang Kurtz, Johannes Keller @@ -1853,11 +1946,21 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) ! ENDDO ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index: i state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) + END DO + end if end subroutine g2l_state_clm @@ -1891,11 +1994,21 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) ! ENDDO ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index i state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) = state_l(i) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) = state_l(i) + END DO + end if end subroutine l2g_state_clm #endif From d1cb4d30f8d15ef51934966b7f25d09c6d312f56 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 11 Dec 2025 13:59:54 +0100 Subject: [PATCH 28/36] bugfix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ``` Error: Symbol ‘levgrnd’ at (1) has no IMPLICIT type; did you mean ‘nlevgrnd’? ``` --- 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 ae3d35e77..1e40acb4d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -1211,7 +1211,7 @@ subroutine update_clm_T(tstartcycle,mype) do lev=1,nlevgrnd t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1+lev)) end do - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,2+levgrnd)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,2+nlevgrnd)) end do endif From afc7202564d8af1d474239c7139b44496c59fc90 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 11 Dec 2025 17:19:28 +0100 Subject: [PATCH 29/36] bugfix: initialize and use begp, endp throughout init_n_domains_clm --- interface/model/eclm/enkf_clm_mod_5.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 1e40acb4d..8ab00c55d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -1708,7 +1708,7 @@ subroutine init_n_domains_clm(n_domains_p) ! TODO: remove unnecessary calls of get_proc_bounds (use clm_begg, ! clm_endg, etc) - call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc) + call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc, begp=begp, endp=endp) if(clmupdate_swc==1) then if(clmstatevec_allcol==1) then @@ -1827,7 +1827,7 @@ subroutine init_n_domains_clm(n_domains_p) ! Columns corresponding to patches do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = patch%column(clm_begp + domain_p - 1) + state_loc2clm_c_p(domain_p) = patch%column(begp + domain_p - 1) end do ! Allocate state_loc2clm_p_p @@ -1839,7 +1839,7 @@ subroutine init_n_domains_clm(n_domains_p) ! Patches do domain_p=1,n_domains_p - state_loc2clm_p_p(domain_p) = clm_begp + domain_p - 1 + state_loc2clm_p_p(domain_p) = begp + domain_p - 1 end do end if From 747dee93671e1756ca5b2204dc25ca79cbb4ccc9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 19 Dec 2025 13:31:59 +0100 Subject: [PATCH 30/36] correct condition for tsoisno_mype.update.* debug output Now the update is written, according to the input `CLM:t_printensemble` from `enkfpf.par`. --- 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 543641f6f..9d59d484d 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -1205,13 +1205,13 @@ subroutine update_clm_T(tstartcycle,mype) endif #ifdef PDAF_DEBUG - IF(clmupdate_T/=0) THEN + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" OPEN(unit=71, file=fn, action="write") WRITE (71,"(es22.15)") t_soisno(:,1) CLOSE(71) - END IF + END IF #endif end subroutine update_clm_T From 38ffe0e739df9d7421e890776469a0a41f68445a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 23 Jan 2026 17:09:55 +0100 Subject: [PATCH 31/36] LSTDA: remove unused variants --- .../running_tsmp_pdaf/input_enkfpf.md | 6 +- interface/framework/init_dim_obs_f_pdaf.F90 | 2 +- interface/framework/init_dim_obs_pdaf.F90 | 2 +- interface/framework/localize_covar_pdaf.F90 | 1 - interface/framework/obs_op_pdaf.F90 | 2 +- interface/model/eclm/enkf_clm_mod_5.F90 | 168 +----------------- 6 files changed, 12 insertions(+), 169 deletions(-) diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index 790b4aa09..6ba12242e 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -484,12 +484,14 @@ CLM (standalone only). `CLM:update_T`: (integer) Flag for updating of ground and vegetation temperature. -Currently only CLM3.5 - - 0: No update of ground and vegetation temperature - 1: Update of ground and vegetation temperature +- 2: Update of soil and vegetation temperature (prognostic variables) + based on comparison of satellite LST with skin temperature + (diagnostic variable) + ### CLM:print_swc ### `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 252bcc359..922ce1691 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -1015,7 +1015,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end do end do - else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then + else if(clmupdate_T==1 .or. clmupdate_T==2) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 4855fe470..f25490c74 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/interface/framework/init_dim_obs_pdaf.F90 @@ -1008,7 +1008,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) end do end do - else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + else if(clmupdate_T==1 .or. clmupdate_T==2) then #ifdef CLMFIVE ! patch loop do g = begg,endg diff --git a/interface/framework/localize_covar_pdaf.F90 b/interface/framework/localize_covar_pdaf.F90 index b0865fb8b..eb93a453c 100644 --- a/interface/framework/localize_covar_pdaf.F90 +++ b/interface/framework/localize_covar_pdaf.F90 @@ -40,7 +40,6 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH) USE shr_kind_mod , only : r8 => shr_kind_r8 USE mod_read_obs, ONLY: clmobs_lon USE mod_read_obs, ONLY: clmobs_lat - USE enkf_clm_mod, ONLY: clmupdate_T USE enkf_clm_mod, ONLY: clm_begc USE enkf_clm_mod, ONLY: clm_endc USE enkf_clm_mod, ONLY: state_pdaf2clm_c_p diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index c58595151..78d3e05cf 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/interface/framework/obs_op_pdaf.F90 @@ -155,7 +155,7 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) endif -if (clmupdate_T==2 .OR. clmupdate_T==3 .OR. clmupdate_T==4 .OR. clmupdate_T==5) then +if (clmupdate_T==2) then lpointobs = .false. diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9d59d484d..63942b2d3 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -439,7 +439,6 @@ subroutine define_clm_statevec_T(mype) endif !end hcp - ! skin temperature state vector if(clmupdate_T==2) then IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) @@ -449,43 +448,6 @@ subroutine define_clm_statevec_T(mype) state_clm2pdaf_p(p,1) = (p - clm_begp + 1) end do - clm_varsize = endp-begp+1 - ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = 3* (endp-begp+1) !TSKIN, then TG and TV - - IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) - allocate(state_pdaf2clm_p_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) - - cc = 0 - - do p=clm_begp,clm_endp - cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TSKIN - state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN - state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TG - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TG - state_pdaf2clm_j_p(cc+clm_varsize) = 1 - state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV - state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 - end do - - endif - - if(clmupdate_T==3) then - - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begp:endp,1)) - - do p=clm_begp,clm_endp - state_clm2pdaf_p(p,1) = (p - clm_begp + 1) - end do - clm_varsize = endp-begp+1 ! clm_paramsize = endp-begp+1 !LAI clm_statevecsize = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and TV @@ -517,77 +479,9 @@ subroutine define_clm_statevec_T(mype) endif - if(clmupdate_T==4) then - - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begp:endp,1)) - - do p=clm_begp,clm_endp - state_clm2pdaf_p(p,1) = (p - clm_begp + 1) - end do - - clm_varsize = endp-begp+1 - ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = 1* (endp-begp+1) !TSOIL - - IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) - allocate(state_pdaf2clm_p_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) - - cc = 0 - - do p=clm_begp,clm_endp - cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TSOIL - state_pdaf2clm_c_p(cc) = patch%column(p) !TSOIL - state_pdaf2clm_j_p(cc) = 1 - end do - - endif - - if(clmupdate_T==5) then - - IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) - allocate(state_clm2pdaf_p(begp:endp,1)) - - do p=clm_begp,clm_endp - state_clm2pdaf_p(p,1) = (p - clm_begp + 1) - end do - - clm_varsize = endp-begp+1 - ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = 2 * (endp-begp+1) !TSKIN and TV - - IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) - allocate(state_pdaf2clm_p_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) - - cc = 0 - - do p=clm_begp,clm_endp - cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TSKIN - state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN - state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_p_p(cc+clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV - state_pdaf2clm_j_p(cc+clm_varsize) = 1 - - ! Check that cc is state_clm2pdaf_p - if (cc /= state_clm2pdaf_p(p,1)) then - print *, "ERROR: Index problem in clmupdate_T==5" - print *, "ERROR: Differences cc,state_clm2pdaf_p(p,1) = ", cc, state_clm2pdaf_p(p,1) - end if - - end do - - endif + if(clmupdate_T < 1 .or. clmupdate_T > 2) then + error stop "Only implemented for clmupdate_T==1 and clmupdate_T==2." + end if end subroutine define_clm_statevec_T @@ -822,18 +716,8 @@ subroutine set_clm_statevec_T(tstartcycle,mype) endif !end hcp LAI - ! skin temperature state vector - if(clmupdate_T==2) then - do cc = 1, clm_varsize - ! t_skin iterated over patches - clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) - clm_statevec(cc+clm_varsize) = t_grnd(state_pdaf2clm_c_p(cc+clm_varsize)) - clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) - end do - endif - ! skin temperature updating state vector with skin, soil and vegetation temperature - if(clmupdate_T==3) then + if(clmupdate_T==2) then do cc = 1, clm_varsize ! t_skin iterated over patches clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) @@ -845,22 +729,6 @@ subroutine set_clm_statevec_T(tstartcycle,mype) end do endif - ! soil temperature updating state vector with soil temperature - if(clmupdate_T==4) then - do cc = 1, clm_varsize - clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) - end do - endif - - ! skin temperature updating state vector with skin and vegetation temperature - if(clmupdate_T==5) then - do cc = 1, clm_varsize - ! t_skin iterated over patches - clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) - clm_statevec(cc+clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+clm_varsize)) - end do - endif - end subroutine set_clm_statevec_T subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") @@ -1166,18 +1034,8 @@ subroutine update_clm_T(tstartcycle,mype) endif ! end hcp TG, TV - ! skin temperature state vector - if(clmupdate_T==2) then - do p = clm_begp, clm_endp - c = patch%column(p) - t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) - t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) - end do - endif - ! skin temperature updating skin, soil and vegetation temperature - if(clmupdate_T==3) then + if(clmupdate_T==2) then do p = clm_begp, clm_endp c = patch%column(p) t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) @@ -1188,22 +1046,6 @@ subroutine update_clm_T(tstartcycle,mype) end do endif - ! soil temperature updating soil temperature - if(clmupdate_T==4) then - do p = clm_begp, clm_endp - c = patch%column(p) - t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) - end do - endif - - ! skin temperature updating skin and vegetation temperature - if(clmupdate_T==5) then - do p = clm_begp, clm_endp - t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) - end do - endif - #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN ! TSMP-PDAF: For debug runs, output the state vector in files From e67133a7d8c3a7fe5c5c55e35ef51e931c3baeaf Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 23 Jan 2026 17:20:04 +0100 Subject: [PATCH 32/36] remove unused index array --- interface/model/eclm/enkf_clm_mod_5.F90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 7360bb4aa..a306ef6ea 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -48,7 +48,6 @@ module enkf_clm_mod integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) - integer,allocatable :: state_pdaf2clm_v_p(:) integer,allocatable :: state_loc2clm_c_p(:) integer,allocatable :: state_loc2clm_p_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model @@ -481,8 +480,6 @@ subroutine define_clm_statevec_T(mype) allocate(state_pdaf2clm_c_p(clm_statevecsize)) IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) allocate(state_pdaf2clm_j_p(clm_statevecsize)) - IF (allocated(state_pdaf2clm_v_p)) deallocate(state_pdaf2clm_v_p) - allocate(state_pdaf2clm_v_p(clm_statevecsize)) cc = 0 @@ -491,18 +488,15 @@ subroutine define_clm_statevec_T(mype) state_pdaf2clm_p_p(cc) = p !TSKIN state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_v_p(cc) = 1 do lev=1,nlevgrnd ! ivar = 2-26: TSOIL state_pdaf2clm_p_p(cc + lev*clm_varsize) = p state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev - state_pdaf2clm_v_p(cc + lev*clm_varsize) = 1+lev end do state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 - state_pdaf2clm_v_p(cc+(1+nlevgrnd)*clm_varsize) = 2+nlevgrnd end do endif From 8eee11054a4c8c3f1573fe27954e4fb7177cefa5 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 23 Jan 2026 17:26:24 +0100 Subject: [PATCH 33/36] small change --- 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 a306ef6ea..9645f96aa 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -459,7 +459,7 @@ subroutine define_clm_statevec_T(mype) cc = (p - clm_begp + 1) ! TSKIN (variable index 1) - state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + state_clm2pdaf_p(p,1) = cc ! TSOIL layers (variable indices 2 to 1+nlevgrnd) do lev=1,nlevgrnd From 249c4ad90dfd424160cedbd73990e2c708bc0f48 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Sat, 24 Jan 2026 23:14:36 +0100 Subject: [PATCH 34/36] remove unused variable `pc` --- interface/framework/init_dim_obs_f_pdaf.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index 922ce1691..ea2d7531f 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/interface/framework/init_dim_obs_f_pdaf.F90 @@ -172,7 +172,6 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) INTEGER :: g ! CLM Gridcell index INTEGER :: cg INTEGER :: pg - INTEGER :: pc INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -1024,7 +1023,6 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) do p = begp,endp pg = patch%gridcell(p) - pc = patch%column(p) if(pg == g) then if(newgridcell) then From 1203bae8c8934358bbe8ba538decf20ab173ae2c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Sat, 24 Jan 2026 23:15:28 +0100 Subject: [PATCH 35/36] remove unused clmupdate_T values --- interface/model/eclm/enkf_clm_mod_5.F90 | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 9645f96aa..a96336701 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -1733,25 +1733,10 @@ subroutine init_dim_l_clm(domain_p, dim_l) endif if(clmupdate_T==2) then - ! TSKIN + TG + TV: 3 temperatures per patch - dim_l = 3 - endif - - if(clmupdate_T==3) then ! TSKIN + TSOIL(nlevgrnd layers) + TV dim_l = 1 + nlevgrnd + 1 endif - if(clmupdate_T==4) then - ! TSOIL (first layer only) - dim_l = 1 - endif - - if(clmupdate_T==5) then - ! TSKIN + TV (first layer only) - dim_l = 2 - endif - end subroutine init_dim_l_clm !> @author Wolfgang Kurtz, Johannes Keller @@ -1791,7 +1776,7 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) END DO end if - if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + if(clmupdate_T==1 .or. clmupdate_T==2) then DO i = 1, dim_l ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P ! Variable index: i @@ -1839,7 +1824,7 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) END DO end if - if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then + if(clmupdate_T==1 .or. clmupdate_T==2) then DO i = 1, dim_l ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P ! Variable index: i From b3d04151c9b366f2c27ecf8d81b33f2ef7b139d9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Sat, 24 Jan 2026 23:18:23 +0100 Subject: [PATCH 36/36] LSTDA: remodel clmupdate_T==2 with gridcell-averages Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- interface/model/eclm/enkf_clm_mod_5.F90 | 357 +++++++++++++++++------- 1 file changed, 253 insertions(+), 104 deletions(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index a96336701..195d71688 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -160,9 +160,9 @@ subroutine define_clm_statevec(mype) ! Allocate statevector-duplicate for saving original column mean ! values used in computing increments during updating the state - ! vector in column-mean-mode. + ! vector in column-mean-mode (SWC) or gridcell-mean-mode (T). IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) - if (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) then + if ((clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_T==2) then allocate(clm_statevec_orig(clm_statevecsize)) end if @@ -372,6 +372,7 @@ end subroutine define_clm_statevec_swc subroutine define_clm_statevec_T(mype) use decompMod , only : get_proc_bounds use clm_varpar , only : nlevgrnd + use clm_varcon , only : ispval use PatchType , only : patch implicit none @@ -425,6 +426,13 @@ subroutine define_clm_statevec_T(mype) IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + cc = 0 do p=clm_begp,clm_endp @@ -442,6 +450,7 @@ subroutine define_clm_statevec_T(mype) if(clmupdate_T==2) then + ! 1) PATCH/GRC: CLM->PDAF ! Allocate with full dimension for all variables/layers ! ! Here, the second index of STATE_CLM2PDAF_P is NOT simply the @@ -454,26 +463,49 @@ subroutine define_clm_statevec_T(mype) ! 1: TSKIN ! 2:(1+nlevgrnd): TSOIL layers ! (2+nlevgrnd): TVEG + do lev=1,(2+nlevgrnd) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do do p=clm_begp,clm_endp - cc = (p - clm_begp + 1) + ! All patches in a gridcell are assigned the index of + ! the gridcell-average in the state vector + cc = (patch%gridcell(p) - clm_begg + 1) ! TSKIN (variable index 1) state_clm2pdaf_p(p,1) = cc ! TSOIL layers (variable indices 2 to 1+nlevgrnd) do lev=1,nlevgrnd - state_clm2pdaf_p(p, 1+lev) = cc + lev*clm_varsize + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) end do ! TVEG (variable index 2+nlevgrnd) - state_clm2pdaf_p(p, 2+nlevgrnd) = cc + (1+nlevgrnd)*clm_varsize + state_clm2pdaf_p(p, 2+nlevgrnd) = cc + (1+nlevgrnd)*(clm_endg - clm_begg + 1) end do - clm_varsize = endp-begp+1 - ! clm_paramsize = endp-begp+1 !LAI - clm_statevecsize = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and TV + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (1 + nlevgrnd + 1)* (endg-begg+1) !TSKIN, then nlevgrnd times TSOIL and TV + ! 3) PATCH/GRC: PDAF->CLM + ! + ! Inverse mapping from state vector index to CLM indices. + ! STATE_PDAF2CLM_*_P is the inverse of STATE_CLM2PDAF_P: + ! state_pdaf2clm_p_p(cc) = p (patch index) + ! state_pdaf2clm_c_p(cc) = c (column index) + ! state_pdaf2clm_j_p(cc) = j (variable index, NOT CLM level) + ! + ! Variable index j layout (same as STATE_CLM2PDAF_P second dimension): + ! j = 1: TSKIN + ! j = 2:(1+nlevgrnd): TSOIL layers (CLM level = j - 1) + ! j = (2+nlevgrnd): TVEG + ! + ! NOTE: To get the CLM level from j for TSOIL, use: clm_level = j - 1 + ! IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) allocate(state_pdaf2clm_p_p(clm_statevecsize)) IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) @@ -481,28 +513,43 @@ subroutine define_clm_statevec_T(mype) IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) allocate(state_pdaf2clm_j_p(clm_statevecsize)) + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + cc = 0 - do p=clm_begp,clm_endp - cc = cc + 1 - state_pdaf2clm_p_p(cc) = p !TSKIN - state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN - state_pdaf2clm_j_p(cc) = 1 - do lev=1,nlevgrnd - ! ivar = 2-26: TSOIL - state_pdaf2clm_p_p(cc + lev*clm_varsize) = p - state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) - state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + do cc=1,clm_statevecsize + + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + ! Set indices for the three temperature variables and then + ! exit loops + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,nlevgrnd + ! variable index 2 to 1+nlevgrnd: TSOIL layers + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 2 + nlevgrnd ! variable index: TVEG + exit + end if end do - state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV - state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV - state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 1 + end do endif if(clmupdate_T < 1 .or. clmupdate_T > 2) then - error stop "Only implemented for clmupdate_T==1 and clmupdate_T==2." + error stop "LST-DA only implemented for clmupdate_T==1 and clmupdate_T==2." end if end subroutine define_clm_statevec_T @@ -681,6 +728,7 @@ subroutine set_clm_statevec_T(tstartcycle,mype) use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevgrnd use PatchType , only : patch + use ColumnType , only : col use shr_kind_mod, only: r8 => shr_kind_r8 implicit none @@ -693,8 +741,8 @@ subroutine set_clm_statevec_T(tstartcycle,mype) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) real(r8), pointer :: tlai(:) - integer :: j,g,cc,c - integer :: n_c + integer :: j,g,cc,c,p + integer :: n_c,n_p integer :: lev character (len = 34) :: fn !TSMP-PDAF: function name for swc output @@ -738,17 +786,74 @@ subroutine set_clm_statevec_T(tstartcycle,mype) endif !end hcp LAI - ! skin temperature updating state vector with skin, soil and vegetation temperature + ! Skin temperature updating state vector with skin, soil and vegetation temperature. + ! Uses gridcell mean: averages over all patches (for TSKIN, TVEG) or + ! columns (for TSOIL) within each gridcell. if(clmupdate_T==2) then + do cc = 1, clm_varsize - ! t_skin iterated over patches - clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + + ! Get gridcell from the reference patch + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- do lev=1,nlevgrnd - clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), & - state_pdaf2clm_j_p(cc+lev*clm_varsize)) + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if end do - clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) + if(n_p > 0) then + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TVEG averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TVEG)" + end if + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + endif end subroutine set_clm_statevec_T @@ -1020,6 +1125,7 @@ subroutine update_clm_T(tstartcycle,mype) use PatchType , only : patch use clm_varpar , only : nlevgrnd use clm_instMod, only : temperature_inst + use IEEE_ARITHMETIC, only: ieee_is_nan implicit none @@ -1030,6 +1136,7 @@ subroutine update_clm_T(tstartcycle,mype) integer :: j integer :: c integer :: p + integer :: cc integer :: lev real(r8), pointer :: t_grnd(:) @@ -1037,6 +1144,9 @@ subroutine update_clm_T(tstartcycle,mype) real(r8), pointer :: t_veg(:) real(r8), pointer :: t_skin(:) + real(r8) :: increment_factor + real(r8) :: t_update + character (len = 31) :: fn !TSMP-PDAF: function name for swc output ! LST @@ -1056,15 +1166,54 @@ subroutine update_clm_T(tstartcycle,mype) endif ! end hcp TG, TV - ! skin temperature updating skin, soil and vegetation temperature + ! Skin temperature updating skin, soil and vegetation temperature. + ! Uses gridcell mean increment factor: applies the ratio of + ! (new gridcell mean / old gridcell mean) to each patch/column value. if(clmupdate_T==2) then do p = clm_begp, clm_endp c = patch%column(p) - t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = t_update + end if + end if + + ! --- TSOIL: update with increment factor for each layer --- do lev=1,nlevgrnd - t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1+lev)) + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = t_update + end if + end if end do - t_veg(p) = clm_statevec(state_clm2pdaf_p(p,2+nlevgrnd)) + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+nlevgrnd) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = t_update + end if + end if + end do endif @@ -1558,10 +1707,11 @@ subroutine init_n_domains_clm(n_domains_p) n_domains_p = endg - begg + 1 end if else if (clmupdate_T/=0) then - ! For LSTDA: patches are domains - ! Each patch is a local domain - ! -> DIM_L: number of layers in patch - n_domains_p = endp - begp + 1 + ! For LSTDA: gridcells are domains + ! Each gridcell is a local domain + ! -> DIM_L: number of temperature variables (each with gridcell + ! -> averages) + n_domains_p = endg - begg + 1 else ! Process-local number of gridcells Default, possibly not tested ! for other updates except SWC @@ -1577,83 +1727,83 @@ subroutine init_n_domains_clm(n_domains_p) ! derived) if(clmupdate_swc==1) then - ! Allocate state_loc2clm_c_p with preliminary n_domains_p - IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) - allocate(state_loc2clm_c_p(n_domains_p)) - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = ispval - end do - - if(clmstatevec_only_active == 1) then - - ! Reset n_domains_p - n_domains_p = 0 - domain_p = 0 - - if(clmstatevec_allcol == 1) then - ! COLUMNS - - ! Each hydrologically active column is a local domain - ! -> DIM_L: number of layers in hydrologically active column - do c=clm_begc,clm_endc - ! Skip state vector loop and directly check if column is - ! hydrologically active - if(col%hydrologically_active(c)) then - domain_p = domain_p + 1 - n_domains_p = n_domains_p + 1 - state_loc2clm_c_p(domain_p) = c - end if - end do + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do - else - ! GRIDCELLS + if(clmstatevec_only_active == 1) then - ! For gridcells - do g = clm_begg,clm_endg + ! Reset n_domains_p + n_domains_p = 0 + domain_p = 0 - ! Search the state vector for col in grc - do cc = 1,clm_statevecsize + if(clmstatevec_allcol == 1) then + ! COLUMNS - if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then - ! Set local domain index + ! Each hydrologically active column is a local domain + ! -> DIM_L: number of layers in hydrologically active column + do c=clm_begc,clm_endc + ! Skip state vector loop and directly check if column is + ! hydrologically active + if(col%hydrologically_active(c)) then domain_p = domain_p + 1 - ! Set new number of local domains n_domains_p = n_domains_p + 1 - ! Set CLM-column-index corresponding to local domain - state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) - ! Exit state vector loop, when fitting column is found - exit + state_loc2clm_c_p(domain_p) = c end if end do - end do + else + ! GRIDCELLS + + ! For gridcells + do g = clm_begg,clm_endg + + ! Search the state vector for col in grc + do cc = 1,clm_statevecsize + + if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then + ! Set local domain index + domain_p = domain_p + 1 + ! Set new number of local domains + n_domains_p = n_domains_p + 1 + ! Set CLM-column-index corresponding to local domain + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) + ! Exit state vector loop, when fitting column is found + exit + end if + end do - end if + end do - else + end if - ! Set state_loc2clm_c_p for non-excluding hydrologically - ! inactive cols/grcs - if(clmstatevec_allcol == 1) then - ! COLUMNS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 - end do else - ! GRIDCELLS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) - end do - end if - end if + ! Set state_loc2clm_c_p for non-excluding hydrologically + ! inactive cols/grcs + if(clmstatevec_allcol == 1) then + ! COLUMNS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 + end do + else + ! GRIDCELLS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + end do + end if + + end if end if ! Possibly: Warning when final n_domains_p actually excludes ! hydrologically inactive gridcells if(clmupdate_T/=0) then - ! For LSTDA: Patches are domains + ! For LSTDA: Gridcells are domains ! Allocate state_loc2clm_c_p with preliminary n_domains_p IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) @@ -1661,22 +1811,21 @@ subroutine init_n_domains_clm(n_domains_p) do domain_p=1,n_domains_p state_loc2clm_c_p(domain_p) = ispval end do - - ! Columns corresponding to patches - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = patch%column(begp + domain_p - 1) - end do - - ! Allocate state_loc2clm_p_p IF (allocated(state_loc2clm_p_p)) deallocate(state_loc2clm_p_p) allocate(state_loc2clm_p_p(n_domains_p)) do domain_p=1,n_domains_p state_loc2clm_p_p(domain_p) = ispval end do - ! Patches + ! Columns from gridcells + ! Patches from gridcells + ! + ! domain_p is equal to cc for the first variable TSKIN + ! Therefore, state_pdaf2clm_c_p / state_pdaf2clm_p_p can be used + ! to derive the corresponding column / patch do domain_p=1,n_domains_p - state_loc2clm_p_p(domain_p) = begp + domain_p - 1 + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + state_loc2clm_p_p(domain_p) = state_pdaf2clm_p_p(domain_p) end do end if