Distinguish between ocean and wetland in fsurdat files

python/ctsm/modify_input_files/modify_fsurdat.py

@@ -351,6 +351,7 @@ def zero_nonveg(self):
         self.setvar_lev0("PCT_WETLAND", 0)
         self.setvar_lev0("PCT_URBAN", 0)
         self.setvar_lev0("PCT_GLACIER", 0)
+        self.setvar_lev0("PCT_OCEAN", 0)
 
     def setvar_lev0(self, var, val):
         """
@@ -421,6 +422,7 @@ def set_idealized(self):
         self.setvar_lev0("PCT_LAKE", pct_not_nat_veg)
         self.setvar_lev0("PCT_URBAN", pct_not_nat_veg)
         self.setvar_lev0("PCT_GLACIER", pct_not_nat_veg)
+        self.setvar_lev0("PCT_OCEAN", pct_not_nat_veg)
         self.setvar_lev0("PCT_NATVEG", pct_nat_veg)
 
         for lev in self.file.nlevsoi:

python/ctsm/site_and_regional/single_point_case.py

@@ -421,6 +421,7 @@ def modify_surfdata_atpoint(self, f_orig):
             f_mod["PCT_WETLAND"][:, :] = 0.0
             f_mod["PCT_URBAN"][:, :, :] = 0.0
             f_mod["PCT_GLACIER"][:, :] = 0.0
+            f_mod["PCT_OCEAN"][:, :] = 0.0
 
             if self.dom_pft is not None:
                 max_dom_pft = max(self.dom_pft)

python/ctsm/test/test_sys_fsurdat_modifier.py

@@ -215,7 +215,7 @@ def test_opt_sections(self):
         np.testing.assert_array_equal(fsurdat_out_data.PCT_CROP, zero0d)
         np.testing.assert_array_equal(fsurdat_out_data.PCT_LAKE, zero0d)
         np.testing.assert_array_equal(fsurdat_out_data.PCT_WETLAND, zero0d)
-        np.testing.assert_array_equal(fsurdat_out_data.PCT_LAKE, zero0d)
+        np.testing.assert_array_equal(fsurdat_out_data.PCT_OCEAN, zero0d)
         np.testing.assert_array_equal(fsurdat_out_data.PCT_GLACIER, zero0d)
         np.testing.assert_array_equal(fsurdat_out_data.PCT_URBAN, pct_urban)
         np.testing.assert_array_equal(fsurdat_out_data.LAKEDEPTH, one0d * 200.0)

python/ctsm/test/test_unit_singlept_data_surfdata.py

@@ -106,6 +106,12 @@ class TestSinglePointCaseSurfaceNoCrop(unittest.TestCase):
                 coords={"lsmlat": lsmlat, "lsmlon": lsmlon},
                 attrs={"long_name": "percent wetland", "units": "unitless"},
             ),
+            "PCT_OCEAN": xr.DataArray(
+                data=np.random.rand(1, 1),
+                dims=["lsmlat", "lsmlon"],
+                coords={"lsmlat": lsmlat, "lsmlon": lsmlon},
+                attrs={"long_name": "percent ocean", "units": "unitless"},
+            ),
             "PCT_URBAN": xr.DataArray(
                 data=np.random.rand(1, 1, 3),
                 dims=["lsmlat", "lsmlon", "numurbl"],
@@ -656,6 +662,12 @@ class TestSinglePointCaseSurfaceCrop(unittest.TestCase):
                 coords={"lsmlat": lsmlat, "lsmlon": lsmlon},
                 attrs={"long_name": "percent wetland", "units": "unitless"},
             ),
+            "PCT_OCEAN": xr.DataArray(
+                data=np.random.rand(1, 1),
+                dims=["lsmlat", "lsmlon"],
+                coords={"lsmlat": lsmlat, "lsmlon": lsmlon},
+                attrs={"long_name": "percent ocean", "units": "unitless"},
+            ),
             "PCT_URBAN": xr.DataArray(
                 data=np.random.rand(1, 1, 3),
                 dims=["lsmlat", "lsmlon", "numurbl"],

python/ctsm/test/testinputs/modify_fsurdat_1x1mexicocity.cfg

@@ -36,6 +36,7 @@ PCT_NATVEG= 0.0
 PCT_GLACIER= 0.0
 PCT_WETLAND= 0.0
 PCT_LAKE = 0.0
+PCT_OCEAN = 0.0
 
 # Section with a list of variables to prcoess
 [modify_fsurdat_variable_list]

python/ctsm/test/testinputs/modify_fsurdat_opt_sections.cfg

@@ -33,6 +33,7 @@ PCT_CROP = 0.0
 PCT_LAKE = 0.0
 PCT_GLACIER = 0.0
 PCT_WETLAND = 0.0
+PCT_OCEAN = 0.0
 # NOTE: PCT_URBAN must be a list of three floats that sum to the total urban area
 PCT_URBAN = 100.0 0.0 0.0
 

src/main/surfrdMod.F90

@@ -390,6 +390,7 @@ subroutine surfrd_special(begg, endg, ncid, ns)
     real(r8),pointer :: pctgla(:)     ! percent of grid cell is glacier
     real(r8),pointer :: pctlak(:)     ! percent of grid cell is lake
     real(r8),pointer :: pctwet(:)     ! percent of grid cell is wetland
+    real(r8),pointer :: pctocn(:)     ! percent of grid cell is ocean
     real(r8),pointer :: pcturb(:,:)   ! percent of grid cell is urbanized
     integer ,pointer :: urban_region_id(:)
     real(r8),pointer :: pcturb_tot(:) ! percent of grid cell is urban (sum over density classes)
@@ -403,6 +404,7 @@ subroutine surfrd_special(begg, endg, ncid, ns)
     allocate(pctgla(begg:endg))
     allocate(pctlak(begg:endg))
     allocate(pctwet(begg:endg))
+    allocate(pctocn(begg:endg))
     allocate(pcturb(begg:endg,numurbl))
     allocate(pcturb_tot(begg:endg))
     allocate(urban_region_id(begg:endg))
@@ -412,6 +414,10 @@ subroutine surfrd_special(begg, endg, ncid, ns)
 
        ! Obtain non-grid surface properties of surface dataset other than percent patch
 
+    call ncd_io(ncid=ncid, varname='PCT_OCEAN', flag='read', data=pctocn, &
+         dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) call endrun( msg=' ERROR: PCT_OCEAN NOT on surfdata file'//errMsg(sourcefile, __LINE__))
+
     call ncd_io(ncid=ncid, varname='PCT_WETLAND', flag='read', data=pctwet, &
          dim1name=grlnd, readvar=readvar)
     if (.not. readvar) call endrun( msg=' ERROR: PCT_WETLAND  NOT on surfdata file'//errMsg(sourcefile, __LINE__))
@@ -475,9 +481,9 @@ subroutine surfrd_special(begg, endg, ncid, ns)
 
     topo_glc_mec(:,:) = max(topo_glc_mec(:,:), 0._r8)
 
-    pctspec = pctwet + pctlak + pcturb_tot + pctgla
+    pctspec = pctwet + pctlak + pcturb_tot + pctgla + pctocn
 
-    ! Error check: glacier, lake, wetland, urban sum must be less than 100
+    ! Error check: sum of glacier, lake, wetland, urban, ocean must be < 100
 
     found = .false.
     do nl = begg,endg
@@ -497,22 +503,26 @@ subroutine surfrd_special(begg, endg, ncid, ns)
 
     do nl = begg,endg
 
-       wt_lunit(nl,istdlak)     = pctlak(nl)/100._r8
+       wt_lunit(nl,istdlak) = pctlak(nl) / 100._r8
 
-       wt_lunit(nl,istwet)      = pctwet(nl)/100._r8
+       ! Until ctsm5.1 we would label ocean points as wetland in fsurdat
+       ! files. Starting with ctsm5.2 we label ocean points as ocean
+       ! (always 100%) and wetland points as wetland. Here we merge them
+       ! again to keep model behavior unchanged for now.
+       wt_lunit(nl,istwet) = (pctocn(nl) + pctwet(nl)) / 100._r8
 
-       wt_lunit(nl,istice)  = pctgla(nl)/100._r8
+       wt_lunit(nl,istice) = pctgla(nl) / 100._r8
 
        do n = isturb_MIN, isturb_MAX
           dens_index = n - isturb_MIN + 1
-          wt_lunit(nl,n)        = pcturb(nl,dens_index) / 100._r8
+          wt_lunit(nl,n) = pcturb(nl,dens_index) / 100._r8
        end do
 
     end do
 
     call CheckUrban(begg, endg, pcturb(begg:endg,:), subname)
 
-    deallocate(pctgla,pctlak,pctwet,pcturb,pcturb_tot,urban_region_id,pctspec)
+    deallocate(pctgla,pctlak,pctwet,pctocn,pcturb,pcturb_tot,urban_region_id,pctspec)
 
   end subroutine surfrd_special
 

tools/contrib/modify_singlept_site

@@ -231,6 +231,7 @@ if create_surfdata:
    f2['PCT_WETLAND'] = 0.
    f2['PCT_URBAN']   = 0.
    f2['PCT_GLACIER'] = 0.
+   f2['PCT_OCEAN'] = 0.
 
    #--  Overwrite global data with raw data  ----------------------------
    f2['LONGXY'] = plon

tools/mksurfdata_esmf/README

@@ -69,10 +69,10 @@ Building the executable (working in tools/mksurfdata_esmf)
 =======================
 Running for a single submission:
 =======================
-# Work in the ctsm_py environment, which requires the following steps:
+# Work in the ctsm_pylib environment, which requires the following steps:
 > module unload python; module load conda
 > cd ../..; ./py_env_create
-> conda activate ctsm_py; cd tools/mksurfdata_esmf
+> conda activate ctsm_pylib; cd tools/mksurfdata_esmf
 # to generate your target namelist:
 > ./gen_mksurfdata_namelist.py --help
 # for example try --res 1.9x2.5 --start-year 1850 --end-year 1850:

tools/mksurfdata_esmf/src/mkfileMod.F90

@@ -456,6 +456,9 @@ subroutine mkfile_define_vars(pioid, dynlanduse)
 
        end if
 
+       call mkpio_def_spatial_var(pioid=pioid, varname='PCT_OCEAN', xtype=xtype, &
+            long_name='percent ocean', units='unitless')
+
        call mkpio_def_spatial_var(pioid=pioid, varname='LAKEDEPTH', xtype=xtype, &
             long_name='lake depth', units='m')
 

tools/mksurfdata_esmf/src/mksurfdata.F90

@@ -162,6 +162,7 @@ program mksurfdata
   real(r8), allocatable           :: pctlak(:)               ! percent of grid cell that is lake
   real(r8), allocatable           :: pctlak_max(:)           ! maximum percent of grid cell that is lake
   real(r8), allocatable           :: pctwet(:)               ! percent of grid cell that is wetland
+  real(r8), allocatable           :: pctocn(:)               ! percent of grid cell that is ocean
   real(r8), allocatable           :: pctgla(:)               ! percent of grid cell that is glacier
   integer , allocatable           :: urban_region(:)         ! urban region ID
   real(r8), allocatable           :: pcturb(:)               ! percent of grid cell that is urbanized (total across all urban classes)
@@ -429,6 +430,13 @@ program mksurfdata
   call mkwetlnd(mksrf_fwetlnd_mesh, mksrf_fwetlnd, mesh_model, pctwet, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkwetlnd')
 
+  ! Initialize pctocn to zero.
+  ! Until ctsm5.1 we set pctwet = 100 at ocean points rather than
+  ! setting a pctocn. Starting with ctsm5.2, we set pctocn = 100 at
+  ! ocean points in subroutine normalize_and_check_landuse.
+  ! No regridding required.
+  allocate ( pctocn(lsize_o)); pctocn(:) = 0._r8
+
   ! -----------------------------------
   ! Make glacier fraction [pctgla] from [fglacier] dataset
   ! -----------------------------------
@@ -648,6 +656,10 @@ program mksurfdata
      call mkfile_output(pioid, mesh_model,  'PCT_WETLAND', pctwet, rc=rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in in mkfile_output for pctwet')
 
+     if (root_task)  write(ndiag, '(a)') trim(subname)//" writing out PCT_OCEAN"
+     call mkfile_output(pioid, mesh_model,  'PCT_OCEAN', pctocn, rc=rc)
+     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in in mkfile_output for pctocn')
+
      if (root_task)  write(ndiag, '(a)') trim(subname)//" writing PCT_NATVEG"
      call get_pct_l2g_array(pctnatpft, pctnatveg)
      call mkfile_output(pioid, mesh_model, 'PCT_NATVEG', pctnatveg, rc=rc)
@@ -1171,7 +1183,7 @@ subroutine normalize_and_check_landuse(ns_o)
 
          if (pct_land < 1.e-6_r8) then
             ! If we have essentially 0 land area, set land area to exactly 0 and put all
-            ! area in wetlands (to simulate ocean). Note that, based on the formulation
+            ! area in pctocn. Note that, based on the formulation
             ! for pct_land above, this should only arise if the non-natveg landunits
             ! already have near-zero area (and the natveg landunit should also have
             ! near-zero area in this case, because its area should be no greater than the
@@ -1184,7 +1196,8 @@ subroutine normalize_and_check_landuse(ns_o)
             pctlak(n) = 0._r8
             pcturb(n) = 0._r8
             pctgla(n) = 0._r8
-            pctwet(n) = 100._r8
+            pctwet(n) = 0._r8
+            pctocn(n) = 100._r8  ! the only asignment of non-zero ocean
          else
             ! Fill the rest of the land area with natveg, then renormalize landunits so
             ! that they are expressed as percent of the land area rather than percent of
@@ -1222,12 +1235,12 @@ subroutine normalize_and_check_landuse(ns_o)
          ! Confirm that we have done the rescaling correctly: now the sum of all landunits
          ! should be 100% within tol_loose
          suma = pctlak(n) + pctwet(n) + pctgla(n) + pcturb(n) + pctcft(n)%get_pct_l2g()
-         suma = suma + pctnatpft(n)%get_pct_l2g()
+         suma = suma + pctnatpft(n)%get_pct_l2g() + pctocn(n)
          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 = '
+            write(6,*) 'n, suma, pctlak, pctwet, pctgla, pcturb, pctnatveg, pctcrop, pctocn = '
             write(6,*) n, suma, pctlak(n), pctwet(n), pctgla(n), pcturb(n), &
-                 pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g()
+                 pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g(), pctocn(n)
             call shr_sys_abort()
          end if
 
@@ -1248,7 +1261,10 @@ subroutine normalize_and_check_landuse(ns_o)
          call pctcft(n)%remove_small_cover(toosmallPFT, nsmall)
          nsmall_tot = nsmall_tot + nsmall
 
-         suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
+         ! Include pctocn in suma but do not include in the
+         ! renormalization. When pctocn /= 0, it is 100, and
+         ! all other terms are 0.
+         suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctocn(n)
          suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
          if ( abs(suma - 100.0_r8) > 2.0*epsilon(suma) )then
             pctlak(n)    = pctlak(n)    * 100._r8/suma
@@ -1299,12 +1315,12 @@ subroutine normalize_and_check_landuse(ns_o)
          end if
 
          suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
-         suma = suma + pctnatpft(n)%get_pct_l2g()
+         suma = suma + pctnatpft(n)%get_pct_l2g() + pctocn(n)
          if ( abs(suma-100._r8) > 1.e-10_r8) then
-            write (6,*) subname, ' error: sum of pctlak, pctwet,', &
+            write (6,*) subname, ' error: sum of pctocn, pctlak, pctwet,', &
                  'pcturb, pctgla, pctnatveg and pctcrop is NOT equal to 100'
-            write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,sum= ', &
-                 n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
+            write (6,*)'n,pctcon,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,sum= ', &
+                 n,pctocn(n),pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(), suma
             call shr_sys_abort()
          end if