From c49edbe9f95e458fde5220e7e664a72cbc2ef459 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Dec 2024 09:37:21 +0100 Subject: [PATCH 01/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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 747dee93671e1756ca5b2204dc25ca79cbb4ccc9 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 19 Dec 2025 13:31:59 +0100 Subject: [PATCH 27/27] 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