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

Change renormalization of landunit areas in mksurfdata_map to be consistent with raw datasets #1564

Closed
wants to merge 10 commits into from
Closed

Change renormalization of landunit areas in mksurfdata_map to be consistent with raw datasets #1564

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
e0e41bf
Change renormalization of landunit areas in mksurfdata_map to be cons…
olyson Nov 29, 2021
972fa29
Remove adjust_total_veg_area
olyson Dec 20, 2021
2f4f0e0
Replace adjust_total_natveg_area subroutine with single equivalent line
olyson Dec 21, 2021
d7dbba1
Remove duplicate suma in normalizencheck_landuse
olyson Dec 21, 2021
9d1d9a6
Remove large threshold check for non-natural-veg landunits (two places)
olyson Dec 21, 2021
d60469a
Rework roundoff error fix
olyson Dec 23, 2021
77d334a
Remove error check (if pctnatveg+pctcrop=0, then sum special landunit…
olyson Dec 23, 2021
e6212ce
Rework final landunit error check
olyson Dec 23, 2021
2cc4b21
Add comment to roundoff error fix
olyson Dec 23, 2021
97cc2f7
Remove some unused variables
olyson Dec 27, 2021
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
42 changes: 1 addition & 41 deletions tools/mksurfdata_map/src/mkpftUtilsMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,7 @@ 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

!
! !PRIVATE MEMBER FUNCTIONS:
Expand Down Expand Up @@ -213,45 +212,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


end module mkpftUtilsMod


200 changes: 57 additions & 143 deletions tools/mksurfdata_map/src/mksurfdat.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -745,21 +745,17 @@ 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)
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)
call abort()
else if (suma > 100._r4) then
! 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 > 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
Expand Down Expand Up @@ -1188,6 +1184,21 @@ 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 > 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)
Expand Down Expand Up @@ -1338,11 +1349,9 @@ 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
implicit none
! !ARGUMENTS:
type(domain_type) :: ldomain
Expand All @@ -1357,11 +1366,7 @@ 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_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
Expand All @@ -1370,7 +1375,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
Expand All @@ -1397,79 +1402,33 @@ 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
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

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
new_total_natveg_pct = max(new_total_natveg_pct, 0._r8)

end if ! pcturb(n) > 0
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%
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 = '
Expand Down Expand Up @@ -1507,29 +1466,29 @@ subroutine normalizencheck_landuse(ldomain)
end if

! Roundoff error fix
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
billsacks marked this conversation as resolved.
Show resolved Hide resolved
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
Expand All @@ -1541,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'
Expand All @@ -1560,45 +1513,6 @@ subroutine normalizencheck_landuse(ldomain)

end do

! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100%

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
Expand Down