Merge branch 'fix_mksurf_soil_color'

Fix rare soil color bug in mksurfdata_map Under rare conditions, mksurfdata_map could put the default soil color in an output cell where there is actually more information. This tag fixes that issue. None of the out-of-the-box surface datasets are impacted by this bug, so I have not recreated any surface datasets. (I checked all out-of-the-box surface datasets except for surfdata_0.125x0.125_mp24_simyr2000_c150114.nc, because it doesn't get remade cleanly out of the box.) Also: - Add some unit tests for the creation of soil color in mksurfdata_map - Point to correct (existing) surface dataset for year-1850 at f05 resolution Resolves #4 (Minor bug in creation of soil color in mksurfdata_map: points can be given the default soil color, when they should have a real color) Resolves #765 (Year-1850 f05 surface dataset missing from inputdata repository)
ESCOMP · Jul 22, 2019 · 842185f · 842185f
2 parents 05fca10 + 37ca8da
commit 842185f
Show file tree

Hide file tree

Showing 11 changed files with 522 additions and 158 deletions.
diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -677,8 +677,6 @@ lnd/clm2/surfdata_map/surfdata_1x1_mexicocityMEX_16pfts_Irrig_CMIP6_simyr2000_c1
 lnd/clm2/surfdata_map/surfdata_1x1_urbanc_alpha_16pfts_Irrig_CMIP6_simyr2000_c171214.nc</fsurdat>
 
 <!-- for pre-industrial simulations - year 1850 without crop -->
-<fsurdat hgrid="0.47x0.63"  sim_year="1850" use_crop=".false." irrigate=".true.">
-lnd/clm2/surfdata_map/surfdata_0.47x0.63_16pfts_Irrig_CMIP6_simyr1850_c180508.nc</fsurdat>
 <fsurdat hgrid="360x720cru"   sim_year="1850" use_crop=".false." irrigate=".true." >
 lnd/clm2/surfdata_map/surfdata_360x720cru_16pfts_Irrig_CMIP6_simyr1850_c170824.nc</fsurdat>
 <fsurdat hgrid="48x96"        sim_year="1850" use_crop=".false." irrigate=".true.">

diff --git a/doc/ChangeLog b/doc/ChangeLog
@@ -1,4 +1,137 @@
 ===============================================================
+Tag name:  ctsm1.0.dev052
+Originator(s):  sacks (Bill Sacks)
+Date: Mon Jul 22 14:02:43 MDT 2019
+One-line Summary: Fix rare soil color bug in mksurfdata_map
+
+Purpose of changes
+------------------
+
+Under rare conditions, mksurfdata_map could put the default soil color
+in an output cell where there is actually more information. This tag
+fixes that issue. None of the out-of-the-box surface datasets are
+impacted by this bug, so I have not recreated any surface datasets. (I
+checked all out-of-the-box surface datasets except for
+surfdata_0.125x0.125_mp24_simyr2000_c150114.nc, because it doesn't get
+remade cleanly out of the box.)
+
+Also:
+
+- Add some unit tests for the creation of soil color in mksurfdata_map
+
+- Point to correct (existing) surface dataset for year-1850 at f05
+  resolution
+
+Bugs fixed or introduced
+------------------------
+
+Issues fixed (include CTSM Issue #):
+- Resolves ESCOMP/ctsm#4 (Minor bug in creation of soil color in
+  mksurfdata_map: points can be given the default soil color, when they
+  should have a real color)
+- Resolves ESCOMP/ctsm#765 (Year-1850 f05 surface dataset missing from
+  inputdata repository)
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+Notes of particular relevance for users
+---------------------------------------
+
+Caveats for users (e.g., need to interpolate initial conditions): none
+
+Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none
+
+Changes made to namelist defaults (e.g., changed parameter values): none
+
+Changes to the datasets (e.g., parameter, surface or initial files):
+- Year-1850 f05 surface dataset now points to a file that actually exists
+
+Substantial timing or memory changes: none
+
+Notes of particular relevance for developers: (including Code reviews and testing)
+---------------------------------------------
+NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide
+
+Caveats for developers (e.g., code that is duplicated that requires double maintenance):
+
+Changes to tests or testing: none
+
+Code reviewed by: self
+
+
+CTSM testing:
+
+ [PASS means all tests PASS and OK means tests PASS other than expected fails.]
+
+  build-namelist tests:
+
+    cheyenne - not run
+
+  tools-tests (test/tools):
+
+    cheyenne - ok
+
+    Most tests pass, including baseline comparisons. The following tests fail, but also failed for me on master:
+
+      019 smiS4 TSMscript_tools.sh ncl_scripts getregional_datasets.pl getregional ....................\c
+       rc=6 FAIL
+      020 bliS4 TBLscript_tools.sh ncl_scripts getregional_datasets.pl getregional ....................\c
+       rc=4 FAIL
+      027 smf84 TSMscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_clm4_5^buildtools ....................\c
+       rc=6 FAIL
+      028 blf84 TBLscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_clm4_5^buildtools ....................\c
+       rc=4 FAIL
+      029 smfc4 TSMscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_Cycle_clm4_5^buildtools ..............\c
+       rc=6 FAIL
+      030 blfc4 TBLscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_Cycle_clm4_5^buildtools ..............\c
+       rc=4 FAIL
+      031 smfg4 TSMscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_Global_clm4_5^buildtools .............\c
+       rc=6 FAIL
+      032 blfg4 TBLscript_tools.sh PTCLM PTCLMmkdata PTCLM_USUMB_Global_clm4_5^buildtools .............\c
+       rc=4 FAIL
+
+  PTCLM testing (tools/shared/PTCLM/test):
+
+    cheyenne - not run
+
+  python testing (see instructions in python/README.md; document testing done):
+
+    (any machine) - not run
+
+  regular tests (aux_clm):
+
+    cheyenne ---- not run
+    hobart ------ not run
+
+If the tag used for baseline comparisons was NOT the previous tag, note that here:
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: NO
+
+Detailed list of changes
+------------------------
+
+List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none
+
+Pull Requests that document the changes (include PR ids): none
+
+===============================================================
+===============================================================
 Tag name:  ctsm1.0.dev051
 Originator(s):  sacks (Bill Sacks)
 Date: Fri Jul 19 13:26:25 MDT 2019

