Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CLM-related index arrays state_clm2pdaf and state_pdaf2clm_* #254

Merged
merged 11 commits into from
Nov 18, 2024
Merged

CLM-related index arrays state_clm2pdaf and state_pdaf2clm_* #254

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
1d07a77
rename `col_index_hydr_act` to `state_clm2pdaf_p`
jjokella Nov 9, 2024
8837136
`cleanup_clm_statevec` routine
jjokella Nov 9, 2024
bba91ac
simplify interface
jjokella Nov 9, 2024
f7b2148
use `state_clm2pdaf_p` throught SWC-update state vector
jjokella Nov 9, 2024
554b949
`state_pdaf2clm_*_p` allocated in `define_clm_statevec`
jjokella Nov 9, 2024
bbc601d
simplify usage of `state_clm2pdaf_p` in `update_clm`
jjokella Nov 9, 2024
2a1c08b
simplify usage of `state_clm2pdaf_p` in `update_clm` II
jjokella Nov 9, 2024
a5a4f39
Merge branch 'master' into dev-state_clm2pdaf
jjokella Nov 14, 2024
a1f69d0
bugfix: no input for `cleanup_clm_statevec`
jjokella Nov 14, 2024
7054b6d
bugfix: add use-statement for `cleanup_clm_statevec`
jjokella Nov 14, 2024
b639ad6
bugfix: remove incorrect de-allocation of `state_pdaf2clm` arrays
jjokella Nov 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
use ColumnType, only : col
! use GetGlobalValuesMod, only: GetGlobalWrite
! use clm_varcon, only: nameg
use enkf_clm_mod, only: col_index_hydr_act
use enkf_clm_mod, only: state_clm2pdaf_p
use enkf_clm_mod, only: clmstatevec_only_active
use enkf_clm_mod, only: clmstatevec_max_layer
#else
Expand Down Expand Up @@ -951,7 +951,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
print *, "clmstatevec_max_layer=", clmstatevec_max_layer
call abort_parallel()
end if
obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i))
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))
Expand Down
4 changes: 2 additions & 2 deletions bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
use ColumnType, only : col
! use GetGlobalValuesMod, only: GetGlobalWrite
! use clm_varcon, only: nameg
use enkf_clm_mod, only: col_index_hydr_act
use enkf_clm_mod, only: state_clm2pdaf_p
use enkf_clm_mod, only: clmstatevec_only_active
use enkf_clm_mod, only: clmstatevec_max_layer
#else
Expand Down Expand Up @@ -944,7 +944,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
print *, "clmstatevec_max_layer=", clmstatevec_max_layer
call abort_parallel()
end if
obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i))
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))
Expand Down
6 changes: 6 additions & 0 deletions bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_5.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,11 +224,17 @@ subroutine clm_finalize() bind(C,name="clm_finalize")

! use ESMF, only : ESMF_Initialize, ESMF_Finalize
use cime_comp_mod, only : cime_final
use enkf_clm_mod, only : cleanup_clm_statevec

implicit none

call cime_final()

#if defined CLMSA
! TSMP-PDAF: Deallocate arrays from `define_clm_statevec`
call cleanup_clm_statevec()
#endif

!--------------------------------------------------------------------------
! Clean-up
!--------------------------------------------------------------------------
Expand Down
215 changes: 120 additions & 95 deletions bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ module enkf_clm_mod
! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model
! LST in LST assimilation (clmupdate_T)
real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI)
integer, allocatable :: col_index_hydr_act(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1)
integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1)
integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc
integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp
integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
Expand Down Expand Up @@ -93,7 +93,10 @@ subroutine define_clm_statevec(mype)
integer,intent(in) :: mype

integer :: i
integer :: j
integer :: jj
integer :: c
integer :: g
integer :: cc
integer :: cccheck

Expand Down Expand Up @@ -123,10 +126,8 @@ subroutine define_clm_statevec(mype)

if(clmstatevec_only_active .eq. 1) then

IF (allocated(col_index_hydr_act)) deallocate(col_index_hydr_act)
allocate(col_index_hydr_act(begc:endc,min(nlevsoi,clmstatevec_max_layer)))

