From 4bda49dee2a4e4c23f29d587f0cfac7efce030b7 Mon Sep 17 00:00:00 2001 From: Mahdi Navari Date: Tue, 18 Jul 2023 23:15:21 -0400 Subject: [PATCH 1/3] SMAP DA soil temeprature bug fix --- .../jules.5.0/da_soilm/jules50_setsoilm.F90 | 23 +++++++++++++------ .../noah.3.6/da_soilm/noah36_setsoilm.F90 | 15 ++++++++---- .../noah.3.9/da_soilm/noah39_setsoilm.F90 | 15 ++++++++---- .../da_soilm/noahmp401_setsoilm.F90 | 17 +++++++++----- 4 files changed, 49 insertions(+), 21 deletions(-) diff --git a/lis/surfacemodels/land/jules.5.0/da_soilm/jules50_setsoilm.F90 b/lis/surfacemodels/land/jules.5.0/da_soilm/jules50_setsoilm.F90 index dfbceb11a..60fadc471 100644 --- a/lis/surfacemodels/land/jules.5.0/da_soilm/jules50_setsoilm.F90 +++ b/lis/surfacemodels/land/jules.5.0/da_soilm/jules50_setsoilm.F90 @@ -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: @@ -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 @@ -159,9 +164,10 @@ 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.) @@ -169,7 +175,8 @@ subroutine jules50_setsoilm(n, LSM_State) 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 @@ -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 @@ -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 diff --git a/lis/surfacemodels/land/noah.3.6/da_soilm/noah36_setsoilm.F90 b/lis/surfacemodels/land/noah.3.6/da_soilm/noah36_setsoilm.F90 index 6101fa196..b487dad67 100644 --- a/lis/surfacemodels/land/noah.3.6/da_soilm/noah36_setsoilm.F90 +++ b/lis/surfacemodels/land/noah.3.6/da_soilm/noah36_setsoilm.F90 @@ -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) @@ -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: @@ -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.) @@ -153,7 +157,8 @@ 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 @@ -161,7 +166,8 @@ subroutine noah36_setsoilm(n, LSM_State) 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 @@ -169,7 +175,8 @@ subroutine noah36_setsoilm(n, LSM_State) 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 diff --git a/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 b/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 index 090b1844f..39781d042 100644 --- a/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 +++ b/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 @@ -13,6 +13,8 @@ ! ! !REVISION HISTORY: ! 21Oct2018: Mahdi Navari; Sujay Kumar ; Initial Specification +! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix +! (add check for frozen soil) ! ! !INTERFACE: subroutine noah39_setsoilm(n, LSM_State) @@ -21,6 +23,7 @@ subroutine noah39_setsoilm(n, LSM_State) use LIS_coreMod use LIS_logMod use noah39_lsmMod + use LIS_constantsMod, only : LIS_CONST_TKFRZ implicit none ! !ARGUMENTS: @@ -142,7 +145,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.) @@ -151,7 +155,8 @@ 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 @@ -159,7 +164,8 @@ subroutine noah39_setsoilm(n, LSM_State) 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 @@ -167,7 +173,8 @@ subroutine noah39_setsoilm(n, LSM_State) 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 diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 index 045bde466..74d117996 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 @@ -27,7 +27,8 @@ ! 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 - +! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix +! (add check for frozen soil) ! !INTERFACE: subroutine NoahMP401_setsoilm(n, LSM_State) @@ -39,7 +40,7 @@ subroutine NoahMP401_setsoilm(n, LSM_State) !use module_sf_noahlsm_36 !MN !use module_sf_noahmpdrv_401, only: parameters use NOAHMP_TABLES_401, ONLY : SMCMAX_TABLE,SMCWLT_TABLE - + use LIS_constantsMod, only : LIS_CONST_TKFRZ implicit none ! !ARGUMENTS: @@ -142,7 +143,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) ! MN: check MIN_THRESHOLD < volumetric liquid soil moisture < threshold if(noahmp401_struc(n)%noahmp401(t)%sh2o(1)+delta1.gt.MIN_THRESHOLD .and.& noahmp401_struc(n)%noahmp401(t)%sh2o(1)+delta1.lt.& - sm_threshold) then + sm_threshold .and.& + noahmp401_struc(n)%noahmp401(t)%tslb(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.) @@ -151,7 +153,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) update_flag_tile(t) = update_flag_tile(t).and.(.false.) endif if(noahmp401_struc(n)%noahmp401(t)%sh2o(2)+delta2.gt.MIN_THRESHOLD .and.& - noahmp401_struc(n)%noahmp401(t)%sh2o(2)+delta2.lt.sm_threshold) then + noahmp401_struc(n)%noahmp401(t)%sh2o(2)+delta2.lt.sm_threshold .and.& + noahmp401_struc(n)%noahmp401(t)%tslb(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 @@ -159,7 +162,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) update_flag_tile(t) = update_flag_tile(t).and.(.false.) endif if(noahmp401_struc(n)%noahmp401(t)%sh2o(3)+delta3.gt.MIN_THRESHOLD .and.& - noahmp401_struc(n)%noahmp401(t)%sh2o(3)+delta3.lt.sm_threshold) then + noahmp401_struc(n)%noahmp401(t)%sh2o(3)+delta3.lt.sm_threshold .and.& + noahmp401_struc(n)%noahmp401(t)%tslb(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 @@ -167,7 +171,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) update_flag_tile(t) = update_flag_tile(t).and.(.false.) endif if(noahmp401_struc(n)%noahmp401(t)%sh2o(4)+delta4.gt.MIN_THRESHOLD .and.& - noahmp401_struc(n)%noahmp401(t)%sh2o(4)+delta4.lt.sm_threshold) then + noahmp401_struc(n)%noahmp401(t)%sh2o(4)+delta4.lt.sm_threshold .and.& + noahmp401_struc(n)%noahmp401(t)%tslb(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 From 17224aec65d6645291b330504b5ad54a5aff56b6 Mon Sep 17 00:00:00 2001 From: Mahdi Navari Date: Wed, 17 Apr 2024 14:36:56 -0400 Subject: [PATCH 2/3] Three bug fixes : 1- Soil temperature bias bug fix (added check for frozen soil) 2- ensemble flag bug fix 3- 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). --- .../noah.3.9/da_soilm/noah39_setsoilm.F90 | 54 +++++++---- .../da_soilm/noahmp401_setsoilm.F90 | 90 +++++++++++-------- 2 files changed, 91 insertions(+), 53 deletions(-) diff --git a/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 b/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 index 39781d042..7c04e0ee3 100644 --- a/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 +++ b/lis/surfacemodels/land/noah.3.9/da_soilm/noah39_setsoilm.F90 @@ -15,7 +15,10 @@ ! 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: @@ -50,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)) @@ -352,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 @@ -372,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 @@ -405,7 +414,8 @@ 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),& @@ -413,7 +423,8 @@ subroutine noah39_setsoilm(n, LSM_State) 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 @@ -427,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 diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 index 74d117996..2f0007002 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 @@ -29,7 +29,11 @@ ! adjust the states ! 8 May 2023: Mahdi Navari; Soil temperature bias bug fix ! (add check for frozen soil) - +! 11 Mar 2024: Mahdi Navari; ensemble flag bug fix +! 10 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 NoahMP401_setsoilm(n, LSM_State) ! !USES: @@ -37,7 +41,7 @@ subroutine NoahMP401_setsoilm(n, LSM_State) use LIS_coreMod use LIS_logMod use NoahMP401_lsmMod - !use module_sf_noahlsm_36 !MN + !use module_sf_noahlsm_36 !use module_sf_noahmpdrv_401, only: parameters use NOAHMP_TABLES_401, ONLY : SMCMAX_TABLE,SMCWLT_TABLE use LIS_constantsMod, only : LIS_CONST_TKFRZ @@ -66,14 +70,13 @@ subroutine NoahMP401_setsoilm(n, LSM_State) real, pointer :: soilm4(:) integer :: t, j,i, gid, m, t_unpert 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)) integer :: SOILTYP ! soil type index [-] real :: SMCMAX , SMCWLT real :: tmp(LIS_rc%nensem(n)), tmp0(LIS_rc%nensem(n)) @@ -90,6 +93,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) real :: smc_rnd, smc_tmp real :: sh2o_tmp, sh2o_rnd INTEGER, DIMENSION (1) :: seed + integer :: row, col + real :: lat, lon call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) call LIS_verify(status,& @@ -199,18 +204,6 @@ subroutine NoahMP401_setsoilm(n, LSM_State) update_flag_new(gid)= update_flag(gid).or.update_flag_ens(gid) ! new flag enddo - ! MN print -#if 0 -if(i.eq.66) then !i=66 ! --> domain's center 1376 - if(LIS_rc%hr.eq.12) then - write(2001,'(I4, 2x, 3(I2,x), 2x, 23(L1,2x))'),& - i, LIS_rc%mo, LIS_rc%da, LIS_rc%hr,update_flag_tile& - ((i-1)*LIS_rc%nensem(n)+1:(i)*LIS_rc%nensem(n)),& - update_flag_ens(i), update_flag_new(i), update_flag(i) - endif !mn - endif -#endif - ! update step ! loop over grid points do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) @@ -218,7 +211,12 @@ subroutine NoahMP401_setsoilm(n, LSM_State) gid = LIS_domain(n)%gindex(& LIS_surface(n,LIS_rc%lsm_index)%tile(i)%col,& LIS_surface(n,LIS_rc%lsm_index)%tile(i)%row) - + + row = LIS_surface(n, LIS_rc%lsm_index)%tile(i)%row + col = LIS_surface(n, LIS_rc%lsm_index)%tile(i)%col + lat = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lat + lon = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lon + !if(update_flag(gid)) then if(update_flag_new(gid)) then !----------------------------------------------------------------------------------------- @@ -328,8 +326,7 @@ subroutine NoahMP401_setsoilm(n, LSM_State) stop endif endif - - + !----------------------------------------------------------------------------------------- ! randomly resample the smc from [MIN_THRESHOLD, Max value from DA @ that tiem step] !----------------------------------------------------------------------------------------- @@ -356,8 +353,7 @@ subroutine NoahMP401_setsoilm(n, LSM_State) smc_tmp = (MaxEnsSM4 - MinEnsSM4)/2 + MinEnsSM4 noahmp401_struc(n)%noahmp401(t)%sh2o(4) = smc_tmp noahmp401_struc(n)%noahmp401(t)%smc(4) = smc_tmp - - + endif ! flag for each tile enddo ! loop over tile @@ -384,15 +380,22 @@ subroutine NoahMP401_setsoilm(n, LSM_State) !t = (i-1)*LIS_rc%nensem(n)+m if(m.ne.LIS_rc%nensem(n)) then + if (noahmp401_struc(n)%noahmp401(t)%tslb(1) .gt. LIS_CONST_TKFRZ) then ! MN delta(j) = delta(j) + & (noahmp401_struc(n)%noahmp401(t)%sh2o(j) - & noahmp401_struc(n)%noahmp401(t_unpert)%sh2o(j)) + else ! MN + delta(j) = delta(j) + & + (noahmp401_struc(n)%noahmp401(t)%smc(j) - & + noahmp401_struc(n)%noahmp401(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 @@ -410,7 +413,8 @@ subroutine NoahMP401_setsoilm(n, LSM_State) noahmp401_struc(n)%noahmp401(t)%smc(j) = & max(noahmp401_struc(n)%noahmp401(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 noahmp401_struc(n)%noahmp401(t)%sh2o(j) = & min(noahmp401_struc(n)%noahmp401(t_unpert)%sh2o(j),& @@ -418,39 +422,53 @@ subroutine NoahMP401_setsoilm(n, LSM_State) noahmp401_struc(n)%noahmp401(t)%smc(j) = & min(noahmp401_struc(n)%noahmp401(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 - + !-------------------------------------------------------------------------- ! Recalculate the deltas and adjust the ensemble !-------------------------------------------------------------------------- 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 - delta(j) = delta(j) + & - (noahmp401_struc(n)%noahmp401(t)%sh2o(j) - & - noahmp401_struc(n)%noahmp401(t_unpert)%sh2o(j)) + + if(m.ne.LIS_rc%nensem(n)) then + if (noahmp401_struc(n)%noahmp401(t)%tslb(1) .gt. LIS_CONST_TKFRZ) then ! MN + delta(j) = delta(j) + & + (noahmp401_struc(n)%noahmp401(t)%sh2o(j) - & + noahmp401_struc(n)%noahmp401(t_unpert)%sh2o(j)) + else ! MN + delta(j) = delta(j) + & + (noahmp401_struc(n)%noahmp401(t)%smc(j) - & + noahmp401_struc(n)%noahmp401(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-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 = noahmp401_struc(n)%noahmp401(t)%sh2o(j) - & delta(j) SOILTYP = noahmp401_struc(n)%noahmp401(t)%soiltype MAX_THRESHOLD = SMCMAX_TABLE(SOILTYP) - if(.not.(tmpval.le.0.0 .or.& tmpval.gt.(MAX_THRESHOLD))) then @@ -474,7 +492,7 @@ subroutine NoahMP401_setsoilm(n, LSM_State) endif enddo enddo - + if(nIter.gt.10.and.bounds_violation) then !-------------------------------------------------------------------------- ! All else fails, set to the bounds @@ -502,17 +520,13 @@ subroutine NoahMP401_setsoilm(n, LSM_State) noahmp401_struc(n)%noahmp401(t)%sh2o(j) = MIN_THRESHOLD noahmp401_struc(n)%noahmp401(t)%smc(j) = MIN_THRESHOLD endif -! print*, i, m -! print*, 'smc',t, noahmp401_struc(n)%noahmp401(t)%smc(:) -! print*, 'sh2o ',t,noahmp401_struc(n)%noahmp401(t)%sh2o(:) -! print*, 'max ',t,MAX_THRESHOLD !noahmp401_struc(n)%noahmp401(t)%smcmax enddo ! call LIS_endrun() enddo endif end do endif - endif + endif enddo From 4bee2d54a691a761676b14a93da1c87e16f68a93 Mon Sep 17 00:00:00 2001 From: "James V. Geiger" Date: Tue, 7 May 2024 13:35:13 -0400 Subject: [PATCH 3/3] Remove unused variables --- .../land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 index 126607bbd..c30bad17c 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_soilm/noahmp401_setsoilm.F90 @@ -93,8 +93,6 @@ subroutine NoahMP401_setsoilm(n, LSM_State) real :: smc_rnd, smc_tmp real :: sh2o_tmp, sh2o_rnd INTEGER, DIMENSION (1) :: seed - integer :: row, col - real :: lat, lon call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) call LIS_verify(status,& @@ -212,11 +210,6 @@ subroutine NoahMP401_setsoilm(n, LSM_State) LIS_surface(n,LIS_rc%lsm_index)%tile(i)%col,& LIS_surface(n,LIS_rc%lsm_index)%tile(i)%row) - row = LIS_surface(n, LIS_rc%lsm_index)%tile(i)%row - col = LIS_surface(n, LIS_rc%lsm_index)%tile(i)%col - lat = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lat - lon = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lon - !if(update_flag(gid)) then if(update_flag_new(gid)) then !-----------------------------------------------------------------------------------------