Skip to content

Commit

Permalink
Merge pull request #1385 from LIS-navari/fix/SMAP_DA_soilTemp_bug_fix
Browse files Browse the repository at this point in the history 
Fix/smap da soil temp bug fix
  • Loading branch information
@jvgeiger
jvgeiger authored May 7, 2024
2 parents f446cb8 + 4bee2d5 commit 0be682f
Show file tree
Hide file tree
Showing 4 changed files with 133 additions and 74 deletions.
23 changes: 16 additions & 7 deletions lis/surfacemodels/land/jules.5.0/da_soilm/jules50_setsoilm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,10 @@
! apply the update for that members and set the other member to the mean
! value of the ensemble (i.e. mean of the members that met the conditions)
! 3- If less then 50% of the ensemble members met the update conditions -->
! adjust the states
! adjust the states
! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix
! (add check for frozen soil)

! !INTERFACE:
subroutine jules50_setsoilm(n, LSM_State)
! !USES:
Expand All @@ -35,7 +38,9 @@ subroutine jules50_setsoilm(n, LSM_State)
!MN: added to use LIS_sfmodel_struc
!use LIS_surfaceModelDataMod, only : LIS_sfmodel_struc
use jules_soil_mod, only: dzsoil !
use jules_surface_mod, only: l_aggregate
use jules_surface_mod, only: l_aggregate
use LIS_constantsMod, only : LIS_CONST_TKFRZ

implicit none
! !ARGUMENTS:
integer, intent(in) :: n
Expand Down Expand Up @@ -159,17 +164,19 @@ subroutine jules50_setsoilm(n, LSM_State)
else
! frozen soil moisture content has been added (Yonghwan Kwon)
if (p_s_sth(1) * sat_p(1) + delta1.gt.MIN_THRESHOLD(1) .and.&
p_s_sth(1) * sat_p(1) + delta1.lt.sm_threshold(1)) then
p_s_sth(1) * sat_p(1) + delta1.lt.sm_threshold(1) .and.&
jules50_struc(n)%jules50(t)%t_soil(1) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
! MN save the flag for each tile (col*row*ens) PILDAS -->(64*44)*20
! MN save the flag for each tile (col*row*ens)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
update_flag(gid) = update_flag(gid).and.(.false.)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif

if (p_s_sth(2) * sat_p(2) + delta2.gt.MIN_THRESHOLD(2) .and.&
p_s_sth(2) * sat_p(2) + delta2.lt.sm_threshold(2)) then
p_s_sth(2) * sat_p(2) + delta2.lt.sm_threshold(2) .and.&
jules50_struc(n)%jules50(t)%t_soil(2) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
Expand All @@ -178,7 +185,8 @@ subroutine jules50_setsoilm(n, LSM_State)
endif

if (p_s_sth(3) * sat_p(3) + delta3.gt.MIN_THRESHOLD(3) .and.&
p_s_sth(3) * sat_p(3) + delta3.lt.sm_threshold(3)) then
p_s_sth(3) * sat_p(3) + delta3.lt.sm_threshold(3) .and.&
jules50_struc(n)%jules50(t)%t_soil(3) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
Expand All @@ -187,7 +195,8 @@ subroutine jules50_setsoilm(n, LSM_State)
endif

if (p_s_sth(4) * sat_p(4) + delta4.gt.MIN_THRESHOLD(4) .and.&
p_s_sth(4) * sat_p(4) + delta4.lt.sm_threshold(4)) then
p_s_sth(4) * sat_p(4) + delta4.lt.sm_threshold(4) .and.&
jules50_struc(n)%jules50(t)%t_soil(4) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
Expand Down
15 changes: 11 additions & 4 deletions lis/surfacemodels/land/noah.3.6/da_soilm/noah36_setsoilm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@
! 27Feb2005: Sujay Kumar; Initial Specification
! 25Jun2006: Sujay Kumar: Updated for the ESMF design
! 28Aug2017: Mahdi Navari; Updated to take into account the latest developments in the SM DA
! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix
! (add check for frozen soil)
!
! !INTERFACE:
subroutine noah36_setsoilm(n, LSM_State)
Expand All @@ -23,6 +25,7 @@ subroutine noah36_setsoilm(n, LSM_State)
use LIS_coreMod
use LIS_logMod
use noah36_lsmMod
use LIS_constantsMod, only : LIS_CONST_TKFRZ