clm_statevecsize = 0
IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p)
allocate(state_clm2pdaf_p(begc:endc,min(nlevsoi,clmstatevec_max_layer)))

cc = 0

Expand All @@ -138,47 +139,119 @@ subroutine define_clm_statevec(mype)
! and layers above bedrock
if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then
cc = cc + 1
col_index_hydr_act(c,i) = cc
state_clm2pdaf_p(c,i) = cc
else
col_index_hydr_act(c,i) = ispval
state_clm2pdaf_p(c,i) = ispval
end if
end if
end do
end do

clm_statevecsize = cc

! ! Check against other method of computation
! cccheck = 0
! do c=clm_begc,clm_endc
! if(col%hydrologically_active(c)) then
! ! Only count non-bedrock layers
! cccheck = cccheck + col%nbedrock(c)
! ! Possible -1 to leave out layer that is partly bedrock
! end if
! end do

! if(clm_statevecsize .ne. cccheck) then
! print *, "clm_statevecsize", clm_statevecsize
! print *, "cccheck", cccheck
! error stop "clm_statevecsize not equal to cccheck"
! end if

! Set `clm_varsize`, even though it is currently not used
! for `clmupdate_swc.eq.1`
clm_varsize = clm_statevecsize
clm_statevecsize = cc

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 i=1,nlevsoi
do c=clm_begc,clm_endc
! Only take into account layers above input maximum layer
if(i<=clmstatevec_max_layer) then
! Only take into account hydrologically active columns
! and layers above bedrock
if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then
cc = cc + 1
state_pdaf2clm_c_p(cc) = c
state_pdaf2clm_j_p(cc) = i
end if
end if
end do
end do
else

IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p)
allocate(state_clm2pdaf_p(begc:endc,nlevsoi))

do i=1,nlevsoi
do c=clm_begc,clm_endc
state_clm2pdaf_p(c,i) = (c - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1)
end do
end do

! #cols values per grid-cell
clm_varsize = (endc-begc+1) * nlevsoi
clm_statevecsize = (endc-begc+1) * nlevsoi

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 i=1,nlevsoi
do c=clm_begc,clm_endc
cc = cc + 1
state_pdaf2clm_c_p(cc) = c
state_pdaf2clm_j_p(cc) = i
end do
end do

end if

else

IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p)
allocate(state_clm2pdaf_p(begc:endc,nlevsoi))

do i=1,nlevsoi
do c=clm_begc,clm_endc
! All columns in a gridcell are assigned the updated
! gridcell-SWC
state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1)
end do
end do

! One value per grid-cell
clm_varsize = (endg-begg+1) * nlevsoi
clm_statevecsize = (endg-begg+1) * nlevsoi

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 i=1,nlevsoi
do j=clm_begg,clm_endg

! SWC from the first column of each gridcell
newgridcell = .true.
do jj=clm_begc,clm_endc
g = col%gridcell(jj)
if (g .eq. j) then
if (newgridcell) then
newgridcell = .false.
! Possibliy: Add state_pdaf2clm_g_p
state_pdaf2clm_c_p(cc) = jj
state_pdaf2clm_j_p(cc) = i
end if
end if
end do

cc = cc + 1
end do
end do



end if
endif

Expand Down Expand Up @@ -207,13 +280,9 @@ subroutine define_clm_statevec(mype)

!write(*,*) 'clm_statevecsize is ',clm_statevecsize
IF (allocated(clm_statevec)) deallocate(clm_statevec)
IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p)
IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p)
if ((clmupdate_swc.ne.0) .or. (clmupdate_T.ne.0) .or. (clmupdate_texture.ne.0)) then
!hcp added condition
allocate(clm_statevec(clm_statevecsize))
allocate(state_pdaf2clm_c_p(clm_statevecsize))
allocate(state_pdaf2clm_j_p(clm_statevecsize))
end if

!write(*,*) 'clm_paramsize is ',clm_paramsize
Expand All @@ -224,6 +293,18 @@ subroutine define_clm_statevec(mype)

end subroutine define_clm_statevec

subroutine cleanup_clm_statevec()

implicit none

! Deallocate arrays from `define_clm_statevec`
IF (allocated(clm_statevec)) deallocate(clm_statevec)
IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p)
IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p)
IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p)

