From e0e41bfc2e62ec121673fb3c8e3aff2de7b5f914 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 29 Nov 2021 13:59:51 -0700 Subject: [PATCH 01/10] Change renormalization of landunit areas in mksurfdata_map to be consistent with raw datasets (Issue #1555) --- tools/mksurfdata_map/src/mkpftUtilsMod.F90 | 38 +++++- tools/mksurfdata_map/src/mksurfdat.F90 | 127 +++++++++------------ 2 files changed, 88 insertions(+), 77 deletions(-) diff --git a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 index 4a9ea12f97..53b8dd04a7 100644 --- a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 +++ b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 @@ -23,8 +23,9 @@ module mkpftUtilsMod ! ! !PUBLIC MEMBER FUNCTIONS: ! - public :: convert_from_p2g ! Convert a p2g array into pct_pft_type objects - public :: adjust_total_veg_area ! Adjust the total vegetated area (natural veg & crop) to a new specified total + public :: convert_from_p2g ! Convert a p2g array into pct_pft_type objects + public :: adjust_total_veg_area ! Adjust the total vegetated area (natural veg & crop) to a new specified total + public :: adjust_total_natveg_area ! Adjust the total natural vegetated area to a new specified total ! ! !PRIVATE MEMBER FUNCTIONS: @@ -251,6 +252,39 @@ subroutine adjust_total_veg_area(new_total_pct, pctnatpft, pctcft) end subroutine adjust_total_veg_area + !----------------------------------------------------------------------- + subroutine adjust_total_natveg_area(new_total_pct, pctnatpft) + ! + ! !DESCRIPTION: + ! Adjust the total natural vegetated area on the grid cell to a new + ! specified total. + ! + ! If the old areas are 0%, then all the new area goes into pctnatpft. + ! + ! !USES: + use mkpctPftTypeMod, only : pct_pft_type + ! + ! !ARGUMENTS: + real(r8), intent(in) :: new_total_pct ! new total % of natural veg landunit + class(pct_pft_type), intent(inout) :: pctnatpft ! natural veg cover information + ! + ! !LOCAL VARIABLES: + real(r8) :: natpft_l2g ! grid cell % cover of natural veg + real(r8) :: old_total ! old total % cover of natural veg + + character(len=*), parameter :: subname = 'adjust_total_natveg_area' + !----------------------------------------------------------------------- + + natpft_l2g = pctnatpft%get_pct_l2g() + old_total = natpft_l2g + if (old_total > 0._r8) then + call pctnatpft%set_pct_l2g(natpft_l2g * new_total_pct / old_total) + else + call pctnatpft%set_pct_l2g(new_total_pct) + end if + + end subroutine adjust_total_natveg_area + end module mkpftUtilsMod diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 9051c57707..706569fde0 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -745,21 +745,23 @@ program mksurfdat pftdata_mask(n) = 1 end if - ! Make sure sum of land cover types does not exceed 100. If it does, - ! subtract excess from most dominant land cover. - - suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + ! Make sure sum of all land cover types except natural vegetation + ! does not exceed 100. If it does, subtract excess from these land cover + ! types proportionally. + + suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() if (suma > 250._r4) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & - 'pcturb and pctgla is greater than 250%' - write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n) + 'pcturb, pctgla, and pctcrop is greater than 250%' + write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctcrop= ', & + n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctcft(n)%get_pct_l2g() call abort() else if (suma > 100._r4) then pctlak(n) = pctlak(n) * 100._r8/suma pctwet(n) = pctwet(n) * 100._r8/suma pcturb(n) = pcturb(n) * 100._r8/suma pctgla(n) = pctgla(n) * 100._r8/suma + call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma) end if end do @@ -1188,6 +1190,27 @@ program mksurfdat call change_landuse(ldomain, dynpft=.true.) + ! Make sure sum of all land cover types except natural vegetation + ! does not exceed 100. If it does, subtract excess from these land cover + ! types proportionally. + + do n = 1,ns_o + suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() + if (suma > 250._r4) then + write (6,*) subname, ' error: sum of pctlak, pctwet,', & + 'pcturb, pctgla, and pctcrop is greater than 250%' + write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctcrop= ', & + n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctcft(n)%get_pct_l2g() + call abort() + else if (suma > 100._r4) then + pctlak(n) = pctlak(n) * 100._r8/suma + pctwet(n) = pctwet(n) * 100._r8/suma + pcturb(n) = pcturb(n) * 100._r8/suma + pctgla(n) = pctgla(n) * 100._r8/suma + call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma) + end if + end do + call normalizencheck_landuse(ldomain) call update_max_array(pctnatpft_max,pctnatpft) @@ -1338,11 +1361,11 @@ subroutine normalizencheck_landuse(ldomain) ! Normalize land use and make sure things add up to 100% as well as ! checking that things are as they should be. ! -! Precondition: pctlak + pctwet + pcturb + pctgla <= 100 (within roundoff) +! Precondition: pctlak + pctwet + pcturb + pctgla + pctcrop <= 100 (within roundoff) ! ! !USES: - use mkpftConstantsMod , only : baregroundindex use mkpftUtilsMod , only : adjust_total_veg_area + use mkpftUtilsMod , only : adjust_total_natveg_area implicit none ! !ARGUMENTS: type(domain_type) :: ldomain @@ -1358,10 +1381,7 @@ subroutine normalizencheck_landuse(ldomain) integer :: nsmall_tot ! total number of small PFT values in all grid cells real(r8) :: suma ! sum for error check real(r8) :: suma2 ! another sum for error check - real(r8) :: new_total_veg_pct ! new % veg (% of grid cell, total of natural veg & crop) - real(r8) :: bare_pct_p2g ! % of bare soil, as % of grid cell - real(r8) :: bare_urb_diff ! difference between bare soil and urban % - real(r8) :: pcturb_excess ! excess urban % not accounted for by bare soil + real(r8) :: new_total_natveg_pct ! new % veg (% of grid cell, natural veg) real(r8) :: sum8, sum8a ! sum for error check real(r4) :: sum4a ! sum for error check real(r8), parameter :: tol_loose = 1.e-4_r8 ! tolerance for some 'loose' error checks @@ -1370,7 +1390,7 @@ subroutine normalizencheck_landuse(ldomain) !----------------------------------------------------------------------- ! ------------------------------------------------------------------------ - ! Normalize vegetated area so that vegetated + special area is 100% + ! Normalize natural vegetated area so that all landunits sum to 100% ! ------------------------------------------------------------------------ ns_o = ldomain%ns @@ -1397,79 +1417,34 @@ subroutine normalizencheck_landuse(ldomain) write(6,*) 'n, pctgla = ', n, pctgla(n) call abort() end if + if ( pctcft(n)%get_pct_l2g() < 0.0_r8 )then + write(6,*) subname, ' ERROR: pctcrop is negative!' + write(6,*) 'n, pctcrop = ', n, pctcft(n)%get_pct_l2g() + call abort() + end if - suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() if (suma > (100._r8 + tol_loose)) then - write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be' + write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla + pctcrop must be' write(6,*) '<= 100% before calling this subroutine' - write(6,*) 'n, pctlak, pctwet, pcturb, pctgla = ', & - n, pctlak(n), pctwet(n), pcturb(n), pctgla(n) + write(6,*) 'n, pctlak, pctwet, pcturb, pctgla, pctcrop = ', & + n, pctlak(n), pctwet(n), pcturb(n), pctgla(n), pctcft(n)%get_pct_l2g() call abort() end if - ! First normalize vegetated (natural veg + crop) cover so that the total of - ! (vegetated + (special excluding urban)) is 100%. We'll deal with urban later. - ! - ! Note that, in practice, the total area of natural veg + crop is typically 100% - ! going into this routine. However, the following code does NOT rely on this, and - ! will work properly regardless of the initial area of natural veg + crop (even if - ! that initial area is 0%). + ! Natural vegetated cover is 100% minus the sum of all other landunits - suma = pctlak(n)+pctwet(n)+pctgla(n) - new_total_veg_pct = 100._r8 - suma + suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)+pctcft(n)%get_pct_l2g() + new_total_natveg_pct = 100._r8 - suma ! correct for rounding error: - new_total_veg_pct = max(new_total_veg_pct, 0._r8) - - call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n)) - - ! Make sure we did the above rescaling correctly + new_total_natveg_pct = max(new_total_natveg_pct, 0._r8) - suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() - if (abs(suma - 100._r8) > tol_loose) then - write(6,*) subname, ' ERROR in rescaling veg based on (special excluding urban' - write(6,*) 'suma = ', suma - call abort() - end if - - ! Now decrease the vegetated area to account for urban area. Urban needs to be - ! handled specially because we replace bare soil preferentially with urban, rather - ! than rescaling all PFTs equally. - - if (pcturb(n) > 0._r8) then - - ! Replace bare soil preferentially with urban - bare_pct_p2g = pctnatpft(n)%get_one_pct_p2g(baregroundindex) - bare_urb_diff = bare_pct_p2g - pcturb(n) - bare_pct_p2g = max(0._r8, bare_urb_diff) - call pctnatpft(n)%set_one_pct_p2g(baregroundindex, bare_pct_p2g) - pcturb_excess = abs(min(0._r8,bare_urb_diff)) - - ! For any urban not accounted for by bare soil, replace other PFTs - ! proportionally - if (pcturb_excess > 0._r8) then - ! Note that, in this case, we will have already reduced bare ground to 0% - - new_total_veg_pct = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() - pcturb_excess - if (new_total_veg_pct < 0._r8) then - if (abs(new_total_veg_pct) < tol_loose) then - ! only slightly less than 0; correct it - new_total_veg_pct = 0._r8 - else - write(6,*) subname, ' ERROR: trying to replace veg with urban,' - write(6,*) 'but pcturb_excess exceeds current vegetation percent' - call abort() - end if - end if - - call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n)) - end if - - end if ! pcturb(n) > 0 + call adjust_total_natveg_area(new_total_natveg_pct, pctnatpft=pctnatpft(n)) ! Confirm that we have done the rescaling correctly: now the sum of all landunits ! should be 100% - suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n) - suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() + suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)+pctcft(n)%get_pct_l2g() + suma = suma + pctnatpft(n)%get_pct_l2g() if (abs(suma - 100._r8) > tol_loose) then write(6,*) subname, ' ERROR: landunits do not sum to 100%' write(6,*) 'n, suma, pctlak, pctwet, pctgla, pcturb, pctnatveg, pctcrop = ' @@ -1507,6 +1482,7 @@ subroutine normalizencheck_landuse(ldomain) end if ! Roundoff error fix + ! KO - Not sure if this needs to be changed or is necessary suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) suma2 = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. & @@ -1561,6 +1537,7 @@ subroutine normalizencheck_landuse(ldomain) end do ! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100% + ! KO - Not sure if this needs to be changed or is necessary if ( .not. outnc_double )then do n = 1,ns_o From 972fa298e02d7f7f619b4258f2fa2e75c3299029 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 20 Dec 2021 16:11:08 -0700 Subject: [PATCH 02/10] Remove adjust_total_veg_area --- tools/mksurfdata_map/src/mkpftUtilsMod.F90 | 39 ---------------------- tools/mksurfdata_map/src/mksurfdat.F90 | 1 - 2 files changed, 40 deletions(-) diff --git a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 index 53b8dd04a7..47d8956bfd 100644 --- a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 +++ b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 @@ -24,7 +24,6 @@ module mkpftUtilsMod ! !PUBLIC MEMBER FUNCTIONS: ! public :: convert_from_p2g ! Convert a p2g array into pct_pft_type objects - public :: adjust_total_veg_area ! Adjust the total vegetated area (natural veg & crop) to a new specified total public :: adjust_total_natveg_area ! Adjust the total natural vegetated area to a new specified total ! @@ -214,44 +213,6 @@ function get_default_cft() result(default_cft) end function get_default_cft - - !----------------------------------------------------------------------- - subroutine adjust_total_veg_area(new_total_pct, pctnatpft, pctcft) - ! - ! !DESCRIPTION: - ! Adjust the total vegetated area on the grid cell (natural veg & crop) to a new - ! specified total. - ! - ! If the old areas are 0%, then all the new area goes into pctnatpft. - ! - ! !USES: - use mkpctPftTypeMod, only : pct_pft_type - ! - ! !ARGUMENTS: - real(r8), intent(in) :: new_total_pct ! new total % of natural veg + crop landunits - class(pct_pft_type), intent(inout) :: pctnatpft ! natural veg cover information - class(pct_pft_type), intent(inout) :: pctcft ! crop cover information - ! - ! !LOCAL VARIABLES: - real(r8) :: natpft_l2g ! grid cell % cover of nat. veg. - real(r8) :: cft_l2g ! grid cell % cover of crop - real(r8) :: old_total ! old total % cover of natural veg + crop landunits - - character(len=*), parameter :: subname = 'adjust_total_veg_area' - !----------------------------------------------------------------------- - - natpft_l2g = pctnatpft%get_pct_l2g() - cft_l2g = pctcft%get_pct_l2g() - old_total = natpft_l2g + cft_l2g - if (old_total > 0._r8) then - call pctnatpft%set_pct_l2g(natpft_l2g * new_total_pct / old_total) - call pctcft%set_pct_l2g(cft_l2g * new_total_pct / old_total) - else - call pctnatpft%set_pct_l2g(new_total_pct) - end if - - end subroutine adjust_total_veg_area - !----------------------------------------------------------------------- subroutine adjust_total_natveg_area(new_total_pct, pctnatpft) ! diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 706569fde0..af1d01128f 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1364,7 +1364,6 @@ subroutine normalizencheck_landuse(ldomain) ! Precondition: pctlak + pctwet + pcturb + pctgla + pctcrop <= 100 (within roundoff) ! ! !USES: - use mkpftUtilsMod , only : adjust_total_veg_area use mkpftUtilsMod , only : adjust_total_natveg_area implicit none ! !ARGUMENTS: From 2f4f0e0cbf6916037e684e9defbba4b782647f3d Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 20 Dec 2021 18:11:04 -0700 Subject: [PATCH 03/10] Replace adjust_total_natveg_area subroutine with single equivalent line --- tools/mksurfdata_map/src/mkpftUtilsMod.F90 | 35 ---------------------- tools/mksurfdata_map/src/mksurfdat.F90 | 3 +- 2 files changed, 1 insertion(+), 37 deletions(-) diff --git a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 index 47d8956bfd..df5eecf6c7 100644 --- a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 +++ b/tools/mksurfdata_map/src/mkpftUtilsMod.F90 @@ -24,7 +24,6 @@ module mkpftUtilsMod ! !PUBLIC MEMBER FUNCTIONS: ! public :: convert_from_p2g ! Convert a p2g array into pct_pft_type objects - public :: adjust_total_natveg_area ! Adjust the total natural vegetated area to a new specified total ! ! !PRIVATE MEMBER FUNCTIONS: @@ -213,40 +212,6 @@ function get_default_cft() result(default_cft) end function get_default_cft - !----------------------------------------------------------------------- - subroutine adjust_total_natveg_area(new_total_pct, pctnatpft) - ! - ! !DESCRIPTION: - ! Adjust the total natural vegetated area on the grid cell to a new - ! specified total. - ! - ! If the old areas are 0%, then all the new area goes into pctnatpft. - ! - ! !USES: - use mkpctPftTypeMod, only : pct_pft_type - ! - ! !ARGUMENTS: - real(r8), intent(in) :: new_total_pct ! new total % of natural veg landunit - class(pct_pft_type), intent(inout) :: pctnatpft ! natural veg cover information - ! - ! !LOCAL VARIABLES: - real(r8) :: natpft_l2g ! grid cell % cover of natural veg - real(r8) :: old_total ! old total % cover of natural veg - - character(len=*), parameter :: subname = 'adjust_total_natveg_area' - !----------------------------------------------------------------------- - - natpft_l2g = pctnatpft%get_pct_l2g() - old_total = natpft_l2g - if (old_total > 0._r8) then - call pctnatpft%set_pct_l2g(natpft_l2g * new_total_pct / old_total) - else - call pctnatpft%set_pct_l2g(new_total_pct) - end if - - end subroutine adjust_total_natveg_area - - end module mkpftUtilsMod diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index af1d01128f..2b6480fe2e 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1364,7 +1364,6 @@ subroutine normalizencheck_landuse(ldomain) ! Precondition: pctlak + pctwet + pcturb + pctgla + pctcrop <= 100 (within roundoff) ! ! !USES: - use mkpftUtilsMod , only : adjust_total_natveg_area implicit none ! !ARGUMENTS: type(domain_type) :: ldomain @@ -1438,7 +1437,7 @@ subroutine normalizencheck_landuse(ldomain) ! correct for rounding error: new_total_natveg_pct = max(new_total_natveg_pct, 0._r8) - call adjust_total_natveg_area(new_total_natveg_pct, pctnatpft=pctnatpft(n)) + call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct) ! Confirm that we have done the rescaling correctly: now the sum of all landunits ! should be 100% From d7dbba18e3fcbff7fcdf1b5c6a4b47806b6d0178 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 20 Dec 2021 18:30:51 -0700 Subject: [PATCH 04/10] Remove duplicate suma in normalizencheck_landuse --- tools/mksurfdata_map/src/mksurfdat.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 2b6480fe2e..95e60415fb 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1432,7 +1432,6 @@ subroutine normalizencheck_landuse(ldomain) ! Natural vegetated cover is 100% minus the sum of all other landunits - suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)+pctcft(n)%get_pct_l2g() new_total_natveg_pct = 100._r8 - suma ! correct for rounding error: new_total_natveg_pct = max(new_total_natveg_pct, 0._r8) From 9d1d9a6d119a5fd7f5f7187fe4765ce10561f461 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Tue, 21 Dec 2021 09:07:10 -0700 Subject: [PATCH 05/10] Remove large threshold check for non-natural-veg landunits (two places) --- tools/mksurfdata_map/src/mksurfdat.F90 | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 95e60415fb..26c0858b87 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -750,13 +750,7 @@ program mksurfdat ! types proportionally. suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() - if (suma > 250._r4) then - write (6,*) subname, ' error: sum of pctlak, pctwet,', & - 'pcturb, pctgla, and pctcrop is greater than 250%' - write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctcrop= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctcft(n)%get_pct_l2g() - call abort() - else if (suma > 100._r4) then + if (suma > 100._r4) then pctlak(n) = pctlak(n) * 100._r8/suma pctwet(n) = pctwet(n) * 100._r8/suma pcturb(n) = pcturb(n) * 100._r8/suma @@ -1196,13 +1190,7 @@ program mksurfdat do n = 1,ns_o suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() - if (suma > 250._r4) then - write (6,*) subname, ' error: sum of pctlak, pctwet,', & - 'pcturb, pctgla, and pctcrop is greater than 250%' - write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctcrop= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctcft(n)%get_pct_l2g() - call abort() - else if (suma > 100._r4) then + if (suma > 100._r4) then pctlak(n) = pctlak(n) * 100._r8/suma pctwet(n) = pctwet(n) * 100._r8/suma pcturb(n) = pcturb(n) * 100._r8/suma From d60469a902c9d0c549c30af91083d475d236909b Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Thu, 23 Dec 2021 08:38:25 -0700 Subject: [PATCH 06/10] Rework roundoff error fix --- tools/mksurfdata_map/src/mksurfdat.F90 | 32 ++++++++++++-------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 26c0858b87..32bd5b44d7 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1366,7 +1366,6 @@ subroutine normalizencheck_landuse(ldomain) integer :: nsmall ! number of small PFT values for a single check integer :: nsmall_tot ! total number of small PFT values in all grid cells real(r8) :: suma ! sum for error check - real(r8) :: suma2 ! another sum for error check real(r8) :: new_total_natveg_pct ! new % veg (% of grid cell, natural veg) real(r8) :: sum8, sum8a ! sum for error check real(r4) :: sum4a ! sum for error check @@ -1467,30 +1466,29 @@ subroutine normalizencheck_landuse(ldomain) end if ! Roundoff error fix - ! KO - Not sure if this needs to be changed or is necessary - suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) - suma2 = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() + suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. & - (suma2 > 0.0_r8 .and. suma2 < 1.e-6_r8) ) then - write (6,*) 'Special land units near 100%, but not quite for n,suma =',n,suma - write (6,*) 'Adjusting special land units to 100%' - if (pctlak(n) >= 25._r8) then - pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n)) - else if (pctwet(n) >= 25._r8) then - pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n)) - else if (pcturb(n) >= 25._r8) then - pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n)) - else if (pctgla(n) >= 25._r8) then - pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n)) + (pctnatpft(n)%get_pct_l2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_l2g() < 1.e-6_r8) ) then + write (6,*) 'Special plus crop land units near 100%, but not quite for n,suma =',n,suma + write (6,*) 'Adjusting special plus crop land units to 100%' + if (pctlak(n) >= 1.0_r8) then + pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()) + else if (pctwet(n) >= 1.0_r8) then + pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()) + else if (pcturb(n) >= 1.0_r8) then + pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n) + pctcft(n)%get_pct_l2g()) + else if (pctgla(n) >= 1.0_r8) then + pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n) + pctcft(n)%get_pct_l2g()) + else if (pctcft(n)%get_pct_l2g() >= 1.0_r8) then + call pctcft(n)%set_pct_l2g(100._r8 - (pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n))) else - write (6,*) subname, 'Error: sum of special land units nearly 100% but none is >= 25% at ', & + write (6,*) subname, 'Error: sum of special plus crop land units nearly 100% but none is >= 1.0% at ', & 'n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma = ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),& pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(),suma call abort() end if call pctnatpft(n)%set_pct_l2g(0._r8) - call pctcft(n)%set_pct_l2g(0._r8) end if if ( any(pctnatpft(n)%get_pct_p2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_p2g() < toosmallPFT ) .or. & any(pctcft(n)%get_pct_p2g() > 0.0_r8 .and. pctcft(n)%get_pct_p2g() < toosmallPFT )) then From 77d334abc84336789b6ba7c570b7c2571f7902db Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Thu, 23 Dec 2021 08:52:08 -0700 Subject: [PATCH 07/10] Remove error check (if pctnatveg+pctcrop=0, then sum special landunits=100) --- tools/mksurfdata_map/src/mksurfdat.F90 | 40 -------------------------- 1 file changed, 40 deletions(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 32bd5b44d7..bd4740dbae 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1519,46 +1519,6 @@ subroutine normalizencheck_landuse(ldomain) end do - ! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100% - ! KO - Not sure if this needs to be changed or is necessary - - if ( .not. outnc_double )then - do n = 1,ns_o - sum8 = real(pctlak(n),r4) - sum8 = sum8 + real(pctwet(n),r4) - sum8 = sum8 + real(pcturb(n),r4) - sum8 = sum8 + real(pctgla(n),r4) - sum4a = real(pctnatpft(n)%get_pct_l2g(),r4) - sum4a = sum4a + real(pctcft(n)%get_pct_l2g(),r4) - if ( sum4a==0.0_r4 .and. sum8 < 100._r4-2._r4*epsilon(sum4a) )then - write (6,*) subname, ' error: sum of pctlak, pctwet,', & - 'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8 - write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n), & - pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g() - call abort() - end if - end do - else - do n = 1,ns_o - sum8 = pctlak(n) - sum8 = sum8 + pctwet(n) - sum8 = sum8 + pcturb(n) - sum8 = sum8 + pctgla(n) - sum8a = pctnatpft(n)%get_pct_l2g() - sum8a = sum8a + pctcft(n)%get_pct_l2g() - if ( sum8a==0._r8 .and. sum8 < (100._r8-4._r8*epsilon(sum8)) )then - write (6,*) subname, ' error: sum of pctlak, pctwet,', & - 'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8 - write (6,*) 'Total error, error/epsilon = ',100._r8-sum8, ((100._r8-sum8)/epsilon(sum8)) - write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,epsilon= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),& - pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(), epsilon(sum8) - call abort() - end if - end do - end if - ! Make sure that there is no vegetation outside the pft mask do n = 1,ns_o if (pftdata_mask(n) == 0 .and. (pctnatpft(n)%get_pct_l2g() > 0 .or. pctcft(n)%get_pct_l2g() > 0)) then From e6212cea39e322c661bf510bcec3d470b238b8ab Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Thu, 23 Dec 2021 09:12:26 -0700 Subject: [PATCH 08/10] Rework final landunit error check --- tools/mksurfdata_map/src/mksurfdat.F90 | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index bd4740dbae..5642011156 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1500,14 +1500,8 @@ subroutine normalizencheck_landuse(ldomain) call abort() end if - suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) - if (suma < 100._r8-epsilon(suma) .and. suma > (100._r8 - 4._r8*epsilon(suma))) then - write (6,*) subname, 'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', & - n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),& - pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g() - call abort() - end if - suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() + suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() + suma = suma + pctnatpft(n)%get_pct_l2g() if ( abs(suma-100._r8) > 1.e-10_r8) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla, pctnatveg and pctcrop is NOT equal to 100' From 2cc4b214c793ebe8df70f6683b9b39d95d924d24 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Thu, 23 Dec 2021 10:19:35 -0700 Subject: [PATCH 09/10] Add comment to roundoff error fix --- tools/mksurfdata_map/src/mksurfdat.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index 5642011156..bda294b2e0 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1465,7 +1465,11 @@ subroutine normalizencheck_landuse(ldomain) call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma) end if - ! Roundoff error fix + ! This roundoff error fix is needed to handle the situation where new_total_natveg_pct + ! ends up near 0 but not exactly 0 due to roundoff issues. In this situation, we set the + ! natveg landunit area to exactly 0 and put the remainder into some other landunit. Since + ! the remainder is very small, it doesn't really matter which other landunit we add it to, + ! so we just pick some landunit that already has at least 1% cover. suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g() if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. & (pctnatpft(n)%get_pct_l2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_l2g() < 1.e-6_r8) ) then From 97cc2f728e83f7c8486454e6955846ad093fba1f Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 27 Dec 2021 13:26:59 -0700 Subject: [PATCH 10/10] Remove some unused variables --- tools/mksurfdata_map/src/mksurfdat.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index bda294b2e0..48e927a274 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -1367,8 +1367,6 @@ subroutine normalizencheck_landuse(ldomain) integer :: nsmall_tot ! total number of small PFT values in all grid cells real(r8) :: suma ! sum for error check real(r8) :: new_total_natveg_pct ! new % veg (% of grid cell, natural veg) - real(r8) :: sum8, sum8a ! sum for error check - real(r4) :: sum4a ! sum for error check real(r8), parameter :: tol_loose = 1.e-4_r8 ! tolerance for some 'loose' error checks real(r8), parameter :: toosmallPFT = 1.e-10_r8 ! tolerance for PFT's to ignore character(len=32) :: subname = 'normalizencheck_landuse' ! subroutine name