implicit none
! !ARGUMENTS:
Expand Down Expand Up @@ -144,7 +147,8 @@ subroutine noah36_setsoilm(n, LSM_State)
! MN: check MIN_THRESHOLD < volumetric liquid soil moisture < threshold
if(noah36_struc(n)%noah(t)%sh2o(1)+delta1.gt.MIN_THRESHOLD .and.&
noah36_struc(n)%noah(t)%sh2o(1)+delta1.lt.&
sm_threshold) then
sm_threshold .and.&
noah36_struc(n)%noah(t)%stc(1) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
! MN save the flag for each tile (col*row*ens) (64*44)*20
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
Expand All @@ -153,23 +157,26 @@ subroutine noah36_setsoilm(n, LSM_State)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah36_struc(n)%noah(t)%sh2o(2)+delta2.gt.MIN_THRESHOLD .and.&
noah36_struc(n)%noah(t)%sh2o(2)+delta2.lt.sm_threshold) then
noah36_struc(n)%noah(t)%sh2o(2)+delta2.lt.sm_threshold .and.&
noah36_struc(n)%noah(t)%stc(2) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
update_flag(gid) = update_flag(gid).and.(.false.)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah36_struc(n)%noah(t)%sh2o(3)+delta3.gt.MIN_THRESHOLD .and.&
noah36_struc(n)%noah(t)%sh2o(3)+delta3.lt.sm_threshold) then
noah36_struc(n)%noah(t)%sh2o(3)+delta3.lt.sm_threshold .and.&
noah36_struc(n)%noah(t)%stc(3) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
update_flag(gid) = update_flag(gid).and.(.false.)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah36_struc(n)%noah(t)%sh2o(4)+delta4.gt.MIN_THRESHOLD .and.&
noah36_struc(n)%noah(t)%sh2o(4)+delta4.lt.sm_threshold) then
noah36_struc(n)%noah(t)%sh2o(4)+delta4.lt.sm_threshold .and.&
noah36_struc(n)%noah(t)%stc(4) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
Expand Down
69 changes: 50 additions & 19 deletions lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,20 @@
!
! !REVISION HISTORY:
! 21Oct2018: Mahdi Navari; Sujay Kumar ; Initial Specification
!
! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix
! (add check for frozen soil)
! 10 Apr 2024: Mahdi Navari; ensemble flag bug fix
! 11 Apr 2024: Mahdi Navari; Fix bug related to computing delta. The bug fix sets
! the code to use total soil moisture (smc) if the ground
! is frozen, rather than the liquid part of soil moisture (sh2o).
! !INTERFACE:
subroutine noah39_setsoilm(n, LSM_State)
! !USES:
use ESMF
use LIS_coreMod
use LIS_logMod
use noah39_lsmMod
use LIS_constantsMod, only : LIS_CONST_TKFRZ

implicit none
! !ARGUMENTS:
Expand All @@ -47,14 +53,13 @@ subroutine noah39_setsoilm(n, LSM_State)
real, pointer :: soilm4(:)
integer :: t, j,i, gid, m, t_unpert, LIS_localP, row , col
integer :: status
real :: delta(4)
real :: delta(4), ens_count(4)
real :: delta1,delta2,delta3,delta4
real :: tmpval
logical :: bounds_violation
integer :: nIter
logical :: update_flag(LIS_rc%ngrid(n))
logical :: ens_flag(LIS_rc%nensem(n))
! mn
logical :: ens_flag(4,LIS_rc%nensem(n))
real :: tmp(LIS_rc%nensem(n)), tmp0(LIS_rc%nensem(n))
real :: tmp1(LIS_rc%nensem(n)),tmp2(LIS_rc%nensem(n)),tmp3(LIS_rc%nensem(n)),tmp4(LIS_rc%nensem(n))
logical :: update_flag_tile(LIS_rc%npatch(n,LIS_rc%lsm_index))
Expand Down Expand Up @@ -142,7 +147,8 @@ subroutine noah39_setsoilm(n, LSM_State)
! MN: check MIN_THRESHOLD < volumetric liquid soil moisture < threshold
if(noah39_struc(n)%noah(t)%sh2o(1)+delta1.gt.MIN_THRESHOLD .and.&
noah39_struc(n)%noah(t)%sh2o(1)+delta1.lt.&
sm_threshold) then
sm_threshold .and.&
noah39_struc(n)%noah(t)%stc(1) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
! MN save the flag for each tile (col*row*ens) (64*44)*20
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
Expand All @@ -151,23 +157,26 @@ subroutine noah39_setsoilm(n, LSM_State)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah39_struc(n)%noah(t)%sh2o(2)+delta2.gt.MIN_THRESHOLD .and.&
noah39_struc(n)%noah(t)%sh2o(2)+delta2.lt.sm_threshold) then
noah39_struc(n)%noah(t)%sh2o(2)+delta2.lt.sm_threshold .and.&
noah39_struc(n)%noah(t)%stc(2) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
update_flag(gid) = update_flag(gid).and.(.false.)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah39_struc(n)%noah(t)%sh2o(3)+delta3.gt.MIN_THRESHOLD .and.&
noah39_struc(n)%noah(t)%sh2o(3)+delta3.lt.sm_threshold) then
noah39_struc(n)%noah(t)%sh2o(3)+delta3.lt.sm_threshold .and.&
noah39_struc(n)%noah(t)%stc(3) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
update_flag(gid) = update_flag(gid).and.(.false.)
update_flag_tile(t) = update_flag_tile(t).and.(.false.)
endif
if(noah39_struc(n)%noah(t)%sh2o(4)+delta4.gt.MIN_THRESHOLD .and.&
noah39_struc(n)%noah(t)%sh2o(4)+delta4.lt.sm_threshold) then
noah39_struc(n)%noah(t)%sh2o(4)+delta4.lt.sm_threshold .and.&
noah39_struc(n)%noah(t)%stc(4) .gt. LIS_CONST_TKFRZ) then
update_flag(gid) = update_flag(gid).and.(.true.)
update_flag_tile(t) = update_flag_tile(t).and.(.true.)
else
Expand Down Expand Up @@ -345,8 +354,8 @@ subroutine noah39_setsoilm(n, LSM_State)
noah39_struc(n)%noah(t)%smc(4) = smc_tmp