end subroutine cleanup_clm_statevec

subroutine set_clm_statevec(tstartcycle, mype)
use clm_instMod, only : soilstate_inst, waterstate_inst
use clm_varpar , only : nlevsoi
Expand Down Expand Up @@ -263,52 +344,10 @@ subroutine set_clm_statevec(tstartcycle, mype)
endif

if(clmupdate_swc.ne.0) then
! write swc values to state vector
cc = 1
do i=1,nlevsoi

if(clmstatevec_allcol.eq.1) then

do jj=clm_begc,clm_endc
! SWC from all columns of each gridcell
if(clmstatevec_only_active.eq.1) then
if(i<=clmstatevec_max_layer .and. col%hydrologically_active(jj) .and. i<=col%nbedrock(jj) ) then
clm_statevec(cc+offset) = swc(jj,i)
state_pdaf2clm_c_p(cc+offset) = jj
state_pdaf2clm_j_p(cc+offset) = i
cc = cc + 1
end if
else
clm_statevec(cc+offset) = swc(jj,i)
state_pdaf2clm_c_p(cc+offset) = jj
state_pdaf2clm_j_p(cc+offset) = i
cc = cc + 1
end if

end do

else

do j=clm_begg,clm_endg
! SWC from the first column of each gridcell
newgridcell = .true.
do jj=clm_begc,clm_endc
g = col%gridcell(jj)
if (g .eq. j) then
if (newgridcell) then
newgridcell = .false.
clm_statevec(cc+offset) = swc(jj,i)
state_pdaf2clm_c_p(cc+offset) = jj
state_pdaf2clm_j_p(cc+offset) = i
endif
endif
end do
cc = cc + 1
end do

end if

end do
! write swc values to state vector
do cc = 1, clm_statevecsize
clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc))
end do
endif

!hcp LAI
Expand Down Expand Up @@ -378,6 +417,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
real(r8) :: rliq,rice
real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm)
real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm)
real(r8) :: swc_update ! updated SWC in loop

integer :: i,j,jj,g,cc=0,offset=0
character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu
Expand Down Expand Up @@ -466,23 +506,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
! do j=clm_begg,clm_endg
do j=clm_begc,clm_endc

! Set cc (the state vector index) from the
! CLM5-grid-index and the `CLM5-layer-index times
! num_gridcells`
if(clmstatevec_allcol.eq.1) then
if(clmstatevec_only_active.eq.1) then
if(i<=clmstatevec_max_layer .and. col%hydrologically_active(j) .and. i<=col%nbedrock(j) ) then
cc = col_index_hydr_act(j,i)
else
cycle
end if
else
cc = (j - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1)
end if
else
cc = (col%gridcell(j) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1)
end if

if(swc(j,i).eq.0.0) then
swc_zero_before_update = .true.

Expand All @@ -499,12 +522,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
!h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice)
end if

if(clm_statevec(cc+offset).le.watmin_check) then
swc_update = clm_statevec(state_clm2pdaf_p(j,i))

if(swc_update.le.watmin_check) then
swc(j,i) = watmin_set
else if(clm_statevec(cc+offset).ge.watsat(j,i)) then
else if(swc_update.ge.watsat(j,i)) then
swc(j,i) = watsat(j,i)
else
swc(j,i) = clm_statevec(cc+offset)
swc(j,i) = swc_update
endif

if (isnan(swc(j,i))) then
Expand Down
6 changes: 6 additions & 0 deletions bldsva/intf_DA/pdaf/model/eclm/enkf_clm_5.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,11 +224,17 @@ subroutine clm_finalize() bind(C,name="clm_finalize")

! use ESMF, only : ESMF_Initialize, ESMF_Finalize
use cime_comp_mod, only : cime_final
use enkf_clm_mod, only : cleanup_clm_statevec

implicit none

call cime_final()

#if defined CLMSA
! TSMP-PDAF: Deallocate arrays from `define_clm_statevec`
call cleanup_clm_statevec()
#endif

!--------------------------------------------------------------------------
! Clean-up
!--------------------------------------------------------------------------
Expand Down
Loading
Loading