From 4f613feb7b19c43aab463fea41f92299e34b44c5 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Thu, 6 Apr 2023 08:34:26 -0400 Subject: [PATCH 1/3] sfc_climo_gen - Output fractions of each vegetation/soil type category (#748) The sfc_climo_gen program was updated to optionally output the fraction of each vegetation and soil type category at a model point. The dominant category is also output. The default is to output just the dominant category. Unit and consistency tests were added to test this new capability. The 'readthedocs' was also updated. Fixes #709. --- docs/source/ufs_utils.rst | 4 + driver_scripts/driver_grid.hera.sh | 15 +- driver_scripts/driver_grid.jet.sh | 6 + driver_scripts/driver_grid.orion.sh | 6 + driver_scripts/driver_grid.wcoss2.sh | 6 + reg_tests/grid_gen/driver.hera.sh | 19 +- reg_tests/grid_gen/driver.jet.sh | 19 +- reg_tests/grid_gen/driver.orion.sh | 19 +- reg_tests/grid_gen/driver.wcoss2.sh | 15 +- reg_tests/grid_gen/esg.regional.pct.cat.sh | 75 ++++++ sorc/sfc_climo_gen.fd/CMakeLists.txt | 3 + sorc/sfc_climo_gen.fd/driver.F90 | 20 +- sorc/sfc_climo_gen.fd/interp.F90 | 59 +++-- sorc/sfc_climo_gen.fd/interp_frac_cats.F90 | 281 +++++++++++++++++++++ sorc/sfc_climo_gen.fd/model_grid.F90 | 71 ++++-- sorc/sfc_climo_gen.fd/output_frac_cats.F90 | 218 ++++++++++++++++ sorc/sfc_climo_gen.fd/program_setup.f90 | 19 +- sorc/sfc_climo_gen.fd/search_frac_cats.f90 | 119 +++++++++ tests/sfc_climo_gen/CMakeLists.txt | 12 +- tests/sfc_climo_gen/data/fort.41 | 21 ++ tests/sfc_climo_gen/ftst_program_setup.F90 | 49 ++++ ush/sfc_climo_gen.sh | 6 + util/sfc_climo_gen/run.wcoss2.sh | 8 +- 23 files changed, 1000 insertions(+), 70 deletions(-) create mode 100755 reg_tests/grid_gen/esg.regional.pct.cat.sh create mode 100644 sorc/sfc_climo_gen.fd/interp_frac_cats.F90 create mode 100644 sorc/sfc_climo_gen.fd/output_frac_cats.F90 create mode 100644 sorc/sfc_climo_gen.fd/search_frac_cats.f90 create mode 100644 tests/sfc_climo_gen/data/fort.41 create mode 100644 tests/sfc_climo_gen/ftst_program_setup.F90 diff --git a/docs/source/ufs_utils.rst b/docs/source/ufs_utils.rst index 48515443d..69acb5417 100644 --- a/docs/source/ufs_utils.rst +++ b/docs/source/ufs_utils.rst @@ -484,10 +484,13 @@ Location of the source code: ./sorc/sfc_climo_gen.fd. Brief description of each * driver.F90 - The main driver routine. * interp.F90 - The interpolation driver routine. Reads the input source data and interpolates it to the model grid. + * interp_frac_cats.F90 - Same as interp.F90, but for computing the fraction of each soil and vegetation type category. (When namelist variable 'vegsoilt_frac' is true). * model_grid.F90 - Defines the ESMF grid object for the model grid. * output.f90 - Writes the output surface data to a NetCDF file. For regional grids, will output separate files with and without the halo. + * output_frac_cats.f90 - Same as output.f90, but for writing fractional soil and vegetation type. (When namelist variable 'vegsoilt_frac' is true). * program_setup.f90 - Reads the namelist and sets up program execution. * search.f90 - Replace undefined values on the model grid with a valid value at a nearby neighbor. Undefined values are typically associated with isolated islands where there is no source data. + * search_frac_cats.f90 - Same as search.f90, but for the fractional soil and vegetation type option. (When namelist variable 'vegsoilt_frac' is true). * source_grid.F90 - Reads the grid specifications and land/sea mask for the source data. Sets up the ESMF grid object for the source grid. * utils.f90 - Contains error handling utility routines. @@ -512,6 +515,7 @@ Program execution is controlled via a namelist. The namelist variables are: * maximum_snow_albedo_method - interpolation method for this field. Bilinear or conservative. Default is bilinear. * snowfree_albedo_method - interpolation method for this field. Bilinear or conservative. Default is bilinear. * vegetation_greenness_method - interpolation method for this field. Bilinear or conservative. Default is bilinear. + * vegsoilt_frac - When 'true', outputs the dominate soil and vegetation type, and the fraction of each category. When 'false', only outputs the dominate categories. Default is 'false'. Program inputs and outputs -------------------------- diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh index 1c9050a7b..3f10f6095 100755 --- a/driver_scripts/driver_grid.hera.sh +++ b/driver_scripts/driver_grid.hera.sh @@ -6,8 +6,9 @@ #SBATCH -o log.fv3_grid_driver #SBATCH -e log.fv3_grid_driver #SBATCH --nodes=1 --ntasks-per-node=24 +##SBATCH --partition=bigmem #SBATCH -q debug -#SBATCH -t 00:30:00 +#SBATCH -t 00:20:00 #----------------------------------------------------------------------- # Driver script to create a cubic-sphere based model grid on Hera. @@ -23,8 +24,9 @@ # Note: The sfc_climo_gen program only runs with an # mpi task count that is a multiple of six. This is # an ESMF library requirement. Large grids may require -# tasks spread across multiple nodes. The orography code -# benefits from threads. +# tasks spread across multiple nodes or to be run on +# 'bigmem' nodes (#SBATCH --partition=bigmem). The +# orography code benefits from threads. # # To run, do the following: # @@ -74,9 +76,16 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg'. + export make_gsl_orog=false # When 'true' will output 'oro' files for # the GSL orographic drag suite. +export vegsoilt_frac='.false.' # When .false., output dominant soil and + # vegetation type category. When .true., + # output fraction of each category and + # the dominant category. A Fortran logical, + # so include the dots. + export veg_type_src="modis.igbp.0.05" # Vegetation type data. # For viirs-based vegetation type data, set to: # 1) "viirs.igbp.0.1" for global 0.10-deg data diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh index 62869e540..2a9ef3e9f 100755 --- a/driver_scripts/driver_grid.jet.sh +++ b/driver_scripts/driver_grid.jet.sh @@ -78,6 +78,12 @@ export gtype=uniform # 'uniform', 'stretch', 'nest', export make_gsl_orog=false # When 'true' will output 'oro' files for # the GSL orographic drag suite. +export vegsoilt_frac='.false.' # When true, outputs percent of each + # soil and veg type category and a + # dominant category. When false, only + # outputs the dominant category. A + # Fortran logical, so include the dots. + export veg_type_src="modis.igbp.0.05" # Vegetation type data. # For viirs-based vegetation type data, set to: # 1) "viirs.igbp.0.1" for global 0.10-deg data diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh index f0a74a33e..3b56094c3 100755 --- a/driver_scripts/driver_grid.orion.sh +++ b/driver_scripts/driver_grid.orion.sh @@ -77,6 +77,12 @@ export gtype=regional_esg # 'uniform', 'stretch', 'nest', export make_gsl_orog=false # When 'true' will output 'oro' files for # the GSL orographic drag suite. +export vegsoilt_frac='.false.' # When true, outputs percent of each + # soil and veg type category and a + # dominant category. When false, only + # outputs the dominant category. A + # Fortran logical, so include the dots. + export veg_type_src="modis.igbp.0.05" # Vegetation type data. # For viirs-based vegetation type data, set to: # 1) "viirs.igbp.0.1" for global 0.10-deg data diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh index a0c290b90..1a9c4bfb3 100755 --- a/driver_scripts/driver_grid.wcoss2.sh +++ b/driver_scripts/driver_grid.wcoss2.sh @@ -76,6 +76,12 @@ export gtype=regional_esg # 'uniform', 'stretch', 'nest', export make_gsl_orog=false # When 'true' will output 'oro' files for # the GSL orographic drag suite. +export vegsoilt_frac='.false.' # When true, outputs percent of each + # soil and veg type category and a + # dominant category. When false, only + # outputs the dominant category. A + # Fortran logical, so include the dots. + export veg_type_src="modis.igbp.0.05" # Vegetation type data. # For viirs-based vegetation type data, set to: # 1) "viirs.igbp.0.1" for global 0.10-deg data diff --git a/reg_tests/grid_gen/driver.hera.sh b/reg_tests/grid_gen/driver.hera.sh index fc7199374..ef3b37b1c 100755 --- a/reg_tests/grid_gen/driver.hera.sh +++ b/reg_tests/grid_gen/driver.hera.sh @@ -86,7 +86,7 @@ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_ -o $LOG_FILE3 -e $LOG_FILE3 ./gfdl.regional.sh) #----------------------------------------------------------------------------- -# esg regional grid +# ESG regional grid (output dominant soil/vegetation type). #----------------------------------------------------------------------------- LOG_FILE4=${LOG_FILE}04 @@ -94,19 +94,28 @@ TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_ -o $LOG_FILE4 -e $LOG_FILE4 ./esg.regional.sh) #----------------------------------------------------------------------------- -# Regional GSL gravity wave drag test. +# ESG regional grid (output percent of each soil and vegetation type and +# the dominant category). #----------------------------------------------------------------------------- LOG_FILE5=${LOG_FILE}05 -TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ - -o $LOG_FILE5 -e $LOG_FILE5 ./regional.gsl.gwd.sh) +TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional.pct.cat \ + -o $LOG_FILE5 -e $LOG_FILE5 ./esg.regional.pct.cat.sh) + +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag test. +#----------------------------------------------------------------------------- + +LOG_FILE6=${LOG_FILE}06 +TEST6=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + -o $LOG_FILE6 -e $LOG_FILE6 ./regional.gsl.gwd.sh) #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - --open-mode=append -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5 << EOF + --open-mode=append -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6 << EOF #!/bin/bash grep -a '<<<' ${LOG_FILE}* > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.jet.sh b/reg_tests/grid_gen/driver.jet.sh index fc5b6aa27..2d3f5fcfc 100755 --- a/reg_tests/grid_gen/driver.jet.sh +++ b/reg_tests/grid_gen/driver.jet.sh @@ -84,7 +84,7 @@ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_ --partition=xjet -o $LOG_FILE3 -e $LOG_FILE3 ./gfdl.regional.sh) #----------------------------------------------------------------------------- -# ESG regional grid +# ESG regional grid (output dominant soil/vegetation type). #----------------------------------------------------------------------------- LOG_FILE4=${LOG_FILE}04 @@ -92,19 +92,28 @@ TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_ --partition=xjet -o $LOG_FILE4 -e $LOG_FILE4 ./esg.regional.sh) #----------------------------------------------------------------------------- -# Regional GSL gravity wave drag. +# ESG regional grid (output percent of each soil and vegetation type and +# the dominant category). #----------------------------------------------------------------------------- LOG_FILE5=${LOG_FILE}05 -TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ - --partition=xjet -o $LOG_FILE5 -e $LOG_FILE5 ./regional.gsl.gwd.sh) +TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional.pct.cat \ + --partition=xjet -o $LOG_FILE5 -e $LOG_FILE5 ./esg.regional.pct.cat.sh) + +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag. +#----------------------------------------------------------------------------- + +LOG_FILE6=${LOG_FILE}06 +TEST6=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:07:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + --partition=xjet -o $LOG_FILE6 -e $LOG_FILE6 ./regional.gsl.gwd.sh) #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --partition=xjet --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5 << EOF + -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6 << EOF #!/bin/bash grep -a '<<<' ${LOG_FILE}* > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.orion.sh b/reg_tests/grid_gen/driver.orion.sh index 08cdd2e1b..a3a60e678 100755 --- a/reg_tests/grid_gen/driver.orion.sh +++ b/reg_tests/grid_gen/driver.orion.sh @@ -82,7 +82,7 @@ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ -o $LOG_FILE3 -e $LOG_FILE3 ./gfdl.regional.sh) #----------------------------------------------------------------------------- -# ESG regional grid +# ESG regional grid (output dominant soil/vegetation type). #----------------------------------------------------------------------------- LOG_FILE4=${LOG_FILE}04 @@ -90,19 +90,28 @@ TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ -o $LOG_FILE4 -e $LOG_FILE4 ./esg.regional.sh) #----------------------------------------------------------------------------- -# Regional grid with GSL gravity wave drag fields. +# ESG regional grid (output percent of each soil and vegetation type and +# the dominant category). #----------------------------------------------------------------------------- LOG_FILE5=${LOG_FILE}05 -TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ - -o $LOG_FILE5 -e $LOG_FILE5 ./regional.gsl.gwd.sh) +TEST5=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional.pct.cat \ + -o $LOG_FILE5 -e $LOG_FILE5 ./esg.regional.pct.cat.sh) + +#----------------------------------------------------------------------------- +# Regional grid with GSL gravity wave drag fields. +#----------------------------------------------------------------------------- + +LOG_FILE6=${LOG_FILE}06 +TEST6=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + -o $LOG_FILE6 -e $LOG_FILE6 ./regional.gsl.gwd.sh) #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5 << EOF + -q $QUEUE -d afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6 << EOF #!/bin/bash grep -a '<<<' ${LOG_FILE}* > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.wcoss2.sh b/reg_tests/grid_gen/driver.wcoss2.sh index efa7a8f44..da8b9d0ed 100755 --- a/reg_tests/grid_gen/driver.wcoss2.sh +++ b/reg_tests/grid_gen/driver.wcoss2.sh @@ -87,7 +87,7 @@ TEST3=$(qsub -V -o $LOG_FILE3 -e $LOG_FILE3 -q $QUEUE -A $PROJECT_CODE -l wallti -N gfdl.regional -l select=1:ncpus=30:mem=40GB $PWD/gfdl.regional.sh) #----------------------------------------------------------------------------- -# esg regional grid +# ESG regional grid (output dominant soil/vegetation type). #----------------------------------------------------------------------------- LOG_FILE4=${LOG_FILE}04 @@ -95,11 +95,20 @@ TEST4=$(qsub -V -o $LOG_FILE4 -e $LOG_FILE4 -q $QUEUE -A $PROJECT_CODE -l wallti -N esg.regional -l select=1:ncpus=30:mem=40GB $PWD/esg.regional.sh) #----------------------------------------------------------------------------- -# Regional GSL gravity wave drag test. +# ESG regional grid (output percent of each soil and vegetation type and +# the dominant category). #----------------------------------------------------------------------------- LOG_FILE5=${LOG_FILE}05 TEST5=$(qsub -V -o $LOG_FILE5 -e $LOG_FILE5 -q $QUEUE -A $PROJECT_CODE -l walltime=00:07:00 \ + -N esg.regional.pct.cat -l select=1:ncpus=30:mem=40GB $PWD/esg.regional.pct.cat.sh) + +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag test. +#----------------------------------------------------------------------------- + +LOG_FILE6=${LOG_FILE}06 +TEST6=$(qsub -V -o $LOG_FILE6 -e $LOG_FILE6 -q $QUEUE -A $PROJECT_CODE -l walltime=00:07:00 \ -N rsg.gsl.gwd -l select=1:ncpus=30:mem=40GB $PWD/regional.gsl.gwd.sh) #----------------------------------------------------------------------------- @@ -107,7 +116,7 @@ TEST5=$(qsub -V -o $LOG_FILE5 -e $LOG_FILE5 -q $QUEUE -A $PROJECT_CODE -l wallti #----------------------------------------------------------------------------- qsub -V -o ${LOG_FILE} -e ${LOG_FILE} -q $QUEUE -A $PROJECT_CODE -l walltime=00:02:00 \ - -N grid_summary -l select=1:ncpus=1:mem=100MB -W depend=afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5 << EOF + -N grid_summary -l select=1:ncpus=1:mem=100MB -W depend=afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6 << EOF #!/bin/bash cd ${this_dir} grep -a '<<<' ${LOG_FILE}* | grep -v echo > $SUM_FILE diff --git a/reg_tests/grid_gen/esg.regional.pct.cat.sh b/reg_tests/grid_gen/esg.regional.pct.cat.sh new file mode 100755 index 000000000..2afbed9e1 --- /dev/null +++ b/reg_tests/grid_gen/esg.regional.pct.cat.sh @@ -0,0 +1,75 @@ +#!/bin/bash + +#----------------------------------------------------------------------- +# Create a regional esg grid. Output dominant soil and vegetation +# categories and well as the percentage of each category. +# Compare output to a set of baseline files using the 'nccmp' +# utility. This script is run by the machine specific driver script. +#----------------------------------------------------------------------- + +set -x + +TEST_NAME="esg.regional.pct.cat" +export TEMP_DIR=${WORK_DIR}/${TEST_NAME}.work +export out_dir=${WORK_DIR}/${TEST_NAME} + +export gtype=regional_esg +export target_lon=-97.5 # Center longitude of the highest resolution tile +export target_lat=35.5 # Center latitude of the highest resolution tile +export idim=1301 # Dimension of grid in 'i' direction +export jdim=600 # Dimension of grid in 'j' direction +export delx=0.0145 # Grid spacing in degrees in 'i' direction +export dely=0.0145 # Grid spacing in degrees in 'j' direction +export halo=3 +export vegsoilt_frac=.true. # Output dominant soil/veg categories as well + # as the percentage of each category. + +NCCMP=${NCCMP:-$(which nccmp)} + +#----------------------------------------------------------------------- +# Start script. +#----------------------------------------------------------------------- + +echo "Starting at: " `date` + +$home_dir/ush/fv3gfs_driver_grid.sh + +iret=$? +if [ $iret -ne 0 ]; then + set +x + echo "<<< ESG REGIONAL PERCENT CATEGORY TEST FAILED. <<<" + exit $iret +fi + +echo "Ending at: " `date` + +#----------------------------------------------------------------------------- +# Compare output to baseline set of data. +#----------------------------------------------------------------------------- + +cd $out_dir/C3113 + +test_failed=0 +for files in *tile*.nc ./fix_sfc/*tile*.nc +do + if [ -f $files ]; then + echo CHECK $files + $NCCMP -dmfqS $files $HOMEreg/${TEST_NAME}/$files + iret=$? + if [ $iret -ne 0 ]; then + test_failed=1 + fi + fi +done + +set +x +if [ $test_failed -ne 0 ]; then + echo "<<< ESG REGIONAL PERCENT CATEGORY TEST FAILED. >>>" + if [ "$UPDATE_BASELINE" = "TRUE" ]; then + $home_dir/reg_tests/update_baseline.sh "${HOMEreg}/.." "${TEST_NAME}" $commit_num + fi +else + echo "<<< ESG REGIONAL PERCENT CATEGORY TEST PASSED. >>>" +fi + +exit 0 diff --git a/sorc/sfc_climo_gen.fd/CMakeLists.txt b/sorc/sfc_climo_gen.fd/CMakeLists.txt index f5e3b9186..a9103106f 100644 --- a/sorc/sfc_climo_gen.fd/CMakeLists.txt +++ b/sorc/sfc_climo_gen.fd/CMakeLists.txt @@ -5,10 +5,13 @@ set(lib_src interp.F90 + interp_frac_cats.F90 model_grid.F90 output.f90 + output_frac_cats.F90 program_setup.f90 search.f90 + search_frac_cats.f90 source_grid.F90 utils.f90) diff --git a/sorc/sfc_climo_gen.fd/driver.F90 b/sorc/sfc_climo_gen.fd/driver.F90 index 60014fda9..520cc99d8 100644 --- a/sorc/sfc_climo_gen.fd/driver.F90 +++ b/sorc/sfc_climo_gen.fd/driver.F90 @@ -80,8 +80,14 @@ program driver !------------------------------------------------------------------------- call define_source_grid(localpet, npets, input_vegetation_type_file) - method=ESMF_REGRIDMETHOD_NEAREST_STOD - call interp(localpet, method, input_vegetation_type_file) + if (fract_vegsoil_type) then + print*,'- WILL OUTPUT VEGETATION TYPE FRACTION AND DOMINANT CATEGORY.' + call interp_frac_cats(localpet, input_vegetation_type_file) + else + print*,'- WILL OUTPUT DOMINANT VEGETATION TYPE.' + method=ESMF_REGRIDMETHOD_NEAREST_STOD + call interp(localpet, method, input_vegetation_type_file) + endif call source_grid_cleanup ! Snow free albedo @@ -136,8 +142,14 @@ program driver if (trim(input_soil_type_file) /= "NULL") then call define_source_grid(localpet, npets, input_soil_type_file) - method=ESMF_REGRIDMETHOD_NEAREST_STOD - call interp(localpet, method, input_soil_type_file) + if (fract_vegsoil_type) then + print*,'- WILL OUTPUT SOIL TYPE FRACTION AND DOMINANT CATEGORY.' + call interp_frac_cats(localpet, input_soil_type_file) + else + print*,'- WILL OUTPUT DOMINANT SOIL TYPE.' + method=ESMF_REGRIDMETHOD_NEAREST_STOD + call interp(localpet, method, input_soil_type_file) + endif call source_grid_cleanup endif diff --git a/sorc/sfc_climo_gen.fd/interp.F90 b/sorc/sfc_climo_gen.fd/interp.F90 index af374478f..877930724 100644 --- a/sorc/sfc_climo_gen.fd/interp.F90 +++ b/sorc/sfc_climo_gen.fd/interp.F90 @@ -15,6 +15,7 @@ subroutine interp(localpet, method, input_file) use esmf use netcdf use model_grid + use program_setup, only : fract_vegsoil_type use source_grid use utils use mpi @@ -185,16 +186,22 @@ subroutine interp(localpet, method, input_file) enddo -! These fields are adjusted at landice. +! Adjust some fields at permanent land ice points. These points are identified by the +! 'permanent ice' vegetation type category. +! +! When outputting the fraction of each vegetation type, land ice points are +! not defined. So don't do this adjustment. - select case (trim(field_names(n))) - case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') - if (localpet == 0) then - allocate(vegt_mdl_one_tile(i_mdl,j_mdl)) - else - allocate(vegt_mdl_one_tile(0,0)) - endif - end select + if (.not. fract_vegsoil_type) then + select case (trim(field_names(n))) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') + if (localpet == 0) then + allocate(vegt_mdl_one_tile(i_mdl,j_mdl)) + else + allocate(vegt_mdl_one_tile(0,0)) + endif + end select + endif OUTPUT_LOOP : do tile = 1, num_tiles @@ -218,30 +225,36 @@ subroutine interp(localpet, method, input_file) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather.", rc) - select case (trim(field_names(n))) - case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') - print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE." - call ESMF_FieldGather(vegt_field_mdl, vegt_mdl_one_tile, rootPet=0, tile=tile, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + if (.not. fract_vegsoil_type) then + select case (trim(field_names(n))) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') + print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE." + call ESMF_FieldGather(vegt_field_mdl, vegt_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather.", rc) - end select + end select + endif if (localpet == 0) then print*,'- CALL SEARCH FOR TILE ',tile call search (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, tile, field_names(n)) - select case (field_names(n)) - case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') - call adjust_for_landice (data_mdl_one_tile, vegt_mdl_one_tile, i_mdl, j_mdl, field_names(n)) - end select + if (.not. fract_vegsoil_type) then + select case (field_names(n)) + case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color') + call adjust_for_landice (data_mdl_one_tile, vegt_mdl_one_tile, i_mdl, j_mdl, field_names(n)) + end select + endif where(mask_mdl_one_tile == 0) data_mdl_one_tile = missing call output (data_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, tile, record, t, n) endif - if (field_names(n) == 'vegetation_type') then - print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE." - call ESMF_FieldScatter(vegt_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + if (.not. fract_vegsoil_type) then + if (field_names(n) == 'vegetation_type') then + print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE." + call ESMF_FieldScatter(vegt_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter.", rc) + endif endif enddo OUTPUT_LOOP diff --git a/sorc/sfc_climo_gen.fd/interp_frac_cats.F90 b/sorc/sfc_climo_gen.fd/interp_frac_cats.F90 new file mode 100644 index 000000000..d82091acf --- /dev/null +++ b/sorc/sfc_climo_gen.fd/interp_frac_cats.F90 @@ -0,0 +1,281 @@ +!> @file +!! @brief Read the input cateogorical source data and +!! interpolate percent of each category to the model grid. +!! +!! @author George Gayno @date 2022 + +!> Read the input source data and interpolate it to the +!! model grid. Outputs the percentage of each category +!! within a model grid box and the dominant category. +!! +!! @param[in] localpet this mpi task +!! @param[in] input_file filename of input source data. +!! @author George Gayno @date 2022 + subroutine interp_frac_cats(localpet, input_file) + + use esmf + use netcdf + use model_grid, only : grid_mdl, i_mdl, j_mdl, & + num_tiles, latitude_field_mdl, & + longitude_field_mdl, mask_field_mdl, & + land_frac_field_mdl + use source_grid + use output_frac_cats, only : output_driver + use utils + use mpi + + implicit none + + character(len=*), intent(in) :: input_file + + integer, intent(in) :: localpet + + integer :: i, j, tile, cat, ncid, status, rc + integer :: varid, water_category, max_cat + integer :: isrctermprocessing + integer :: category, num_categories + + integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:) + integer(esmf_kind_i4), pointer :: unmapped_ptr(:) + + real(esmf_kind_r4), allocatable :: data_src_global(:,:) + real(esmf_kind_r4), allocatable :: data_src_global2(:,:,:) + real(esmf_kind_r4), allocatable :: data_mdl_one_tile(:,:,:) + real(esmf_kind_r4), allocatable :: dom_cat_mdl_one_tile(:,:) + real(esmf_kind_r4), allocatable :: lat_mdl_one_tile(:,:) + real(esmf_kind_r4), allocatable :: sum_mdl_one_tile(:,:) + real(esmf_kind_r4), allocatable :: lon_mdl_one_tile(:,:) + real(esmf_kind_r4), allocatable :: land_frac_mdl_one_tile(:,:) + real(esmf_kind_r4) :: max_frac + + type(esmf_regridmethod_flag) :: method + type(esmf_field) :: data_field_src + type(esmf_field) :: data_field_mdl + type(esmf_routehandle) :: regrid_data + type(esmf_polemethod_flag) :: pole + + if (localpet == 0) then + allocate(data_src_global(i_src,j_src)) + else + allocate(data_src_global(0,0)) + endif + + if (localpet == 0) then + print*,'- OPEN SOURCE FILE ', trim(input_file) + status = nf90_open(input_file, nf90_nowrite, ncid) + call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE") + status = nf90_inq_varid(ncid, field_names(1), varid) + call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID") + status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,1/), count=(/i_src,j_src,1/)) + call netcdf_err(status, "IN ROUTINE INTERP READING FIELD") + print*,'number of cats ',maxval(data_src_global) + num_categories = nint(maxval(data_src_global)) + status = nf90_get_att(ncid, varid, 'water_category', water_category) + call netcdf_err(status, "IN ROUTINE INTERP READING water_category") + print*,'water cat ',water_category + endif + + call mpi_bcast(num_categories,1,MPI_INTEGER,0,MPI_COMM_WORLD,status) + + print*,"- CALL FieldCreate FOR SOURCE GRID DATA." + data_field_src = ESMF_FieldCreate(grid_src, & + typekind=ESMF_TYPEKIND_R4, & + indexflag=ESMF_INDEX_GLOBAL, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/num_categories/), & + name="source data", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR model GRID veg DATA." + data_field_mdl = ESMF_FieldCreate(grid_mdl, & + typekind=ESMF_TYPEKIND_R4, & + indexflag=ESMF_INDEX_GLOBAL, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/num_categories/), & + name="mdl data", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + if (localpet == 0) then + allocate(data_src_global2(i_src,j_src,num_categories)) + allocate(data_mdl_one_tile(i_mdl,j_mdl,num_categories)) + allocate(dom_cat_mdl_one_tile(i_mdl,j_mdl)) + allocate(mask_mdl_one_tile(i_mdl,j_mdl)) + allocate(land_frac_mdl_one_tile(i_mdl,j_mdl)) + allocate(lat_mdl_one_tile(i_mdl,j_mdl)) + allocate(sum_mdl_one_tile(i_mdl,j_mdl)) + allocate(lon_mdl_one_tile(i_mdl,j_mdl)) + else + allocate(data_src_global2(0,0,0)) + allocate(data_mdl_one_tile(0,0,0)) + allocate(dom_cat_mdl_one_tile(0,0)) + allocate(mask_mdl_one_tile(0,0)) + allocate(land_frac_mdl_one_tile(0,0)) + allocate(lat_mdl_one_tile(0,0)) + allocate(sum_mdl_one_tile(0,0)) + allocate(lon_mdl_one_tile(0,0)) + endif + + if (localpet == 0) then + data_src_global2 = 0.0 + do j = 1, j_src + do i = 1, i_src + category = nint(data_src_global(i,j)) + if (category < 1) cycle + data_src_global2(i,j,category) = 1.0 + enddo + enddo + endif + + deallocate(data_src_global) + + print*,"- CALL FieldScatter FOR SOURCE GRID DATA." + call ESMF_FieldScatter(data_field_src, data_src_global2, rootPet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter.", rc) + + deallocate(data_src_global2) + + isrctermprocessing = 1 + + method = ESMF_REGRIDMETHOD_CONSERVE + pole = ESMF_POLEMETHOD_NONE + + print*,"- CALL FieldRegridStore." + nullify(unmapped_ptr) + call ESMF_FieldRegridStore(data_field_src, & + data_field_mdl, & + srcmaskvalues=(/0/), & + dstmaskvalues=(/0/), & + polemethod=pole, & + unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & + normtype=ESMF_NORMTYPE_FRACAREA, & + srctermprocessing=isrctermprocessing, & + routehandle=regrid_data, & + regridmethod=method, & + unmappedDstList=unmapped_ptr, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldRegridStore.", rc) + + print*,"- CALL Field_Regrid." + call ESMF_FieldRegrid(data_field_src, & + data_field_mdl, & + routehandle=regrid_data, & + termorderflag=ESMF_TERMORDER_SRCSEQ, & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldRegrid.", rc) + + print*,"- CALL FieldRegridRelease." + call ESMF_FieldRegridRelease(routehandle=regrid_data, rc=rc) + + print*,"- CALL FieldDestroy." + call ESMF_FieldDestroy(data_field_src, rc=rc) + + OUTPUT_LOOP : do tile = 1, num_tiles + + print*,"- CALL FieldGather FOR MODEL LATITUDE." + call ESMF_FieldGather(latitude_field_mdl, lat_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + print*,"- CALL FieldGather FOR MODEL LONGITUDE." + call ESMF_FieldGather(longitude_field_mdl, lon_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + print*,"- CALL FieldGather FOR MODEL GRID DATA." + call ESMF_FieldGather(data_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + print*,"- CALL FieldGather FOR MODEL GRID LAND MASK." + call ESMF_FieldGather(mask_field_mdl, mask_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + print*,"- CALL FieldGather FOR MODEL GRID LAND FRACTION." + call ESMF_FieldGather(land_frac_field_mdl, land_frac_mdl_one_tile, rootPet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather.", rc) + + if (localpet == 0) then + print*,'- CALL SEARCH FOR TILE ',tile + +! Where sum is zero, the regridding did not find any input data for the model point +! (ex. and isolated island). Call the search routine at these points. + sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3) + do j = 1, j_mdl + do i = 1, i_mdl + if (mask_mdl_one_tile(i,j) == 1 .and. sum_mdl_one_tile(i,j) == 0.0) then + data_mdl_one_tile(i,j,:) = -9999.9 ! flag to tell search routine to search. + endif + enddo + enddo + call search_frac_cats (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, field_names(1)) + print*,'after regrid ',data_mdl_one_tile(i_mdl/2,j_mdl/2,:) + +! These points are all non-land. Set to 100% of the water category. + do j = 1, j_mdl + do i = 1, i_mdl + if (mask_mdl_one_tile(i,j) == 0) then + data_mdl_one_tile(i,j,water_category) = 1.0 + endif + enddo + enddo + +! For fractional grids, need to rescale the category percentages by the +! fraction of land in the model grid cell. + +! When running with fractional grids, the land_frac_mdl_one_tile array will +! contain a fraction between 0 and 1. When not running with fractional +! grids, this array will contain negative fill values. + + if (maxval(land_frac_mdl_one_tile) > 0.0) then + print*,'before rescale ',tile,land_frac_mdl_one_tile(42,82),data_mdl_one_tile(42,82,:) + do j = 1, j_mdl + do i = 1, i_mdl + if (mask_mdl_one_tile(i,j) == 1) then + data_mdl_one_tile(i,j,:) = data_mdl_one_tile(i,j,:) * land_frac_mdl_one_tile(i,j) + data_mdl_one_tile(i,j,water_category) = 1.0 - land_frac_mdl_one_tile(i,j) + max_frac = 0.0 + max_cat = -9 + do cat = 1, num_categories + if (cat == water_category) cycle + if(data_mdl_one_tile(i,j,cat) > max_frac) then + max_frac = data_mdl_one_tile(i,j,cat) + max_cat = cat + endif + enddo + dom_cat_mdl_one_tile(i,j) = max_cat + else + dom_cat_mdl_one_tile(i,j) = water_category + endif + enddo + enddo + else ! non-fractonal grids. + dom_cat_mdl_one_tile = 0.0 + dom_cat_mdl_one_tile = maxloc(data_mdl_one_tile,dim=3) + endif + call output_driver (data_mdl_one_tile, dom_cat_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, num_categories, tile) + endif + + enddo OUTPUT_LOOP + + status=nf90_close(ncid) + + deallocate(data_mdl_one_tile, dom_cat_mdl_one_tile, mask_mdl_one_tile) + deallocate(lat_mdl_one_tile, lon_mdl_one_tile, sum_mdl_one_tile, land_frac_mdl_one_tile) + + print*,"- CALL FieldDestroy." + call ESMF_FieldDestroy(data_field_mdl, rc=rc) + + call mpi_barrier(mpi_comm_world, rc) + + end subroutine interp_frac_cats diff --git a/sorc/sfc_climo_gen.fd/model_grid.F90 b/sorc/sfc_climo_gen.fd/model_grid.F90 index 74763b1a1..abc6b4b37 100644 --- a/sorc/sfc_climo_gen.fd/model_grid.F90 +++ b/sorc/sfc_climo_gen.fd/model_grid.F90 @@ -4,7 +4,7 @@ !> This module defines the model grid. !! -!! Variables named with 'mdl' refer to the model grid. +!! Variables named with '_mdl' refer to the model grid. !! !! @author George Gayno @date 2018 module model_grid @@ -28,8 +28,16 @@ module model_grid type(esmf_grid), public :: grid_mdl !< ESMF grid object for the model grid. type(esmf_field), public :: data_field_mdl !< ESMF field object that holds the !! data interpolated to model grid. + type(esmf_field), public :: land_frac_field_mdl !< ESMF field object that holds the + !! model land fraction. When running + !! with fractional grids, will be + !! between zero and one. For non- + !! fractional grids, will contain a + !! fill value. type(esmf_field), public :: mask_field_mdl !< ESMF field object that holds the - !! model land mask. + !! model land mask. Equal to '1' if + !! point is partial or all land. Equal + !! to zero is point is all non-land. type(esmf_field), public :: latitude_field_mdl !< ESMF field object that holds the !! model grid latitude. type(esmf_field), public :: longitude_field_mdl !< ESMF field object that holds the @@ -75,6 +83,7 @@ subroutine define_model_grid(localpet, npets) real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:) real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:) + real(esmf_kind_r4), allocatable :: land_frac_one_tile(:,:) !----------------------------------------------------------------------- ! Get the number of tiles from the mosaic file. @@ -164,14 +173,16 @@ subroutine define_model_grid(localpet, npets) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate", rc) - print*,"- CALL FieldCreate FOR VEGETATION TYPE INTERPOLATED TO MODEL GRID." - vegt_field_mdl = ESMF_FieldCreate(grid_mdl, & + if (.not. fract_vegsoil_type) then + print*,"- CALL FieldCreate FOR VEGETATION TYPE INTERPOLATED TO MODEL GRID." + vegt_field_mdl = ESMF_FieldCreate(grid_mdl, & typekind=ESMF_TYPEKIND_R4, & staggerloc=ESMF_STAGGERLOC_CENTER, & name="veg type on model grid", & rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate", rc) + endif print*,"- CALL FieldCreate FOR MODEL GRID LATITUDE." latitude_field_mdl = ESMF_FieldCreate(grid_mdl, & @@ -212,12 +223,23 @@ subroutine define_model_grid(localpet, npets) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) + print*,"- CALL FieldCreate FOR MODEL GRID LAND FRACTION." + land_frac_field_mdl = ESMF_FieldCreate(grid_mdl, & + typekind=ESMF_TYPEKIND_R4, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="model grid land fraction", & + rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + if (localpet == 0) then allocate(mask_mdl_one_tile(i_mdl,j_mdl)) + allocate(land_frac_one_tile(i_mdl,j_mdl)) allocate(latitude_one_tile(i_mdl,j_mdl)) allocate(longitude_one_tile(i_mdl,j_mdl)) else allocate(mask_mdl_one_tile(0,0)) + allocate(land_frac_one_tile(0,0)) allocate(latitude_one_tile(0,0)) allocate(longitude_one_tile(0,0)) endif @@ -225,7 +247,7 @@ subroutine define_model_grid(localpet, npets) do tile = 1, num_tiles if (localpet == 0) then the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(tile)) - call get_model_info(trim(the_file), mask_mdl_one_tile, & + call get_model_info(trim(the_file), mask_mdl_one_tile, land_frac_one_tile, & latitude_one_tile, longitude_one_tile, i_mdl, j_mdl) endif @@ -234,6 +256,11 @@ subroutine define_model_grid(localpet, npets) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) + print*,"- CALL FieldScatter FOR MODEL GRID LAND FRACTION. TILE IS: ", tile + call ESMF_FieldScatter(land_frac_field_mdl, land_frac_one_tile, rootpet=0, tile=tile, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + print*,"- CALL FieldScatter FOR MODEL LATITUDE. TILE IS: ", tile call ESMF_FieldScatter(latitude_field_mdl, latitude_one_tile, rootpet=0, tile=tile, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -246,7 +273,7 @@ subroutine define_model_grid(localpet, npets) enddo - deallocate(mask_mdl_one_tile, latitude_one_tile, longitude_one_tile) + deallocate(mask_mdl_one_tile, latitude_one_tile, longitude_one_tile, land_frac_one_tile) print*,"- CALL GridAddItem FOR MODEL GRID." call ESMF_GridAddItem(grid_mdl, & @@ -271,16 +298,17 @@ end subroutine define_model_grid !> Get model information !! -!! Read model land/sea mask and lat/lon from the orography file. +!! Read model land/sea mask, land fraction and lat/lon from the orography file. !! !! @param[in] orog_file the orography file -!! @param[out] mask land/sea mask +!! @param[out] mask land/sea mask 0-all non-land; 1-some land. +!! @param[out] land_frac land fraction between 0 and 1. !! @param[out] lat2d latitude !! @param[out] lon2d longitude !! @param[in] idim i dimension of the model tile !! @param[in] jdim j dimension of the model tile !! @author George Gayno @date 2018 - subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) + subroutine get_model_info(orog_file, mask, land_frac, lat2d, lon2d, idim, jdim) use esmf use netcdf @@ -295,6 +323,7 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) real(esmf_kind_r4), intent(out) :: lat2d(idim,jdim) real(esmf_kind_r4), intent(out) :: lon2d(idim,jdim) + real(esmf_kind_r4), intent(out) :: land_frac(idim,jdim) integer :: error, lat, lon, i, j integer :: ncid, id_dim, id_var @@ -328,11 +357,17 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) allocate(dummy(idim,jdim)) !----------------------------------------------------------------------- -! If the lake maker was used, there will be a 'lake_frac' record. +! If the lake maker was used, we are running with a fractional +! land/non-land grid and there will be a 'lake_frac' record. ! In that case, land/non-land is determined by 'land_frac'. ! ! If the lake maker was not used, use 'slmsk', which is defined ! as the nint(land_frac). +! +! In summary, if 'mask' is one, the point is all land or +! partial land and surface data will be mapped to it. Otherwise, +! when 'mask' is zero, then the point is all non-land and +! surface data will not be mapped to it. !----------------------------------------------------------------------- error=nf90_inq_varid(ncid, 'lake_frac', id_var) @@ -343,16 +378,17 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) error=nf90_get_var(ncid, id_var, dummy) call netcdf_err(error, "READING SLMSK") mask = nint(dummy) + land_frac = -999. else print*,"- READ LAND FRACTION" error=nf90_inq_varid(ncid, 'land_frac', id_var) call netcdf_err(error, "READING LAND_FRAC ID") - error=nf90_get_var(ncid, id_var, dummy) + error=nf90_get_var(ncid, id_var, land_frac) call netcdf_err(error, "READING LAND_FRAC") mask = 0 do j = 1, lat do i = 1, lon - if (dummy(i,j) > 0.0) then + if (land_frac(i,j) > 0.0) then mask(i,j) = 1 endif enddo @@ -396,11 +432,16 @@ subroutine model_grid_cleanup print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK." call ESMF_FieldDestroy(mask_field_mdl,rc=rc) + print*,"- CALL FieldDestroy FOR MODEL GRID LAND FRACTION." + call ESMF_FieldDestroy(land_frac_field_mdl,rc=rc) + print*,"- CALL FieldDestroy FOR MODEL GRID DATA FIELD." call ESMF_FieldDestroy(data_field_mdl,rc=rc) - print*,"- CALL FieldDestroy FOR MODEL GRID VEGETATION TYPE." - call ESMF_FieldDestroy(vegt_field_mdl,rc=rc) + if (ESMF_FieldIsCreated(vegt_field_mdl)) then + print*,"- CALL FieldDestroy FOR MODEL GRID VEGETATION TYPE." + call ESMF_FieldDestroy(vegt_field_mdl,rc=rc) + endif print*,"- CALL FieldDestroy FOR MODEL GRID LATITUDE." call ESMF_FieldDestroy(latitude_field_mdl,rc=rc) diff --git a/sorc/sfc_climo_gen.fd/output_frac_cats.F90 b/sorc/sfc_climo_gen.fd/output_frac_cats.F90 new file mode 100644 index 000000000..732eced58 --- /dev/null +++ b/sorc/sfc_climo_gen.fd/output_frac_cats.F90 @@ -0,0 +1,218 @@ +!> @file +!! @brief Write model categorical data for a single tile. +!! @author George Gayno NCEP/EMC @date 2022 + +!> Output categorical data such as vegetation type. Include +!! percentage of each category within a model grid box and +!! the dominant category. +!! +!! @author George Gayno NCEP/EMC @date 2022 + module output_frac_cats + + implicit none + + private + + public :: output_driver + + contains + +!> Driver routine to output model categorical data. +!! +!! @param[in] data_one_tile The percentage of each category within a model grid cell. +!! @param[in] dom_cat_one_tile The dominant category within a model grid cell. +!! @param[in] lat_one_tile Latitude of each model grid cell. +!! @param[in] lon_one_tile Longitude of each model grid cell. +!! @param[in] i_mdl i dimension of model grid. +!! @param[in] j_mdl j dimension of model grid. +!! @param[in] num_categories Number of categories. +!! @param[in] tile Tile number. +!! @author George Gayno @date 2022 + subroutine output_driver(data_one_tile, dom_cat_one_tile, lat_one_tile, lon_one_tile, & + i_mdl, j_mdl, num_categories, tile) + + use mpi + use esmf + use source_grid, only : field_names + use model_grid, only : grid_tiles + use program_setup, only : halo + + implicit none + + integer, intent(in) :: i_mdl, j_mdl, tile, num_categories + + real(esmf_kind_r4), intent(in) :: data_one_tile(i_mdl,j_mdl,num_categories) + real(esmf_kind_r4), intent(in) :: dom_cat_one_tile(i_mdl,j_mdl) + real(esmf_kind_r4), intent(in) :: lat_one_tile(i_mdl,j_mdl) + real(esmf_kind_r4), intent(in) :: lon_one_tile(i_mdl,j_mdl) + + character(len=200) :: out_file + character(len=200) :: out_file_with_halo + + integer :: field_idx + integer :: ierr + integer :: i_out, j_out + integer :: i_start, i_end, j_start, j_end + + field_idx = 1 + + select case (field_names(field_idx)) + case ('soil_type') + out_file = "./soil_type." // grid_tiles(tile) // ".nc" + out_file_with_halo = "./soil_type." // grid_tiles(tile) // ".halo.nc" + case ('vegetation_type') + out_file = "./vegetation_type." // grid_tiles(tile) // ".nc" + out_file_with_halo = "./vegetation_type." // grid_tiles(tile) // ".halo.nc" + case default + print*,'- FATAL ERROR IN ROUTINE OUTPUT. UNIDENTIFIED FIELD : ', field_names(field_idx) + call mpi_abort(mpi_comm_world, 67, ierr) + end select + +!---------------------------------------------------------------------- +! If user specified a halo (for running stand-alone regional grid), +! remove it. +!---------------------------------------------------------------------- + + if (halo > 0) then + print*,"- WILL WRITE WITHOUT HALO REGION OF ", halo, " ROWS/COLS." + i_start = 1 + halo + i_end = i_mdl - halo + j_start = 1 + halo + j_end = j_mdl - halo + i_out = i_end - i_start + 1 + j_out = j_end - j_start + 1 + call writeit(out_file, i_out, j_out, num_categories, & + lat_one_tile(i_start:i_end,j_start:j_end), & + lon_one_tile(i_start:i_end,j_start:j_end), & + data_one_tile(i_start:i_end,j_start:j_end,:), & + dom_cat_one_tile(i_start:i_end,j_start:j_end) ) + print*,"- WILL WRITE FULL DOMAIN INCLUDING HALO." + call writeit(out_file_with_halo, i_mdl, j_mdl, num_categories, & + lat_one_tile, lon_one_tile, data_one_tile, dom_cat_one_tile) + else + print*,"- WILL WRITE DATA." + call writeit(out_file, i_mdl, j_mdl, num_categories, & + lat_one_tile, lon_one_tile, data_one_tile, dom_cat_one_tile) + endif + + return + + end subroutine output_driver + +!> Write data to a netcdf file. +!! +!! @param[in] out_file Output file name. +!! @param[in] iout i-dimension of data. +!! @param[in] jout j-dimension of data. +!! @param[in] num_categories Number of categories. +!! @param[in] latitude Latitude of data. +!! @param[in] longitude Longitude of data. +!! @param[in] data_pct Percentage of each category in each model grid cell. +!! @param[in] dominant_cat Dominant category in each model grid cell. + subroutine writeit(out_file, iout, jout, num_categories, & + latitude, longitude, data_pct, dominant_cat) + + use esmf + use netcdf + use utils + use source_grid, only : day_of_rec, source, field_names, num_time_recs + + implicit none + + character(len=*), intent(in) :: out_file + + integer, intent(in) :: iout, jout, num_categories + + real(esmf_kind_r4), intent(in) :: latitude(iout,jout) + real(esmf_kind_r4), intent(in) :: longitude(iout,jout) + real(esmf_kind_r4), intent(in) :: data_pct(iout,jout,num_categories) + real(esmf_kind_r4), intent(in) :: dominant_cat(iout,jout) + + character(len=200) :: field_names_pct + integer :: header_buffer_val = 16384 + integer :: ncid, dim_x, dim_y, dim_z, dim_time + integer :: id_times, id_lat, id_lon, id_data_pct + integer :: id_data_dom_cat +!integer :: id_sum + integer :: error + +!real :: sum_all(iout,jout) + + print*,"- OPEN AND WRITE: ",trim(out_file) + error = nf90_create(out_file, NF90_NETCDF4, ncid) + call netcdf_err(error, 'ERROR IN NF90_CREATE' ) + error = nf90_def_dim(ncid, 'nx', iout, dim_x) + call netcdf_err(error, 'DEFINING NX DIMENSION' ) + error = nf90_def_dim(ncid, 'ny', jout, dim_y) + call netcdf_err(error, 'DEFINING NY DIMENSION' ) + error = nf90_def_dim(ncid, 'num_categories', num_categories, dim_z) + call netcdf_err(error, 'DEFINING NZ DIMENSION' ) + error = nf90_def_dim(ncid, 'time', num_time_recs, dim_time) + call netcdf_err(error, 'DEFINING TIME DIMENSION' ) + + error = nf90_def_var(ncid, 'time', NF90_FLOAT, dim_time, id_times) + call netcdf_err(error, 'DEFINING TIME VARIABLE' ) + error = nf90_put_att(ncid, id_times, "units", "days since 2015-1-1") + call netcdf_err(error, 'DEFINING TIME ATTRIBUTE' ) + if (len_trim(source) > 0) then + error = nf90_put_att(ncid, nf90_global, 'source', source) + call netcdf_err(error, 'DEFINING GLOBAL SOURCE ATTRIBUTE' ) + endif + + error = nf90_def_var(ncid, 'geolat', NF90_FLOAT, (/dim_x,dim_y/), id_lat) + call netcdf_err(error, 'DEFINING GEOLAT FIELD' ) + error = nf90_put_att(ncid, id_lat, "long_name", "Latitude") + call netcdf_err(error, 'DEFINING GEOLAT NAME ATTRIBUTE' ) + error = nf90_put_att(ncid, id_lat, "units", "degrees_north") + call netcdf_err(error, 'DEFINING GEOLAT UNIT ATTRIBUTE' ) + error = nf90_def_var(ncid, 'geolon', NF90_FLOAT, (/dim_x,dim_y/), id_lon) + call netcdf_err(error, 'DEFINING GEOLON FIELD' ) + error = nf90_put_att(ncid, id_lon, "long_name", "Longitude") + call netcdf_err(error, 'DEFINING GEOLON NAME ATTRIBUTE' ) + error = nf90_put_att(ncid, id_lon, "units", "degrees_east") + call netcdf_err(error, 'DEFINING GEOLON UNIT ATTRIBUTE' ) + + field_names_pct = trim(field_names(1)) // "_pct" + error = nf90_def_var(ncid, trim(field_names_pct), NF90_FLOAT, (/dim_x,dim_y,dim_z,dim_time/), id_data_pct) + call netcdf_err(error, 'DEFINING FIELD' ) + error = nf90_put_att(ncid, id_data_pct, "units", "percent coverage each category") + call netcdf_err(error, 'DEFINING FIELD ATTRIBUTE' ) + error = nf90_put_att(ncid, id_data_pct, "coordinates", "geolon geolat") + call netcdf_err(error, 'DEFINING COORD ATTRIBUTE' ) + + error = nf90_def_var(ncid, trim(field_names(1)), NF90_FLOAT, (/dim_x,dim_y,dim_time/), id_data_dom_cat) + call netcdf_err(error, 'DEFINING FIELD' ) + error = nf90_put_att(ncid, id_data_dom_cat, "units", "dominant category") + call netcdf_err(error, 'DEFINING FIELD ATTRIBUTE' ) + error = nf90_put_att(ncid, id_data_dom_cat, "coordinates", "geolon geolat") + call netcdf_err(error, 'DEFINING COORD ATTRIBUTE' ) + +!error = nf90_def_var(ncid, 'sum', NF90_FLOAT, (/dim_x,dim_y,dim_time/), id_sum) +!call netcdf_err(error, 'DEFINING FIELD' ) + + error = nf90_enddef(ncid, header_buffer_val,4,0,4) + + error = nf90_put_var( ncid, id_times, day_of_rec) + call netcdf_err(error, 'WRITING TIME FIELD' ) + + error = nf90_put_var( ncid, id_lat, latitude) + call netcdf_err(error, 'IN NF90_PUT_VAR FOR GEOLAT' ) + + error = nf90_put_var( ncid, id_lon, longitude) + call netcdf_err(error, 'IN NF90_PUT_VAR FOR GEOLON' ) + + error = nf90_put_var( ncid, id_data_pct, data_pct) + call netcdf_err(error, 'IN NF90_PUT_VAR' ) + + error = nf90_put_var( ncid, id_data_dom_cat, dominant_cat) + call netcdf_err(error, 'IN NF90_PUT_VAR' ) + +! Temporary output of sum of %. +!sum_all = sum(data_pct, dim=3) +!error = nf90_put_var( ncid, id_sum, sum_all) + + error = nf90_close(ncid) + + end subroutine writeit + + end module output_frac_cats diff --git a/sorc/sfc_climo_gen.fd/program_setup.f90 b/sorc/sfc_climo_gen.fd/program_setup.f90 index 608befd03..6c1a9068c 100644 --- a/sorc/sfc_climo_gen.fd/program_setup.f90 +++ b/sorc/sfc_climo_gen.fd/program_setup.f90 @@ -36,14 +36,22 @@ module program_setup character(len=500), public :: orog_dir_mdl = "NULL" !< Directory containing the model grid orography files. character(len=500), public :: orog_files_mdl(6) = "NULL" !< Model grid orography filenames. - character(len=50), public :: leaf_area_index_method='bilinear' !< Interpolation method for leaf area index. Conservative or bilinear (default). - character(len=50), public :: maximum_snow_albedo_method='bilinear' !< Interpolation method for max snow albedo. Conservative or bilinear (default). - character(len=50), public :: snowfree_albedo_method='bilinear' !< Interpolation method for snowfree albedo. Conservative or bilinear (default). - character(len=50), public :: vegetation_greenness_method='bilinear' !< Interpolation method for vegetation greenness. Conservative or bilinear (default). + character(len=50), public :: leaf_area_index_method='bilinear' !< Interpolation method for leaf area index. + !! Conservative or bilinear (default). + character(len=50), public :: maximum_snow_albedo_method='bilinear' !< Interpolation method for max snow albedo. + !! Conservative or bilinear (default). + character(len=50), public :: snowfree_albedo_method='bilinear' !< Interpolation method for snowfree albedo. + !! Conservative or bilinear (default). + character(len=50), public :: vegetation_greenness_method='bilinear' !< Interpolation method for vegetation greenness. + !! Conservative or bilinear (default). integer, public :: halo = 0 !< Number of row/cols defining the lateral !! boundary halo. Used for regional nests. + logical, public :: fract_vegsoil_type = .false. !< When true, output the percentage of each soil + !! and vegetation type category, and the dominant + !! category within a model grid box. + public :: read_setup_namelist contains @@ -69,7 +77,8 @@ subroutine read_setup_namelist(localpet) input_vegetation_greenness_file, mosaic_file_mdl, & orog_dir_mdl, orog_files_mdl, halo, & vegetation_greenness_method, leaf_area_index_method, & - maximum_snow_albedo_method, snowfree_albedo_method + maximum_snow_albedo_method, snowfree_albedo_method, & + fract_vegsoil_type print*,"- READ SETUP NAMELIST, LOCALPET: ", localpet diff --git a/sorc/sfc_climo_gen.fd/search_frac_cats.f90 b/sorc/sfc_climo_gen.fd/search_frac_cats.f90 new file mode 100644 index 000000000..09ead8929 --- /dev/null +++ b/sorc/sfc_climo_gen.fd/search_frac_cats.f90 @@ -0,0 +1,119 @@ +!> @file +!! @brief Replace undefined values on the model grid. +!! @author George Gayno @date 2022 + +!> Replace undefined values on the model grid with valid +!! values at a nearby neighbor. Undefined values are typically +!! associated with isolated islands where there is no source data. +!! Routine searches a neighborhood with a radius of 100 grid points. +!! If no valid values are found, a default value is used. This +!! routine works for one tile of a cubed sphere grid. It does +!! not consider valid values at adjacent faces. This routine +!! works for fractional categorical fields, such as soil +!! type. +!! +!! @param[inout] field - input: Field before missing values are replaced. +!! - output: Field after missing values are replaced. +!! @param[in] mask Field bitmap. Field defined where mask=1. +!! @param[in] idim i dimension of tile. +!! @param[in] jdim j dimension of tile. +!! @param[in] num_categories Number of veg/soil categories. +!! @param[in] tile Tile number. +!! @param[in] field_name Field name. +!! @author George Gayno @date 2022 + subroutine search_frac_cats (field, mask, idim, jdim, num_categories, tile, field_name) + + use mpi + use esmf + + implicit none + + integer, intent(in) :: idim, jdim, tile, num_categories + integer(esmf_kind_i4), intent(in) :: mask(idim,jdim) + + real(esmf_kind_r4), intent(inout) :: field(idim,jdim,num_categories) + + character(len=*) :: field_name + + integer :: i, j, krad, ii, jj + integer :: istart, iend + integer :: jstart, jend + integer :: ierr + integer :: default_category + + real(esmf_kind_r4), allocatable :: field_save(:,:,:) + +!----------------------------------------------------------------------- +! Set default category. +!----------------------------------------------------------------------- + + select case (field_name) + case ('soil_type') ! soil type + default_category = 3 + case ('vegetation_type') ! vegetation type + default_category = 3 + case default + print*,'- FATAL ERROR IN ROUTINE SEARCH. UNIDENTIFIED FIELD : ', field + call mpi_abort(mpi_comm_world, 77, ierr) + end select + +!----------------------------------------------------------------------- +! Perform search and replace. +!----------------------------------------------------------------------- + + allocate (field_save(idim,jdim,num_categories)) + field_save = field + + J_LOOP : do j = 1, jdim + I_LOOP : do i = 1, idim + + if (mask(i,j) == 1 .and. field_save(i,j,1) < -9999.0) then + + KRAD_LOOP : do krad = 1, 100 + + istart = i - krad + iend = i + krad + jstart = j - krad + jend = j + krad + + JJ_LOOP : do jj = jstart, jend + II_LOOP : do ii = istart, iend + +!----------------------------------------------------------------------- +! Search only along outer square. +!----------------------------------------------------------------------- + + if ((jj == jstart) .or. (jj == jend) .or. & + (ii == istart) .or. (ii == iend)) then + + if (jj < 1 .or. jj > jdim) cycle JJ_LOOP + if (ii < 1 .or. ii > idim) cycle II_LOOP + + if (mask(ii,jj) == 1 .and. maxval(field_save(ii,jj,:)) > 0.0) then + field(i,j,:) = field_save(ii,jj,:) + write(6,100) tile,i,j,ii,jj + cycle I_LOOP + endif + + endif + + enddo II_LOOP + enddo JJ_LOOP + + enddo KRAD_LOOP + + field(i,j,:) = 0.0 + field(i,j,default_category) = 1.0 ! Search failed. Use 100% of default category. + + write(6,101) tile,i,j,default_category + + endif + enddo I_LOOP + enddo J_LOOP + + deallocate(field_save) + + 100 format(1x,"- MISSING POINT TILE: ",i2," I/J: ",i5,i5," SET TO VALUE AT: ",i5,i5) + 101 format(1x,"- MISSING POINT TILE: ",i2," I/J: ",i5,i5," SET TO DEFAULT VALUE OF: ",i3) + + end subroutine search_frac_cats diff --git a/tests/sfc_climo_gen/CMakeLists.txt b/tests/sfc_climo_gen/CMakeLists.txt index 82e90207f..0067900d7 100644 --- a/tests/sfc_climo_gen/CMakeLists.txt +++ b/tests/sfc_climo_gen/CMakeLists.txt @@ -33,10 +33,20 @@ include_directories(${PROJECT_SOURCE_DIR}) execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/LSanSuppress.supp ${CMAKE_CURRENT_BINARY_DIR}/LSanSuppress.supp) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/fort.41 ${CMAKE_CURRENT_BINARY_DIR}/fort.41) + add_executable(ftst_interp ftst_interp.F90) target_link_libraries(ftst_interp sfc_climo_gen_lib) -# Cause test to be run with MPI. +add_executable(ftst_pgm_setup ftst_program_setup.F90) +target_link_libraries(ftst_pgm_setup sfc_climo_gen_lib) + +add_mpi_test(sfc_climo_gen-ftst_pgm_setup + EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_pgm_setup + NUMPROCS 1 + TIMEOUT 60) + add_mpi_test(sfc_climo_gen-ftst_interp EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_interp NUMPROCS 1 diff --git a/tests/sfc_climo_gen/data/fort.41 b/tests/sfc_climo_gen/data/fort.41 new file mode 100644 index 000000000..c42cf7df9 --- /dev/null +++ b/tests/sfc_climo_gen/data/fort.41 @@ -0,0 +1,21 @@ + &config + input_leaf_area_index_file="leaf" + input_facsf_file="facsf" + input_substrate_temperature_file="substrate_temp" + input_maximum_snow_albedo_file="maxsnow_alb" + input_snowfree_albedo_file="snowfree_alb" + input_slope_type_file="slope" + input_soil_type_file="soil_type" + input_soil_color_file="soil_color" + input_vegetation_type_file="veg_type" + input_vegetation_greenness_file="greenness" + mosaic_file_mdl="mosaic.nc" + orog_dir_mdl="./orog" + orog_files_mdl="oro.tile1.nc","oro.tile2.nc","oro.tile3.nc","oro.tile4.nc","oro.tile5.nc","oro.tile6.nc" + leaf_area_index_method="bilinear" + maximum_snow_albedo_method="conserve" + snowfree_albedo_method="bilinear" + vegetation_greenness_method="conserve" + halo=4 + fract_vegsoil_type=.true. + / diff --git a/tests/sfc_climo_gen/ftst_program_setup.F90 b/tests/sfc_climo_gen/ftst_program_setup.F90 new file mode 100644 index 000000000..e63115e8e --- /dev/null +++ b/tests/sfc_climo_gen/ftst_program_setup.F90 @@ -0,0 +1,49 @@ +! Unit test for sfc_climo_gen utility, program_setup.f90. + +program ftst_pgm_setup + use program_setup + use mpi + implicit none + + integer :: my_rank, ierr + + call mpi_init(ierr) + call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr) + + if (my_rank .eq. 0) print*,"Starting test of program_setup." + if (my_rank .eq. 0) print*,"Call routine read_setup_namelist." + + call read_setup_namelist(my_rank) + + if (trim(input_leaf_area_index_file) /= "leaf") stop 2 + if (trim(input_facsf_file) /= "facsf") stop 4 + if (trim(input_substrate_temperature_file) /= "substrate_temp") stop 6 + if (trim(input_maximum_snow_albedo_file) /= "maxsnow_alb") stop 8 + if (trim(input_snowfree_albedo_file) /= "snowfree_alb") stop 10 + if (trim(input_slope_type_file) /= "slope") stop 12 + if (trim(input_soil_type_file) /= "soil_type") stop 14 + if (trim(input_soil_color_file) /= "soil_color") stop 16 + if (trim(input_vegetation_type_file) /= "veg_type") stop 18 + if (trim(input_vegetation_greenness_file) /= "greenness") stop 20 + if (trim(mosaic_file_mdl) /= "mosaic.nc") stop 22 + if (trim(orog_dir_mdl) /= "./orog") stop 24 + if (trim(orog_files_mdl(1)) /= "oro.tile1.nc") stop 26 + if (trim(orog_files_mdl(2)) /= "oro.tile2.nc") stop 28 + if (trim(orog_files_mdl(3)) /= "oro.tile3.nc") stop 30 + if (trim(orog_files_mdl(4)) /= "oro.tile4.nc") stop 32 + if (trim(orog_files_mdl(5)) /= "oro.tile5.nc") stop 34 + if (trim(orog_files_mdl(6)) /= "oro.tile6.nc") stop 36 + if (trim(leaf_area_index_method) /= "bilinear") stop 38 + if (trim(maximum_snow_albedo_method) /= "conserve") stop 40 + if (trim(snowfree_albedo_method) /= "bilinear") stop 42 + if (trim(vegetation_greenness_method) /= "conserve") stop 44 + if (.not. fract_vegsoil_type) stop 46 + if (halo /= 4) stop 46 + + if (my_rank .eq. 0) print*, "OK" + + if (my_rank .eq. 0) print*, "SUCCESS!" + + call mpi_finalize(ierr) + +end program ftst_pgm_setup diff --git a/ush/sfc_climo_gen.sh b/ush/sfc_climo_gen.sh index 43d601c32..020ec0d13 100755 --- a/ush/sfc_climo_gen.sh +++ b/ush/sfc_climo_gen.sh @@ -25,6 +25,10 @@ # WORK_DIR Temporary working directory # SOIL_TYPE_FILE Path/name of input soil type data. # VEG_TYPE_FILE Path/name of input vegetation type data. +# vegsoilt_frac When true, outputs dominant soil and +# vegetation type category and the +# fractional value of each category. +# When false, outputs dominant category. #------------------------------------------------------------------------- set -eux @@ -39,6 +43,7 @@ FIX_FV3=${FIX_FV3:-/scratch4/NCEPDEV/global/save/glopara/git/fv3gfs/fix/fix_fv3_ input_sfc_climo_dir=${input_sfc_climo_dir:?} mosaic_file=${mosaic_file:-$FIX_FV3/C${res}_mosaic.nc} HALO=${HALO:-0} +vegsoilt_frac=${vegsoilt_frac:-.false.} veg_type_src=${veg_type_src:-"modis.igbp.0.05"} VEG_TYPE_FILE=${VEG_TYPE_FILE:-${input_sfc_climo_dir}/vegetation_type.${veg_type_src}.nc} soil_type_src=${soil_type_src:-"statsgo.0.05"} @@ -80,6 +85,7 @@ halo=$HALO maximum_snow_albedo_method="bilinear" snowfree_albedo_method="bilinear" vegetation_greenness_method="bilinear" +fract_vegsoil_type=${vegsoilt_frac} / EOF diff --git a/util/sfc_climo_gen/run.wcoss2.sh b/util/sfc_climo_gen/run.wcoss2.sh index f8ecc9588..a96eceeec 100755 --- a/util/sfc_climo_gen/run.wcoss2.sh +++ b/util/sfc_climo_gen/run.wcoss2.sh @@ -13,7 +13,7 @@ #PBS -A GFS-DEV #PBS -N grid_fv3 #PBS -l walltime=00:10:00 -#PBS -l select=1:ncpus=24:mem=75GB +#PBS -l select=1:ncpus=24:mem=250GB set -x @@ -79,6 +79,12 @@ export veg_type_src="modis.igbp.0.05" export soil_type_src="statsgo.0.05" +export vegsoilt_frac='.false.' # When true, outputs percent of each + # soil and veg type category and a + # dominate category. When false, only + # outputs the dominate category. A + # Fortran logical, so include the dots. + #------------------------------------- # Set working directory and directory where output files will be saved. #------------------------------------- From 5101237fb568ee68d8ee3fe900456ac86054d509 Mon Sep 17 00:00:00 2001 From: Larissa Reames <52886575+LarissaReames-NOAA@users.noreply.github.com> Date: Wed, 12 Apr 2023 09:19:57 -0500 Subject: [PATCH 2/3] Fix bug in chgres_cube subroutine search_many (#808) The routine is called from all MPI tasks, but some arrays passed to the routine were not allocated on all tasks. This problem was only seen when compiling/running under 'Debug' mode (all regression tests crashed at the call to search_many). Fixes #797. --- sorc/chgres_cube.fd/surface.F90 | 84 +++++++++++-------- .../chgres_cube/ftst_surface_search_many.F90 | 19 ++--- 2 files changed, 59 insertions(+), 44 deletions(-) diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index 26066e362..b77f83f80 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -851,11 +851,12 @@ subroutine interp(localpet) if (localpet == 0) then where(mask_target_one_tile == 1) mask_target_one_tile = 0 where(mask_target_one_tile == 2) mask_target_one_tile = 1 + call search_many(num_fields,bundle_seaice_target,tile, search_nums,localpet, & + mask=mask_target_one_tile) + else + call search_many(num_fields,bundle_seaice_target, tile,search_nums,localpet) endif - - call search_many(num_fields,bundle_seaice_target,data_one_tile, mask_target_one_tile,tile,search_nums,localpet, & - field_data_3d=data_one_tile_3d) enddo deallocate(search_nums) @@ -977,10 +978,12 @@ subroutine interp(localpet) allocate(water_target_one_tile(i_target,j_target)) water_target_one_tile = 0 where(mask_target_one_tile == 0) water_target_one_tile = 1 - endif - call search_many(num_fields,bundle_water_target,data_one_tile, water_target_one_tile,& - tile,search_nums,localpet,latitude=latitude_one_tile) + call search_many(num_fields,bundle_water_target, tile,search_nums,localpet, & + latitude=latitude_one_tile,mask=water_target_one_tile) + else + call search_many(num_fields,bundle_water_target, tile,search_nums,localpet) + endif if (localpet == 0) deallocate(water_target_one_tile) @@ -1068,10 +1071,12 @@ subroutine interp(localpet) allocate(land_target_one_tile(i_target,j_target)) land_target_one_tile = 0 where(mask_target_one_tile == 1) land_target_one_tile = 1 - endif - call search_many(num_fields,bundle_allland_target,data_one_tile, land_target_one_tile,& - tile,search_nums,localpet) + call search_many(num_fields,bundle_allland_target, & + tile,search_nums,localpet, mask=land_target_one_tile) + else + call search_many(num_fields,bundle_allland_target, tile,search_nums,localpet) + endif if (localpet == 0) deallocate(land_target_one_tile) enddo @@ -1202,8 +1207,12 @@ subroutine interp(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather", rc) - call search_many(num_fields,bundle_landice_target,data_one_tile, land_target_one_tile,& - tile,search_nums,localpet,terrain_land=data_one_tile2,field_data_3d=data_one_tile_3d) + if (localpet==0) then + call search_many(num_fields,bundle_landice_target,tile,search_nums,localpet,& + terrain_land=data_one_tile2,mask=land_target_one_tile) + else + call search_many(num_fields,bundle_landice_target,tile,search_nums,localpet) + endif enddo deallocate (veg_type_target_one_tile) @@ -1416,9 +1425,12 @@ subroutine interp(localpet) call ESMF_FieldGather(soil_type_target_grid, data_one_tile2, rootPet=0,tile=tile, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldGather", rc) - - call search_many(num_fields,bundle_nolandice_target,data_one_tile, mask_target_one_tile,& - tile,search_nums,localpet,soilt_climo=data_one_tile2, field_data_3d=data_one_tile_3d) + if (localpet==0) then + call search_many(num_fields,bundle_nolandice_target,tile,search_nums,localpet, & + soilt_climo=data_one_tile2, mask=mask_target_one_tile) + else + call search_many(num_fields,bundle_nolandice_target, tile,search_nums,localpet) + endif print*,"- CALL FieldGather FOR TARGET GRID TOTAL SOIL MOISTURE, TILE: ", tile call ESMF_FieldGather(soilm_tot_target_grid, data_one_tile_3d, rootPet=0, tile=tile, rc=rc) @@ -3290,20 +3302,18 @@ end subroutine regrid_many !! !! @param[in] num_field Number of fields to process. !! @param[inout] bundle_target ESMF FieldBundle holding target fields to search -!! @param[inout] field_data_2d A real array of size i_target,j_target to temporarily hold data for searching -!! @param[inout] mask An integer array of size i_target,j_target that holds masked (0) and unmasked (1) -!! values indicating where to execute search (only at unmasked points). !! @param[in] tile Current cubed sphere tile. !! @param[inout] search_nums Array length num_field holding search field numbers corresponding to each field provided for searching. !! @param[in] localpet ESMF local persistent execution thread. !! @param[in] latitude (optional) A real array size i_target,j_target of latitude on the target grid !! @param[in] terrain_land (optional) A real array size i_target,j_target of terrain height (m) on the target grid !! @param[in] soilt_climo (optional) A real array size i_target,j_target of climatological soil type on the target grid -!! @param[in] field_data_3d (optional) An empty real array of size i_target,j_target,lsoil_target to temporarily hold soil data for searching +!! @param[inout] mask (optional) An integer array of size i_target,j_target that holds masked (0) and unmasked (1) +!! values indicating where to execute search (only at +!unmasked points). !! @author Larissa Reames, OU CIMMS/NOAA/NSSL - subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & - search_nums,localpet,latitude,terrain_land,soilt_climo,& - field_data_3d) + subroutine search_many(num_field,bundle_target,tile,search_nums,localpet,latitude, & + terrain_land,soilt_climo, mask) use model_grid, only : i_target,j_target, lsoil_target use program_setup, only : external_model, input_type @@ -3313,14 +3323,14 @@ subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & integer, intent(in) :: num_field type(esmf_fieldbundle), intent(inout) :: bundle_target - real(esmf_kind_r8), intent(inout) :: field_data_2d(i_target,j_target) - real(esmf_kind_r8), intent(inout), optional :: field_data_3d(i_target,j_target,lsoil_target) + real(esmf_kind_r8), intent(inout), optional :: latitude(i_target,j_target) real(esmf_kind_r8), intent(inout), optional :: terrain_land(i_target,j_target) real(esmf_kind_r8), intent(inout), optional :: soilt_climo(i_target,j_target) - integer(esmf_kind_i8), intent(inout) :: mask(i_target,j_target) + integer(esmf_kind_i8), intent(inout), optional :: mask(i_target,j_target) - + real(esmf_kind_r8), allocatable :: field_data_2d(:,:) + real(esmf_kind_r8), allocatable :: field_data_3d(:,:,:) integer, intent(in) :: tile,localpet integer, intent(inout) :: search_nums(num_field) @@ -3331,6 +3341,7 @@ subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & integer, parameter :: TERRAIN_FIELD_NUM= 7 integer :: j,k, rc, ndims + do k = 1,num_field call ESMF_FieldBundleGet(bundle_target,k,temp_field, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& @@ -3338,39 +3349,37 @@ subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & call ESMF_FieldGet(temp_field, name=fname, dimcount=ndims,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldGet", rc) + if (localpet==0) then + allocate(field_data_2d(i_target,j_target)) + else + allocate(field_data_2d(0,0)) + endif if (ndims .eq. 2) then - print*, "processing 2d field ", trim(fname) - print*, "FieldGather" call ESMF_FieldGather(temp_field,field_data_2d,rootPet=0,tile=tile, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldGather", rc) if (localpet == 0) then if (present(latitude) .and. search_nums(k).eq.SST_FIELD_NUM) then ! Sea surface temperatures; pass latitude field to search - print*, "search1" call search(field_data_2d, mask, i_target, j_target, tile,search_nums(k),latitude=latitude) elseif (present(terrain_land) .and. search_nums(k) .eq. TERRAIN_FIELD_NUM) then ! Terrain height; pass optional climo terrain array to search - print*, "search2" call search(field_data_2d, mask, i_target, j_target, tile,search_nums(k),terrain_land=terrain_land) elseif (search_nums(k) .eq. SOTYP_LAND_FIELD_NUM) then ! Soil type over land if (fname .eq. "soil_type_target_grid") then ! Soil type over land when interpolating input data to target grid ! *with* the intention of retaining interpolated data in output - print*, "search3" call search(field_data_2d, mask, i_target, j_target, tile,search_nums(k),soilt_climo=soilt_climo) elseif (present(soilt_climo)) then if (maxval(field_data_2d) > 0 .and. (trim(external_model) .ne. "GFS" .or. trim(input_type) .ne. "grib2")) then ! Soil type over land when interpolating input data to target grid ! *without* the intention of retaining data in output file - print*, "search4" call search(field_data_2d, mask, i_target, j_target, tile, search_nums(k)) else ! If no soil type field exists in input data (e.g., GFS grib2) then don't search ! but simply set data to the climo field. This may result in ! somewhat inaccurate soil moistures as no scaling will occur - print*, "search5" field_data_2d = soilt_climo endif !check field value endif !sotype from target grid @@ -3384,12 +3393,17 @@ subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldScatter", rc) else + if (localpet==0) then + allocate(field_data_3d(i_target,j_target,lsoil_target)) + else + allocate(field_data_3d(0,0,0)) + endif + ! Process 3d fields soil temperature, moisture, and liquid - print*, "FieldGather" call ESMF_FieldGather(temp_field,field_data_3d,rootPet=0,tile=tile,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldGather", rc) - print*, "processing 3d field ", trim(fname) + if (localpet==0) then do j = 1, lsoil_target field_data_2d = field_data_3d(:,:,j) @@ -3400,7 +3414,9 @@ subroutine search_many(num_field,bundle_target,field_data_2d,mask, tile, & call ESMF_FieldScatter(temp_field, field_data_3d, rootPet=0, tile=tile,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldScatter", rc) + deallocate(field_data_3d) endif !ndims + deallocate(field_data_2d) end do !fields end subroutine search_many diff --git a/tests/chgres_cube/ftst_surface_search_many.F90 b/tests/chgres_cube/ftst_surface_search_many.F90 index f8700ba9f..2c7520ca0 100644 --- a/tests/chgres_cube/ftst_surface_search_many.F90 +++ b/tests/chgres_cube/ftst_surface_search_many.F90 @@ -295,8 +295,8 @@ program surface_interp input_type="restart" !Call the search many routine to test search and replace - call search_many(num_fields,bundle_search1,dummy_2d,mask_target_search,1,field_nums,localpet, & - soilt_climo=soilt_climo) + call search_many(num_fields,bundle_search1,1,field_nums,localpet, & + soilt_climo=soilt_climo,mask=mask_target_search) call ESMF_FieldBundleDestroy(bundle_search1,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -342,8 +342,8 @@ program surface_interp external_model="HRRR" !Call the search many routine to test search and replace - call search_many(num_fields,bundle_search2,dummy_2d,mask_target_search,1,field_nums,localpet, & - soilt_climo=soilt_climo) + call search_many(num_fields,bundle_search2,1,field_nums,localpet, & + soilt_climo=soilt_climo,mask=mask_target_search) call ESMF_FieldBundleDestroy(bundle_search2,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -385,8 +385,8 @@ program surface_interp allocate(field_nums(num_fields)) field_nums = (/11,7,224/) !Call the search many routine to test some branches of default behavior - call search_many(num_fields,bundle_default1,dummy_2d,mask_default,1,field_nums,localpet, & - latitude=latitude_default,terrain_land=terrain_land,soilt_climo=soilt_climo) + call search_many(num_fields,bundle_default1,1,field_nums,localpet, & + latitude=latitude_default,terrain_land=terrain_land,soilt_climo=soilt_climo,mask=mask_default) print*,"Check results for bundle_default1." @@ -441,8 +441,8 @@ program surface_interp input_type="grib2" external_model="GFS" !Call the search many routine to test behavior for GFS grib2 soil type - call search_many(num_fields,bundle_default2,dummy_2d,mask_default,1,field_nums,localpet,& - soilt_climo=soilt_climo) + call search_many(num_fields,bundle_default2,1,field_nums,localpet,& + soilt_climo=soilt_climo,mask=mask_default) call ESMF_FieldBundleDestroy(bundle_default2,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& @@ -478,8 +478,7 @@ program surface_interp field_nums(:) = (/21/) !Call the search many routine to test behavior for GFS grib2 soil type - call search_many(num_fields,bundle_3d_search,dummy_2d,mask_target_search,1,field_nums,localpet,& - field_data_3d=dummy_3d) + call search_many(num_fields,bundle_3d_search,1,field_nums,localpet,mask=mask_target_search) call ESMF_FieldBundleDestroy(bundle_3d_search,rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& From ebbe14b0da399ef3b9ec14b03c824569ec2dff91 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Wed, 12 Apr 2023 11:00:33 -0400 Subject: [PATCH 3/3] Run sfc_climo_gen utility on Hera (#791) Add script to run the utility on Hera. Move common logic and configuration variables to a new script - sfc_gen.sh. Fixes #788. --- util/sfc_climo_gen/readme.md | 6 +- util/sfc_climo_gen/run.hera.sh | 30 ++++++++ util/sfc_climo_gen/run.wcoss2.sh | 96 ++--------------------- util/sfc_climo_gen/sfc_gen.sh | 127 +++++++++++++++++++++++++++++++ 4 files changed, 167 insertions(+), 92 deletions(-) create mode 100755 util/sfc_climo_gen/run.hera.sh create mode 100755 util/sfc_climo_gen/sfc_gen.sh diff --git a/util/sfc_climo_gen/readme.md b/util/sfc_climo_gen/readme.md index 95745a1ec..1606199f9 100644 --- a/util/sfc_climo_gen/readme.md +++ b/util/sfc_climo_gen/readme.md @@ -1,4 +1,8 @@ Run the sfc_climo_gen program stand-alone on WCOSS2 using -pre-exiting 'grid' and 'orography' files. +pre-exiting 'grid' and 'orography' files. See the +sfc_gen.sh script for details. Outputs surface fields such as soil and vegetation type. + +Set the configuration variables in sfc_gen.sh. Then +run using the machine specific driver script. diff --git a/util/sfc_climo_gen/run.hera.sh b/util/sfc_climo_gen/run.hera.sh new file mode 100755 index 000000000..697d0f350 --- /dev/null +++ b/util/sfc_climo_gen/run.hera.sh @@ -0,0 +1,30 @@ +#!/bin/bash + +#------------------------------------------------------------ +# Run the sfc_climo_gen program stand-alone on Hera using +# pre-exiting 'grid' and 'orography' files. See the +# sfc_gen.sh script for details. +# +# Set the configuration variables in sfc_gen.sh. Then +# run this script as follows: 'sbatch $script' +#------------------------------------------------------------ + +#SBATCH -J sfc_climo_gen +#SBATCH -A fv3-cpu +#SBATCH --open-mode=truncate +#SBATCH -o log +#SBATCH -e log +#SBATCH --nodes=1 --ntasks-per-node=24 +#SBATCH --partition=bigmem +#SBATCH -q debug +#SBATCH -t 00:10:00 + +set -x + +export APRUN_SFC="srun" + +export BASE_DIR=$SLURM_SUBMIT_DIR/../.. + +$SLURM_SUBMIT_DIR/sfc_gen.sh + +exit diff --git a/util/sfc_climo_gen/run.wcoss2.sh b/util/sfc_climo_gen/run.wcoss2.sh index a96eceeec..ad767015d 100755 --- a/util/sfc_climo_gen/run.wcoss2.sh +++ b/util/sfc_climo_gen/run.wcoss2.sh @@ -2,9 +2,11 @@ #------------------------------------------------------------ # Run the sfc_climo_gen program stand-alone on WCOSS2 using -# pre-exiting 'grid' and 'orography files. +# pre-exiting 'grid' and 'orography' files. See the +# sfc_gen.sh script for details. # -# To run, type: 'qsub $script' +# Set the configuration variables in sfc_gen.sh. Then, +# run this script as follows: 'qsub $script' #------------------------------------------------------------ #PBS -o log @@ -22,94 +24,6 @@ export APRUN_SFC="mpiexec -n 24 -ppn 24 -cpu-bind core" export BASE_DIR=$PBS_O_WORKDIR/../.. -source ${BASE_DIR}/sorc/machine-setup.sh > /dev/null 2>&1 -module use ${BASE_DIR}/modulefiles -module load build.$target.intel -module list - -#------------------------------------- -# Set model resolution. -#------------------------------------- - -export res=384 - -#------------------------------------- -# Where the model "grid", "mosaic" and "oro" files reside. -#------------------------------------- - -export FIX_FV3=${BASE_DIR}/fix/orog/C${res} - -#------------------------------------- -# Uncomment for regional grids. -#------------------------------------- - -##HALO=3 -##export GRIDTYPE=regional - -#------------------------------------------------------------- -# Choose which soil type and vegetation type data to use. -# -# For viirs-based vegetation type data, set to: -# 1) "viirs.igbp.0.1" for global 0.10-deg data -# 2) "viirs.igbp.0.05" for global 0.05-deg data -# 3) "viirs.igbp.0.03" for global 0.03-deg data -# 4) "viirs.igbp.conus.30s" for CONUS 30s data -# 4) "viirs.igbp.nh.30s" for NH 30s data -# 4) "viirs.igbp.30s" for global 30s data -# -# For the modis-based vegetation data, set to: -# 1) "modis.igbp.0.05" for global 0.05-deg data -# 2) "modis.igbp.0.03" for global 0.03-deg data -# 3) "modis.igbp.conus.30s" for CONUS 30s data -# 4) "modis.igbp.nh.30s" for NH 30s data -# 5) "modis.igbp.30s" for global 30s data -# -# For STATSGO soil type data -# 1) "statsgo.0.05" for global 0.05-deg data -# 2) "statsgo.0.03" for global 0.03-deg data -# 3) "statsgo.conus.30s" for CONUS 30s data -# 4) "statsgo.nh.30s" for NH 30s data -# 5) "statsgo.30s" for global 30s data -# -# For Beijing Norm. Univ. soil type data -# 1) "bnu.30s" for global 30s data -#------------------------------------------------------------- - -export veg_type_src="modis.igbp.0.05" - -export soil_type_src="statsgo.0.05" - -export vegsoilt_frac='.false.' # When true, outputs percent of each - # soil and veg type category and a - # dominate category. When false, only - # outputs the dominate category. A - # Fortran logical, so include the dots. - -#------------------------------------- -# Set working directory and directory where output files will be saved. -#------------------------------------- - -export WORK_DIR=/lfs/h2/emc/stmp/$LOGNAME/work.sfc -export SAVE_DIR=/lfs/h2/emc/stmp/$LOGNAME/sfc.C${res} - -#------------------------------------- -# Should not have to touch anything below here. -#------------------------------------- - -if [[ $GRIDTYPE = "regional" ]]; then - HALO=$(( $HALO + 1 )) - export HALO - ln -fs $FIX_FV3/C${res}_grid.tile7.halo${HALO}.nc $FIX_FV3/C${res}_grid.tile7.nc - ln -fs $FIX_FV3/C${res}_oro_data.tile7.halo${HALO}.nc $FIX_FV3/C${res}_oro_data.tile7.nc -fi - -export input_sfc_climo_dir=${BASE_DIR}/fix/sfc_climo - -ulimit -a -ulimit -s unlimited - -rm -fr $WORK_DIR $SAVE_DIR - -${BASE_DIR}/ush/sfc_climo_gen.sh +$PBS_O_WORKDIR/sfc_gen.sh exit diff --git a/util/sfc_climo_gen/sfc_gen.sh b/util/sfc_climo_gen/sfc_gen.sh new file mode 100755 index 000000000..9ba7b747c --- /dev/null +++ b/util/sfc_climo_gen/sfc_gen.sh @@ -0,0 +1,127 @@ +#!/bin/bash + +#----------------------------------------------------------------------- +# +# This script is run by the machine specific driver script. +# +# Set the following variables: +# +# res - Grid resolution. Example: 384 or 384.mx025. +# +# FIX_FV3 - Location of the pre-existing 'grid' and 'orography' +# files. Defaults to ${BASE_DIR}/fix/orog/C${res}, where +# BASE_DIR is the location of the checked out repository. +# +# The required files are: +# +# 'mosaic' file - C${res}_mosaic.nc (Note: 'res' without +# the 'mx' extension.) +# +# 'grid' files - C${res}_grid.tile7.halo${HALO}.nc (regional grids). +# C${res}_grid.tile[1-6].nc (global grids). +# Note: 'res' without the 'mx' extension. +# +# 'orog' files - C${res}_oro_data.tile7.halo${HALO}.nc (regional grids). +# C${res}_oro_data.tile[1-6].nc (global grids). +# +# GRIDTYPE - set to 'regional' for regional grids. Otherwise, +# comment out. +# +# FIX_REG - For regional grids. Hold links to the 'grid' and 'orog' files +# with names expected by the program. +# +# HALO - The number of halo rows/cols. Only for regional grids. +# Otherwise, comment out. +# +# WORK_DIR - Working directory. +# +# SAVE_DIR - Directory where the surface files will be saved. +# +# veg_type_src - Input vegetation type data. Choices are: +# For viirs-based vegetation type data, set to: +# - "viirs.igbp.0.1" for global 0.10-deg data +# - "viirs.igbp.0.05" for global 0.05-deg data +# - "viirs.igbp.0.03" for global 0.03-deg data +# - "viirs.igbp.conus.30s" for CONUS 30s data +# - "viirs.igbp.nh.30s" for NH 30s data +# - "viirs.igbp.30s" for global 30s data +# For the modis-based vegetation data, set to: +# - "modis.igbp.0.05" for global 0.05-deg data +# - "modis.igbp.0.03" for global 0.03-deg data +# - "modis.igbp.conus.30s" for CONUS 30s data +# - "modis.igbp.nh.30s" for NH 30s data +# - "modis.igbp.30s" for global 30s data +# +# soil_type_src - Input soil type data. Choices are: +# For STATSGO soil type data +# - "statsgo.0.05" for global 0.05-deg data +# - "statsgo.0.03" for global 0.03-deg data +# - "statsgo.conus.30s" for CONUS 30s data +# - "statsgo.nh.30s" for NH 30s data +# - "statsgo.30s" for global 30s data +# For Beijing Norm. Univ. soil type data +# - "bnu.30s" for global 30s data +# +# vegsoilt_frac - When .true., output the fraction of each +# vegetation and soil type and the dominant +# category. When .false., output dominant +# category only. +#----------------------------------------------------------------------- + +set -x + +#export res=96 +export res=96.mx100 + +#HALO=4 +#export GRIDTYPE=regional +#FIX_REG=/lfs/h2/emc/stmp/$LOGNAME/fix.reg + +export veg_type_src="modis.igbp.0.05" + +export soil_type_src="statsgo.0.05" + +export WORK_DIR=/lfs/h2/emc/stmp/$LOGNAME/work.sfc +export SAVE_DIR=/lfs/h2/emc/stmp/$LOGNAME/sfc.C${res} + +export FIX_FV3=${BASE_DIR}/fix/orog/C${res} + +export vegsoilt_frac=.true. + +#------------------------------------------------------------------------ +#------------------------------------------------------------------------ +# Should not have to touch anything below here. +#------------------------------------------------------------------------ +#------------------------------------------------------------------------ + +if [[ "$GRIDTYPE" = "regional" ]]; then + mkdir -p $FIX_REG + ln -fs $FIX_FV3/C${res}_grid.tile7.halo${HALO}.nc $FIX_REG/C${res}_grid.tile7.halo${HALO}.nc + ln -fs $FIX_FV3/C${res}_oro_data.tile7.halo${HALO}.nc $FIX_REG/C${res}_oro_data.tile7.nc + ln -fs $FIX_FV3/C${res}_mosaic.nc $FIX_REG/C${res}_mosaic.nc + export mosaic_file=$FIX_REG/C${res}_mosaic.nc + export FIX_FV3=$FIX_REG + HALO=$(( $HALO + 1 )) + export HALO +else + res2=${res//".mx"*} + export mosaic_file=$FIX_FV3/C${res2}_mosaic.nc +fi + +export input_sfc_climo_dir=${BASE_DIR}/fix/sfc_climo + +ulimit -a +ulimit -s unlimited + +source ${BASE_DIR}/sorc/machine-setup.sh > /dev/null 2>&1 +module use ${BASE_DIR}/modulefiles +module load build.$target.intel +module list + +rm -fr $WORK_DIR $SAVE_DIR + +export APRUN + +${BASE_DIR}/ush/sfc_climo_gen.sh + +exit