diff --git a/doc/ChangeSum b/doc/ChangeSum
@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+  ctsm1.0.dev052    sacks 07/22/2019 Fix rare soil color bug in mksurfdata_map
   ctsm1.0.dev051    sacks 07/19/2019 Update water tracers for remainder of first stage of hydrology
   ctsm1.0.dev050   slevis 07/15/2019 dz --> dz_lake bug-fix in LakeTemperatureMod.F90 line 960
   ctsm1.0.dev049     erik 06/23/2019 Update mosart and intel to intel-19 on cheyenne

diff --git a/tools/mksurfdata_map/src/CMakeLists.txt b/tools/mksurfdata_map/src/CMakeLists.txt
@@ -17,6 +17,7 @@ list(APPEND mksurfdat_sources
   mkpctPftTypeMod.F90
   mkutilsMod.F90
   mkpftUtilsMod.F90
+  mksoilUtilsMod.F90
   mkvarctl.F90
   mkvarpar.F90
   shr_const_mod.F90

diff --git a/tools/mksurfdata_map/src/Srcfiles b/tools/mksurfdata_map/src/Srcfiles
@@ -12,6 +12,7 @@ mkpftConstantsMod.F90
 mkpftUtilsMod.F90
 mkpctPftTypeMod.F90
 mkpftMod.F90
+mksoilUtilsMod.F90
 mksoilMod.F90
 mkvocefMod.F90
 mksurfdat.F90

diff --git a/tools/mksurfdata_map/src/mkpftMod.F90 b/tools/mksurfdata_map/src/mkpftMod.F90
@@ -524,7 +524,8 @@ subroutine mkpft(ldomain, mapfname, fpft, ndiag, allow_no_crops, &
         call gridmap_areaave(tgridmap, pctcrop_i,   pctcrop_o,   nodata=0._r8)
 
         do m = 0, num_natpft
-           call gridmap_areaave_scs(tgridmap, pct_nat_pft_i(:,m), pct_nat_pft_o(:,m), nodata=0._r8,src_wt=pctnatveg_i*0.01_r8,dst_wt=pctnatveg_o*0.01_r8)
+           call gridmap_areaave_scs(tgridmap, pct_nat_pft_i(:,m), pct_nat_pft_o(:,m), &
+                nodata=0._r8,src_wt=pctnatveg_i*0.01_r8,dst_wt=pctnatveg_o*0.01_r8)
            do no = 1,ns_o
               if (pctlnd_o(no) < 1.0e-6 .or. pctnatveg_o(no) < 1.0e-6) then
                  if (m == 0) then
@@ -536,7 +537,8 @@ subroutine mkpft(ldomain, mapfname, fpft, ndiag, allow_no_crops, &
            enddo
         end do
         do m = 1, num_cft
-           call gridmap_areaave_scs(tgridmap, pct_cft_i(:,m), pct_cft_o(:,m), nodata=0._r8,src_wt=pctcrop_i*0.01_r8,dst_wt=pctcrop_o*0.01_r8)
+           call gridmap_areaave_scs(tgridmap, pct_cft_i(:,m), pct_cft_o(:,m), &
+                nodata=0._r8,src_wt=pctcrop_i*0.01_r8,dst_wt=pctcrop_o*0.01_r8)
            do no = 1,ns_o
               if (pctlnd_o(no) < 1.0e-6 .or. pctcrop_o(no) < 1.0e-6) then
                  if (m == 1) then

diff --git a/tools/mksurfdata_map/src/mksoilMod.F90 b/tools/mksurfdata_map/src/mksoilMod.F90
@@ -15,6 +15,7 @@ module mksoilMod
   use shr_kind_mod, only : r8 => shr_kind_r8, r4=>shr_kind_r4
   use shr_sys_mod , only : shr_sys_flush
   use mkdomainMod , only : domain_checksame
+  use mksoilUtilsMod, only : mkrank, dominant_soil_color
   implicit none
 
   SAVE
@@ -43,7 +44,6 @@ module mksoilMod
 ! !PRIVATE DATA MEMBERS:
 !
 ! !PRIVATE MEMBER FUNCTIONS:
-  private :: mkrank
   private :: mksoiltexInit  ! Soil texture Initialization
   private :: mksoilcolInit  ! Soil color Initialization
   private :: mksoilfmaxInit ! Soil fmax Initialization
@@ -597,21 +597,15 @@ subroutine mksoilcol(ldomain, mapfname, datfname, ndiag, &
 !EOP
   type(gridmap_type)    :: tgridmap
   type(domain_type)    :: tdomain          ! local domain
-  integer, parameter :: num=2               ! set soil mapunit number
-  integer  :: wsti(num)                     ! index to 1st and 2nd largest wst
-  real(r8), allocatable :: wst(:,:)         ! overlap weights, by surface type
   real(r8), allocatable :: gast_i(:)        ! global area, by surface type
   real(r8), allocatable :: gast_o(:)        ! global area, by surface type
   integer , allocatable :: soil_color_i(:)  ! input grid: BATS soil color
-  integer , allocatable :: color(:)         ! 0: none; 1: some
-  real(r8) :: wt                            ! map overlap weight
   real(r8) :: sum_fldi                      ! global sum of dummy input fld
   real(r8) :: sum_fldo                      ! global sum of dummy output fld
   character(len=35), allocatable :: col(:)  ! name of each color
-  integer  :: k,l,n,m,ni,no,ns_i,ns_o       ! indices
+  integer  :: k,l,m,ni,no,ns_i,ns_o       ! indices
   integer  :: ncid,dimid,varid              ! input netCDF id's
   integer  :: ier                           ! error status
-  integer  :: miss = 99999                  ! missing data indicator
   real(r8) :: relerr = 0.00001              ! max error: sum overlap wts ne 1
   character(len=32) :: subname = 'mksoilcol'
 !-----------------------------------------------------------------------
@@ -703,62 +697,14 @@ subroutine mksoilcol(ldomain, mapfname, datfname, ndiag, &
 
      call domain_checksame( tdomain, ldomain, tgridmap )
 
-     ! find area of overlap for each soil color for each no
-
-     allocate(wst(0:nsoicol,ns_o))
-     wst(0:nsoicol,:) = 0
-     allocate(color(ns_o))
-     color(:) = 0
-
-     ! TODO: need to do a loop to determine
-     ! the maximum number of over lap cells throughout the grid 
-     ! first get an array that is novr(ns_o) and fill this in - then set
-     ! maxovr - to max(novr) - then allocate the array wst to be size of
-     ! maxovr,ns_o or 0:nsoilcol,ns_o
-
-     do n = 1,tgridmap%ns
-        ni = tgridmap%src_indx(n)
-        no = tgridmap%dst_indx(n)
-        wt = tgridmap%wovr(n)
-        k  = soil_color_i(ni) * tdomain%mask(ni)
-        wst(k,no) = wst(k,no) + wt
-        if (k>0 .and. wst(k,no)>0.) then
-           color(no) = 1
-           wst(0,no) = 0.0
-        end if
-     enddo
-
-     soil_color_o(:) = 0
-     do no = 1,ns_o
-
-        ! Rank non-zero weights by color type. wsti(1) is the most extensive
-        ! color type. 
-
-        if (color(no) == 1) then
-           call mkrank (nsoicol, wst(0:nsoicol,no), miss, wsti, num)
-           soil_color_o(no) = wsti(1)
-        end if
-
-        ! If land but no color, set color to 15 (in older dataset generic 
-        ! soil color 4)
-
-        if (nsoicol == 8) then
-           if (soil_color_o(no)==0) soil_color_o(no) = 4
-        else if (nsoicol == 20) then
-           if (soil_color_o(no)==0) soil_color_o(no) = 15
-        end if
-
-        ! Error checks
+     ! Determine dominant soil color for each output cell
 
-        if (soil_color_o(no) < 0 .or. soil_color_o(no) > nsoicol) then
-           write (6,*) 'MKSOILCOL error: land model soil color = ', &
-                soil_color_o(no),' is not valid for lon,lat = ',no
-           call abort()
-        end if
-
-     enddo
-     deallocate (wst)
-     deallocate (color)
+     call dominant_soil_color( &
+          tgridmap = tgridmap, &
+          mask_i = tdomain%mask, &
+          soil_color_i = soil_color_i, &
+          nsoicol = nsoicol, &
+          soil_color_o = soil_color_o)
 
      ! Global sum of output field 
 
@@ -983,96 +929,6 @@ subroutine mkorganic(ldomain, mapfname, datfname, ndiag, organic_o)
 
 end subroutine mkorganic
 
-!-----------------------------------------------------------------------
-!BOP
-!
-! !ROUTINE: mkrank
-!
-! !INTERFACE:
-subroutine mkrank (n, a, miss, iv, num)
-!
-! !DESCRIPTION:
-! Return indices of largest [num] values in array [a]. Private method
-! only used for soil color and soil texture.
-!
-! !USES:
-!
-! !ARGUMENTS:
-  implicit none
-  integer , intent(in) :: n        !array length
-  real(r8), intent(in) :: a(0:n)   !array to be ranked
-  integer , intent(in) :: miss     !missing data value
-  integer , intent(in) :: num      !number of largest values requested
-  integer , intent(out):: iv(num)  !index to [num] largest values in array [a]
-!
-! !CALLED FROM:
-! subroutine mksoilcol 
-! subroutine mksoiltex
-!
-! !REVISION HISTORY:
-! Author: Gordon Bonan
-!
-!
-! !LOCAL VARIABLES:
-!EOP
-  real(r8) a_max       !maximum value in array
-  integer i            !array index
-  real(r8) delmax      !tolerance for finding if larger value
-  integer m            !do loop index
-  integer k            !do loop index
-  logical exclude      !true if data value has already been chosen
-!-----------------------------------------------------------------------
-
-  delmax = 1.e-06
-
-  ! Find index of largest non-zero number
-
-  iv(1) = miss
-  a_max = -9999.
-
-  do i = 0, n
-     if (a(i)>0. .and. (a(i)-a_max)>delmax) then
-        a_max = a(i)
-        iv(1)  = i
-     end if
-  end do
-
-  ! iv(1) = miss indicates no values > 0. this is an error
-
-  if (iv(1) == miss) then
-     write (6,*) 'MKRANK error: iv(1) = missing'
-     call abort()
-  end if
-
-  ! Find indices of the next [num]-1 largest non-zero number.
-  ! iv(m) = miss if there are no more values > 0
-
-  do m = 2, num
-     iv(m) = miss
-     a_max = -9999.
-     do i = 0, n
-
-        ! exclude if data value has already been chosen
-
-        exclude = .false.
-        do k = 1, m-1
-           if (i == iv(k)) exclude = .true.
-        end do
-
-        ! if not already chosen, see if it is the largest of
-        ! the remaining values
-
-        if (.not. exclude) then
-           if (a(i)>0. .and. (a(i)-a_max)>delmax) then
-              a_max = a(i)
-              iv(m)  = i
-           end if
-        end if
-     end do
-  end do
-
-end subroutine mkrank
-
 !-----------------------------------------------------------------------
 !BOP
 !