From d0802b6434cb906a2f2674dda6addb5cf1e2af82 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 3 Dec 2024 17:08:22 +0100 Subject: [PATCH 01/21] PDAF: LST-update for eCLM-PDAF --- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 72 +++++++++++++++++-- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 72 +++++++++++++++++-- 2 files changed, 130 insertions(+), 14 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index f480577b4..5cf4682d5 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -96,6 +96,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: p integer :: c integer :: g integer :: cc @@ -273,7 +274,33 @@ 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_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_c_p(cc) = p !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_c_p(cc+clm_varsize) = p !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + endif !end hcp @@ -292,7 +319,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 @@ -310,11 +337,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 @@ -332,6 +363,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 @@ -360,7 +397,16 @@ 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 + clm_statevec(cc) = t_grnd(patch%column(state_pdaf2clm_c_p(cc))) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_c_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 @@ -403,8 +449,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 @@ -428,7 +477,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 @@ -460,6 +509,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 @@ -606,7 +660,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index f480577b4..5cf4682d5 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -96,6 +96,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: p integer :: c integer :: g integer :: cc @@ -273,7 +274,33 @@ 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_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_c_p(cc) = p !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_c_p(cc+clm_varsize) = p !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + endif !end hcp @@ -292,7 +319,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 @@ -310,11 +337,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 @@ -332,6 +363,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 @@ -360,7 +397,16 @@ 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 + clm_statevec(cc) = t_grnd(patch%column(state_pdaf2clm_c_p(cc))) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_c_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 @@ -403,8 +449,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 @@ -428,7 +477,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 @@ -460,6 +509,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 @@ -606,7 +660,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 bc1ad979cc3c062a10af2b7c77a92af52819c524 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 4 Dec 2024 11:47:28 +0100 Subject: [PATCH 02/21] PDAF: LST update, first implementation init_dim_obs / obs_op --- .../pdaf/framework/init_dim_obs_f_pdaf.F90 | 159 +++++++++++++----- .../pdaf/framework/init_dim_obs_pdaf.F90 | 159 +++++++++++++----- bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 | 14 +- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 16 +- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 16 +- 5 files changed, 270 insertions(+), 94 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index d3f1484d3..9292d4b52 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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) + + !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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index e1a6ee4e0..a0ecf4e81 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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) + + !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/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index 97becfc66..84088f217 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 @@ -133,13 +133,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 & @@ -151,6 +152,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/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 5cf4682d5..c2692c173 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/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) @@ -88,6 +89,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 @@ -286,6 +288,8 @@ 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) + 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) @@ -295,9 +299,11 @@ subroutine define_clm_statevec(mype) do p=clm_begp,clm_endp cc = cc + 1 - state_pdaf2clm_c_p(cc) = p !TG + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_c_p(cc+clm_varsize) = p !TV + 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 @@ -398,8 +404,10 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then do cc = 1, clm_varsize - clm_statevec(cc) = t_grnd(patch%column(state_pdaf2clm_c_p(cc))) - clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_c_p(cc+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 diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 5cf4682d5..c2692c173 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/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) @@ -88,6 +89,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 @@ -286,6 +288,8 @@ 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) + 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) @@ -295,9 +299,11 @@ subroutine define_clm_statevec(mype) do p=clm_begp,clm_endp cc = cc + 1 - state_pdaf2clm_c_p(cc) = p !TG + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG state_pdaf2clm_j_p(cc) = 1 - state_pdaf2clm_c_p(cc+clm_varsize) = p !TV + 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 @@ -398,8 +404,10 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then do cc = 1, clm_varsize - clm_statevec(cc) = t_grnd(patch%column(state_pdaf2clm_c_p(cc))) - clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_c_p(cc+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 From e60d7481f35bc59eae2b6b2ccf697febf4615cfb Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 4 Dec 2024 13:53:42 +0100 Subject: [PATCH 03/21] bugfix: unbalanced parentheses clm5_0/enkf_clm_mod_5.F90(404): error #5276: Unbalanced parentheses --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index c2692c173..ff13fea88 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -406,7 +406,7 @@ subroutine set_clm_statevec(tstartcycle, mype) 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) = t_grnd(state_pdaf2clm_c_p(cc)) clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) end do diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index c2692c173..ff13fea88 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -406,7 +406,7 @@ subroutine set_clm_statevec(tstartcycle, mype) 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) = t_grnd(state_pdaf2clm_c_p(cc)) clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) end do From 6cb013648a0c1cb005948d279e9bbba88470298b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 4 Dec 2024 16:24:45 +0100 Subject: [PATCH 04/21] bugfix: missing pointer declarations ``` clm5_0/enkf_clm_mod_5.F90(372): error #6404: This name does not have a type, and must have an explicit type. [T_GRND] ``` --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 6 ++++++ bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index ff13fea88..8157d513b 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -360,6 +360,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 @@ -477,6 +480,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(:,:) diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index ff13fea88..8157d513b 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -360,6 +360,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 @@ -477,6 +480,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(:,:) From f10dba9affae217be101f6cab48901e1b4ab2c6d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 4 Dec 2024 17:28:48 +0100 Subject: [PATCH 05/21] bugfix: incorrect number of subscripts ``` init_dim_obs_f_pdaf.F90(1022): error #6351: The number of subscripts is incorrect. [STATE_CLM2PDAF_P] ``` --- bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 | 2 +- bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 9292d4b52..60c2e75c3 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 @@ -1019,7 +1019,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) ! Set index in state vector, LST will be computed ! for first patch appearing here - obs_index_p(cnt) = state_clm2pdaf_p(p) + 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) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index a0ecf4e81..4a5635a23 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 @@ -1012,7 +1012,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) ! Set index in state vector, LST will be computed ! for first patch appearing here - obs_index_p(cnt) = state_clm2pdaf_p(p) + 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) From b97ed6e3a065bb1671d99e5d91b696e3608e7abb Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Dec 2024 09:52:26 +0100 Subject: [PATCH 06/21] bugfix: re-introduce default for general update gridcell loop with layer information --- .../pdaf/framework/init_dim_obs_f_pdaf.F90 | 24 +++++++++++++++++-- .../pdaf/framework/init_dim_obs_pdaf.F90 | 24 +++++++++++++++++-- 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 60c2e75c3..589aa4085 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 4a5635a23..90af36521 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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 363a2df435332cc213d08247786fe5242158b93a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 6 Dec 2024 17:35:03 +0100 Subject: [PATCH 07/21] doc: restriction to CLM3.5 no longer needed --- doc/content/setup_tsmp/input_enkfpf.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index 4d1082a1a..ccdd64c9b 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -469,10 +469,8 @@ CLM (standalone only). ### CLM:update_T ### -`CLM:update_T`: (integer) Flag for updating of ground and vegetation -temperature. - -Currently only CLM3.5 +`CLM:update_T`: (integer) Flag for updating of CLM's ground and +vegetation temperature. - 0: No update of ground and vegetation temperature From 4d2db0fc728ea685bc8e7f283d51610598ab70e6 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 13:49:21 +0100 Subject: [PATCH 08/21] bugfix: lai is a patch-array --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 8157d513b..0725aedfe 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -416,7 +416,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 8157d513b..0725aedfe 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -416,7 +416,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 27a019aedf528b4b2ebfbe5b934138c1d5847762 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 13:56:54 +0100 Subject: [PATCH 09/21] bugfix: typo --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 2 +- bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 0725aedfe..1b4b72775 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -288,7 +288,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 0725aedfe..1b4b72775 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -288,7 +288,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 6938117b744b6100c453dab75752694199e5c0c7 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 12 Dec 2024 14:17:48 +0100 Subject: [PATCH 10/21] dev: skin temperature state vector, first implementation --- bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 | 11 ++++ .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 56 ++++++++++++++++++- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 56 ++++++++++++++++++- 3 files changed, 121 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index 84088f217..315f1dc7c 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 @@ -153,6 +153,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 diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 1b4b72775..2bc8cf1c5 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -310,6 +310,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 @@ -324,7 +360,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 @@ -375,6 +411,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 @@ -421,6 +458,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" @@ -526,6 +571,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 @@ -682,6 +728,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 1b4b72775..2bc8cf1c5 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -310,6 +310,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 @@ -324,7 +360,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 @@ -375,6 +411,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 @@ -421,6 +458,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" @@ -526,6 +571,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 @@ -682,6 +728,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 dddec953b949be9ec599deb71a306b79fc0b4c35 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 13 Dec 2024 09:43:34 +0100 Subject: [PATCH 11/21] 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] ``` --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 2 ++ bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 2bc8cf1c5..943b7baf2 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -398,6 +398,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 @@ -527,6 +528,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 2bc8cf1c5..943b7baf2 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -398,6 +398,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 @@ -527,6 +528,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 eb32d882f8256bb97512a922049cd7c354bdd721 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 08:20:29 +0100 Subject: [PATCH 12/21] 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. --- bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 | 4 ++-- bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 589aa4085..9581edfd8 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 90af36521..6ab7fc2cc 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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 3970ec9540638618146a840fd5016bcd5258955e Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 08:47:50 +0100 Subject: [PATCH 13/21] LST-DA with TSKIN: use same obs_index_p setting as for TG/TV --- bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 | 2 +- bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 9581edfd8..98d5abb3b 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 6ab7fc2cc..9220d106c 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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 572e60450e1fe6125fbbfa38a66ceabb716bba56 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 09:13:32 +0100 Subject: [PATCH 14/21] LST-DA: state-vector with TSKIN, plus TG and TV --- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 21 ++++++++++++------- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 21 ++++++++++++------- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 943b7baf2..f69b035d1 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -322,7 +322,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)) @@ -335,12 +335,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 @@ -461,9 +464,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 @@ -735,6 +740,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 943b7baf2..f69b035d1 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -322,7 +322,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)) @@ -335,12 +335,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 @@ -461,9 +464,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 @@ -735,6 +740,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 e19b7aafbb731c902d2a4ddb9158533edf04be96 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 09:36:44 +0100 Subject: [PATCH 15/21] LST-DA: state-vector with TSKIN, plus TSOIL and TV `clmupdate_T.eq.3` --- .../pdaf/framework/init_dim_obs_f_pdaf.F90 | 2 +- .../pdaf/framework/init_dim_obs_pdaf.F90 | 2 +- bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 | 2 +- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 63 ++++++++++++++++++- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 63 ++++++++++++++++++- 5 files changed, 127 insertions(+), 5 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 98d5abb3b..070482313 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 9220d106c..0fe391a9e 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index 315f1dc7c..abc9f0d57 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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) then +if (clmupdate_T.EQ.2 .OR. clmupdate_T.EQ.3) then lpointobs = .false. diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index f69b035d1..f9bea2758 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -348,6 +348,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 @@ -400,6 +437,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(:) @@ -412,10 +450,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 @@ -472,6 +511,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" @@ -532,6 +581,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(:) @@ -577,6 +627,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 @@ -745,6 +796,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index f69b035d1..f9bea2758 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -348,6 +348,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 @@ -400,6 +437,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(:) @@ -412,10 +450,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 @@ -472,6 +511,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" @@ -532,6 +581,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(:) @@ -577,6 +627,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 @@ -745,6 +796,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 730bc99f6e78d71cbf3907480b7a7d2ffd34cabb Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 10:28:56 +0100 Subject: [PATCH 16/21] doc: update `CLM:update_T` documentation --- doc/content/setup_tsmp/input_enkfpf.md | 41 ++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index ccdd64c9b..bfec14040 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -472,9 +472,46 @@ CLM (standalone only). `CLM:update_T`: (integer) Flag for updating of CLM's ground and vegetation temperature. -- 0: No update of ground and vegetation temperature +`0` +: No update of ground and vegetation temperature + +`1` +: Update of ground and vegetation temperature, based on + Kustas2009-observation operator + +`2` +: Update of ground and vegetation temperature based on TSKIN in + observation operator + +`3` +: Update of first-layer-soil-temperature and vegetation temperature + based on TSKIN in observation operator + +#### Kustas2009 observation operator + +Variable names in this subsection reflect the conventions from +Kustas2009, not from CLM model. + +Eq(7) from Kustas2009 +() + +{attribution="Kustas2009"} +>\begin{align*} +>T_R(\phi) \approx \left[ f_C(\phi) T_{C}^{4} + (1 - f_C(\phi)) T_{S}^{4} \right]^{\frac{1}{4}} +>\end{align*} +> +> where $T_C$ is canopy temperature, $T_{S}$ is soil temperature, and +> $f_C(\phi)$ is the fractional vegetation cover observed at the +> radiometer view angle $\phi$. For a canopy with a spherical leaf +> angle distribution and leaf area index $LAI$, +> +>\begin{align*} +> f_C(\phi) = 1 - \exp\left(\frac{-0.5 \, \Omega \, LAI}{\cos\phi}\right) +>\end{align*} +> +> where the factor $\Omega$ indicates the degree to which vegetation +> is clumped as in row crops or sparsely vegetated shrubland canopies. -- 1: Update of ground and vegetation temperature ### CLM:print_swc ### From 454d8a615922478eecb97062e0df2647a76b940c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 23 Jan 2025 10:40:37 +0100 Subject: [PATCH 17/21] doc: update `CLM:update_T` doc --- doc/content/setup_tsmp/input_enkfpf.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index bfec14040..0831ac59e 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -472,18 +472,18 @@ CLM (standalone only). `CLM:update_T`: (integer) Flag for updating of CLM's ground and vegetation temperature. -`0` +`CLM:update_T=0` : No update of ground and vegetation temperature -`1` +`CLM:update_T=1` : Update of ground and vegetation temperature, based on Kustas2009-observation operator -`2` +`CLM:update_T=2` : Update of ground and vegetation temperature based on TSKIN in observation operator -`3` +`CLM:update_T=3` : Update of first-layer-soil-temperature and vegetation temperature based on TSKIN in observation operator @@ -492,10 +492,10 @@ vegetation temperature. Variable names in this subsection reflect the conventions from Kustas2009, not from CLM model. -Eq(7) from Kustas2009 -() +LST is compared to the surface temperature $T_R$ defined by +Kustas2009: -{attribution="Kustas2009"} +{attribution="Kustas2009 (), Eq(7)"} >\begin{align*} >T_R(\phi) \approx \left[ f_C(\phi) T_{C}^{4} + (1 - f_C(\phi)) T_{S}^{4} \right]^{\frac{1}{4}} >\end{align*} From 1ce57b074655587341f272d97613a8c1561f345b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 24 Jan 2025 14:50:58 +0100 Subject: [PATCH 18/21] LST-DA: add debug output for first layer `t_soisno` --- .../pdaf/model/clm3_5/enkf_clm_mod.F90 | 67 +++++++++++------- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 68 ++++++++++++------- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 68 ++++++++++++------- 3 files changed, 125 insertions(+), 78 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 index 5a8114de0..58aa5e9f9 100755 --- a/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 +++ b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 @@ -214,6 +214,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 @@ -473,32 +481,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 @@ -534,6 +516,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/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index f9bea2758..bd7d66560 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -469,6 +469,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 @@ -748,32 +756,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 @@ -831,6 +813,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index f9bea2758..bd7d66560 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -469,6 +469,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 @@ -748,32 +756,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 @@ -831,6 +813,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 4ccc5f3076b87e2323b1cec0567801c5a868880f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 28 Jan 2025 20:46:46 +0100 Subject: [PATCH 19/21] LST-DA: all temperature layers updated --- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 36 ++++++++++++------- .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 36 ++++++++++++------- 2 files changed, 48 insertions(+), 24 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index bd7d66560..1235bf7c8 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -87,6 +87,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 @@ -98,6 +99,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev integer :: p integer :: c integer :: g @@ -359,7 +361,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)) @@ -375,13 +377,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 @@ -424,6 +428,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 @@ -442,6 +447,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 @@ -524,8 +530,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 @@ -566,6 +574,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 @@ -602,6 +611,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 @@ -783,8 +793,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index bd7d66560..1235bf7c8 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -87,6 +87,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 @@ -98,6 +99,7 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev integer :: p integer :: c integer :: g @@ -359,7 +361,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)) @@ -375,13 +377,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 @@ -424,6 +428,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 @@ -442,6 +447,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 @@ -524,8 +530,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 @@ -566,6 +574,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 @@ -602,6 +611,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 @@ -783,8 +793,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 fb7e8fd13b725b8f3f47725c5d80172d411984f3 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 29 Jan 2025 09:23:48 +0100 Subject: [PATCH 20/21] bugfix: LST-DA missing `end do` --- bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 1 + bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 | 1 + 2 files changed, 2 insertions(+) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 1235bf7c8..5e89a725f 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -386,6 +386,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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 1235bf7c8..5e89a725f 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -386,6 +386,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 1bca0d54096cb4606c70ba2d83e288b297ce264a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 25 Apr 2025 09:54:59 +0200 Subject: [PATCH 21/21] introduce `clmupdate_T.eq.4` --- .../pdaf/framework/init_dim_obs_f_pdaf.F90 | 2 +- .../pdaf/framework/init_dim_obs_pdaf.F90 | 2 +- bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 | 2 +- .../pdaf/model/clm5_0/enkf_clm_mod_5.F90 | 46 +++++++++++++++++++ .../pdaf/model/eclm/enkf_clm_mod_5.F90 | 46 +++++++++++++++++++ 5 files changed, 95 insertions(+), 3 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index 070482313..98b44597c 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index 0fe391a9e..c17dc28f8 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index abc9f0d57..62ef12e5a 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index a58ecce74..bdb7df336 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/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/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index a58ecce74..bdb7df336 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/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