diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 7c5859e0ad..1aafae5659 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2645,6 +2645,13 @@ sub setup_logic_do_transient_lakes { my $var = 'do_transient_lakes'; + # Start by assuming a default value of '.true.'. Then check a number of + # conditions under which do_transient_lakes cannot be true. Under these + # conditions: (1) set default value to '.false.'; (2) make sure that the + # value is indeed false (e.g., that the user didn't try to set it to true). + + my $default_val = ".true."; + # cannot_be_true will be set to a non-empty string in any case where # do_transient_lakes should not be true; if it turns out that # do_transient_lakes IS true in any of these cases, a fatal error will be @@ -2668,7 +2675,7 @@ sub setup_logic_do_transient_lakes { # Note that, if the variable cannot be true, we don't call add_default # - so that we don't clutter up the namelist with variables that don't # matter for this case - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, val=>$default_val); } # Make sure the value is false when it needs to be false - i.e., that the @@ -2708,6 +2715,13 @@ sub setup_logic_do_transient_urban { my $var = 'do_transient_urban'; + # Start by assuming a default value of '.true.'. Then check a number of + # conditions under which do_transient_urban cannot be true. Under these + # conditions: (1) set default value to '.false.'; (2) make sure that the + # value is indeed false (e.g., that the user didn't try to set it to true). + + my $default_val = ".true."; + # cannot_be_true will be set to a non-empty string in any case where # do_transient_urban should not be true; if it turns out that # do_transient_urban IS true in any of these cases, a fatal error will be @@ -2731,7 +2745,7 @@ sub setup_logic_do_transient_urban { # Note that, if the variable cannot be true, we don't call add_default # - so that we don't clutter up the namelist with variables that don't # matter for this case - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, val=>$default_val); } # Make sure the value is false when it needs to be false - i.e., that the diff --git a/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm index 923abcdaec..aede147461 100644 --- a/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/smallville_dynlakes_monthly/user_nl_clm @@ -1,7 +1,7 @@ do_transient_lakes = .true. ! This file was created with the following command: -! ncap2 -s 'PCT_LAKE=array(0.0,0.0,PCT_CROP); PCT_LAKE={0.,50.,25.,25.,25.,25.}; HASLAKE=array(1.,1.,AREA); PCT_CROP=array(0.0,0.0,PCT_LAKE); PCT_CROP={0.,25.,12.,12.,12.,12.}' landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_c160127.nc landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200928.nc +! ncap2 -s 'PCT_LAKE=array(0.0,0.0,PCT_CROP); PCT_LAKE={0.,50.,25.,25.,25.,25.}; PCT_LAKE_MAX=array(1.,1.,AREA); PCT_CROP=array(0.0,0.0,PCT_LAKE); PCT_CROP={0.,25.,12.,12.,12.,12.}' landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_c160127.nc landuse.timeseries_1x1_smallvilleIA_hist_78pfts_simyr1850-1855_dynLakes_c200928.nc ! Key points are that lake area starts as 0, increases after the first year, then decreases after the second year. ! PCT_CROP is also changed so that PCT_LAKE + PCT_CROP <= 100. (Here, PCT_CROP increases and decreases at the same time as PCT_LAKE in order to exercise the simultaneous increase or decrease of two landunits, but that isn't a critical part of this test.) ! Note that the use of this file means that this testmod can only be used with the 1x1_smallvilleIA grid. diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index b01f2aa599..f1034d3510 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -16,7 +16,7 @@ module clm_initializeMod use clm_varctl , only : use_lch4, use_cn, use_cndv, use_c13, use_c14, use_fates use clm_varctl , only : use_soil_moisture_streams use clm_instur , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft - use clm_instur , only : irrig_method, wt_glc_mec, topo_glc_mec, haslake, pct_urban_max + use clm_instur , only : irrig_method, wt_glc_mec, topo_glc_mec, pct_lake_max, pct_urban_max use perf_mod , only : t_startf, t_stopf use readParamsMod , only : readParameters use ncdio_pio , only : file_desc_t @@ -218,7 +218,7 @@ subroutine initialize2(ni,nj) allocate (irrig_method (begg:endg, cft_lb:cft_ub )) allocate (wt_glc_mec (begg:endg, maxpatch_glc )) allocate (topo_glc_mec (begg:endg, maxpatch_glc )) - allocate (haslake (begg:endg )) + allocate (pct_lake_max (begg:endg )) allocate (pct_urban_max(begg:endg, numurbl )) ! Read list of Patches and their corresponding parameter values @@ -295,7 +295,7 @@ subroutine initialize2(ni,nj) ! Some things are kept until the end of initialize2; urban_valid is kept through the ! end of the run for error checking, pct_urban_max is kept through the end of the run ! for reweighting in subgridWeights. - deallocate (wt_lunit, wt_cft, wt_glc_mec, haslake) + deallocate (wt_lunit, wt_cft, wt_glc_mec, pct_lake_max) ! Determine processor bounds and clumps for this processor call get_proc_bounds(bounds_proc) diff --git a/src/main/clm_varsur.F90 b/src/main/clm_varsur.F90 index d360941d23..8c50d9810e 100644 --- a/src/main/clm_varsur.F90 +++ b/src/main/clm_varsur.F90 @@ -47,7 +47,7 @@ module clm_instur real(r8), pointer :: topo_glc_mec(:,:) ! whether we have lake to initialise in each grid cell - logical , pointer :: haslake(:) + real(r8), pointer :: pct_lake_max(:) ! whether we have urban to initialize in each grid cell ! (second dimension goes 1:numurbl) diff --git a/src/main/subgridMod.F90 b/src/main/subgridMod.F90 index 645d02a603..7987d85cd5 100644 --- a/src/main/subgridMod.F90 +++ b/src/main/subgridMod.F90 @@ -568,11 +568,11 @@ function lake_landunit_exists(gi) result(exists) ! ! !DESCRIPTION: ! Returns true if a land unit for lakes should be created in memory - ! which is defined for gridcells which will grow lake, given by haslake + ! which is defined for gridcells which will grow lake, given by pct_lake_max ! ! !USES: use dynSubgridControlMod , only : get_do_transient_lakes - use clm_instur , only : haslake + use clm_instur , only : pct_lake_max ! ! !ARGUMENTS: logical :: exists ! function result @@ -584,10 +584,10 @@ function lake_landunit_exists(gi) result(exists) !----------------------------------------------------------------------- if (get_do_transient_lakes()) then - ! To support dynamic landunits, we initialise a lake land unit in each grid cell in which there are lakes. - ! This is defined by the haslake variable + ! To support dynamic landunits, we initialise a lake land unit in + ! each grid cell in which there are lakes as defined by pct_lake_max - if (haslake(gi)) then + if (pct_lake_max(gi) > 0._r8) then exists = .true. else exists = .false. diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index 057a199e25..deb8a3d93f 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -800,7 +800,7 @@ subroutine surfrd_lakemask(begg, endg) ! Necessary for the initialisation of the lake land units ! ! !USES: - use clm_instur , only : haslake + use clm_instur , only : pct_lake_max use dynSubgridControlMod , only : get_flanduse_timeseries use clm_varctl , only : fname_len use fileutils , only : getfil @@ -836,9 +836,9 @@ subroutine surfrd_lakemask(begg, endg) call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0) ! read the lakemask - call ncd_io(ncid=ncid_dynuse, varname='HASLAKE' , flag='read', data=haslake, & + call ncd_io(ncid=ncid_dynuse, varname='PCT_LAKE_MAX' , flag='read', data=pct_lake_max, & dim1name=grlnd, readvar=readvar) - if (.not. readvar) call endrun( msg=' ERROR: HASLAKE is not on landuse.timeseries file'//errMsg(sourcefile, __LINE__)) + if (.not. readvar) call endrun( msg=' ERROR: PCT_LAKE_MAX is not on landuse.timeseries file'//errMsg(sourcefile, __LINE__)) ! close landuse_timeseries file again call ncd_pio_closefile(ncid_dynuse) diff --git a/tools/mksurfdata_esmf/README b/tools/mksurfdata_esmf/README index 309b57d611..7e1b49311c 100644 --- a/tools/mksurfdata_esmf/README +++ b/tools/mksurfdata_esmf/README @@ -13,10 +13,9 @@ Build Requirements: =================== mksurfdata_esmf is a distributed memory parallel program (using Message Passing -Interface -- MPI) that utilizes -both ESMF (Earth System Modelling Framework) for regridding as well as Parallel -I/O (PIO) and NetCDF output. As -such libraries must be built for the following: +Interface -- MPI) that utilizes both ESMF (Earth System Modelling Framework) +for regridding as well as PIO (Parallel I/O) and NetCDF output. As +such, libraries must be built for the following: 1) MPI 2) NetCDF @@ -25,22 +24,22 @@ such libraries must be built for the following: In addition for the build: python, bash-shell, CMake and GNU-Make are required -These libraries need to be built such that they can all work together in the same executable. Hence, the above -order may be required in building them. +These libraries need to be built such that they can all work together in the +same executable. Hence, the above order may be required in building them. ========================================= Use cime to manage the build requirements ========================================= -For users working on cime machines you can use the build -script to build the tool for you. On other machines you'll need to do a port to cime -and tell how to build for that machine. That's talked about in the cime documentation. +For users working on cime machines you can use the build script to build the +tool. On other machines you'll need to do a port to cime and tell how to build +for that machine. That's talked about in the cime documentation. And you'll have to make some modifications to the build script. https://github.com/ESMCI/cime/wiki/Porting-Overview -Machines that already run CTSM or CESM have been ported to cime. So if you can run the model -on your machine you will be able to build the tool there. +Machines that already run CTSM or CESM have been ported to cime. So if you can +run the model on your machine, you will be able to build the tool there. To get a list of the machines that have been ported to cime: @@ -48,12 +47,13 @@ cd ../../cime/scripts # assumes we are in tools/mksurfdata_esmf ./query_config --machines NOTE: -In addition to having a port to cime, the machine also needs to have PIO built and able -to be referenced with the env variable PIO which will need to be in the porting instructions -for the machine. Currently an independent PIO library is not available on cime ported machines. +In addition to having a port to cime, the machine also needs to have PIO built +and able to be referenced with the env variable PIO which will need to be in +the porting instructions for the machine. Currently an independent PIO library +is not available on cime ported machines. ======================= -building the executable (working in tools/mksurfdata_esmf) +Building the executable (working in tools/mksurfdata_esmf) ======================= > ./gen_mksurfdata_build.sh # For machines with a cime build @@ -63,7 +63,7 @@ building the executable (working in tools/mksurfdata_esmf) # a default value does get set for other machines. ======================= -running for a single submission: +Running for a single submission: ======================= # to generate your target namelist: > ./gen_mksurfdata_namelist.py --help @@ -83,7 +83,7 @@ running for a single submission: # Read note about regional grids below. ======================= -running for the generation of multiple datasets +Running for the generation of multiple datasets ======================= # Notes: # - gen_mksurfdata_jobscript_multi.py runs ./gen_mksurfdata_namelist.py for you @@ -102,6 +102,7 @@ slevis HAS RUN THESE CASES and HAS LISTED ISSUES ENCOUNTERED ------------------------------------------------------------ REMEMBER TO compare against existing fsurdat files in /glade/p/cesmdata/cseg/inputdata/lnd/clm2/surfdata_map/release-clm5.0.18 +0) New 30-sec raw data for soil texture fails. Try requesting more mem. 1) Soil color & texture and ag fire peak month outputs too high in .log TODO? Change frac_o from always 1. 2) Pct lake has chged in the .log bc the old diagnostic omitted mask_i frac_o diff --git a/tools/mksurfdata_esmf/gen_mksurfdata_namelist.py b/tools/mksurfdata_esmf/gen_mksurfdata_namelist.py index 921028d607..f9226cac2c 100755 --- a/tools/mksurfdata_esmf/gen_mksurfdata_namelist.py +++ b/tools/mksurfdata_esmf/gen_mksurfdata_namelist.py @@ -196,7 +196,19 @@ def get_parser(): [default: %(default)s] """, action="store_true", - dest="hres_flag", + dest="hres_pft", + default=False, + ) + parser.add_argument( + "--hires_soitex", + help=""" + If you want to use the high-resolution soil texture dataset rather + than the default lower resolution dataset. + (Low resolution is 5x5min, high resolution 30-second) + [default: %(default)s] + """, + action="store_true", + dest="hres_soitex", default=False, ) parser.add_argument( @@ -230,19 +242,6 @@ def get_parser(): dest="potveg_flag", default=False, ) - parser.add_argument( - "--merge_gis", - help=""" - If you want to use the glacier dataset that merges in - the Greenland Ice Sheet data that CISM uses (typically - used only if consistency with CISM is important) - [default: %(default)s] - """, - action="store", - dest="merge_gis", - choices=["on","off"], - default="off", - ) return parser def main (): @@ -267,8 +266,8 @@ def main (): glc_flag = args.glc_flag potveg = args.potveg_flag glc_nec = args.glc_nec - merge_gis = args.merge_gis - if args.hres_flag: + + if args.hres_pft: if (start_year == 1850 and end_year == 1850) or \ (start_year == 2005 and end_year == 2005): hires_pft = 'on' @@ -279,6 +278,13 @@ def main (): else: hires_pft = 'off' + if args.hres_soitex: + hires_soitex = 'on' + else: + hires_soitex = 'off' + + verbose = args.verbose + if force_model_mesh_file is not None: # open mesh_file to read element_count and, if available, orig_grid_dims mesh_file = netCDF4.Dataset(force_model_mesh_file, 'r') @@ -369,10 +375,10 @@ def main (): # create attribute list for parsing xml file attribute_list = {'hires_pft':hires_pft, + 'hires_soitex':hires_soitex, 'pft_years':pft_years, 'pft_years_ssp':pft_years_ssp, 'ssp_rcp':ssp_rcp, - 'mergeGIS':merge_gis, 'res':res} # create dictionary for raw data files names @@ -431,9 +437,10 @@ def main (): print('WARNING: run ./download_input_data to try TO ' \ 'OBTAIN MISSING FILES') _must_run_download_input_data = True - elif 'urban_properties' in rawdata_files[child1.tag]: - # Time-slice cases pull urban_properties from the transient - # urban_properties data files + elif 'urban_properties' in rawdata_files[child1.tag] or \ + 'lake' in rawdata_files[child1.tag]: + # Time-slice cases pull urban_properties and %lake from the + # corresponding transient files rawdata_files[child1.tag] = rawdata_files[child1.tag]. \ replace("%y",str(start_year)) @@ -455,6 +462,10 @@ def main (): new_key = f"{child1.tag}_urban" rawdata_files[new_key] = os.path.join(input_path, item.text) + if item.tag == 'lookup_filename': + new_key = f"{child1.tag}_lookup" + rawdata_files[new_key] = os.path.join(input_path, item.text) + # determine output mesh xml_path = os.path.join(tool_path, '../../ccs_config/component_grids_nuopc.xml') tree2 = ET.parse(xml_path) diff --git a/tools/mksurfdata_esmf/gen_mksurfdata_namelist.xml b/tools/mksurfdata_esmf/gen_mksurfdata_namelist.xml index fd86bfd364..acbf408b72 100644 --- a/tools/mksurfdata_esmf/gen_mksurfdata_namelist.xml +++ b/tools/mksurfdata_esmf/gen_mksurfdata_namelist.xml @@ -20,16 +20,9 @@ - + - - - lnd/clm2/rawdata/mksrf_organic_10level_5x5min_ISRIC-WISE-NCSCD_nlev7_c120830.nc - lnd/clm2/mappingdata/grids/UNSTRUCTgrid_5x5min_nomask_cdf5_c200129.nc - - - @@ -50,9 +43,15 @@ - lnd/clm2/rawdata/mksrf_soitex.10level.c010119.nc + lnd/clm2/rawdata/mksrf_soil_mapunits_5x5min_WISE.c220330.nc + lnd/clm2/rawdata/mksrf_soil_lookup.10level.WISE.c220330.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_5x5min_nomask_cdf5_c200129.nc + + lnd/clm2/rawdata/mksrf_soil_mapunits_30sec_WISE.c220330.nc + lnd/clm2/rawdata/mksrf_soil_lookup.10level.WISE.c220330.nc + lnd/clm2/rawdata/mksrf_soil_mapunits_30sec_WISE.c220330.nc + @@ -66,12 +65,21 @@ - + - lnd/clm2/rawdata/mksrf_LakePnDepth_3x3min_simyr2004_csplk_c151015.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc - + + + + + + + lnd/clm2/rawdata/mksrf_LakePnDepth_3x3min_simyr2017_MODISgrid.cdf5.c200305.nc + lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc + + @@ -97,11 +105,11 @@ - lnd/clm2/rawdata/gao_oneill_urban/historical/urban_properties_GaoOneil_05deg_ThreeClass_%y_c20220910.nc + lnd/clm2/rawdata/gao_oneill_urban/historical/urban_properties_GaoOneil_05deg_ThreeClass_%y_cdf5_c20220910.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/pot_veg/mksrf_urban_0.05x0.05_zerourbanpct.cdf5.c181014.nc + lnd/clm2/rawdata/mksrf_urban_0.05x0.05_zerourbanpct.cdf5.c181014.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc @@ -134,11 +142,7 @@ version of the raw dataset will probably go away. - lnd/clm2/rawdata/mksrf_glacier_3x3min_simyr2000.cdf5.c120926.nc - lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc - - - lnd/clm2/rawdata/mksrf_glacier_3x3min_simyr2000_mergeGreenland.c120921.nc + lnd/clm2/rawdata/mksrf_glacier_3x3min_simyr2000.c20210708.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_3x3min_nomask_cdf5_c200129.nc @@ -160,7 +164,7 @@ version of the raw dataset will probably go away. lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.5x0.5_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/pot_veg/mksrf_gdp_0.5x0_zerogdp.cdf5.c200413.nc + lnd/clm2/rawdata/mksrf_gdp_0.5x0.5_zerogdp.cdf5.c200413.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.5x0.5_nomask_cdf5_c200129.nc @@ -185,7 +189,7 @@ version of the raw dataset will probably go away. lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.5x0.5_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/pot_veg/mksrf_abm_0.5x0.5_missingabm.cdf5.c200413.nc + lnd/clm2/rawdata/mksrf_abm_0.5x0.5_missingabm.cdf5.c200413.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.5x0.5_nomask_cdf5_c200129.nc @@ -219,21 +223,21 @@ version of the raw dataset will probably go away. - /glade/p/cgd/tss/people/slevis/pot_veg/mksrf_landuse_potvegclm50_LUH2.cdf5.c181106.nc + lnd/clm2/rawdata/mksrf_landuse_potvegclm50_LUH2.cdf5.c181106.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.LUH2.histsimyr1850-2015.cdf5.c170629/mksrf_landuse_histclm50_LUH2_1850.cdf5.c170629.nc + /glade/p/cesm/sdwg_dev/thesis/data/cesm_tools/surfacedata/ctsm52finalrawdata/globalctsm52nrhistLUH2deg025_220918/mksrf_landuse_ctsm52_nrhistLUH2_1850.c220918.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.LUH2.histsimyr1850-2015.cdf5.c170629/mksrf_landuse_histclm50_LUH2_2000.cdf5.c170629.nc + /glade/p/cesm/sdwg_dev/thesis/data/cesm_tools/surfacedata/ctsm52finalrawdata/globalctsm52nrhistLUH2deg025_220918/mksrf_landuse_ctsm52_nrhistLUH2_2000.c220918.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.LUH2.histsimyr1850-2015.cdf5.c170629/mksrf_landuse_histclm50_LUH2_2005.cdf5.c170629.nc + /glade/p/cesm/sdwg_dev/thesis/data/cesm_tools/surfacedata/ctsm52finalrawdata/globalctsm52nrhistLUH2deg025_220918/mksrf_landuse_ctsm52_nrhistLUH2_2005.c220918.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc @@ -249,10 +253,10 @@ version of the raw dataset will probably go away. - lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.LUH2.histsimyr1850-2015.c20220305/mksrf_landuse_ctsm52_histLUH2_%y.c220305.nc + /glade/p/cesm/sdwg_dev/thesis/data/cesm_tools/surfacedata/ctsm52finalrawdata/globalctsm52nrhistLUH2deg025_220918/mksrf_landuse_ctsm52_nrhistLUH2_%y.c220918.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/historical/urban_properties_GaoOneil_05deg_ThreeClass_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/historical/urban_properties_GaoOneil_05deg_ThreeClass_%y_cdf5_c20220910.nc @@ -265,57 +269,57 @@ version of the raw dataset will probably go away. lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP1-2.6.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP1RCP26_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp1/urban_properties_GaoOneil_05deg_ThreeClass_ssp1_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp1/urban_properties_GaoOneil_05deg_ThreeClass_ssp1_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP1-1.9.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP1RCP19_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp1/urban_properties_GaoOneil_05deg_ThreeClass_ssp1_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp1/urban_properties_GaoOneil_05deg_ThreeClass_ssp1_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP2-4.5.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP2RCP45_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp2/urban_properties_GaoOneil_05deg_ThreeClass_ssp2_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp2/urban_properties_GaoOneil_05deg_ThreeClass_ssp2_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP3-7.0.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP3RCP70_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp3/urban_properties_GaoOneil_05deg_ThreeClass_ssp3_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp3/urban_properties_GaoOneil_05deg_ThreeClass_ssp3_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP4-3.4.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP4RCP34_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp4/urban_properties_GaoOneil_05deg_ThreeClass_ssp4_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp4/urban_properties_GaoOneil_05deg_ThreeClass_ssp4_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP4-6.0.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP4RCP60_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp4/urban_properties_GaoOneil_05deg_ThreeClass_ssp4_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp4/urban_properties_GaoOneil_05deg_ThreeClass_ssp4_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP5-8.5.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP5RCP85_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp5/urban_properties_GaoOneil_05deg_ThreeClass_ssp5_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp5/urban_properties_GaoOneil_05deg_ThreeClass_ssp5_%y_cdf5_c20220910.nc lnd/clm2/rawdata/pftcftdynharv.0.25x0.25.SSP5-3.4.simyr2015-2100.c20220305/mksrf_landuse_ctsm52_SSP5RCP34_%y.c220305.nc lnd/clm2/mappingdata/grids/UNSTRUCTgrid_0.25x0.25_nomask_cdf5_c200129.nc - /glade/p/cgd/tss/people/slevis/transient_lake_data/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc - lnd/clm2/rawdata/gao_oneill_urban/ssp5/urban_properties_GaoOneil_05deg_ThreeClass_ssp5_%y_c20220910.nc + lnd/clm2/rawdata/lake_area/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_%y.cdf5.c20220325.nc + lnd/clm2/rawdata/gao_oneill_urban/ssp5/urban_properties_GaoOneil_05deg_ThreeClass_ssp5_%y_cdf5_c20220910.nc diff --git a/tools/mksurfdata_esmf/src/CMakeLists.txt b/tools/mksurfdata_esmf/src/CMakeLists.txt index 8720f3188a..7d845d27bf 100644 --- a/tools/mksurfdata_esmf/src/CMakeLists.txt +++ b/tools/mksurfdata_esmf/src/CMakeLists.txt @@ -19,7 +19,6 @@ set(SRCFILES mkagfirepkmonthMod.F90 mkinputMod.F90 mklaiMod.F90 mklanwatMod.F90 - mkorganicMod.F90 mkpeatMod.F90 mkpioMod.F90 mkpftMod.F90 diff --git a/tools/mksurfdata_esmf/src/mkfileMod.F90 b/tools/mksurfdata_esmf/src/mkfileMod.F90 index 4582e584b8..09f406bc9e 100644 --- a/tools/mksurfdata_esmf/src/mkfileMod.F90 +++ b/tools/mksurfdata_esmf/src/mkfileMod.F90 @@ -135,8 +135,10 @@ subroutine mkfile_define_atts(pioid, dynlanduse) ! Raw data file names str = get_filename(mksrf_fgrid_mesh) rcode = pio_put_att(pioid, pio_global, 'Input_grid_dataset', trim(str)) - str = get_filename(mksrf_flakwat) - rcode = pio_put_att(pioid, pio_global, 'Inland_lake_raw_data_file_name', trim(str)) + str = get_filename(mksrf_fpctlak) + rcode = pio_put_att(pioid, pio_global, 'Percent_lake_raw_data_file_name', trim(str)) + str = get_filename(mksrf_flakdep) + rcode = pio_put_att(pioid, pio_global, 'Lake_depth_raw_data_file_name', trim(str)) str = get_filename(mksrf_fwetlnd) rcode = pio_put_att(pioid, pio_global, 'Inland_wetland_raw_data_file_name', trim(str)) str = get_filename(mksrf_fglacier) @@ -167,8 +169,10 @@ subroutine mkfile_define_atts(pioid, dynlanduse) ! Mesh file names str = get_filename(mksrf_fvegtyp_mesh) rcode = pio_put_att(pioid, pio_global, 'mesh_pft_file_name', trim(str)) - str = get_filename(mksrf_flakwat_mesh) - rcode = pio_put_att(pioid, pio_global, 'mesh_lakwat_file', trim(str)) + str = get_filename(mksrf_fpctlak_mesh) + rcode = pio_put_att(pioid, pio_global, 'mesh_pctlak_file', trim(str)) + str = get_filename(mksrf_flakdep_mesh) + rcode = pio_put_att(pioid, pio_global, 'mesh_lakdep_file', trim(str)) str = get_filename(mksrf_fwetlnd_mesh) rcode = pio_put_att(pioid, pio_global, 'mesh_wetlnd_file', trim(str)) str = get_filename(mksrf_fglacier_mesh) @@ -179,8 +183,6 @@ subroutine mkfile_define_atts(pioid, dynlanduse) rcode = pio_put_att(pioid, pio_global, 'mesh_soil_texture_file', trim(str)) str = get_filename(mksrf_fsoicol_mesh) rcode = pio_put_att(pioid, pio_global, 'mesh_soil_color_file', trim(str)) - str = get_filename(mksrf_forganic_mesh) - rcode = pio_put_att(pioid, pio_global, 'mesh_soil_organic_file', trim(str)) str = get_filename(mksrf_furban_mesh) rcode = pio_put_att(pioid, pio_global, 'mesh_urban_file', trim(str)) str = get_filename(mksrf_fmax_mesh) @@ -216,11 +218,11 @@ subroutine mkfile_define_atts(pioid, dynlanduse) str = get_filename(mksrf_fsoicol) rcode = pio_put_att(pioid, pio_global, 'soil_color_raw_data_file_name', trim(str)) str = get_filename(mksrf_fsoitex) - rcode = pio_put_att(pioid, pio_global, 'soil_texture_raw_data_file_name', trim(str)) + rcode = pio_put_att(pioid, pio_global, 'soil_texture_mapunit_raw_data_file_name', trim(str)) + str = get_filename(mksrf_fsoitex_lookup) + rcode = pio_put_att(pioid, pio_global, 'soil_texture_lookup_raw_data_file_name', trim(str)) str = get_filename(mksrf_fmax) rcode = pio_put_att(pioid, pio_global, 'fmax_raw_data_file_name', trim(str)) - str = get_filename(mksrf_forganic) - rcode = pio_put_att(pioid, pio_global, 'organic_matter_raw_data_file_name', trim(str)) str = get_filename(mksrf_fvocef) rcode = pio_put_att(pioid, pio_global, 'VOC_EF_raw_data_file_name', trim(str)) end if @@ -261,21 +263,36 @@ subroutine mkfile_define_vars(pioid, dynlanduse) call mkpio_def_spatial_var(pioid=pioid, varname='SOIL_COLOR', xtype=PIO_INT, & long_name='soil color', units='unitless') - call mkpio_def_spatial_var(pioid=pioid, varname='PCT_SAND', xtype=xtype, & + call mkpio_def_spatial_var(pioid=pioid, varname='PCT_SAND', xtype=PIO_REAL, & lev1name='nlevsoi', & long_name='percent sand', units='unitless') - call mkpio_def_spatial_var(pioid=pioid, varname='PCT_CLAY', xtype=xtype, & + call mkpio_def_spatial_var(pioid=pioid, varname='PCT_CLAY', xtype=PIO_REAL, & lev1name='nlevsoi', & long_name='percent clay', units='unitless') call mkpio_def_spatial_var(pioid=pioid, varname='mapunits', xtype=PIO_INT, & long_name='soil texture map units', units='unitless') - call mkpio_def_spatial_var(pioid=pioid, varname='ORGANIC', xtype=xtype, & + call mkpio_def_spatial_var(pioid=pioid, varname='ORGANIC', xtype=PIO_REAL, & lev1name='nlevsoi', & long_name='organic matter density at soil levels', & units='kg/m3 (assumed carbon content 0.58 gC per gOM)') + call mkpio_def_spatial_var(pioid=pioid, varname='ORGC', xtype=PIO_REAL, & + lev1name='nlevsoi', & + long_name='soil organic carbon', units='gC/kg soil') + + call mkpio_def_spatial_var(pioid=pioid, varname='BULK', xtype=PIO_REAL, & + lev1name='nlevsoi', & + long_name='bulk density', units='g cm-3') + + call mkpio_def_spatial_var(pioid=pioid, varname='CFRAG', xtype=PIO_REAL, & + lev1name='nlevsoi', & + long_name='coarse fragments', units='vol% > 2 mm') + + call mkpio_def_spatial_var(pioid=pioid, varname='PHAQ', xtype=PIO_REAL, & + lev1name='nlevsoi', & + long_name='pH measured in water', units='unitless') call mkpio_def_spatial_var(pioid=pioid, varname='FMAX', xtype=xtype, & long_name='maximum fractional saturated area', units='unitless') diff --git a/tools/mksurfdata_esmf/src/mkinputMod.F90 b/tools/mksurfdata_esmf/src/mkinputMod.F90 index 1d5c3a38a9..465d653313 100644 --- a/tools/mksurfdata_esmf/src/mkinputMod.F90 +++ b/tools/mksurfdata_esmf/src/mkinputMod.F90 @@ -44,13 +44,11 @@ module mkinputMod character(CX) , public :: mksrf_fhrvtyp = ' ' ! harvest data file name character(CX) , public :: mksrf_fhrvtyp_mesh = ' ' ! harvest mesh file name - character(CX) , public :: mksrf_forganic = ' ' ! organic matter data file name - character(CX) , public :: mksrf_forganic_mesh = ' ' ! organic matter mesh file name - character(CX) , public :: mksrf_fsoicol = ' ' ! soil color data file name character(CX) , public :: mksrf_fsoicol_mesh = ' ' ! soil color mesh file name - character(CX) , public :: mksrf_fsoitex = ' ' ! soil texture data file name + character(CX) , public :: mksrf_fsoitex = ' ' ! soil texture mapunit data file name + character(CX) , public :: mksrf_fsoitex_lookup = ' ' ! soil texture lookup data file name character(CX) , public :: mksrf_fsoitex_mesh = ' ' ! soil texture mesh file name character(CX) , public :: mksrf_fmax = ' ' ! fmax data file name @@ -68,8 +66,10 @@ module mkinputMod character(CX) , public :: mksrf_fgdp = ' ' ! gdp data file names character(CX) , public :: mksrf_fgdp_mesh = ' ' ! gdp mesh file names - character(CX) , public :: mksrf_flakwat = ' ' ! inland lake data file name - character(CX) , public :: mksrf_flakwat_mesh = ' ' ! inland lake mesh file name + character(CX) , public :: mksrf_fpctlak = ' ' ! percent lake data file name + character(CX) , public :: mksrf_fpctlak_mesh = ' ' ! percent lake file name + character(CX) , public :: mksrf_flakdep = ' ' ! lake depth data file name + character(CX) , public :: mksrf_flakdep_mesh = ' ' ! lake depth file name character(CX) , public :: mksrf_fwetlnd = ' ' ! inland wetlands data file name character(CX) , public :: mksrf_fwetlnd_mesh = ' ' ! inland wetlands mesh file name @@ -136,15 +136,16 @@ subroutine read_namelist_input() mksrf_fhrvtyp, & mksrf_fhrvtyp_mesh, & mksrf_fsoitex, & + mksrf_fsoitex_lookup, & mksrf_fsoitex_mesh, & - mksrf_forganic, & - mksrf_forganic_mesh, & mksrf_fsoicol, & mksrf_fsoicol_mesh, & mksrf_fvocef, & mksrf_fvocef_mesh, & - mksrf_flakwat, & - mksrf_flakwat_mesh, & + mksrf_fpctlak, & + mksrf_fpctlak_mesh, & + mksrf_flakdep, & + mksrf_flakdep_mesh, & mksrf_fwetlnd, & mksrf_fwetlnd_mesh, & mksrf_fglacier, & @@ -224,14 +225,12 @@ subroutine read_namelist_input() call mpi_bcast (mksrf_fhrvtyp_mesh, len(mksrf_fhrvtyp_mesh), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fsoitex, len(mksrf_fsoitex), MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (mksrf_fsoitex_lookup, len(mksrf_fsoitex_lookup), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fsoitex_mesh, len(mksrf_fsoitex_mesh), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fmax, len(mksrf_fmax), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fmax_mesh, len(mksrf_fmax_mesh), MPI_CHARACTER, 0, mpicom, ier) - call mpi_bcast (mksrf_forganic, len(mksrf_forganic), MPI_CHARACTER, 0, mpicom, ier) - call mpi_bcast (mksrf_forganic_mesh, len(mksrf_forganic_mesh), MPI_CHARACTER, 0, mpicom, ier) - call mpi_bcast (mksrf_fsoicol, len(mksrf_fsoicol), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fsoicol_mesh, len(mksrf_fsoicol_mesh), MPI_CHARACTER, 0, mpicom, ier) @@ -247,8 +246,10 @@ subroutine read_namelist_input() call mpi_bcast (mksrf_fgdp, len(mksrf_fgdp), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fgdp_mesh, len(mksrf_fgdp_mesh), MPI_CHARACTER, 0, mpicom, ier) - call mpi_bcast (mksrf_flakwat, len(mksrf_flakwat), MPI_CHARACTER, 0, mpicom, ier) - call mpi_bcast (mksrf_flakwat_mesh, len(mksrf_flakwat_mesh), MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (mksrf_fpctlak, len(mksrf_fpctlak), MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (mksrf_fpctlak_mesh, len(mksrf_fpctlak_mesh), MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (mksrf_flakdep, len(mksrf_flakdep), MPI_CHARACTER, 0, mpicom, ier) + call mpi_bcast (mksrf_flakdep_mesh, len(mksrf_flakdep_mesh), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fwetlnd, len(mksrf_fwetlnd), MPI_CHARACTER, 0, mpicom, ier) call mpi_bcast (mksrf_fwetlnd_mesh, len(mksrf_fwetlnd_mesh), MPI_CHARACTER, 0, mpicom, ier) @@ -341,18 +342,18 @@ subroutine write_namelist_input() write(ndiag,'(a)')' PFTs from: '//trim(mksrf_fvegtyp) write(ndiag,'(a)')' mesh for pft '//trim(mksrf_fvegtyp_mesh) write(ndiag,*) - write(ndiag,'(a)')' inland lake from: '//trim(mksrf_flakwat) - write(ndiag,'(a)')' mesh for lake water '//trim(mksrf_flakwat_mesh) + write(ndiag,'(a)')' percent lake from: '//trim(mksrf_fpctlak) + write(ndiag,'(a)')' mesh for percent lake '//trim(mksrf_fpctlak_mesh) + write(ndiag,'(a)')' lake depth from: '//trim(mksrf_flakdep) + write(ndiag,'(a)')' mesh for lake depth '//trim(mksrf_flakdep_mesh) write(ndiag,*) write(ndiag,'(a)')' inland wetland from: '//trim(mksrf_fwetlnd) write(ndiag,'(a)')' mesh for wetland '//trim(mksrf_fwetlnd_mesh) write(ndiag,*) - write(ndiag,'(a)')' soil texture from: '//trim(mksrf_fsoitex) + write(ndiag,'(a)')' soil texture mapunits from: '//trim(mksrf_fsoitex) + write(ndiag,'(a)')' soil texture (sand/clay, orgc) lookup: '//trim(mksrf_fsoitex_lookup) write(ndiag,'(a)')' mesh for soil texture '//trim(mksrf_fsoitex_mesh) write(ndiag,*) - write(ndiag,'(a)')' soil organic from: '//trim(mksrf_forganic) - write(ndiag,'(a)')' mesh for soil organic '//trim(mksrf_forganic_mesh) - write(ndiag,*) write(ndiag,'(a)')' soil color from: '//trim(mksrf_fsoicol) write(ndiag,'(a)')' mesh for soil color '//trim(mksrf_fsoicol_mesh) write(ndiag,*) diff --git a/tools/mksurfdata_esmf/src/mklanwatMod.F90 b/tools/mksurfdata_esmf/src/mklanwatMod.F90 index a01ca38515..90c8cf2b0a 100644 --- a/tools/mksurfdata_esmf/src/mklanwatMod.F90 +++ b/tools/mksurfdata_esmf/src/mklanwatMod.F90 @@ -20,7 +20,8 @@ module mklanwatMod implicit none private - public :: mklakwat ! make %lake and lake dpeth + public :: mkpctlak ! make %lake + public :: mklakdep ! make lake depth public :: mkwetlnd ! make %wetland public :: update_max_array_lake ! Update the maximum lake percent @@ -34,13 +35,11 @@ module mklanwatMod contains !=============================================================== - subroutine mklakwat(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, & - fsurdat, rc, do_depth) + subroutine mkpctlak(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, rc) ! ------------------- - ! make %lake and lake depth - ! PCT_LAKE is written out to fsurdat in mksurfdata after adjustments are made - ! LAKE_DEPTH is written out to fsurdat here + ! make %lake + ! PCT_LAKE is written to fsurdat in mksurfdata after adjustments are made ! ------------------- ! uses @@ -52,23 +51,18 @@ subroutine mklakwat(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, & type(ESMF_Mesh) , intent(in) :: mesh_o type(file_desc_t) , intent(inout) :: pioid_o real(r8) , intent(out) :: lake_o(:) ! output grid: %lake - character(len=*) , intent(in) :: fsurdat integer , intent(out) :: rc - logical, optional , intent(in) :: do_depth ! do the depth part of the subroutine ! local variables type(ESMF_Mesh) :: mesh_i type(file_desc_t) :: pioid_i - integer , allocatable :: mask_i(:) + integer , allocatable :: mask_i(:) real(r8), allocatable :: rmask_i(:) real(r8), allocatable :: lake_i(:) ! input grid: percent lake - real(r8), allocatable :: lakedepth_i(:) ! iput grid: lake depth (m) - real(r8), allocatable :: lakedepth_o(:) ! output grid: lake depth (m) integer :: ni,no,k ! indices integer :: ns_i,ns_o ! local sizes integer :: ier,rcode ! error status - real(r8), parameter :: min_valid_lakedepth = 0._r8 - character(len=*), parameter :: subname = ' mklakwat ' + character(len=*), parameter :: subname = ' mkpctlak ' !----------------------------------------------------------------------- rc = ESMF_SUCCESS @@ -77,15 +71,13 @@ subroutine mklakwat(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, & write(ndiag,*) write(ndiag,'(1x,80a1)') ('=',k=1,80) write(ndiag,*) - write(ndiag,'(a)')'Attempting to make lake data' + write(ndiag,'(a)')'Attempting to make %lake' write(ndiag,'(a)') ' Input file is '//trim(file_data_i) write(ndiag,'(a)') ' Input mesh file is '//trim(file_mesh_i) end if ! ---------------------------------------- ! Read i input data and input mesh and create route handle - ! ** ASSUME ** for now that have only 1 input data file and 1 input mesh - ! for both %lake and lakedepth ! ---------------------------------------- ! Open raw data file @@ -158,73 +150,176 @@ subroutine mklakwat(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, & deallocate (lake_i) + ! ---------------------------------------- + ! Wrap things up + ! ---------------------------------------- + + call pio_closefile(pioid_i) + + ! Release memory + if (mksrf_fdynuse == ' ') then ! ...else we will reuse it + deallocate(frac_o_mklak) + call ESMF_RouteHandleDestroy(routehandle_mklak, nogarbage = .true., rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() + end if + call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() + + if (root_task) then + write (ndiag,*) 'Successfully made %lake' + end if + call ESMF_VMLogMemInfo("At end of "//trim(subname)) + + end subroutine mkpctlak + +!=============================================================== + subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) + + ! ------------------- + ! make lake depth + ! LAKE_DEPTH is written out to fsurdat here + ! ------------------- + + ! uses + + ! input/output variables + character(len=*) , intent(in) :: file_mesh_i ! input mesh file name + character(len=*) , intent(in) :: file_data_i ! input data file name + type(ESMF_Mesh) , intent(in) :: mesh_o + type(file_desc_t) , intent(inout) :: pioid_o + character(len=*) , intent(in) :: fsurdat + integer , intent(out) :: rc + + ! local variables + type(ESMF_Mesh) :: mesh_i + type(file_desc_t) :: pioid_i + integer , allocatable :: mask_i(:) + real(r8), allocatable :: rmask_i(:) + real(r8), allocatable :: lakedepth_i(:) ! iput grid: lake depth (m) + real(r8), allocatable :: lakedepth_o(:) ! output grid: lake depth (m) + integer :: ni,no,k ! indices + integer :: ns_i,ns_o ! local sizes + integer :: ier,rcode ! error status + real(r8), parameter :: min_valid_lakedepth = 0._r8 + character(len=*), parameter :: subname = ' mklakdep ' + !----------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + if (root_task) then + write(ndiag,*) + write(ndiag,'(1x,80a1)') ('=',k=1,80) + write(ndiag,*) + write(ndiag,'(a)')'Attempting to make lake depth' + write(ndiag,'(a)') ' Input file is '//trim(file_data_i) + write(ndiag,'(a)') ' Input mesh file is '//trim(file_mesh_i) + end if + + ! ---------------------------------------- + ! Read i input data and input mesh and create route handle + ! ---------------------------------------- + + ! Open raw data file + rcode = pio_openfile(pio_iosystem, pioid_i, pio_iotype, trim(file_data_i), pio_nowrite) + + ! Read in input mesh + mesh_i = ESMF_MeshCreate(filename=trim(file_mesh_i), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After create mesh_i in "//trim(subname)) + + ! Determine ns_i + call ESMF_MeshGet(mesh_i, numOwnedElements=ns_i, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Determine ns_o + call ESMF_MeshGet(mesh_o, numOwnedElements=ns_o, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Get the landmask from the input data file and reset the mesh mask based on that + allocate(rmask_i(ns_i), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(mask_i(ns_i), stat=ier) + if (ier/=0) call shr_sys_abort() + call mkpio_get_rawdata(pioid_i, 'LANDMASK', mesh_i, rmask_i, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + do ni = 1,ns_i + if (rmask_i(ni) > 0._r8) then + mask_i(ni) = 1 + else + mask_i(ni) = 0 + end if + end do + call ESMF_MeshSet(mesh_i, elementMask=mask_i, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Create a route handle between the input and output mesh + if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak)) then + allocate(frac_o_mklak(ns_o)) + call create_routehandle_r8(mesh_i, mesh_o, routehandle_mklak, frac_o=frac_o_mklak, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) + end if + ! ---------------------------------------- ! Create lake parameter (lakedepth) ! ---------------------------------------- - if (present(do_depth)) then - if (do_depth) then - if (root_task) then - write (ndiag,*) 'Attempting to make lake depth .....' - end if - - ! lakedepth - allocate(lakedepth_i(ns_i), stat=rcode) - if (rcode/=0) call shr_sys_abort() - call mkpio_get_rawdata(pioid_i, 'LAKEDEPTH', mesh_i, lakedepth_i, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! regrid lakedepth_i to lakedepth_o - this also returns lakedepth_i to be used in the global sums below - allocate (lakedepth_o(ns_o)); lakedepth_o(:) = spval - call regrid_rawdata(mesh_i, mesh_o, routehandle_mklak, lakedepth_i, lakedepth_o, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After regrid_rawdata for lakedepth in "//trim(subname)) - do no = 1,ns_o - if (frac_o_mklak(no) == 0._r8) then - lakedepth_o(no) = 10._r8 - end if - enddo - if (min_bad(lakedepth_o, min_valid_lakedepth, 'lakedepth')) then - call shr_sys_abort() - end if - - if (fsurdat /= ' ') then - if (root_task) write(ndiag, '(a)') trim(subname)//" writing out lakedepth" - call mkfile_output(pioid_o, mesh_o, 'LAKEDEPTH', lakedepth_o, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') - call pio_syncfile(pioid_o) - end if - - ! Check global areas for lake depth - call output_diagnostics_continuous(mesh_i, mesh_o, & - lakedepth_i, lakedepth_o, "lake depth", "m", & - ndiag=ndiag, rc=rc, mask_i=mask_i, frac_o=frac_o_mklak) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (root_task) then + write (ndiag,*) 'Attempting to make lake depth .....' + end if + + ! lakedepth + allocate(lakedepth_i(ns_i), stat=rcode) + if (rcode/=0) call shr_sys_abort() + call mkpio_get_rawdata(pioid_i, 'LAKEDEPTH', mesh_i, lakedepth_i, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! regrid lakedepth_i to lakedepth_o - this also returns lakedepth_i to be used in the global sums below + allocate (lakedepth_o(ns_o)); lakedepth_o(:) = spval + call regrid_rawdata(mesh_i, mesh_o, routehandle_mklak, lakedepth_i, lakedepth_o, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After regrid_rawdata for lakedepth in "//trim(subname)) + do no = 1,ns_o + if (frac_o_mklak(no) == 0._r8) then + lakedepth_o(no) = 10._r8 end if + enddo + if (min_bad(lakedepth_o, min_valid_lakedepth, 'lakedepth')) then + call shr_sys_abort() + end if + + if (fsurdat /= ' ') then + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out lakedepth" + call mkfile_output(pioid_o, mesh_o, 'LAKEDEPTH', lakedepth_o, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') + call pio_syncfile(pioid_o) end if + ! Check global areas for lake depth + call output_diagnostics_continuous(mesh_i, mesh_o, & + lakedepth_i, lakedepth_o, "lake depth", "m", & + ndiag=ndiag, rc=rc, mask_i=mask_i, frac_o=frac_o_mklak) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ---------------------------------------- ! Wrap things up ! ---------------------------------------- - ! Close the single input data file call pio_closefile(pioid_i) ! Release memory - if (mksrf_fdynuse == ' ') then ! ...else we will reuse it - deallocate(frac_o_mklak) - call ESMF_RouteHandleDestroy(routehandle_mklak, nogarbage = .true., rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() - end if + deallocate(frac_o_mklak) + call ESMF_RouteHandleDestroy(routehandle_mklak, nogarbage = .true., rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() if (root_task) then - write (ndiag,*) 'Successfully made lake data' + write (ndiag,*) 'Successfully made lake depth' end if call ESMF_VMLogMemInfo("At end of "//trim(subname)) - end subroutine mklakwat + end subroutine mklakdep !=============================================================== subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) @@ -252,7 +347,7 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) integer :: ni,no,k ! indices integer :: ns_i,ns_o ! local sizes integer :: ier,rcode ! error status - character(len=*), parameter :: subname = ' mklakwat ' + character(len=*), parameter :: subname = ' mkwetlnd ' !----------------------------------------------------------------------- rc = ESMF_SUCCESS diff --git a/tools/mksurfdata_esmf/src/mkorganicMod.F90 b/tools/mksurfdata_esmf/src/mkorganicMod.F90 deleted file mode 100644 index fce5dc0e92..0000000000 --- a/tools/mksurfdata_esmf/src/mkorganicMod.F90 +++ /dev/null @@ -1,195 +0,0 @@ -module mkorganicMod - - ! make organic matter dataset - - use ESMF - use pio - use shr_kind_mod , only : r8 => shr_kind_r8, r4 => shr_kind_r4 - use shr_sys_mod , only : shr_sys_abort - use mkpioMod , only : mkpio_get_rawdata, mkpio_get_dimlengths - use mkpioMod , only : pio_iotype, pio_ioformat, pio_iosystem - use mkesmfMod , only : regrid_rawdata, create_routehandle_r8 - use mkutilsMod , only : chkerr - use mkvarctl , only : root_task, ndiag, spval - use mkvarpar , only : nlevsoi - use mkfileMod , only : mkfile_output - - implicit none - private - - public :: mkorganic ! Set organic soil - - character(len=*) , parameter :: u_FILE_u = & - __FILE__ - -!================================================================================= -contains -!================================================================================= - - subroutine mkorganic(file_mesh_i, file_data_i, mesh_o, pctlnd_pft_o, lat_o, pioid_o, rc) - - ! input/output variables - character(len=*) , intent(in) :: file_mesh_i ! input mesh file name - character(len=*) , intent(in) :: file_data_i ! input data file name - type(ESMF_Mesh) , intent(in) :: mesh_o - real(r8) , intent(in) :: pctlnd_pft_o(:) - real(r8) , intent(in) :: lat_o(:) - type(file_desc_t) , intent(inout) :: pioid_o - integer , intent(out) :: rc - - ! local variables - type(ESMF_RouteHandle) :: routehandle - type(ESMF_Mesh) :: mesh_i - type(file_desc_t) :: pioid_i - type(var_desc_t) :: pio_varid - integer :: ni,no - integer :: ns_i, ns_o - integer :: nlay - integer :: n,l,k ! indices - integer :: rcode, ier ! error status - integer :: ndims - integer , allocatable :: dimlengths(:) - integer , allocatable :: mask_i(:) - real(r8), allocatable :: frac_i(:) - real(r8), allocatable :: frac_o(:) - real(r8), allocatable :: area_i(:) - real(r8), allocatable :: area_o(:) - real(r8), allocatable :: data_i(:,:) - real(r8), allocatable :: data_o(:,:) - real(r8), allocatable :: organic_o(:,:) ! output grid: %lake - character(len=*), parameter :: subname = 'mkorganic' - !----------------------------------------------------------------------- - - rc = ESMF_SUCCESS - - if (root_task) then - write(ndiag,*) - write(ndiag,'(1x,80a1)') ('=',k=1,80) - write(ndiag,*) - write(ndiag,'(a)')'Attempting to make organic mater dataset .....' - write(ndiag,'(a)') ' Input file is '//trim(file_data_i) - write(ndiag,'(a)') ' Input mesh file is '//trim(file_mesh_i) - end if - - ! Open input data file - rcode = pio_openfile(pio_iosystem, pioid_i, pio_iotype, trim(file_data_i), pio_nowrite) - call ESMF_VMLogMemInfo("After pio_openfile "//trim(file_data_i)) - - ! Get dimensions of raw data. - ! - raw data is dimensions (lon,lat,lev) - ! - input read from pio has dimensions(n,lev) - ! - esmf field dataptr has dimensions (lev,n) - allocate(dimlengths(3)) - call mkpio_get_dimlengths(pioid_i, 'ORGANIC', ndims, dimlengths) - nlay = dimlengths(3) - if (nlay /= nlevsoi) then - write(6,*)'nlay, nlevsoi= ',nlay,nlevsoi,' do not match' - call shr_sys_abort() - end if - deallocate(dimlengths) - - ! Read in input mesh - call ESMF_VMLogMemInfo("Before create mesh_i in "//trim(subname)) - mesh_i = ESMF_MeshCreate(filename=trim(file_mesh_i), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After create mesh_i in "//trim(subname)) - - ! Determine ns_i - call ESMF_MeshGet(mesh_i, numOwnedElements=ns_i, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! Determine ns_o and allocate output data - call ESMF_MeshGet(mesh_o, numOwnedElements=ns_o, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - allocate ( organic_o(ns_o,nlevsoi)); organic_o(:,:) = spval - - ! Get the landmask from the file and reset the mesh mask based on that - allocate(frac_i(ns_i), stat=ier) - if (ier/=0) call shr_sys_abort() - allocate(mask_i(ns_i), stat=ier) - if (ier/=0) call shr_sys_abort() - call mkpio_get_rawdata(pioid_i, 'LANDMASK', mesh_i, frac_i, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - do ni = 1,ns_i - if (frac_i(ni) > 0._r8) then - mask_i(ni) = 1 - else - mask_i(ni) = 0 - end if - end do - call ESMF_MeshSet(mesh_i, elementMask=mask_i, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - - ! Create a route handle between the input and output mesh and get frac_o - allocate(frac_o(ns_o),stat=ier) - if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) - - ! Read in input data - ! - levels are the innermost dimension for esmf fields - ! - levels are the outermost dimension in pio reads - ! Input data is read into (ns_i,nlay) array and then transferred to data_i(nlay,ns_i) - allocate(data_i(nlay,ns_i),stat=ier) - if (ier/=0) call shr_sys_abort() - call mkpio_get_rawdata(pioid_i, 'ORGANIC', mesh_i, data_i, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After mkpio_getrawdata in "//trim(subname)) - - ! Regrid data_i to data_o - allocate(data_o(nlay,ns_o),stat=ier) - if (ier/=0) call shr_sys_abort() - call regrid_rawdata(mesh_i, mesh_o, routehandle, data_i, data_o, 1, nlay, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After regrid organic_i in "//trim(subname)) - - ! Set organic_o - do l = 1,nlevsoi - do no = 1,ns_o - organic_o(no,l) = data_o(l,no) - end do - end do - do l = 1,nlevsoi - do no = 1,ns_o - if ((organic_o(no,l)) > 130.000001_r8) then - write (6,*) trim(subname)//' error: organic = ',organic_o(no,l), & - ' greater than 130.000001 for n,lev ',no,l - call shr_sys_abort() - end if - enddo - end do - do no = 1,ns_o - ! If have very small values of pctlnd - set organic_o to zero - if (pctlnd_pft_o(no) < 1.e-6_r8) then - organic_o(no,:) = 0._r8 - end if - ! If have pole points on grid - set south pole to glacier and north pole is not land - if (abs((lat_o(no) - 90._r8)) < 1.e-6_r8) then - organic_o(no,:) = 0._r8 - end if - end do - - ! Write out the output data - if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil organic matter density" - call mkfile_output(pioid_o, mesh_o, 'ORGANIC', organic_o, lev1name='nlevsoi', rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') - - ! Close the file - call pio_closefile(pioid_i) - call ESMF_VMLogMemInfo("After pio_closefile in "//trim(subname)) - - ! Release memory - call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() - call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() - call ESMF_VMLogMemInfo("After destroy operations in "//trim(subname)) - - if (root_task) then - write (ndiag,'(a)') 'Successfully made organic matter ' - end if - - end subroutine mkorganic - -end module mkorganicMod diff --git a/tools/mksurfdata_esmf/src/mksoiltexMod.F90 b/tools/mksurfdata_esmf/src/mksoiltexMod.F90 index 04e68c60d0..295217230f 100644 --- a/tools/mksurfdata_esmf/src/mksoiltexMod.F90 +++ b/tools/mksurfdata_esmf/src/mksoiltexMod.F90 @@ -33,25 +33,26 @@ module mksoiltexMod contains !================================================================================= - subroutine mksoiltex(file_mesh_i, file_data_i, mesh_o, pioid_o, pctlnd_pft_o, rc) + subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, pctlnd_pft_o, rc) ! - ! make %sand and %clay from IGBP soil data, which includes - ! igbp soil 'mapunits' and their corresponding textures + ! make %sand, %clay, organic carbon content, coarse fragments, bulk density, + ! and pH measured in H2O ! ! input/output variables - character(len=*) , intent(in) :: file_mesh_i ! input mesh file name - character(len=*) , intent(in) :: file_data_i ! input data file name - type(ESMF_Mesh) , intent(in) :: mesh_o ! output mesh + character(len=*) , intent(in) :: file_mesh_i ! input mesh/grid file name + character(len=*) , intent(in) :: file_mapunit_i ! input mapunit file name + character(len=*) , intent(in) :: file_lookup_i ! input data file name + type(ESMF_Mesh) , intent(in) :: mesh_o ! output mesh type(file_desc_t) , intent(inout) :: pioid_o real(r8) , intent(in) :: pctlnd_pft_o(:) ! PFT data: % of gridcell for PFTs integer , intent(out) :: rc ! local variables type(ESMF_RouteHandle) :: routehandle + type(ESMF_Grid) :: grid_i type(ESMF_Mesh) :: mesh_i type(ESMF_Field) :: field_i type(ESMF_Field) :: field_o - type(ESMF_Field) :: field_dstfrac type(file_desc_t) :: pioid_i type(var_desc_t) :: pio_varid integer :: pio_vartype @@ -59,27 +60,42 @@ subroutine mksoiltex(file_mesh_i, file_data_i, mesh_o, pioid_o, pctlnd_pft_o, rc integer :: ni,no integer :: ns_i, ns_o integer :: k,l,m,n - integer :: nlay ! number of soil layers + integer :: nlay ! number of soil layers + integer :: n_scid integer , allocatable :: mask_i(:) - real(r4), allocatable :: rmask_i(:) - real(r8), allocatable :: frac_o(:) - real(r4), allocatable :: sand_i(:,:) ! input grid: percent sand - real(r4), allocatable :: clay_i(:,:) ! input grid: percent clay - real(r4), allocatable :: mapunit_i(:) ! input grid: igbp soil mapunits - real(r8), allocatable :: sand_o(:,:) ! % sand (output grid) - real(r8), allocatable :: clay_o(:,:) ! % clay (output grid) - integer , allocatable :: mapunit_o(:) real(r4), pointer :: dataptr(:) - real(r8), pointer :: dataptr_r8(:) - real(r4) :: sumtex - integer :: mapunit ! temporary igbp soil mapunit - integer, allocatable :: soil_i(:) - integer, allocatable :: soil_o(:) - integer :: rcode, ier ! error status + integer :: mapunit ! temporary igbp soil mapunit + integer, allocatable :: sand_i(:,:,:) ! input grid: percent sand + integer, allocatable :: clay_i(:,:,:) ! input grid: percent clay + integer, allocatable :: cfrag_i(:,:,:) ! input grid: coarse fragments (vol% > 2 mm) + real(r4), allocatable :: bulk_i(:,:,:) ! input grid: bulk density (g cm-3) + real(r4), allocatable :: orgc_i(:,:,:) ! input grid: organic carbon content (gC kg-1) + real(r4), allocatable :: phaq_i(:,:,:) ! input grid: soil pH measured in H2O (unitless) + real(r4), allocatable :: sand_o(:,:) ! output grid: % sand + real(r4), allocatable :: clay_o(:,:) ! output grid: % clay + real(r4), allocatable :: orgc_o(:,:) ! output grid: organic carbon content (gC kg-1) + real(r4), allocatable :: cfrag_o(:,:) ! output grid: coarse fragments (vol% > 2 mm) + real(r4), allocatable :: bulk_o(:,:) ! output grid: bulk density (g cm-3) + real(r4), allocatable :: phaq_o(:,:) ! output grid: soil pH measured in H2O (unitless) + real(r4), allocatable :: organic_o(:,:) ! output grid: organic matter (kg m-3) + integer :: n_mapunits + integer :: lookup_index + integer :: SCID + real(r4), allocatable :: mapunit_i(:) ! input grid: igbp soil mapunits + integer , allocatable :: mapunit_o(:) ! output grid: igbp soil mapunits + integer , allocatable :: MapUnits(:) + integer , allocatable :: mapunit_lookup(:) + type(var_desc_t) :: pio_varid_sand + type(var_desc_t) :: pio_varid_clay + type(var_desc_t) :: pio_varid_orgc + type(var_desc_t) :: pio_varid_cfrag + type(var_desc_t) :: pio_varid_bulk + type(var_desc_t) :: pio_varid_phaq + type(var_desc_t) :: pio_varid_organic + integer :: starts(3) ! starting indices for reading lookup table + integer :: counts(3) ! dimension counts for reading lookup table integer :: srcTermProcessing_Value = 0 - integer, parameter :: nlsm=4 ! number of soil textures - character(len=38) :: soil(0:nlsm) ! name of each soil texture - character(len=38) :: typ ! soil texture based on ... + integer :: rcode, ier ! error status character(len=*), parameter :: subname = 'mksoiltex' !----------------------------------------------------------------------- @@ -87,9 +103,10 @@ subroutine mksoiltex(file_mesh_i, file_data_i, mesh_o, pioid_o, pctlnd_pft_o, rc write(ndiag,*) write(ndiag,'(1x,80a1)') ('=',k=1,80) write(ndiag,*) - write(ndiag,'(a)') 'Attempting to make %sand and %clay .....' - write(ndiag,'(a)') ' Input file is '//trim(file_data_i) - write(ndiag,'(a)') ' Input mesh file is '//trim(file_mesh_i) + write(ndiag,'(a)') 'Attempting to make %sand, %clay, orgc, cfrag, bulk, phaq .....' + write(ndiag,'(a)') ' Input mapunit file is '//trim(file_mapunit_i) + write(ndiag,'(a)') ' Input lookup table file is '//trim(file_lookup_i) + write(ndiag,'(a)') ' Input mesh/grid file is '//trim(file_mesh_i) end if ! Determine ns_o and allocate output data @@ -98,237 +115,436 @@ subroutine mksoiltex(file_mesh_i, file_data_i, mesh_o, pioid_o, pctlnd_pft_o, rc allocate(mapunit_o(ns_o)) ; mapunit_o(:) = 0 allocate(sand_o(ns_o,nlevsoi)) ; sand_o(:,:) = spval allocate(clay_o(ns_o,nlevsoi)) ; clay_o(:,:) = spval + allocate(orgc_o(ns_o,nlevsoi)) ; orgc_o(:,:) = spval + allocate(cfrag_o(ns_o,nlevsoi)) ; cfrag_o(:,:) = spval + allocate(bulk_o(ns_o,nlevsoi)) ; bulk_o(:,:) = spval + allocate(phaq_o(ns_o,nlevsoi)) ; phaq_o(:,:) = spval + allocate(organic_o(ns_o,nlevsoi)) ; organic_o(:,:) = spval + + !--------------------------------- + ! Determine mapunits on output grid + !--------------------------------- + + ! Determine input mesh + if (trim(file_mesh_i) == trim(file_mapunit_i)) then + ! input format is GRIDSPEC and read in grid and then create mesh + if (root_task) write(ndiag,*)"reading grid_i and then creating mesh_i in "//trim(subname) + call ESMF_VMLogMemInfo("Before create read in grid_i in "//trim(subname)) + grid_i = ESMF_GridCreate(filename=trim(file_mesh_i), & + fileformat=ESMF_FILEFORMAT_GRIDSPEC, addCornerStagger=.true., addmask=.true., varname='MU', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("Before create mesh_i in "//trim(subname)) + mesh_i = esmf_meshcreate(grid_i, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After create mesh_i in "//trim(subname)) + else + ! Read in mesh directly + if (root_task) write(ndiag,*)"reading mesh_i directly in "//trim(subname) + mesh_i = ESMF_MeshCreate(filename=trim(file_mesh_i), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! Open input data file - call ESMF_VMLogMemInfo("Before pio_openfile for "//trim(file_data_i)) - rcode = pio_openfile(pio_iosystem, pioid_i, pio_iotype, trim(file_data_i), pio_nowrite) - - ! Read in input mesh - call ESMF_VMLogMemInfo("Before create mesh_i in "//trim(subname)) - mesh_i = ESMF_MeshCreate(filename=trim(file_mesh_i), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After create mesh_i in "//trim(subname)) - - ! Determine ns_i + ! Determine ns_i (use the distgrid to the number of elements) call ESMF_MeshGet(mesh_i, numOwnedElements=ns_i, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Get the landmask from the file and reset the mesh mask based on that - allocate(rmask_i(ns_i), stat=ier) + ! Read in mapunit data + if (root_task) write(ndiag,*)"Reading in mapunit data in "//trim(subname) + call ESMF_VMLogMemInfo("Before pio_openfile for "//trim(file_mapunit_i)) + rcode = pio_openfile(pio_iosystem, pioid_i, pio_iotype, trim(file_mapunit_i), pio_nowrite) + allocate(mapunit_i(ns_i), stat=ier) if (ier/=0) call shr_sys_abort() + call mkpio_get_rawdata(pioid_i, 'MU', mesh_i, mapunit_i, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After mkpio_getrawdata in "//trim(subname)) + call pio_closefile(pioid_i) + + ! Set mesh mask to zero where the mapunit values are 0 + if (root_task) write(ndiag,*)"Setting mask in mesh where mapunit data is 0 "//trim(subname) allocate(mask_i(ns_i), stat=ier) if (ier/=0) call shr_sys_abort() - call mkpio_get_rawdata(pioid_i, 'LANDMASK', mesh_i, rmask_i, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return do ni = 1,ns_i - if (rmask_i(ni) > 0._r4) then - mask_i(ni) = 1 - else + if (mapunit_i(ni) == 0.) then mask_i(ni) = 0 end if end do call ESMF_MeshSet(mesh_i, elementMask=mask_i, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - ! Read in mapunit data - allocate(mapunit_i(ns_i), stat=ier) - if (ier/=0) call shr_sys_abort() - call mkpio_get_rawdata(pioid_i, 'MAPUNITS', mesh_i, mapunit_i, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After mkpio_getrawdata in "//trim(subname)) - - ! Set mapunit values to zero where the input mask is 0 - do ni = 1,ns_i - mapunit_i(ni) = mapunit_i(ni) * mask_i(ni) - end do - - ! Determine mapunit_value_max (set it as a module variable so that it can be - ! accessible to gen_dominant_mapunit) - rcode = pio_inq_dimid (pioid_i, 'max_value_mapunit', dimid) - rcode = pio_inq_dimlen (pioid_i, dimid, mapunit_value_max) - ! Create ESMF fields that will be used below field_i = ESMF_FieldCreate(mesh_i, ESMF_TYPEKIND_R4, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return field_o = ESMF_FieldCreate(mesh_o, ESMF_TYPEKIND_R4, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - field_dstfrac = ESMF_FieldCreate(mesh_o, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle + if (root_task) write(ndiag,*)" before route handle creation "//trim(subname) call ESMF_FieldRegridStore(field_i, field_o, routehandle=routehandle, & regridmethod=ESMF_REGRIDMETHOD_CONSERVE, srcTermProcessing=srcTermProcessing_Value, & - ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & - dstFracField= field_dstfrac, rc=rc) + ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After regridstore in "//trim(subname)) - - ! Determin frac_o - call ESMF_FieldGet(field_dstfrac, farrayptr=dataptr_r8, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - allocate(frac_o(ns_o)) - frac_o(:) = dataptr_r8(:) + if (root_task) write(ndiag,*)" after route handle creation "//trim(subname) ! Create a dynamic mask object ! The dynamic mask object further holds a pointer to the routine that will be called in order to ! handle dynamically masked elements - in this case its DynMaskProc (see below) + if (root_task) write(ndiag,*)" before call to dynamic mask set creation "//trim(subname) call ESMF_DynamicMaskSetR4R8R4(dynamicMask, dynamicMaskRoutine=get_dominant_mapunit, & handleAllElements=.true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (root_task) write(ndiag,*)" after call to dynamic mask set creation "//trim(subname) - ! Determine dominant soil color in the field regrid call below + ! Determine values in field_i call ESMF_FieldGet(field_i, farrayptr=dataptr, rc=rc) dataptr(:) = real(mapunit_i(:), kind=r4) + + ! Determine mapunit_value_max (set it as a module variable so that it can be + ! accessible to gen_dominant_mapunit) - this is needed in the dynamic mask routine + mapunit_value_max = maxval(dataptr) + + ! Determine values in field_o call ESMF_FieldGet(field_o, farrayptr=dataptr, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return dataptr(:) = 0._r4 + ! Determine mapunit_o call ESMF_FieldRegrid(field_i, field_o, routehandle=routehandle, dynamicMask=dynamicMask, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldGet(field_o, farrayptr=dataptr, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return do no = 1,ns_o mapunit_o(no) = int(dataptr(no)) end do - ! Get dimensions from input file and allocate memory for sand_i and clay_i - rcode = pio_inq_dimid (pioid_i, 'number_of_layers', dimid) + do no = 1,ns_o + if (mapunit_o(no) > mapunit_value_max) then + write(6,*)'mapunit_o is out of bounds ',mapunit_o(no) + ! call shr_sys_abort("mapunit_o is out of bounds") + end if + end do + + ! Write out mapunit_o + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out mapunits " + call mkfile_output(pioid_o, mesh_o, 'mapunits', mapunit_o, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output for mapunits') + call pio_syncfile(pioid_o) + + !--------------------------------- + ! Determine %sand, %clay, orgc, cfrag, bulk, phaq on output grid using + ! mapunits + !--------------------------------- + + rcode = pio_openfile(pio_iosystem, pioid_i, pio_iotype, trim(file_lookup_i), pio_nowrite) + + rcode = pio_inq_dimid (pioid_i, 'MapUnit', dimid) + rcode = pio_inq_dimlen (pioid_i, dimid, n_mapunits) + + rcode = pio_inq_dimid (pioid_i, 'soil_layer', dimid) rcode = pio_inq_dimlen (pioid_i, dimid, nlay) - allocate(sand_i(mapunit_value_max,nlay), stat=ier) - if (ier/=0) call shr_sys_abort() - allocate(clay_i(mapunit_value_max,nlay), stat=ier) + + rcode = pio_inq_dimid (pioid_i, 'SCID', dimid) + rcode = pio_inq_dimlen (pioid_i, dimid, n_scid) + + ! Read In MapUnits from the input file + allocate(MapUnits(n_mapunits), stat=ier) if (ier/=0) call shr_sys_abort() + rcode = pio_inq_varid(pioid_i, 'MapUnit', pio_varid) + rcode = pio_get_var(pioid_i, pio_varid, MapUnits) + + ! Determine the mapunit lookup index from the value of the MapUnit variable + mapunit_value_max = maxval(MapUnits) + allocate(mapunit_lookup(mapunit_value_max)) + mapunit_lookup(:) = -999 + do n = 1,size(MapUnits) + mapunit_lookup(MapUnits(n)) = n + end do - ! read in sand_i and clay_i (they will read in total on all processors) - rcode = pio_inq_varid(pioid_i, 'PCT_SAND', pio_varid) - rcode = pio_get_var(pioid_i, pio_varid, sand_i) - rcode = pio_inq_varid(pioid_i, 'PCT_CLAY', pio_varid) - rcode = pio_get_var(pioid_i, pio_varid, clay_i) + allocate(sand_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(clay_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(orgc_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(cfrag_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(bulk_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() + allocate(phaq_i(nlay,n_scid,n_mapunits), stat=ier) + if (ier/=0) call shr_sys_abort() - ! Set soil texture as follows: - ! a. Use dominant igbp soil mapunit based on area of overlap unless 'no data' is dominant - ! b. If this has no data, use loam for soil texture + ! Get dimensions from input file and allocate memory for sand_i, clay_i, + ! organic carbon content, coarse fragments, bulk density, pH measured in H2O + rcode = pio_inq_varid(pioid_i, 'PCT_SAND', pio_varid_sand) + rcode = pio_inq_varid(pioid_i, 'PCT_CLAY', pio_varid_clay) + rcode = pio_inq_varid(pioid_i, 'ORGC', pio_varid_orgc) + rcode = pio_inq_varid(pioid_i, 'CFRAG', pio_varid_cfrag) + rcode = pio_inq_varid(pioid_i, 'BULK', pio_varid_bulk) + rcode = pio_inq_varid(pioid_i, 'PHAQ', pio_varid_phaq) + + rcode = pio_get_var(pioid_i, pio_varid_sand, sand_i) + rcode = pio_get_var(pioid_i, pio_varid_clay, clay_i) + rcode = pio_get_var(pioid_i, pio_varid_orgc, orgc_i) + rcode = pio_get_var(pioid_i, pio_varid_cfrag, cfrag_i) + rcode = pio_get_var(pioid_i, pio_varid_bulk, bulk_i) + rcode = pio_get_var(pioid_i, pio_varid_phaq, phaq_i) do no = 1,ns_o - if (mapunit_o(no) > 0) then - ! valid value is obtained - if (mapunit_o(no) > mapunit_value_max) then - write(6,*)'mapunit_o is out of bounds ',mapunit_o(no) - ! call shr_sys_abort("mapunit_o is out of bounds") + + if (pctlnd_pft_o(no) < 1.e-6_r8 .or. mapunit_o(no) == 0) then + + ! Set sand and clay to loam if pctlnd_pft is < 1.e-6 or mapunit is 0 + sand_o(no,:) = 43._r4 + clay_o(no,:) = 18._r4 + orgc_o(no,:) = 0._r4 + cfrag_o(no,:) = 0._r4 + bulk_o(no,:) = 1.5_r4 ! TODO Ok as a fill value? + phaq_o(no,:) = 7._r4 + organic_o(no,:) = 0._r4 ! TODO Rm bef merging PR #1732 + + else + + ! Determine lookup_index + lookup_index = mapunit_lookup(mapunit_o(no)) + + ! Determine top soil layer sand_o + ! If less than 0 search within the SCID array for the first index + ! that gives a value greater than or equal to 0 + ! Then determine the other soil layers + sand_o(no,1) = float(sand_i(1,1,lookup_index)) + if (sand_o(no,1) < 0._r4) then + do l = 2,n_scid + if (float(sand_i(1,l,lookup_index)) >= 0._r4) then + sand_o(no,1) = float(sand_i(1,l,lookup_index)) + exit + end if + end do end if - do l = 1, nlay - sand_o(no,l) = sand_i(mapunit_o(no),l) - clay_o(no,l) = clay_i(mapunit_o(no),l) + if (sand_o(no,1) < 0._r4) then + if (int(sand_o(no,1)) == -4) then + write(6,'(a,i8)')'WARNING: changing sand_o from -4 to 99% at no = ',no + sand_o(no,:) = 99._r4 + else + write(6,'(a,i8,a,i8)')'WARNING: changing sand_o from ',int(sand_o(no,1)),' to 43 at no = ',no + sand_o(no,:) = 43._r4 + end if + end if + do l = 2,nlay + sand_o(no,l) = float(sand_i(l,1,lookup_index)) + if (sand_o(no,l) < 0._r4) then + sand_o(no,l) = sand_o(no,l-1) + end if end do - else - ! use loam - do l = 1, nlay - sand_o(no,l) = 43. - clay_o(no,l) = 18. + + ! Same algorithm for clay_o as for sand_o + clay_o(no,1) = float(clay_i(1,1,lookup_index)) + if (clay_o(no,1) < 0._r4) then + do l = 2,n_scid + if (float(clay_i(1,l,lookup_index)) >= 0._r4) then + clay_o(no,1) = float(clay_i(1,l,lookup_index)) + exit + end if + end do + end if + if (clay_o(no,1) < 0._r4) then + if (int(clay_o(no,1)) == -4) then + write(6,'(a,i8)')'WARNING: changing clay_o from -4 to 1% at no = ',no + clay_o(no,:) = 1._r4 + else + write(6,'(a,i8,a,i8)')'WARNING: changing clay_o from ',int(clay_o(no,1)),' to 18 at no = ',no + clay_o(no,:) = 18._r4 + end if + end if + if (clay_o(no,1) < 0._r4) then + write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index + call shr_sys_abort('could not find a value >= 0 for clay_i') + end if + do l = 2,nlay + clay_o(no,l) = float(clay_i(l,1,lookup_index)) + if (clay_o(no,l) < 0._r4) then + clay_o(no,l) = clay_o(no,l-1) + end if end do - end if - end do - ! Define the model surface types: 0:4 - soil(0) = 'no soil: ocean, glacier, lake, no data' - soil(1) = 'clays ' - soil(2) = 'sands ' - soil(3) = 'loams ' - soil(4) = 'silts ' - - ! input grid: determine soil_i for global area comparison - allocate(soil_i(ns_i)) - do l = 1, nlay - do ni = 1,ns_i - mapunit = nint(mapunit_i(ni)) - if (mapunit==0) then - typ = 'no soil: ocean, glacier, lake, no data' - else if (clay_i(mapunit,l) >= 40.) then - typ = 'clays' - else if (sand_i(mapunit,l) >= 50.) then - typ = 'sands' - else if (clay_i(mapunit,l)+sand_i(mapunit,l) < 50.) then - if (rmask_i(ni) /= 0.) then - typ = 'silts' + ! Same algorithm for orgc_o as for sand_o + ! organic_o OPTION 2 (commented out) + ! Calculate from multiple input variables to get the output variable + orgc_o(no,1) = orgc_i(1,1,lookup_index) +! organic_o(no,1) = orgc_i(1,1,lookup_index) * bulk_i(1,1,lookup_index) * float(100 - cfrag_i(1,1,lookup_index)) * 0.01_r4 / 0.58_r4 + if (orgc_o(no,1) < 0._r4) then + do l = 2,n_scid + if (orgc_i(1,l,lookup_index) >= 0._r4) then + orgc_o(no,1) = orgc_i(1,l,lookup_index) +! organic_o(no,1) = orgc_i(1,l,lookup_index) * bulk_i(1,l,lookup_index) * float(100 - cfrag_i(1,l,lookup_index)) * 0.01_r4 / 0.58_r4 + exit + end if + end do + end if + if (orgc_o(no,1) < 0._r4) then + if (int(orgc_o(no,1)) == -4) then ! sand dunes + write(6,'(a,i8)')'WARNING: changing orgc_o from -4 to 1 at no = ', no + orgc_o(no,:) = 1._r4 +! write(6,'(a,i8)')'WARNING: changing organic_o from -4 to 1 at no = ', no +! organic_o(no,:) = 1._r4 else - typ = 'no soil: ocean, glacier, lake, no data' + write(6,'(a,i8,a,i8)')'WARNING: changing orgc_o from ',int(orgc_o(no,1)),' to 0 at no = ',no + orgc_o(no,:) = 0._r4 +! write(6,'(a,i8,a,i8)')'WARNING: changing organic_o from ',int(organic_o(no,1)),' to 0 at no = ', no +! organic_o(no,:) = 0._r4 end if - else - typ = 'loams' end if - do m = 0, nlsm - if (typ == soil(m)) go to 101 + if (orgc_o(no,1) < 0._r4) then + write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index + call shr_sys_abort('could not find a value >= 0 for orgc_i') + end if + do l = 2,nlay + orgc_o(no,l) = orgc_i(l,1,lookup_index) +! organic_o(no,l) = orgc_i(l,1,lookup_index) * bulk_i(l,1,lookup_index) * float(100 - cfrag_i(l,1,lookup_index)) * 0.01_r4 / 0.58_r4 + if (orgc_o(no,l) < 0._r4) then + orgc_o(no,l) = orgc_o(no,l-1) +! organic_o(no,l) = organic_o(no,l-1) + end if end do - write (6,*) 'MKSOILTEX error: sand = ',sand_i(mapunit,l), & - ' clay = ',clay_i(mapunit,l),' not assigned to soil type for input grid lon,lat,layer = ',ni,l - call shr_sys_abort() -101 continue - soil_i(ni) = m ! allocate memory for this - end do - end do - ! output grid: determine soil_o for global area comparison - allocate(soil_o(ns_o)) - do l = 1, nlay - do no = 1,ns_o - if (clay_o(no,l)==0. .and. sand_o(no,l)==0.) then - typ = 'no soil: ocean, glacier, lake, no data' - else if (clay_o(no,l) >= 40.) then - typ = 'clays' - else if (sand_o(no,l) >= 50.) then - typ = 'sands' - else if (clay_o(no,l)+sand_o(no,l) < 50.) then - typ = 'silts' - else - typ = 'loams' + ! Same algorithm for cfrag_o as for sand_o + cfrag_o(no,1) = float(cfrag_i(1,1,lookup_index)) + if (cfrag_o(no,1) < 0._r4) then + do l = 2,n_scid + if (float(cfrag_i(1,l,lookup_index)) >= 0._r4) then + cfrag_o(no,1) = float(cfrag_i(1,l,lookup_index)) + exit + end if + end do end if - do m = 0, nlsm - if (typ == soil(m)) go to 102 + if (cfrag_o(no,1) < 0._r4) then + if (int(cfrag_o(no,1)) == -4) then ! sand dunes + write(6,'(a,i8)')'WARNING: changing cfrag_o from -4 to 1 at no = ',no + cfrag_o(no,:) = 1._r4 + else + write(6,'(a,i8,a,i8)')'WARNING: changing cfrag_o from ',int(cfrag_o(no,1)),' to 0 at no = ',no + cfrag_o(no,:) = 0._r4 + end if + end if + if (cfrag_o(no,1) < 0._r4) then + write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index + call shr_sys_abort('could not find a value >= 0 for cfrag_i') + end if + do l = 2,nlay + cfrag_o(no,l) = float(cfrag_i(l,1,lookup_index)) + if (cfrag_o(no,l) < 0._r4) then + cfrag_o(no,l) = cfrag_o(no,l-1) + end if end do - write (6,*) 'MKSOILTEX error: sand = ',sand_o(no,l), & - ' clay = ',clay_o(no,l), ' not assigned to soil type for output grid lon,lat,layer = ',no,l - call shr_sys_abort() -102 continue - soil_o(no) = m - end do - end do - ! Compare global areas - if (root_task) then - write (ndiag,'(a)') 'The following table of soil texture classes is for comparison only.' - write (ndiag,'(a)') 'The actual data is continuous %sand, %silt and %clay not textural classes' - end if - call output_diagnostics_index(mesh_i, mesh_o, mask_i, frac_o, & - 0, nlsm, soil_i, soil_o, 'soil texture class', ndiag, rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + ! Same algorithm for bulk_o as for sand_o + bulk_o(no,1) = bulk_i(1,1,lookup_index) + if (bulk_o(no,1) < 0._r4) then + do l = 2,n_scid + if (bulk_i(1,l,lookup_index) >= 0._r4) then + bulk_o(no,1) = bulk_i(1,l,lookup_index) + exit + end if + end do + end if + if (bulk_o(no,1) < 0._r4) then + if (int(bulk_o(no,1)) == -4) then ! sand dunes + write(6,'(a,i8)')'WARNING: changing bulk_o from -4 to 1 at no = ',no + bulk_o(no,:) = 1.5_r4 ! TODO Ok for sand dunes? + else + write(6,'(a,i8,a,i8)')'WARNING: changing bulk_o from ',int(bulk_o(no,1)),' to 0 at no = ',no + bulk_o(no,:) = 1.5_r4 ! TODO Ok for -7? + end if + end if + if (bulk_o(no,1) < 0._r4) then + write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index + call shr_sys_abort('could not find a value >= 0 for bulk_i') + end if + do l = 2,nlay + bulk_o(no,l) = bulk_i(l,1,lookup_index) + if (bulk_o(no,l) < 0._r4) then + bulk_o(no,l) = bulk_o(no,l-1) + end if + end do + + ! Same algorithm for phaq_o as for sand_o + phaq_o(no,1) = phaq_i(1,1,lookup_index) + if (phaq_o(no,1) < 0._r4) then + do l = 2,n_scid + if (phaq_i(1,l,lookup_index) >= 0._r4) then + phaq_o(no,1) = phaq_i(1,l,lookup_index) + exit + end if + end do + end if + if (phaq_o(no,1) < 0._r4) then + if (int(phaq_o(no,1)) == -4) then ! sand dunes + write(6,'(a,i8)')'WARNING: changing phaq_o from -4 to 1 at no = ',no + phaq_o(no,:) = 7._r4 + else + write(6,'(a,i8,a,i8)')'WARNING: changing phaq_o from ',int(phaq_o(no,1)),' to 0 at no = ',no + phaq_o(no,:) = 7._r4 + end if + end if + if (phaq_o(no,1) < 0._r4) then + write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index + call shr_sys_abort('could not find a value >= 0 for phaq_i') + end if + do l = 2,nlay + phaq_o(no,l) = phaq_i(l,1,lookup_index) + if (phaq_o(no,l) < 0._r4) then + phaq_o(no,l) = phaq_o(no,l-1) + end if + end do + + ! --------------------------------------------------------------- + ! organic_o OPTION 1, as we plan to calculate organic in the CTSM + ! --------------------------------------------------------------- + ! Calculate organic from orgc_o, cfrag_o, and bulk_o, i.e. after + ! these terms have been regridded. The plan is to move this + ! calculation step from here to the CTSM. This approach keeps + ! ORGC the same in fsurdat as in the raw data. + ! Alternative approach considered but not selected: Regrid organic_i + ! (calculated from orgc_i, cfrag_i, and bulk_i) to organic_o. This + ! approach first calculates organic_i and then regrids to organic_o + ! rather than regridding all the terms first and then calculating + ! organic_o. Commented out code above belongs to that option. + do l = 1, nlay + organic_o(no,l) = orgc_o(no,l) * bulk_o(no,l) * & + (100._r4 - cfrag_o(no,l)) * 0.01_r4 / 0.58_r4 + end do - ! Adjust pct sand and pct clay to be nearest integers and to be loam if pctlnd_pft is < 1.e-6 - ! Truncate all percentage fields on output grid. This is needed to insure that wt is zero - ! (not a very small number such as 1e-16) where it really should be zero - do no = 1,ns_o - do k = 1,nlevsoi - sand_o(no,k) = float(nint(sand_o(no,k))) - clay_o(no,k) = float(nint(clay_o(no,k))) - end do - if (pctlnd_pft_o(no) < 1.e-6_r8) then - sand_o(no,:) = 43._r8 - clay_o(no,:) = 18._r8 end if - end do - ! Write out fields - if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil percent sand" - call mkfile_output(pioid_o, mesh_o, 'mapunits', mapunit_o, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for mapunits') + end do if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil percent sand" call mkfile_output(pioid_o, mesh_o, 'PCT_SAND', sand_o, lev1name='nlevsoi', rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil percent clay" call mkfile_output(pioid_o, mesh_o, 'PCT_CLAY', clay_o, lev1name='nlevsoi', rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil organic matter" + call mkfile_output(pioid_o, mesh_o, 'ORGANIC', organic_o, lev1name='nlevsoi', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil organic carbon content" + call mkfile_output(pioid_o, mesh_o, 'ORGC', orgc_o, lev1name='nlevsoi', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out coarse fragments in soil" + call mkfile_output(pioid_o, mesh_o, 'CFRAG', cfrag_o, lev1name='nlevsoi', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil bulk density" + call mkfile_output(pioid_o, mesh_o, 'BULK', bulk_o, lev1name='nlevsoi', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil pH measured in H2O" + call mkfile_output(pioid_o, mesh_o, 'PHAQ', phaq_o, lev1name='nlevsoi', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error calling mkfile_output') + call pio_syncfile(pioid_o) ! Release memory @@ -339,7 +555,7 @@ subroutine mksoiltex(file_mesh_i, file_data_i, mesh_o, pioid_o, pctlnd_pft_o, rc call ESMF_VMLogMemInfo("After destroy operations in "//trim(subname)) if (root_task) then - write (ndiag,'(a)') 'Successfully made %sand and %clay' + write(ndiag,'(a)') 'Successfully made %sand, %clay, orgc, cfrag, bulk, phaq .....' end if end subroutine mksoiltex diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index d3228b06c1..a7f021e9b4 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -22,17 +22,18 @@ program mksurfdata ! mksrf_fglacierregion_mesh - Mesh for mksrf_fglacierregion ! mksrf_flai - Leaf Area Index dataset ! mksrf_flai_mesh - Mesh for mksrf_flai - ! mksrf_flakwat - Lake water dataset - ! mksrf_flakwat_mesh - Mesh for mksrf_flakwat + ! mksrf_fpctlak - Percent lake dataset + ! mksrf_fpctlak_mesh - Mesh for mksrf_fpctlak + ! mksrf_flakdep - Lake depth dataset + ! mksrf_flakdep_mesh - Mesh for mksrf_flakdep ! mksrf_fwetlnd - Wetland water dataset ! mksrf_fwetlnd_mesh - Mesh for mksrf_fwetlnd - ! mksrf_forganic - Organic soil carbon dataset - ! mksrf_forganic_mesh - Mesh for mksrf_forganic ! mksrf_fmax - Max fractional saturated area dataset ! mksrf_fmax_mesh - Mesh for mksrf_fmax ! mksrf_fsoicol - Soil color dataset ! mksrf_fsoicol_mesh - Mesh for mksrf_fsoicol - ! mksrf_fsoitex - Soil texture dataset + ! mksrf_fsoitex - Soil texture dataset in mapunits + ! mksrf_fsoitex_lookup - Soil texture lookup for converting mapunits to sand/silt/clay and organic carbon content ! mksrf_fsoitex_mesh - Mesh for mksrf_fsoitex ! mksrf_furbtopo - Topography dataset (for limiting urban areas) ! mksrf_furbtopo_mesh - Mesh for mksrf_furbtopo @@ -107,8 +108,7 @@ program mksurfdata use mksoildepthMod , only : mksoildepth use mksoilcolMod , only : mksoilcol use mkurbanparMod , only : mkurbanInit, mkurban, mkurbanpar, mkurban_topo, numurbl, update_max_array_urban - use mklanwatMod , only : mklakwat, mkwetlnd, update_max_array_lake - use mkorganicMod , only : mkorganic + use mklanwatMod , only : mkpctlak, mklakdep, mkwetlnd, update_max_array_lake use mkutilsMod , only : normalize_classes_by_gcell, chkerr use mkfileMod , only : mkfile_define_dims, mkfile_define_atts, mkfile_define_vars use mkfileMod , only : mkfile_output @@ -135,7 +135,7 @@ program mksurfdata integer :: ntim ! time sample for dynamic land use integer :: year ! year for dynamic land use integer :: year2 ! year for dynamic land use for harvest file - real(r8) :: suma ! local sum for error check + real(r8) :: suma ! local sum for error check real(r8) :: loc_suma, glob_suma ! local and global sum for error check with mpi_allreduce ! model grid @@ -428,15 +428,18 @@ program mksurfdata end if ! ----------------------------------- - ! Make inland water [pctlak, pctwet] [flakwat] [fwetlnd] + ! Make inland water [pctlak, pctwet] ! ----------------------------------- ! LAKEDEPTH is written out in the subroutine ! Need to keep pctlak and pctwet external for use below allocate ( pctlak(lsize_o)) ; pctlak(:) = spval allocate ( pctlak_max(lsize_o)) ; pctlak_max(:) = spval - call mklakwat(mksrf_flakwat_mesh, mksrf_flakwat, mesh_model, pctlak, pioid, & - fsurdat, rc=rc, do_depth=.true.) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mklatwat') + call mkpctlak(mksrf_fpctlak_mesh, mksrf_fpctlak, mesh_model, pctlak, pioid, & + rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkpctlak') + call mklakdep(mksrf_flakdep_mesh, mksrf_flakdep, mesh_model, pioid, fsurdat, & + rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mklakdep') allocate ( pctwet(lsize_o)) ; pctwet(:) = spval allocate ( pctwet_orig(lsize_o)) ; pctwet_orig(:) = spval @@ -461,11 +464,11 @@ program mksurfdata end if ! ----------------------------------- - ! Make soil texture [pctsand, pctclay] + ! Make soil texture and organic carbon content [pctsand, pctclay, organic] ! ----------------------------------- if (fsurdat /= ' ') then - ! mapunits, PCT_SAND and PCT_CLAY are written out in the subroutine - call mksoiltex( mksrf_fsoitex_mesh, mksrf_fsoitex, mesh_model, pioid, pctlnd_pft, rc=rc) + call mksoiltex( mksrf_fsoitex_mesh, file_mapunit_i=mksrf_fsoitex, file_lookup_i=mksrf_fsoitex_lookup, & + mesh_o=mesh_model, pioid_o=pioid, pctlnd_pft_o=pctlnd_pft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mksoiltex') end if @@ -568,18 +571,10 @@ program mksurfdata if (fsurdat /= ' ') then if (outnc_vic) then call mkVICparams ( mksrf_fvic_mesh, mksrf_fvic, mesh_model, pioid, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkorganic') + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkVICparams') end if end if - ! ----------------------------------- - ! Make organic matter density [organic] [forganic] - ! ----------------------------------- - if (fsurdat /= ' ') then - call mkorganic( mksrf_forganic_mesh, mksrf_forganic, mesh_model, pctlnd_pft, lat, pioid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkorganic') - end if - ! ----------------------------------- ! Make VOC emission factors for isoprene [ef1_btr,ef1_fet,ef1_fdt,ef1_shr,ef1_grs,ef1_crp] ! ----------------------------------- @@ -604,7 +599,7 @@ program mksurfdata end do - ! Save special land unit areas of surface dataset + ! Save special land unit areas of surface dataset pctwet_orig(:) = pctwet(:) pctgla_orig(:) = pctgla(:) @@ -700,7 +695,7 @@ program mksurfdata end if call mkfile_output(pioid, mesh_model, 'PCT_NAT_PFT', pct_nat_pft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_NAT_PFT') - + if (num_cft > 0) then if (root_task) write(ndiag, '(a)') trim(subname)//" writing PCT_CFT" if (lsize_o /= 0) then @@ -932,9 +927,9 @@ program mksurfdata call pio_syncfile(pioid) ! Create pctlak data at model resolution (use original mapping file from lake data) - call mklakwat(mksrf_flakwat_mesh, flakname, mesh_model, pctlak, pioid, & - fsurdat, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mklakwat') + call mkpctlak(mksrf_fpctlak_mesh, flakname, mesh_model, pctlak, pioid, & + rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkpctlak') call pio_syncfile(pioid) call mkurban(mksrf_furban_mesh, furbname, mesh_model, pcturb, &