!MN 4 print
tmp0(m) = noah39_struc(n)%noah(t)%smc(1)
!MN 4 print
tmp0(m) = noah39_struc(n)%noah(t)%smc(1)

endif ! flag for each tile
enddo ! loop over tile
Expand All @@ -365,23 +374,30 @@ subroutine noah39_setsoilm(n, LSM_State)
do while(bounds_violation)
niter = niter + 1
!t_unpert = i*LIS_rc%nensem(n)
t_unpert = i+LIS_rc%nensem(n)-1
t_unpert = i+LIS_rc%nensem(n)-1
do j=1,4
delta(j) = 0.0
do m=1,LIS_rc%nensem(n)-1
t = i+m-1
!t = (i-1)*LIS_rc%nensem(n)+m

if(m.ne.LIS_rc%nensem(n)) then
if(m.ne.LIS_rc%nensem(n)) then
if (noah39_struc(n)%noah(t)%stc(1) .gt. LIS_CONST_TKFRZ) then ! MN
delta(j) = delta(j) + &
(noah39_struc(n)%noah(t)%sh2o(j) - &
noah39_struc(n)%noah(t_unpert)%sh2o(j))
else ! MN
delta(j) = delta(j) + &
(noah39_struc(n)%noah(t)%smc(j) - &
noah39_struc(n)%noah(t_unpert)%smc(j))
endif
endif

enddo
enddo

do j=1,4
ens_count(j) = 0.0
delta(j) =delta(j)/(LIS_rc%nensem(n)-1)
do m=1,LIS_rc%nensem(n)-1
t = i+m-1
Expand All @@ -398,15 +414,17 @@ subroutine noah39_setsoilm(n, LSM_State)
noah39_struc(n)%noah(t)%smc(j) = &
max(noah39_struc(n)%noah(t_unpert)%smc(j),&
MIN_THRESHOLD)
ens_flag(m) = .false.
ens_flag(j,m) = .false.
ens_count(j) = ens_count(j) + 1
elseif(tmpval.ge.sm_threshold) then
noah39_struc(n)%noah(t)%sh2o(j) = &
min(noah39_struc(n)%noah(t_unpert)%sh2o(j),&
sm_threshold)
noah39_struc(n)%noah(t)%smc(j) = &
min(noah39_struc(n)%noah(t_unpert)%smc(j),&
sm_threshold)
ens_flag(m) = .false.
ens_flag(j,m) = .false.
ens_count(j) = ens_count(j) + 1
endif
enddo
enddo
Expand All @@ -420,20 +438,33 @@ subroutine noah39_setsoilm(n, LSM_State)
t = i+m-1
!t = (i-1)*LIS_rc%nensem(n)+m
if(m.ne.LIS_rc%nensem(n)) then
delta(j) = delta(j) + &
(noah39_struc(n)%noah(t)%sh2o(j) - &
noah39_struc(n)%noah(t_unpert)%sh2o(j))
if (noah39_struc(n)%noah(t)%stc(1) .gt. LIS_CONST_TKFRZ) then ! MN
delta(j) = delta(j) + &
(noah39_struc(n)%noah(t)%sh2o(j) - &
noah39_struc(n)%noah(t_unpert)%sh2o(j))
else ! MN
delta(j) = delta(j) + &
(noah39_struc(n)%noah(t)%smc(j) - &
noah39_struc(n)%noah(t_unpert)%smc(j))
endif
endif
enddo
enddo

do j=1,4
delta(j) =delta(j)/(LIS_rc%nensem(n)-1)
if (ens_count(j).lt.(LIS_rc%nensem(n)-1)) then
!delta(j) =delta(j)/(LIS_rc%nensem(n)-1)
delta(j) =delta(j)/(LIS_rc%nensem(n)-1-ens_count(j))
else
delta(j) =delta(j)/(LIS_rc%nensem(n)-1) !We apply this
!to remove the residual bias after setting all ensemble members
!to the bond values (MIN_THRESHOLD, sm_threshold, or unperturbed)
endif
do m=1,LIS_rc%nensem(n)-1
t = i+m-1
!t = (i-1)*LIS_rc%nensem(n)+m

if(ens_flag(m)) then
if(ens_flag(j,m)) then
tmpval = noah39_struc(n)%noah(t)%sh2o(j) - &
delta(j)
MAX_THRESHOLD = noah39_struc(n)%noah(t)%smcmax
Expand Down
Loading

0 comments on commit 0be682f

Please sign in to comment.