From f8dc806e710bbe380d6894f1412eeddaa28c8548 Mon Sep 17 00:00:00 2001 From: JeffBeck-NOAA <55201531+JeffBeck-NOAA@users.noreply.github.com> Date: Thu, 15 Oct 2020 09:31:11 -0600 Subject: [PATCH 1/2] Update chgres_cube to ingest HRRR, RAP, and NAM grib2 files Update chgres_cube to ingest HRRR, RAP and NAM grib2 files. Add UFS_UTILS build modules for machines Cheyenne, Stampede2, and Odin. In the build script and all utility/regression test scripts, do a 'module load' of the machine-specific build modules instead of sourcing them. This is more portable. --- build_all.sh | 20 +- driver_scripts/driver_grid.cray.sh | 3 +- driver_scripts/driver_grid.dell.sh | 3 +- driver_scripts/driver_grid.hera.sh | 3 +- driver_scripts/driver_grid.jet.sh | 3 +- driver_scripts/driver_grid.orion.sh | 3 +- modulefiles/build.cheyenne | 30 + modulefiles/build.jet | 13 +- modulefiles/build.odin | 41 + modulefiles/build.orion | 4 +- modulefiles/build.stampede | 41 + modulefiles/build.wcoss_cray | 13 +- modulefiles/build.wcoss_dell_p3 | 3 +- parm/varmap_tables/GSDphys_var_map.txt | 27 + reg_tests/chgres_cube/c96.regional.sh | 3 + reg_tests/chgres_cube/driver.cray.sh | 3 +- reg_tests/chgres_cube/driver.dell.sh | 4 +- reg_tests/chgres_cube/driver.hera.sh | 4 +- reg_tests/chgres_cube/driver.jet.sh | 4 +- reg_tests/chgres_cube/driver.orion.sh | 3 +- reg_tests/global_cycle/driver.cray.sh | 3 +- reg_tests/global_cycle/driver.dell.sh | 4 +- reg_tests/global_cycle/driver.hera.sh | 4 +- reg_tests/global_cycle/driver.jet.sh | 4 +- reg_tests/global_cycle/driver.orion.sh | 4 +- reg_tests/grid_gen/driver.cray.sh | 3 +- reg_tests/grid_gen/driver.dell.sh | 4 +- reg_tests/grid_gen/driver.hera.sh | 4 +- reg_tests/grid_gen/driver.jet.sh | 4 +- reg_tests/grid_gen/driver.orion.sh | 4 +- reg_tests/ice_blend/driver.cray.sh | 3 +- reg_tests/ice_blend/driver.dell.sh | 4 +- reg_tests/ice_blend/driver.hera.sh | 4 +- reg_tests/ice_blend/driver.jet.sh | 4 +- reg_tests/ice_blend/driver.orion.sh | 4 +- reg_tests/snow2mdl/driver.cray.sh | 3 +- reg_tests/snow2mdl/driver.dell.sh | 4 +- reg_tests/snow2mdl/driver.hera.sh | 4 +- reg_tests/snow2mdl/driver.jet.sh | 4 +- reg_tests/snow2mdl/driver.orion.sh | 4 +- sorc/chgres_cube.fd/atmosphere.F90 | 30 +- sorc/chgres_cube.fd/grib2_util.F90 | 21 + sorc/chgres_cube.fd/input_data.F90 | 973 +++++++++++++++++---- sorc/chgres_cube.fd/model_grid.F90 | 407 ++++++++- sorc/chgres_cube.fd/program_setup.f90 | 145 +++- sorc/chgres_cube.fd/search_util.f90 | 50 +- sorc/chgres_cube.fd/surface.F90 | 1088 +++++++++++++++++++++++- sorc/chgres_cube.fd/write_data.F90 | 49 +- sorc/emcsfc_snow2mdl.fd/.gitignore | 3 - sorc/global_chgres.fd/.gitignore | 1 - sorc/machine-setup.sh | 9 +- sorc/nst_tf_chg.fd/.gitignore | 3 - sorc/sfc_climo_gen.fd/.gitignore | 3 - util/gdas_init/driver.cray.sh | 3 +- util/gdas_init/driver.dell.sh | 4 +- util/gdas_init/driver.hera.sh | 4 +- util/vcoord_gen/run.cray.sh | 3 +- 57 files changed, 2791 insertions(+), 309 deletions(-) create mode 100644 modulefiles/build.cheyenne create mode 100644 modulefiles/build.odin create mode 100644 modulefiles/build.stampede create mode 100644 parm/varmap_tables/GSDphys_var_map.txt delete mode 100644 sorc/emcsfc_snow2mdl.fd/.gitignore delete mode 100644 sorc/global_chgres.fd/.gitignore delete mode 100644 sorc/nst_tf_chg.fd/.gitignore delete mode 100644 sorc/sfc_climo_gen.fd/.gitignore diff --git a/build_all.sh b/build_all.sh index 64d159c1d..273b40f1e 100755 --- a/build_all.sh +++ b/build_all.sh @@ -3,20 +3,22 @@ set -eux target=${target:-"NULL"} -if [[ "$target" == "linux.gnu" || "$target" == "linux.intel" ]]; then +export MOD_PATH + +if [[ "$target" == "linux.*" || "$target" == "macosx.*" ]]; then unset -f module + set +x + source ./modulefiles/build.$target > /dev/null 2>&1 + set -x else set +x - source ./sorc/machine-setup.sh > /dev/null 2>&1 + source ./sorc/machine-setup.sh + module use ./modulefiles + module load build.$target > /dev/null 2>&1 + module list set -x fi -export MOD_PATH -set +x -source ./modulefiles/build.$target > /dev/null 2>&1 -module list -set -x - # --- Build all programs. # @@ -26,7 +28,7 @@ cd ./build CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=../ -DEMC_EXEC_DIR=ON" -if [[ "$target" != "wcoss_cray" ]]; then +if [[ "$target" != "wcoss_cray" && "$target" != "odin" ]]; then CMAKE_FLAGS+=" -DCMAKE_Fortran_COMPILER=ifort -DCMAKE_C_COMPILER=icc" fi diff --git a/driver_scripts/driver_grid.cray.sh b/driver_scripts/driver_grid.cray.sh index f2faa4a78..1aa1088bb 100755 --- a/driver_scripts/driver_grid.cray.sh +++ b/driver_scripts/driver_grid.cray.sh @@ -57,7 +57,8 @@ #----------------------------------------------------------------------- source ../sorc/machine-setup.sh > /dev/null 2>&1 -source ../modulefiles/build.$target +module use ../modulefiles +module load build.$target module list #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.dell.sh b/driver_scripts/driver_grid.dell.sh index d2d6d08d0..760c44151 100755 --- a/driver_scripts/driver_grid.dell.sh +++ b/driver_scripts/driver_grid.dell.sh @@ -59,7 +59,8 @@ #----------------------------------------------------------------------- source ../sorc/machine-setup.sh > /dev/null 2>&1 -source ../modulefiles/build.$target +module use ../modulefiles +module load build.$target module list #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh index e00c7de1e..6fbe90c37 100755 --- a/driver_scripts/driver_grid.hera.sh +++ b/driver_scripts/driver_grid.hera.sh @@ -59,7 +59,8 @@ set -x source ../sorc/machine-setup.sh > /dev/null 2>&1 -source ../modulefiles/build.$target +module use ../modulefiles +module load build.$target module list #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh index 2df3f01a5..88e7e7cda 100755 --- a/driver_scripts/driver_grid.jet.sh +++ b/driver_scripts/driver_grid.jet.sh @@ -60,7 +60,8 @@ set -x source ../sorc/machine-setup.sh > /dev/null 2>&1 -source ../modulefiles/build.$target +module use ../modulefiles +module load build.$target module list #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh index 9eebd32a7..97dce226c 100755 --- a/driver_scripts/driver_grid.orion.sh +++ b/driver_scripts/driver_grid.orion.sh @@ -59,7 +59,8 @@ set -x source ../sorc/machine-setup.sh > /dev/null 2>&1 -source ../modulefiles/build.$target +module use ../modulefiles +module load build.$target module list #----------------------------------------------------------------------- diff --git a/modulefiles/build.cheyenne b/modulefiles/build.cheyenne new file mode 100644 index 000000000..474981369 --- /dev/null +++ b/modulefiles/build.cheyenne @@ -0,0 +1,30 @@ +#%Module##################################################### +## Build and run module for Cheyenne +############################################################# + +module purge +module load ncarenv/1.3 +module load intel/19.1.1 +module load mpt/2.19 +module load ncarcompilers/0.5.0 +module load cmake/3.16.4 + +module use -a /glade/p/ral/jntp/GMTB/tools/NCEPLIBS-ufs-v2.0.0/intel-19.1.1/mpt-2.19/modules + +module load bacio/2.4.1 +module load g2/3.4.1 +module load ip/3.3.3 +module load nemsio/2.5.2 +module load sp/2.3.3 +module load w3emc/2.7.3 +module load w3nco/2.4.1 +module load sigio/2.3.2 + +module load sfcio/1.4.1 +module load gfsio/1.4.1 +module load nemsiogfs/2.5.3 +module load landsfcutil/2.4.1 +module load wgrib2/2.0.8 +module load netcdf/4.7.4 + +setenv ESMFMKFILE /glade/p/ral/jntp/GMTB/tools/NCEPLIBS-ufs-v2.0.0/intel-19.1.1/mpt-2.19/lib64/esmf.mk diff --git a/modulefiles/build.jet b/modulefiles/build.jet index b81a52cd0..e0156505f 100644 --- a/modulefiles/build.jet +++ b/modulefiles/build.jet @@ -8,12 +8,12 @@ module load szip/2.1 module load hdf5/1.10.5 module load netcdf/4.7.0 -export CC=icc -export FC=ifort -export CXX=icpc +setenv CC icc +setenv FC ifort +setenv CXX icpc -export ESMFMKFILE=/lfs4/HFIP/hfv3gfs/software/NCEPLIBS-ufs-v2.0.0beta01/intel-18.0.5.274/impi-2018.4.274/lib64/esmf.mk -export Jasper_ROOT=/lfs4/HFIP/hfv3gfs/software/NCEPLIBS-ufs-v2.0.0beta01/intel-18.0.5.274/impi-2018.4.274 +setenv ESMFMKFILE /lfs4/HFIP/hfv3gfs/software/NCEPLIBS-ufs-v2.0.0beta01/intel-18.0.5.274/impi-2018.4.274/lib64/esmf.mk +setenv Jasper_ROOT /lfs4/HFIP/hfv3gfs/software/NCEPLIBS-ufs-v2.0.0beta01/intel-18.0.5.274/impi-2018.4.274 module use /lfs4/HFIP/hfv3gfs/software/NCEPLIBS-ufs-v2.0.0beta01/intel-18.0.5.274/impi-2018.4.274/modules module load w3nco/2.4.1 @@ -29,6 +29,3 @@ module load gfsio/1.4.1 module load landsfcutil/2.4.1 module load g2/3.4.1 module load wgrib2/2.0.8 - - - diff --git a/modulefiles/build.odin b/modulefiles/build.odin new file mode 100644 index 000000000..5be328674 --- /dev/null +++ b/modulefiles/build.odin @@ -0,0 +1,41 @@ +#%Module##################################################### +## Build and run module for Odin +############################################################# + +module load craype/2.6.2 +module load craype-ivybridge +module load PrgEnv-intel +module swap intel/19.0.5.281 +module load cray-mpich/7.7.10 +module load cray-libsci +module load cray-netcdf-hdf5parallel +module load cray-parallel-netcdf +module load cray-hdf5-parallel + +setenv NETCDF /opt/cray/pe/netcdf-hdf5parallel/4.6.3.2/INTEL/19.0 + +#module use -a /oldscratch/ywang/external/modulefiles +#module use /oldscratch/ywang/external/NCEPLIBS_SRW/modules +module use /oldscratch/ywang/external/NCEPLIBS_SRWv2.0/modules +module load w3nco +module load w3emc +module load sp +module load ip +module load bacio +module load sigio +module load sfcio +module load nemsio +module load nemsiogfs +module load gfsio +module load landsfcutil +module load g2 +module load wgrib2 + +#module load esmf/8.0.0 +#setenv ESMFMKFILE /oldscratch/ywang/external/NCEPLIBS_SRW/lib64/esmf.mk +setenv ESMFMKFILE /oldscratch/ywang/external/NCEPLIBS_SRWv2.0/lib64/esmf.mk + +setenv CMAKE_Fortran_COMPILER ftn +setenv CMAKE_C_COMPILER cc + +#setenv WGRIB2_ROOT /oldscratch/ywang/external/NCEPLIBS_SRWv2.0/wgrib2-2.0.8 diff --git a/modulefiles/build.orion b/modulefiles/build.orion index 5c623989a..c2c30d6f7 100644 --- a/modulefiles/build.orion +++ b/modulefiles/build.orion @@ -1,4 +1,4 @@ -############################################################# +#%Module##################################################### ## Build module for Orion ############################################################# @@ -21,7 +21,7 @@ module load nemsiogfs/2.5.3 module load landsfcutil/2.4.1 module load wgrib2/2.0.8 -export Jasper_ROOT="/apps/jasper-1.900.1" +setenv Jasper_ROOT /apps/jasper-1.900.1 module use -a /apps/contrib/NCEPLIBS/lib/modulefiles module load netcdfp/4.7.4.release diff --git a/modulefiles/build.stampede b/modulefiles/build.stampede new file mode 100644 index 000000000..a26bfae82 --- /dev/null +++ b/modulefiles/build.stampede @@ -0,0 +1,41 @@ +#%Module##################################################### +## Build and run module for Stampede +############################################################# + +module purge + +module load libfabric/1.7.0 +module load git/2.24.1 +module load autotools/1.1 +module load xalt/2.8 +module load TACC + +module load python3/3.7.0 +module load intel/18.0.2 +module load cmake/3.16.1 +module load impi/18.0.2 +module load pnetcdf/1.11.0 +module load netcdf/4.6.2 +module li + +setenv NETCDF /opt/apps/intel18/netcdf/4.6.2/x86_64 + +#module use /work/00315/tg455890/stampede2/regional_fv3/ufs_testing/INSTALL/modules +module use /work/00315/tg455890/stampede2/regional_fv3/NCEPLIBS_SRWv2.0/modules +module load esmf/8.0.0 + +module load w3nco +module load w3emc +module load sp +module load ip +module load bacio +module load sigio +module load sfcio +module load nemsio +module load nemsiogfs +module load gfsio +module load landsfcutil +module load g2 +module load wgrib2 + +#setenv ESMFMKFILE /work/00315/tg455890/stampede2/regional_fv3/NCEPLIBS_SRWv2.0/lib64/esmf.mk diff --git a/modulefiles/build.wcoss_cray b/modulefiles/build.wcoss_cray index b6ce5aa57..c54daea65 100644 --- a/modulefiles/build.wcoss_cray +++ b/modulefiles/build.wcoss_cray @@ -31,14 +31,13 @@ module load nemsiogfs/2.5.3 module load landsfcutil/2.4.1 module load wgrib2/2.0.8 -export ZLIB_ROOT=/usrx/local/prod/zlib/1.2.7/intel/haswell -export PNG_ROOT=/usrx/local/prod//png/1.2.49/intel/haswell -export Jasper_ROOT=/usrx/local/prod/jasper/1.900.1/intel/haswell +setenv ZLIB_ROOT /usrx/local/prod/zlib/1.2.7/intel/haswell +setenv PNG_ROOT /usrx/local/prod/png/1.2.49/intel/haswell +setenv Jasper_ROOT /usrx/local/prod/jasper/1.900.1/intel/haswell module use /gpfs/hps3/emc/nems/noscrub/emc.nemspara/soft/modulefiles -module load esmf/8.0.0 -export ESMFMKFILE=/gpfs/hps3/emc/nems/noscrub/emc.nemspara/soft/esmf/8.0.0/lib/esmf.mk -export NETCDF=/opt/cray/netcdf/4.3.3.1/INTEL/14.0 +#module load esmf/8.0.0 +setenv ESMFMKFILE /gpfs/hps3/emc/nems/noscrub/emc.nemspara/soft/esmf/8.0.0/lib/esmf.mk +setenv NETCDF /opt/cray/netcdf/4.3.3.1/INTEL/14.0 module rm gcc module load gcc/6.3.0 - diff --git a/modulefiles/build.wcoss_dell_p3 b/modulefiles/build.wcoss_dell_p3 index ae9d55e2c..bb9093b03 100644 --- a/modulefiles/build.wcoss_dell_p3 +++ b/modulefiles/build.wcoss_dell_p3 @@ -8,9 +8,8 @@ module load cmake/3.16.2 module load ips/18.0.1.163 module load impi/18.0.1 - module load jasper/1.900.29 -export Jasper_ROOT="/usrx/local/prod/packages/gnu/4.8.5/jasper/1.900.29" +setenv Jasper_ROOT /usrx/local/prod/packages/gnu/4.8.5/jasper/1.900.29 module use /usrx/local/nceplibs/dev/NCEPLIBS/modulefiles module load netcdf_parallel/4.7.4 diff --git a/parm/varmap_tables/GSDphys_var_map.txt b/parm/varmap_tables/GSDphys_var_map.txt new file mode 100644 index 000000000..22f6b47f5 --- /dev/null +++ b/parm/varmap_tables/GSDphys_var_map.txt @@ -0,0 +1,27 @@ +dzdt dzdt set_to_fill 0 D +delta_p delp skip 0 D +sphum sphum set_to_fill 1E-7 T +liq_wat liq_wat set_to_fill 0 T +o3mr o3mr set_to_fill 1E-7 T +rainwat rainwat set_to_fill 0 T +ice_wat ice_wat set_to_fill 0 T +snowwat snowwat set_to_fill 0 T +graupel graupel set_to_fill 0 T +ice_nc ice_nc set_to_fill -1.0 T +rain_nc rain_nc set_to_fill -1.0 T +water_nc water_nc set_to_fill -1.0 T +liq_aero liq_aero set_to_fill 0 T +ice_aero ice_aero set_to_fill 0 T +sgs_tke sgs_tke set_to_fill 0 T +vtype vtype skip 0 S +sotype stype skip 0 S +vfrac vfrac skip 0 S +fricv uustar skip 0 S +sfcr zorl set_to_fill 0.01 S +soilw smc stop 0 S +soilt stc stop 0 S +cnwat cnwat set_to_fill 0.0 S +icetk icetk set_to_fill 265.0 S +weasd weasd set_to_fill 0.0 S +snod snod set_to_fill 0.0 S +tprcp tprcp set_to_fill 0.0 S diff --git a/reg_tests/chgres_cube/c96.regional.sh b/reg_tests/chgres_cube/c96.regional.sh index cecb199cc..e67df060b 100755 --- a/reg_tests/chgres_cube/c96.regional.sh +++ b/reg_tests/chgres_cube/c96.regional.sh @@ -60,6 +60,9 @@ fi cd $DATA +mv out.sfc.tile7.nc out.sfc.tile1.nc +mv out.atm.tile7.nc out.atm.tile1.nc + test_failed=0 for files in *.nc do diff --git a/reg_tests/chgres_cube/driver.cray.sh b/reg_tests/chgres_cube/driver.cray.sh index cae83b182..ee3f4d424 100755 --- a/reg_tests/chgres_cube/driver.cray.sh +++ b/reg_tests/chgres_cube/driver.cray.sh @@ -23,7 +23,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list export OUTDIR=/gpfs/hps3/stmp/$LOGNAME/chgres_reg_tests diff --git a/reg_tests/chgres_cube/driver.dell.sh b/reg_tests/chgres_cube/driver.dell.sh index 4cbd48361..97247cbfb 100755 --- a/reg_tests/chgres_cube/driver.dell.sh +++ b/reg_tests/chgres_cube/driver.dell.sh @@ -23,7 +23,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export OUTDIR=/gpfs/dell1/stmp/$LOGNAME/chgres_reg_tests QUEUE="debug" diff --git a/reg_tests/chgres_cube/driver.hera.sh b/reg_tests/chgres_cube/driver.hera.sh index c61b72c87..b6f78c461 100755 --- a/reg_tests/chgres_cube/driver.hera.sh +++ b/reg_tests/chgres_cube/driver.hera.sh @@ -25,7 +25,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export OUTDIR=/scratch2/NCEPDEV/stmp1/$LOGNAME/chgres_reg_tests PROJECT_CODE="fv3-cpu" diff --git a/reg_tests/chgres_cube/driver.jet.sh b/reg_tests/chgres_cube/driver.jet.sh index a9ee38707..c093eb7cb 100755 --- a/reg_tests/chgres_cube/driver.jet.sh +++ b/reg_tests/chgres_cube/driver.jet.sh @@ -25,7 +25,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export OUTDIR=/lfs4/HFIP/emcda/$LOGNAME/stmp/chgres_reg_tests PROJECT_CODE="hfv3gfs" diff --git a/reg_tests/chgres_cube/driver.orion.sh b/reg_tests/chgres_cube/driver.orion.sh index 479eedb96..0886871ad 100755 --- a/reg_tests/chgres_cube/driver.orion.sh +++ b/reg_tests/chgres_cube/driver.orion.sh @@ -25,7 +25,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list export OUTDIR=/work/noaa/stmp/$LOGNAME/chgres_reg_tests diff --git a/reg_tests/global_cycle/driver.cray.sh b/reg_tests/global_cycle/driver.cray.sh index ffd5048dd..288c38edf 100755 --- a/reg_tests/global_cycle/driver.cray.sh +++ b/reg_tests/global_cycle/driver.cray.sh @@ -28,7 +28,8 @@ #BSUB -extsched 'CRAYLINUX[]' source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list export DATA=/gpfs/hps3/stmp/$LOGNAME/reg_tests.cycle diff --git a/reg_tests/global_cycle/driver.dell.sh b/reg_tests/global_cycle/driver.dell.sh index 383124240..178355667 100755 --- a/reg_tests/global_cycle/driver.dell.sh +++ b/reg_tests/global_cycle/driver.dell.sh @@ -33,7 +33,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA=/gpfs/dell1/stmp/$LOGNAME/reg_tests.cycle diff --git a/reg_tests/global_cycle/driver.hera.sh b/reg_tests/global_cycle/driver.hera.sh index ecc05cd9b..9ee1cb64a 100755 --- a/reg_tests/global_cycle/driver.hera.sh +++ b/reg_tests/global_cycle/driver.hera.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA=/scratch2/NCEPDEV/stmp1/$LOGNAME/reg_tests.cycle diff --git a/reg_tests/global_cycle/driver.jet.sh b/reg_tests/global_cycle/driver.jet.sh index 74183c565..aa2d7285a 100755 --- a/reg_tests/global_cycle/driver.jet.sh +++ b/reg_tests/global_cycle/driver.jet.sh @@ -31,7 +31,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA=/lfs4/HFIP/emcda/$LOGNAME/stmp/reg_tests.cycle diff --git a/reg_tests/global_cycle/driver.orion.sh b/reg_tests/global_cycle/driver.orion.sh index 63c51f9f7..5731b9063 100755 --- a/reg_tests/global_cycle/driver.orion.sh +++ b/reg_tests/global_cycle/driver.orion.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA=/work/noaa/stmp/$LOGNAME/reg_tests.cycle diff --git a/reg_tests/grid_gen/driver.cray.sh b/reg_tests/grid_gen/driver.cray.sh index 535ce28f4..e2f8e5927 100755 --- a/reg_tests/grid_gen/driver.cray.sh +++ b/reg_tests/grid_gen/driver.cray.sh @@ -20,7 +20,8 @@ #----------------------------------------------------------------------------- source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list set -x diff --git a/reg_tests/grid_gen/driver.dell.sh b/reg_tests/grid_gen/driver.dell.sh index 8027b84e3..5c2ae6a7f 100755 --- a/reg_tests/grid_gen/driver.dell.sh +++ b/reg_tests/grid_gen/driver.dell.sh @@ -20,7 +20,9 @@ #----------------------------------------------------------------------------- source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/grid_gen/driver.hera.sh b/reg_tests/grid_gen/driver.hera.sh index f2a8536cf..0dc66996c 100755 --- a/reg_tests/grid_gen/driver.hera.sh +++ b/reg_tests/grid_gen/driver.hera.sh @@ -22,7 +22,9 @@ #----------------------------------------------------------------------------- source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/grid_gen/driver.jet.sh b/reg_tests/grid_gen/driver.jet.sh index 5c1facea3..d38c08f8a 100755 --- a/reg_tests/grid_gen/driver.jet.sh +++ b/reg_tests/grid_gen/driver.jet.sh @@ -22,7 +22,9 @@ #----------------------------------------------------------------------------- source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/grid_gen/driver.orion.sh b/reg_tests/grid_gen/driver.orion.sh index 12e56d21b..cdad89801 100755 --- a/reg_tests/grid_gen/driver.orion.sh +++ b/reg_tests/grid_gen/driver.orion.sh @@ -22,7 +22,9 @@ #----------------------------------------------------------------------------- source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/ice_blend/driver.cray.sh b/reg_tests/ice_blend/driver.cray.sh index f2904b0df..c34cb3bca 100755 --- a/reg_tests/ice_blend/driver.cray.sh +++ b/reg_tests/ice_blend/driver.cray.sh @@ -29,7 +29,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list export DATA=/gpfs/hps3/stmp/$LOGNAME/reg_tests.ice_blend diff --git a/reg_tests/ice_blend/driver.dell.sh b/reg_tests/ice_blend/driver.dell.sh index 9ca96d131..1296c7436 100755 --- a/reg_tests/ice_blend/driver.dell.sh +++ b/reg_tests/ice_blend/driver.dell.sh @@ -27,7 +27,9 @@ #BSUB -P GFS-DEV source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/ice_blend/driver.hera.sh b/reg_tests/ice_blend/driver.hera.sh index 4866e7974..63967711a 100755 --- a/reg_tests/ice_blend/driver.hera.sh +++ b/reg_tests/ice_blend/driver.hera.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/scratch2/NCEPDEV/stmp1/$LOGNAME/reg_test.ice_blend" diff --git a/reg_tests/ice_blend/driver.jet.sh b/reg_tests/ice_blend/driver.jet.sh index bc55aa272..160b21a94 100755 --- a/reg_tests/ice_blend/driver.jet.sh +++ b/reg_tests/ice_blend/driver.jet.sh @@ -29,7 +29,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/lfs4/HFIP/emcda/$LOGNAME/stmp/reg_test.ice_blend" diff --git a/reg_tests/ice_blend/driver.orion.sh b/reg_tests/ice_blend/driver.orion.sh index 188ede2bb..df2b108e1 100755 --- a/reg_tests/ice_blend/driver.orion.sh +++ b/reg_tests/ice_blend/driver.orion.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/work/noaa/stmp/$LOGNAME/reg_test.ice_blend" diff --git a/reg_tests/snow2mdl/driver.cray.sh b/reg_tests/snow2mdl/driver.cray.sh index f39c5b3ce..064b8a7d8 100755 --- a/reg_tests/snow2mdl/driver.cray.sh +++ b/reg_tests/snow2mdl/driver.cray.sh @@ -29,7 +29,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list export DATA=/gpfs/hps3/stmp/$LOGNAME/reg_tests.snow2mdl diff --git a/reg_tests/snow2mdl/driver.dell.sh b/reg_tests/snow2mdl/driver.dell.sh index 610f726d8..2b3e603d8 100755 --- a/reg_tests/snow2mdl/driver.dell.sh +++ b/reg_tests/snow2mdl/driver.dell.sh @@ -27,7 +27,9 @@ #BSUB -P GFS-DEV source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list set -x diff --git a/reg_tests/snow2mdl/driver.hera.sh b/reg_tests/snow2mdl/driver.hera.sh index cd0457104..8c074f66c 100755 --- a/reg_tests/snow2mdl/driver.hera.sh +++ b/reg_tests/snow2mdl/driver.hera.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/scratch2/NCEPDEV/stmp1/$LOGNAME/reg_tests.snow2mdl" diff --git a/reg_tests/snow2mdl/driver.jet.sh b/reg_tests/snow2mdl/driver.jet.sh index c4679031e..e12a69790 100755 --- a/reg_tests/snow2mdl/driver.jet.sh +++ b/reg_tests/snow2mdl/driver.jet.sh @@ -29,7 +29,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/lfs4/HFIP/emcda/$LOGNAME/stmp/reg_tests.snow2mdl" diff --git a/reg_tests/snow2mdl/driver.orion.sh b/reg_tests/snow2mdl/driver.orion.sh index fbc331510..be991b359 100755 --- a/reg_tests/snow2mdl/driver.orion.sh +++ b/reg_tests/snow2mdl/driver.orion.sh @@ -30,7 +30,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list export DATA="/work/noaa/stmp/$LOGNAME/reg_tests.snow2mdl" diff --git a/sorc/chgres_cube.fd/atmosphere.F90 b/sorc/chgres_cube.fd/atmosphere.F90 index a76ea36b7..f372d5dbb 100644 --- a/sorc/chgres_cube.fd/atmosphere.F90 +++ b/sorc/chgres_cube.fd/atmosphere.F90 @@ -164,7 +164,7 @@ subroutine atmosphere_driver(localpet) real(esmf_kind_r8), parameter :: exponent = rd*lapse/grav real(esmf_kind_r8), parameter :: one_over_exponent = 1.0 / exponent - real(esmf_kind_r8), pointer :: psptr(:,:) + real(esmf_kind_r8), pointer :: psptr(:,:), tempptr(:,:,:) !----------------------------------------------------------------------------------- ! Read atmospheric fields on the input grid. @@ -212,7 +212,6 @@ subroutine atmosphere_driver(localpet) temp_b4adj_target_grid, & polemethod=ESMF_POLEMETHOD_ALLAVG, & srctermprocessing=isrctermprocessing, & - extrapmethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & routehandle=regrid_bl, & regridmethod=method, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -256,7 +255,25 @@ subroutine atmosphere_driver(localpet) 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) + + nullify(tempptr) + print*,"- CALL FieldGet FOR INPUT GRID VERTICAL VEL." + call ESMF_FieldGet(dzdt_input_grid, & + farrayPtr=tempptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*, "MIN MAX W INPUT = ", minval(tempptr), maxval(tempptr) + nullify(tempptr) + print*,"- CALL FieldGet FOR VERTICAL VEL B4ADJ." + call ESMF_FieldGet(dzdt_b4adj_target_grid, & + farrayPtr=tempptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*, "MIN MAX W B4ADJ = ", minval(tempptr), maxval(tempptr) + nullify(psptr) print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE." call ESMF_FieldGet(ps_input_grid, & @@ -364,8 +381,8 @@ subroutine atmosphere_driver(localpet) wind_w_target_grid, & polemethod=ESMF_POLEMETHOD_ALLAVG, & srctermprocessing=isrctermprocessing, & - extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & routehandle=regrid_bl, & + extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & regridmethod=method, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldRegridStore", rc) @@ -391,8 +408,8 @@ subroutine atmosphere_driver(localpet) wind_s_target_grid, & polemethod=ESMF_POLEMETHOD_ALLAVG, & srctermprocessing=isrctermprocessing, & - extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & routehandle=regrid_bl, & + extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & regridmethod=method, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldRegridStore", rc) @@ -722,6 +739,7 @@ subroutine convert_winds enddo enddo enddo + print*,"- CALL FieldGet FOR 3-D WIND_W." call ESMF_FieldGet(wind_w_target_grid, & @@ -1221,7 +1239,6 @@ subroutine horiz_interp_thomp_mp_climo qnifa_climo_b4adj_target_grid, & polemethod=ESMF_POLEMETHOD_ALLAVG, & srctermprocessing=isrctermprocessing, & - extrapmethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, & routehandle=regrid_bl, & regridmethod=method, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -1478,7 +1495,7 @@ SUBROUTINE VINTG C1(:,:,:,1) = WIND1PTR(:,:,:,1) C1(:,:,:,2) = WIND1PTR(:,:,:,2) C1(:,:,:,3) = WIND1PTR(:,:,:,3) - + print*,"- CALL FieldGet FOR VERTICAL VELOCITY." call ESMF_FieldGet(dzdt_b4adj_target_grid, & farrayPtr=DZDT1PTR, rc=rc) @@ -1486,6 +1503,7 @@ SUBROUTINE VINTG call error_handler("IN FieldGet", rc) C1(:,:,:,4) = DZDT1PTR(:,:,:) + print*,"MIN MAX W TARGETB4 IN VINTG = ", minval(DZDT1PTR(:,:,:)), maxval(DZDT1PTR(:,:,:)) print*,"- CALL FieldGet FOR 3-D TEMP." call ESMF_FieldGet(temp_b4adj_target_grid, & diff --git a/sorc/chgres_cube.fd/grib2_util.F90 b/sorc/chgres_cube.fd/grib2_util.F90 index 64a1cdb71..195d5c71d 100644 --- a/sorc/chgres_cube.fd/grib2_util.F90 +++ b/sorc/chgres_cube.fd/grib2_util.F90 @@ -74,4 +74,25 @@ subroutine convert_omega(omega,p,t,q,clb,cub) end subroutine convert_omega +function to_upper(strIn) result(strOut) +! Adapted from http://www.star.le.ac.uk/~cgp/fortran.html (25 May 2012) +! Original author: Clive Page + + implicit none + + character(len=*), intent(in) :: strIn + character(len=len(strIn)) :: strOut + integer :: i,j + + do i = 1, len(strIn) + j = iachar(strIn(i:i)) + if (j>= iachar("a") .and. j<=iachar("z") ) then + strOut(i:i) = achar(iachar(strIn(i:i))-32) + else + strOut(i:i) = strIn(i:i) + end if + end do + +end function to_upper + end module grib2_util diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index f5847fe6d..9011a1e1b 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -39,7 +39,12 @@ module input_data orog_files_input_grid, & tracers_input, num_tracers, & input_type, tracers, & - get_var_cond, read_from_input + get_var_cond, read_from_input, & + geogrid_file_input_grid, & + external_model, & + vgfrc_from_climo, & + minmax_vgfrc_from_climo, & + lai_from_climo use model_grid, only : input_grid, & i_input, j_input, & @@ -47,7 +52,7 @@ module input_data num_tiles_input_grid, & latitude_input_grid, & longitude_input_grid, & - inv_file + inv_file!, the_file_hrrr implicit none @@ -98,9 +103,13 @@ module input_data type(esmf_field), public :: ustar_input_grid ! fric velocity type(esmf_field), public :: veg_type_input_grid ! vegetation type type(esmf_field), public :: z0_input_grid ! roughness length + type(esmf_field), public :: veg_greenness_input_grid ! vegetation fraction + type(esmf_field), public :: lai_input_grid ! leaf area index + type(esmf_field), public :: max_veg_greenness_input_grid ! shdmax + type(esmf_field), public :: min_veg_greenness_input_grid ! shdmin - integer, public :: lsoil_input=4 ! # of soil layers, - ! # hardwire for now + integer, public :: lsoil_input=4 ! # of soil layers, no longer hardwired to allow + ! # for 7 layers of soil for the RUC LSM character(len=50), private, allocatable :: slevs(:) @@ -537,6 +546,42 @@ subroutine read_input_sfc_data(localpet) ungriddedUBound=(/lsoil_input/), rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate", rc) + + + + if (.not. vgfrc_from_climo) then + print*,"- CALL FieldCreate FOR INPUT VEGETATION GREENNESS." + veg_greenness_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + endif + + if (.not. minmax_vgfrc_from_climo) then + print*,"- CALL FieldCreate FOR INPUT MIN VEGETATION GREENNESS." + min_veg_greenness_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, 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 INPUT MAX VEGETATION GREENNESS." + max_veg_greenness_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldCreate", rc) + endif + + if (.not. lai_from_climo) then + print*,"- CALL FieldCreate FOR INPUT LEAF AREA INDEX." + lai_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldCreate", rc) + endif !------------------------------------------------------------------------------- ! Read the tiled 'warm' restart files. @@ -2432,11 +2477,13 @@ subroutine read_input_atm_grib2_file(localpet) logical :: lret logical :: conv_omega=.false., & - hasspfh=.true. + hasspfh=.true., & + isnative=.false. real(esmf_kind_r8), allocatable :: rlevs(:) real(esmf_kind_r4), allocatable :: dummy2d(:,:) - real(esmf_kind_r8), allocatable :: dummy3d(:,:,:), dummy2d_8(:,:) + real(esmf_kind_r8), allocatable :: dummy3d(:,:,:), dummy2d_8(:,:),& + u_tmp_3d(:,:,:), v_tmp_3d(:,:,:) real(esmf_kind_r8), pointer :: presptr(:,:,:), psptr(:,:),tptr(:,:,:), & qptr(:,:,:), wptr(:,:,:), & uptr(:,:,:), vptr(:,:,:) @@ -2472,17 +2519,24 @@ subroutine read_input_atm_grib2_file(localpet) if (.not.lret) call error_handler("OPENING GRIB2 ATM FILE.", iret) print*,"- READ VERTICAL COORDINATE." - iret = grb2_inq(the_file,inv_file,":var_0_2","_0_0:"," hybrid level:") + iret = grb2_inq(the_file,inv_file,":var0_2","_0_0:"," hybrid ") if (iret <= 0) then lvl_str = "mb:" lvl_str_space = " mb:" lvl_str_space_len = 4 + isnative = .false. iret = grb2_inq(the_file,inv_file,":UGRD:",lvl_str_space) lev_input=iret if (localpet == 0) print*,"- DATA IS ON ", lev_input, " ISOBARIC LEVELS." else - call error_handler("HYBRID VERTICAL COORD DATA NOT SUPPORTED", -1) + lvl_str = " level:" + lvl_str_space = " hybrid " + lvl_str_space_len = 7 + isnative = .true. + iret = grb2_inq(the_file,inv_file,":UGRD:",lvl_str_space, " level:") + if (iret < 0) call error_handler("READING VERTICAL LEVEL TYPE.", iret) + lev_input=iret endif allocate(slevs(lev_input)) @@ -2500,8 +2554,8 @@ subroutine read_input_atm_grib2_file(localpet) read(metadata(j:k),*) rlevs(i) - slevs(i) = metadata(j-1:k) - rlevs(i) = rlevs(i) * 100.0 + slevs(i) = metadata(j-1:k) + if (.not. isnative) rlevs(i) = rlevs(i) * 100.0 if (localpet==0) print*, "- LEVEL = ", slevs(i) enddo @@ -2509,28 +2563,30 @@ subroutine read_input_atm_grib2_file(localpet) call quicksort(rlevs,1,lev_input) - do i = 1,lev_input - write(slevs(i),"(F20.10)") rlevs(i)/100.0 - len_str = len_trim(slevs(i)) + if (.not. isnative) then + do i = 1,lev_input + write(slevs(i),"(F20.10)") rlevs(i)/100.0 + len_str = len_trim(slevs(i)) - do while (slevs(i)(len_str:len_str) .eq. '0') - slevs(i) = slevs(i)(:len_str-1) - len_str = len_str - 1 - end do + do while (slevs(i)(len_str:len_str) .eq. '0') + slevs(i) = slevs(i)(:len_str-1) + len_str = len_str - 1 + end do - if (slevs(i)(len_str:len_str) .eq. '.') then + if (slevs(i)(len_str:len_str) .eq. '.') then slevs(i) = slevs(i)(:len_str-1) len_str = len_str - 1 - end if + end if - slevs(i) = trim(slevs(i)) - - slevs(i) = ":"//trim(adjustl(slevs(i)))//" mb:" - if (localpet==0) print*, "- LEVEL AFTER SORT = ",slevs(i) - enddo + slevs(i) = trim(slevs(i)) + slevs(i) = ":"//trim(adjustl(slevs(i)))//" mb:" + if (localpet==0) print*, "- LEVEL AFTER SORT = ",slevs(i) + enddo + endif + if (localpet == 0) print*,"- FIND SPFH OR RH IN FILE" - iret = grb2_inq(the_file,inv_file,trac_names_grib_1(1),trac_names_grib_2(1),lvl_str_space) + iret = grb2_inq(the_file,inv_file,trim(trac_names_grib_1(1)),trac_names_grib_2(1),lvl_str_space) if (iret <= 0) then iret = grb2_inq(the_file,inv_file, ':var0_2','_1_1:',lvl_str_space) @@ -2541,8 +2597,53 @@ subroutine read_input_atm_grib2_file(localpet) else if (localpet == 0) print*,"- FILE CONTAINS SPFH." endif + + if (localpet == 0) print*,"- FIND ICMR, SCLIWC, OR CICE IN FILE" + iret = grb2_inq(the_file,inv_file,trac_names_grib_1(4),trac_names_grib_2(4),lvl_str_space) + + if (iret <= 0) then + vname = trac_names_vmap(4) + print*, "vname = ", vname + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + this_field_var_name=tmpstr,loc=varnum) + iret = grb2_inq(the_file,inv_file, ':var0_2','_1_84:',lvl_str_space) + if (iret <= 0) then + iret = grb2_inq(the_file,inv_file, ':var0_2','_6_0:',lvl_str_space) + if (iret <= 0 ) then + call handle_grib_error(vname, slevs(1),method,value,varnum,rc,var=dummy2d) + else + trac_names_grib_2(4) = '_6_0' + if (localpet == 0) print*,"- FILE CONTAINS CICE." + endif + else + trac_names_grib_2(4)='_1_84:' + if (localpet == 0) print*,"- FILE CONTAINS SCLIWC." + endif + else + if (localpet == 0) print*,"- FILE CONTAINS ICMR." + endif + + if (localpet == 0) print*,"- FIND CLWMR or SCLLWC IN FILE" + iret = grb2_inq(the_file,inv_file,trac_names_grib_1(5),trac_names_grib_2(5),lvl_str_space) + + if (iret <= 0) then + vname = trac_names_vmap(5) + print*, "vname = ", vname + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + this_field_var_name=tmpstr,loc=varnum) + iret = grb2_inq(the_file,inv_file, ':var0_2','_1_83:',lvl_str_space) + if (iret <= 0) then + call handle_grib_error(vname, slevs(1),method,value,varnum,rc,var=dummy2d) + elseif (iret <=0 .and. rc .ne. 1) then + call error_handler("READING CLOUD WATER VARIABLE.", iret) + else + trac_names_grib_2(4)='_1_83:' + if (localpet == 0) print*,"- FILE CONTAINS SCLLWC." + endif + else + if (localpet == 0) print*,"- FILE CONTAINS CLWMR." + endif - print*,"- COUNT NUMBER OF TRACERS TO BE READ IN BASED ON PHYSICS SUITE TABLE" do n = 1, num_tracers vname = tracers_input(n) @@ -2556,7 +2657,7 @@ subroutine read_input_atm_grib2_file(localpet) enddo - if (localpet==0) print*, "- NUMBER OF TRACERS IN FILE = ", num_tracers + if (localpet==0) print*, "- NUMBER OF TRACERS TO BE PROCESSED = ", num_tracers !--------------------------------------------------------------------------- ! Initialize esmf atmospheric fields. @@ -2574,10 +2675,11 @@ subroutine read_input_atm_grib2_file(localpet) allocate(dummy3d(0,0,0)) endif -!----------------------------------------------------------------------- -! Fields in non-native files read in from top to bottom. We will -! flip indices later. This program expects bottom to top. -!----------------------------------------------------------------------- +!---------------------------------------------------------------------------------- +! This program expects field levels from bottom to top. Fields in non-native +! files read in from top to bottom. We will flip indices later. Fields on +! native vertical coordinates read from bottom to top so those need no adjustments. +!---------------------------------------------------------------------------------- if (localpet == 0) then print*,"- READ TEMPERATURE." @@ -2609,6 +2711,7 @@ subroutine read_input_atm_grib2_file(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) endif + if (localpet == 0) then vname = trim(tracers_input_grib_1(n)) vname2 = trim(tracers_input_grib_2(n)) @@ -2623,8 +2726,6 @@ subroutine read_input_atm_grib2_file(localpet) trim(vname2) == ":14:192:") then call error_handler("READING IN "//trim(vname)//" AT LEVEL "//trim(slevs(vlev))& //". SET A FILL VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret) - else - exit endif endif endif @@ -2645,51 +2746,20 @@ subroutine read_input_atm_grib2_file(localpet) enddo - if (localpet==0) then - do vlev = 1, lev_input - - vname = ":var0_2" - vname2 = "_2_2:" - iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) - if (iret<=0) then - call error_handler("READING UWIND AT LEVEL "//trim(slevs(vlev)),iret) - endif - - print*, 'max, min U ', minval(dummy2d), maxval(dummy2d) - dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) - - enddo - endif +call read_winds(the_file,inv_file,u_tmp_3d,v_tmp_3d, localpet) if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT U-WIND." - call ESMF_FieldScatter(u_input_grid, dummy3d, rootpet=0, rc=rc) + call ESMF_FieldScatter(u_input_grid, u_tmp_3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) - if (localpet==0) then - do vlev = 1, lev_input - - vname = ":var0_2" - vname2 = "_2_3:" - iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) - if (iret<=0) then - call error_handler("READING VWIND AT LEVEL "//trim(slevs(vlev)),iret) - endif - - print*, 'max, min V ', minval(dummy2d), maxval(dummy2d) - dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) - - enddo - endif - if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT V-WIND." - call ESMF_FieldScatter(v_input_grid, dummy3d, rootpet=0, rc=rc) + call ESMF_FieldScatter(v_input_grid, v_tmp_3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) if (localpet == 0) then print*,"- READ SURFACE PRESSURE." - !vname = ":PRES:" vname = ":var0_2" vname2 = "_3_0:" vlevtyp = ":surface:" @@ -2708,18 +2778,14 @@ subroutine read_input_atm_grib2_file(localpet) vname = "dzdt" call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & loc=varnum) - !vname = ":DZDT:" vname = ":var0_2" vname2 = "_2_9:" do vlev = 1, lev_input iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) if (iret <= 0 ) then print*,"DZDT not available at level ", trim(slevs(vlev)), " so checking for VVEL" - !vname = ":VVEL:" vname2 = "_2_8:" iret = grb2_inq(the_file,inv_file,vname,vname2,slevs(vlev),data2=dummy2d) - - if (iret <= 0) then call handle_grib_error(vname, slevs(vlev),method,value,varnum,iret,var=dummy2d) if (iret==1) then ! missing_var_method == skip @@ -2742,7 +2808,6 @@ subroutine read_input_atm_grib2_file(localpet) if (localpet == 0) then print*,"- READ TERRAIN." - !vname = ":HGT:" vname = ":var0_2" vname2 = "_3_5:" vlevtyp = ":surface:" @@ -2756,72 +2821,73 @@ subroutine read_input_atm_grib2_file(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) - deallocate(dummy2d, dummy3d, dummy2d_8) + deallocate(dummy2d, dummy2d_8) -!--------------------------------------------------------------------------- -! Flip 'z' indices to all 3-d variables. Data is read in from model -! top to surface. This program expects surface to model top. -!--------------------------------------------------------------------------- - - if (localpet == 0) print*,"- CALL FieldGet FOR SURFACE PRESSURE." - nullify(psptr) - call ESMF_FieldGet(ps_input_grid, & - farrayPtr=psptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +if (.not. isnative) then + !--------------------------------------------------------------------------- + ! Flip 'z' indices to all 3-d variables. Data is read in from model + ! top to surface. This program expects surface to model top. + !--------------------------------------------------------------------------- + + if (localpet == 0) print*,"- CALL FieldGet FOR SURFACE PRESSURE." + nullify(psptr) + call ESMF_FieldGet(ps_input_grid, & + farrayPtr=psptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - - nullify(presptr) - if (localpet == 0) print*,"- CALL FieldGet FOR 3-D PRESSURE." - call ESMF_FieldGet(pres_input_grid, & - computationalLBound=clb, & - computationalUBound=cub, & - farrayPtr=presptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + + nullify(presptr) + if (localpet == 0) print*,"- CALL FieldGet FOR 3-D PRESSURE." + call ESMF_FieldGet(pres_input_grid, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=presptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - nullify(tptr) - if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE." - call ESMF_FieldGet(temp_input_grid, & - farrayPtr=tptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + nullify(tptr) + if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE." + call ESMF_FieldGet(temp_input_grid, & + farrayPtr=tptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - nullify(uptr) - if (localpet == 0) print*,"- CALL FieldGet FOR U" - call ESMF_FieldGet(u_input_grid, & - farrayPtr=uptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + nullify(uptr) + if (localpet == 0) print*,"- CALL FieldGet FOR U" + call ESMF_FieldGet(u_input_grid, & + farrayPtr=uptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - - nullify(vptr) - if (localpet == 0) print*,"- CALL FieldGet FOR V" - call ESMF_FieldGet(v_input_grid, & - farrayPtr=vptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + + nullify(vptr) + if (localpet == 0) print*,"- CALL FieldGet FOR V" + call ESMF_FieldGet(v_input_grid, & + farrayPtr=vptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - nullify(wptr) - if (localpet == 0) print*,"- CALL FieldGet FOR W" - call ESMF_FieldGet(dzdt_input_grid, & - farrayPtr=wptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + nullify(wptr) + if (localpet == 0) print*,"- CALL FieldGet FOR W" + call ESMF_FieldGet(dzdt_input_grid, & + farrayPtr=wptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - if (localpet == 0) print*,"- CALL FieldGet FOR TRACERS." - do n=1,num_tracers + if (localpet == 0) print*,"- CALL FieldGet FOR TRACERS." + do n=1,num_tracers nullify(qptr) call ESMF_FieldGet(tracers_input_grid(n), & - farrayPtr=qptr, rc=rc) + farrayPtr=qptr, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) do i = clb(1),cub(1) do j = clb(2),cub(2) - qptr(i,j,:) = qptr(i,j,lev_input:1:-1) + qptr(i,j,:) = qptr(i,j,lev_input:1:-1) end do end do - end do - - do i = clb(1),cub(1) + end do + + do i = clb(1),cub(1) do j = clb(2),cub(2) presptr(i,j,:) = rlevs(lev_input:1:-1) tptr(i,j,:) = tptr(i,j,lev_input:1:-1) @@ -2829,17 +2895,39 @@ subroutine read_input_atm_grib2_file(localpet) vptr(i,j,:) = vptr(i,j,lev_input:1:-1) wptr(i,j,:) = wptr(i,j,lev_input:1:-1) end do - end do + end do - if (localpet == 0) then - print*,'psfc is ',clb(1),clb(2),psptr(clb(1),clb(2)) - print*,'pres is ',cub(1),cub(2),presptr(cub(1),cub(2),:) - - print*,'pres check 1',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2),1)), & - minval(presptr(clb(1):cub(1),clb(2):cub(2),1)) - print*,'pres check lev',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2), & - lev_input)),minval(presptr(clb(1):cub(1),clb(2):cub(2),lev_input)) + if (localpet == 0) then + print*,'psfc is ',clb(1),clb(2),psptr(clb(1),clb(2)) + print*,'pres is ',cub(1),cub(2),presptr(cub(1),cub(2),:) + + print*,'pres check 1',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2),1)), & + minval(presptr(clb(1):cub(1),clb(2):cub(2),1)) + print*,'pres check lev',localpet,maxval(presptr(clb(1):cub(1),clb(2):cub(2), & + lev_input)),minval(presptr(clb(1):cub(1),clb(2):cub(2),lev_input)) + endif + +else + ! For native files, read in pressure field directly from file but don't flip levels + if (localpet == 0) then + print*,"- READ PRESSURE." + vname = ":PRES:" + do vlev = 1, lev_input + iret = grb2_inq(the_file,inv_file,vname,slevs(vlev),data2=dummy2d) + if (iret<=0) then + call error_handler("READING IN PRESSURE AT LEVEL "//trim(slevs(vlev)),iret) + endif + dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8) + print*,'pres check after read ',vlev, dummy3d(1,1,vlev) + enddo + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID PRESSURE." + call ESMF_FieldScatter(pres_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) endif + deallocate(dummy3d) !--------------------------------------------------------------------------- ! Convert from 2-d to 3-d component winds. @@ -4498,31 +4586,44 @@ end subroutine read_input_sfc_netcdf_file subroutine read_input_sfc_grib2_file(localpet) use wgrib2api + use grib2_util, only : to_upper + use program_setup, only : vgtyp_from_climo, sotyp_from_climo + use model_grid, only : input_grid_type + use search_util + implicit none integer, intent(in) :: localpet character(len=250) :: the_file + character(len=250) :: geo_file character(len=20) :: vname, vname_file,slev character(len=50) :: method integer :: rc, varnum, iret, i, j,k + integer :: ncid2d, varid integer, parameter :: icet_default = 265.0 - logical :: exist + logical :: exist, rap_latlon real(esmf_kind_r4) :: value real(esmf_kind_r4), allocatable :: dummy2d(:,:),tsk_save(:,:),icec_save(:,:) - real(esmf_kind_r8), allocatable :: dummy2d_8(:,:) - real(esmf_kind_r8), allocatable :: dummy3d(:,:,:) + real(esmf_kind_r4), allocatable :: dummy1d(:) + real(esmf_kind_r8), allocatable :: dummy2d_8(:,:),dummy2d_82(:,:) + real(esmf_kind_r8), allocatable :: dummy3d(:,:,:), dummy3d_stype(:,:,:) integer(esmf_kind_i4), allocatable :: slmsk_save(:,:) + integer(esmf_kind_i8), allocatable :: dummy2d_i(:,:) + + rap_latlon = trim(to_upper(external_model))=="RAP" .and. trim(input_grid_type) == "rotated_latlon" the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid) - + geo_file = trim(geogrid_file_input_grid) + + print*,"- READ SFC DATA FROM GRIB2 FILE: ", trim(the_file) inquire(file=the_file,exist=exist) if (.not.exist) then @@ -4533,17 +4634,58 @@ subroutine read_input_sfc_grib2_file(localpet) lsoil_input = grb2_inq(the_file, inv_file, ':TSOIL:',' below ground:') print*, "- FILE HAS ", lsoil_input, " SOIL LEVELS" if (lsoil_input <= 0) call error_handler("COUNTING SOIL LEVELS.", rc) + + !We need to recreate the soil fields if we have something other than 4 levels + if (lsoil_input /= 4) then + + call ESMF_FieldDestroy(soil_temp_input_grid, rc=rc) + call ESMF_FieldDestroy(soilm_tot_input_grid, rc=rc) + call ESMF_FieldDestroy(soilm_liq_input_grid, rc=rc) + + print*,"- CALL FieldCreate FOR INPUT SOIL TEMPERATURE." + soil_temp_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_input/), 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 INPUT TOTAL SOIL MOISTURE." + soilm_tot_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_input/), 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 INPUT LIQUID SOIL MOISTURE." + soilm_liq_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + endif if (localpet == 0) then allocate(dummy2d(i_input,j_input)) allocate(slmsk_save(i_input,j_input)) + allocate(dummy2d_i(i_input,j_input)) allocate(tsk_save(i_input,j_input)) allocate(icec_save(i_input,j_input)) allocate(dummy2d_8(i_input,j_input)) + allocate(dummy2d_82(i_input,j_input)) allocate(dummy3d(i_input,j_input,lsoil_input)) + allocate(dummy3d_stype(i_input,j_input,16)) + allocate(dummy1d(16)) else allocate(dummy3d(0,0,0)) allocate(dummy2d_8(0,0)) + allocate(dummy2d_82(0,0)) allocate(dummy2d(0,0)) endif @@ -4582,8 +4724,8 @@ subroutine read_input_sfc_grib2_file(localpet) !---------------------------------------------------------------------------------- ! GFS v14 and v15.2 grib data has two land masks. LANDN is created by ! nearest neighbor interpolation. LAND is created by bilinear interpolation. -! LANDN matches the bitmap. So use it first. For other GFS versions, use LAND. -! Mask in grib file is '1' (land), '0' (not land). Add sea/lake ice category +! LANDN matches the bitmap. So use it first. For other GFS versions or other models, +! use LAND. Mask in grib file is '1' (land), '0' (not land). Add sea/lake ice category ! '2' based on ice concentration. !---------------------------------------------------------------------------------- @@ -4647,7 +4789,7 @@ subroutine read_input_sfc_grib2_file(localpet) if(dummy2d(i,j) == grb2_UNDEFINED) dummy2d(i,j) = 0.0_esmf_kind_r4 enddo enddo -! print*,'weasd ',maxval(dummy2d),minval(dummy2d) + print*,'weasd ',maxval(dummy2d),minval(dummy2d) endif print*,"- CALL FieldScatter FOR INPUT GRID SNOW LIQUID EQUIVALENT." @@ -4662,7 +4804,7 @@ subroutine read_input_sfc_grib2_file(localpet) where(dummy2d == grb2_UNDEFINED) dummy2d = 0.0_esmf_kind_r4 dummy2d = dummy2d*1000.0 ! Grib2 files have snow depth in (m), fv3 expects it in mm where(slmsk_save == 0) dummy2d = 0.0_esmf_kind_r4 -! print*,'snod ',maxval(dummy2d),minval(dummy2d) + print*,'snod ',maxval(dummy2d),minval(dummy2d) endif print*,"- CALL FieldScatter FOR INPUT GRID SNOW DEPTH." @@ -4727,9 +4869,100 @@ subroutine read_input_sfc_grib2_file(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldScatter", rc) -! Soil type is not available. Set to a large negative fill value. + if (localpet == 0) then + print*,"- READ SOIL TYPE." + slev=":surface:" + vname=":SOTYP:" + rc = grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + !failed => rc = 0 + if (rc <= 0 .and. (trim(to_upper(external_model))=="HRRR" .or. rap_latlon)) then + ! Some HRRR and RAP files don't have dominant soil type in the output, but the geogrid files + ! do, so this gives users the option to provide the geogrid file and use input soil + ! type + print*, "OPEN GEOGRID FILE ", trim(geo_file) + rc = nf90_open(geo_file,NF90_NOWRITE,ncid2d) + ! failed => rc < 0 + if (rc == 0) then + print*, "INQUIRE ABOUT SOIL TYPE FROM GEOGRID FILE" + rc = nf90_inq_varid(ncid2d,"SCT_DOM",varid) + ! failed => rc < 0 + if (rc<0) print*, "ERROR FINDING SCT_DOM IN GEOGRID FILE" + if (rc == 0) then + print*, "READ SOIL TYPE FROM GEOGRID FILE " + rc = nf90_get_var(ncid2d,varid,dummy2d) + ! failed => rc < 0 + if (rc<0) print*, "ERROR READING SCT_DOM FROM FILE" + print*, "min max dummy2d = ", minval(dummy2d), maxval(dummy2d) + endif + print*, "INQUIRE ABOUT SOIL TYPE FRACTIONS FROM GEOGRID FILE" + rc = nf90_inq_varid(ncid2d,"SOILCTOP",varid) + ! failed => rc < 0 + if (rc<0) print*, "ERROR FINDING SOILCTOP IN GEOGRID FILE" + if (rc == 0) then + print*, "READ SOIL TYPE FRACTIONS FROM GEOGRID FILE " + rc = nf90_get_var(ncid2d,varid,dummy3d_stype) + ! failed => rc < 0 + if (rc<0) print*, "ERROR READING SCT_DOM FROM FILE" + print*, "min max dummy3d_stype = ", minval(dummy3d_stype), maxval(dummy3d_stype) + endif - dummy2d_8 = -99999.0_esmf_kind_r8 + print*, "CLOSE GEOGRID FILE " + iret = nf90_close(ncid2d) + endif + + ! There's an issue with the geogrid file containing soil type water at land points. + ! This correction replaces the soil type at these points with the soil type with + ! the next highest fractional coverage. + do j = 1, j_input + do i = 1, i_input + if(dummy2d(i,j) == 14.0_esmf_kind_r4 .and. slmsk_save(i,j) == 1) then + dummy1d(:) = dummy3d_stype(i,j,:) + dummy1d(14) = 0.0_esmf_kind_r4 + dummy2d(i,j) = real(MAXLOC(dummy1d, 1),esmf_kind_r4) + endif + enddo + enddo + endif + + if ((rc <= 0 .and. trim(to_upper(external_model)) /= "HRRR" .and. .not. rap_latlon) & + .or. (rc < 0 .and. (trim(to_upper(external_model)) == "HRRR" .or. rap_latlon))) then + if (.not. sotyp_from_climo) then + call error_handler("COULD NOT FIND SOIL TYPE IN FILE. PLEASE SET SOTYP_FROM_CLIMO=.TRUE. . EXITING", rc) + else + vname = "sotyp" + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + loc=varnum) + call handle_grib_error(vname, slev ,method,value,varnum,rc, var= dummy2d) + if (rc == 1) then ! missing_var_method == skip or no entry in varmap table + print*, "WARNING: "//trim(vname)//" NOT AVAILABLE IN FILE. WILL NOT "//& + "SCALE SOIL MOISTURE FOR DIFFERENCES IN SOIL TYPE. " + dummy2d(:,:) = -99999.0_esmf_kind_r4 + endif + endif + endif + + ! In the event that the soil type on the input grid still contains mismatches between + ! soil type and landmask, this correction is a last-ditch effort to replace these points + ! with soil type from a nearby land point. + if (.not. sotyp_from_climo) then + do j = 1, j_input + do i = 1, i_input + if(dummy2d(i,j) == 14.0_esmf_kind_r4 .and. slmsk_save(i,j) == 1) dummy2d(i,j) = -99999.9 + enddo + enddo + + dummy2d_8 = real(dummy2d,esmf_kind_r8) + dummy2d_i(:,:) = 0 + where(slmsk_save == 1) dummy2d_i = 1 + + call search(dummy2d_8,dummy2d_i,i_input,j_input,1,230) + endif + + print*,'sotype ',maxval(dummy2d_8),minval(dummy2d_8) + deallocate(dummy2d_i) + deallocate(dummy3d_stype) + endif + print*,"- CALL FieldScatter FOR INPUT GRID SOIL TYPE." call ESMF_FieldScatter(soil_type_input_grid,dummy2d_8, rootpet=0, rc=rc) @@ -4738,10 +4971,127 @@ subroutine read_input_sfc_grib2_file(localpet) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Begin variables whose presence in grib2 files varies, but no climatological - ! data is - ! available, so we have to account for values in the varmap table + ! data is available, so we have to account for values in the varmap table !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (.not. vgfrc_from_climo) then + if (localpet == 0) then + print*,"- READ VEG FRACTION." + vname="vfrac" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + loc=varnum) + !! Changing these for GSD internal runs using new HRRR files + vname=":VEG:" + rc= grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d) + + if (rc > 1) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1105:', data2=dummy2d) + if (rc <= 0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1101:', data2=dummy2d) + if (rc <= 0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1151:', data2=dummy2d) + if (rc <= 0) call error_handler("COULD NOT DETERMINE VEGETATION FRACTION IN FILE. & + RECORD NUMBERS MAY HAVE CHANGED. PLEASE SET VGFRC_FROM_CLIMO=.TRUE. EXITING", rc) + endif + endif + elseif (rc <= 0) then + call error_handler("COULD NOT FIND VEGETATION FRACTION IN FILE. & + PLEASE SET VGFRC_FROM_CLIMO=.TRUE. EXITING", rc) + endif + if(maxval(dummy2d) > 2.0) dummy2d = dummy2d / 100.0_esmf_kind_r4 + print*,'vfrac ',maxval(dummy2d),minval(dummy2d) + endif + + + print*,"- CALL FieldScatter FOR INPUT GRID VEG GREENNESS." + call ESMF_FieldScatter(veg_greenness_input_grid,real(dummy2d,esmf_kind_r8), rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + endif + + if (.not. minmax_vgfrc_from_climo) then + if (localpet == 0) then + print*,"- READ MIN VEG FRACTION." + vname="vfrac_min" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":VEG:" + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1106:',data2=dummy2d) + + if (rc <= 0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1102:',data2=dummy2d) + if (rc <= 0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1152:',data2=dummy2d) + if (rc<=0) call error_handler("COULD NOT FIND MIN VEGETATION FRACTION IN FILE. & + PLEASE SET MINMAX_VGFRC_FROM_CLIMO=.TRUE. . EXITING",rc) + endif + endif + if(maxval(dummy2d) > 2.0) dummy2d = dummy2d / 100.0_esmf_kind_r4 + print*,'vfrac min',maxval(dummy2d),minval(dummy2d) + + endif + + print*,"- CALL FieldScatter FOR INPUT GRID MIN VEG GREENNESS." + call ESMF_FieldScatter(min_veg_greenness_input_grid,real(dummy2d,esmf_kind_r8), rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) then + print*,"- READ MAX VEG FRACTION." + vname="vfrac_max" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + + vname=":VEG:" + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1107:',data2=dummy2d) + if (rc <=0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1103:',data2=dummy2d) + if (rc <=0) then + rc= grb2_inq(the_file, inv_file, vname,slev,'n=1153:',data2=dummy2d) + if (rc <= 0) call error_handler("COULD NOT FIND MAX VEGETATION FRACTION IN FILE. & + PLEASE SET MINMAX_VGFRC_FROM_CLIMO=.TRUE. . EXITING",rc) + endif + endif + if(maxval(dummy2d) > 2.0) dummy2d = dummy2d / 100.0_esmf_kind_r4 + print*,'vfrac max',maxval(dummy2d),minval(dummy2d) + + endif !localpet==0 + + print*,"- CALL FieldScatter FOR INPUT GRID MAX VEG GREENNESS." + call ESMF_FieldScatter(max_veg_greenness_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + endif !minmax_vgfrc_from_climo + + if (.not. lai_from_climo) then + if (localpet == 0) then + print*,"- READ LAI." + vname="lai" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & + loc=varnum) + vname=":var0_7_198:" + rc= grb2_inq(the_file, inv_file, vname,slev,':n=1108:',data2=dummy2d) + if (rc <=0) then + rc= grb2_inq(the_file, inv_file, vname,slev,':n=1104:',data2=dummy2d) + if (rc <=0) then + rc= grb2_inq(the_file, inv_file, vname,slev,':n=1154:',data2=dummy2d) + if (rc <= 0) call error_handler("COULD NOT FIND LAI IN FILE. & + PLEASE SET LAI_FROM_CLIMO=.TRUE. . EXITING",rc) + endif + endif + print*,'lai',maxval(dummy2d),minval(dummy2d) + endif !localpet==0 + + print*,"- CALL FieldScatter FOR INPUT GRID LAI." + call ESMF_FieldScatter(lai_input_grid,real(dummy2d,esmf_kind_r8),rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + + endif if (localpet == 0) then print*,"- READ SEAICE DEPTH." vname="hice" @@ -4929,7 +5279,6 @@ subroutine read_input_sfc_grib2_file(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldScatter", rc) - deallocate(dummy2d) if (localpet == 0) then print*,"- READ LIQUID SOIL MOISTURE." @@ -4953,38 +5302,91 @@ subroutine read_input_sfc_grib2_file(localpet) call read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) print*,'soilm ',maxval(dummy3d),minval(dummy3d) endif + + print*,"- CALL FieldScatter FOR INPUT TOTAL SOIL MOISTURE." + call ESMF_FieldScatter(soilm_tot_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldScatter", rc) + +!---------------------------------------------------------------------------------------- +! Vegetation type is not available in some files. However, it is needed to identify +! permanent land ice points. At land ice, the total soil moisture is a flag value of +! '1'. Use this flag as a temporary solution. +!---------------------------------------------------------------------------------------- -!----------------------------------------------------------------------- -! Vegetation type is not available. However, it is needed to identify -! permanent land ice points. At land ice, the total soil moisture -! is a flag value of '1'. Use this flag as a temporary solution. -!----------------------------------------------------------------------- - + print*, "- CALL FieldGather for INPUT SOIL TYPE." + call ESMF_FieldGather(soil_type_input_grid, dummy2d_82, rootPet=0, tile=1, 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 - dummy2d_8(:,:) = 0.0_esmf_kind_r8 + print*,"- READ VEG TYPE." + vname="vtype" + slev=":surface:" + call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, & + loc=varnum) + !Note: sometimes the grib files don't have this one named. Searching for this string + ! ensures that the data is found when it exists + + vname="var2_2" + rc= grb2_inq(the_file, inv_file, vname,"_0_198:",slev,' hour fcst:', data2=dummy2d) + if (rc <= 0) then + rc= grb2_inq(the_file, inv_file, vname,"_0_198:",slev,':anl:', data2=dummy2d) + if (rc <= 0) then + if (.not. vgtyp_from_climo) then + call error_handler("COULD NOT FIND VEGETATION TYPE IN FILE. PLEASE SET VGTYP_FROM_CLIMO=.TRUE. . EXITING", rc) + else + do j = 1, j_input + do i = 1, i_input + dummy2d(i,j) = 0.0_esmf_kind_r4 + if(slmsk_save(i,j) == 1 .and. dummy3d(i,j,1) > 0.99) & + dummy2d(i,j) = real(veg_type_landice_input,esmf_kind_r4) + enddo + enddo + endif ! replace_vgtyp + endif !not find :anl: + endif !not find hour fcst: + + if (trim(external_model) .ne. "GFS") then do j = 1, j_input - do i = 1, i_input - if(slmsk_save(i,j) == 1_esmf_kind_i4 .and. dummy3d(i,j,1) > 0.99) & - dummy2d_8(i,j) = real(veg_type_landice_input,esmf_kind_r8) - enddo + do i = 1,i_input + if (dummy2d(i,j) == 15.0_esmf_kind_r4 .and. slmsk_save(i,j) == 1) then + if (dummy3d(i,j,1) < 0.6) then + dummy2d(i,j) = real(veg_type_landice_input,esmf_kind_r4) + elseif (dummy3d(i,j,1) > 0.99) then + slmsk_save(i,j) = 0 + dummy2d(i,j) = 0.0_esmf_kind_r4 + dummy2d_82(i,j) = 0.0_esmf_kind_r8 + endif + elseif (dummy2d(i,j) == 17.0_esmf_kind_r4 .and. slmsk_save(i,j)==0) then + dummy2d(i,j) = 0.0_esmf_kind_r4 + endif + enddo enddo - endif - + endif + dummy2d_8= real(dummy2d,esmf_kind_r8) + print*,'vgtyp ',maxval(dummy2d),minval(dummy2d) + endif !localpet + deallocate(dummy2d) print*,"- CALL FieldScatter FOR INPUT VEG TYPE." call ESMF_FieldScatter(veg_type_input_grid, dummy2d_8, rootpet=0, 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 INPUT TOTAL SOIL MOISTURE." - call ESMF_FieldScatter(soilm_tot_input_grid, dummy3d, rootpet=0, rc=rc) + print*,"- CALL FieldScatter FOR INPUT VEG TYPE." + call ESMF_FieldScatter(soil_type_input_grid, dummy2d_82, rootpet=0, 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 INPUT LANDSEA MASK." + call ESMF_FieldScatter(landsea_mask_input_grid,real(slmsk_save,esmf_kind_r8),rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& call error_handler("IN FieldScatter", rc) !--------------------------------------------------------------------------------- ! At open water (slmsk==0), the soil temperature array is not used and set ! to the filler value of SST. At lake/sea ice points (slmsk=2), the soil -! temperature array holds ice column temperature. That field is not available -! in GFS grib data, so set to a default value. +! temperature array holds ice column temperature. This field is not available +! in the grib data, so set to a default value. !--------------------------------------------------------------------------------- if (localpet == 0) then @@ -5604,6 +6006,166 @@ SUBROUTINE READ_FV3_GRID_DATA_NETCDF(FIELD,TILE_NUM,IMO,JMO,LMO, & ERROR = NF90_CLOSE(NCID) END SUBROUTINE READ_FV3_GRID_DATA_NETCDF + + !--------------------------------------------------------------------------- +! Read winds from a grib2 file +!--------------------------------------------------------------------------- + + subroutine read_winds(file,inv,u,v,localpet) + + use wgrib2api + use netcdf + use program_setup, only : get_var_cond, fix_dir_input_grid + use model_grid, only : input_grid_type + implicit none + + character(len=250), intent(in) :: file + character(len=10), intent(in) :: inv + integer, intent(in) :: localpet + real(esmf_kind_r8), intent(inout), allocatable :: u(:,:,:),v(:,:,:) + + real(esmf_kind_r4), dimension(i_input,j_input) :: alpha + real(esmf_kind_r8), dimension(i_input,j_input) :: lon, lat + real(esmf_kind_r4), allocatable :: u_tmp(:,:),v_tmp(:,:) + real(esmf_kind_r4), dimension(i_input,j_input) :: ws,wd + real(esmf_kind_r4) :: value_u, value_v,lov,latin1,latin2 + real(esmf_kind_r8) :: d2r + + integer :: varnum_u, varnum_v, vlev, & !ncid, id_var, & + error, iret, i,istr + + character(len=20) :: vname + character(len=50) :: method_u, method_v + character(len=250) :: file_coord, cmdline_msg + character(len=10000) :: temp_msg + + d2r=acos(-1.0_esmf_kind_r8) / 180.0_esmf_kind_r8 + if (localpet==0) then + allocate(u(i_input,j_input,lev_input)) + allocate(v(i_input,j_input,lev_input)) + else + allocate(u(0,0,0)) + allocate(v(0,0,0)) + endif + + file_coord = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc" + + vname = "u" + call get_var_cond(vname,this_miss_var_method=method_u, this_miss_var_value=value_u, & + loc=varnum_u) + vname = "v" + call get_var_cond(vname,this_miss_var_method=method_v, this_miss_var_value=value_v, & + loc=varnum_v) + + if (trim(input_grid_type)=="rotated_latlon") then + print*,"- CALL FieldGather FOR INPUT GRID LONGITUDE" + call ESMF_FieldGather(longitude_input_grid, lon, rootPet=0, tile=1, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + print*,"- CALL FieldGather FOR INPUT GRID LATITUDE" + call ESMF_FieldGather(latitude_input_grid, lat, rootPet=0, tile=1, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + + if (localpet==0) then + print*,"- CALCULATE ROTATION ANGLE FOR ROTATED_LATLON INPUT GRID" + error = grb2_inq(file, inv,grid_desc=temp_msg) + !1:0:grid_template=32769:winds(grid): + ! I am not an Arakawa E-grid. + ! I am rotated but have no rotation angle. + ! I am staggered. What am I? + ! (953 x 834) units 1e-06 input WE:SN output WE:SN res 56 + ! lat0 -10.590603 lat-center 54.000000 dlat 121.813000 + ! lon0 220.914154 lon-center 254.000000 dlon 121.813000 #points=794802 + + istr = index(temp_msg, "lat-center ") + len("lat_center ") + read(temp_msg(istr:istr+9),"(F8.5)") latin1 + istr = index(temp_msg, "lon-center ") + len("lon-center ") + read(temp_msg(istr:istr+10),"(F9.6)") lov + + print*, "- CALL CALCALPHA_ROTLATLON with center lat,lon = ",latin1,lov + call calcalpha_rotlatlon(lat,lon,latin1,lov,alpha) + print*, " alpha min/max = ",MINVAL(alpha),MAXVAL(alpha) + endif + elseif (trim(input_grid_type) == "lambert") then + !# NG this has been edited to correctly calculate gridrot for Lambert grids + ! Previously was incorrectly using polar-stereographic formation + print*,"- CALL FieldGather FOR INPUT GRID LONGITUDE" + call ESMF_FieldGather(longitude_input_grid, lon, rootPet=0, tile=1, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + + if (localpet==0) then + error = grb2_inq(file, inv,grid_desc=temp_msg) + !1:0:grid_template=30:winds(grid): + ! Lambert Conformal: (1799 x 1059) input WE:SN output WE:SN res 8 + ! Lat1 21.138123 Lon1 237.280472 LoV 262.500000 + ! LatD 38.500000 Latin1 38.500000 Latin2 38.500000 + ! LatSP 0.000000 LonSP 0.000000 + ! North Pole (1799 x 1059) Dx 3000.000000 m Dy 3000.000000 m mode 8 + + istr = index(temp_msg, "LoV ") + len("LoV ") + read(temp_msg(istr:istr+10),"(F9.6)") lov + istr = index(temp_msg, "Latin1 ") + len("Latin1 ") + read(temp_msg(istr:istr+9),"(F8.5)") latin1 + istr = index(temp_msg, "Latin2 ") + len("Latin2 ") + read(temp_msg(istr:istr+9),"(F8.5)") latin2 + + print*, "- CALL GRIDROT for LC grid with lov,latin1/2 = ",lov,latin1,latin2 + call gridrot(lov,latin1,latin2,lon,alpha) + print*, " alpha min/max = ",MINVAL(alpha),MAXVAL(alpha) + endif + endif + + if (localpet==0) then + do vlev = 1, lev_input + + vname = ":UGRD:" + iret = grb2_inq(file,inv,vname,slevs(vlev),data2=u_tmp) + if (iret <= 0) then + call handle_grib_error(vname, slevs(vlev),method_u,value_u,varnum_u,iret,var=u_tmp) + if (iret==1) then ! missing_var_method == skip + call error_handler("READING IN U AT LEVEL "//trim(slevs(vlev))//". SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret) + endif + endif + + vname = ":VGRD:" + iret = grb2_inq(file,inv,vname,slevs(vlev),data2=v_tmp) + if (iret <= 0) then + call handle_grib_error(vname, slevs(vlev),method_v,value_v,varnum_v,iret,var=v_tmp) + if (iret==1) then ! missing_var_method == skip + call error_handler("READING IN V AT LEVEL "//trim(slevs(vlev))//". SET A FILL "// & + "VALUE IN THE VARMAP TABLE IF THIS ERROR IS NOT DESIRABLE.",iret) + endif + endif + + if (trim(input_grid_type) == "latlon") then + if (external_model == 'UKMET') then + u(:,:,vlev) = u_tmp + v(:,:,vlev) = (v_tmp(:,2:jp1_input) + v_tmp(:,1:j_input))/2 + else + u(:,:,vlev) = u_tmp + v(:,:,vlev) = v_tmp + endif + else if (trim(input_grid_type) == "rotated_latlon") then + ws = sqrt(u_tmp**2 + v_tmp**2) + wd = atan2(-u_tmp,-v_tmp) / d2r ! calculate grid-relative wind direction + wd = wd + alpha + 180.0 ! Rotate from grid- to earth-relative direction + wd = 270.0 - wd ! Convert from meteorological (true N) to mathematical direction + u(:,:,vlev) = -ws*cos(wd*d2r) + v(:,:,vlev) = -ws*sin(wd*d2r) + else + u(:,:,vlev) = real(u_tmp * cos(alpha) + v_tmp * sin(alpha),esmf_kind_r8) + v(:,:,vlev) = real(v_tmp * cos(alpha) - u_tmp * sin(alpha),esmf_kind_r8) + endif + + print*, 'max, min U ', minval(u(:,:,vlev)), maxval(u(:,:,vlev)) + print*, 'max, min V ', minval(v(:,:,vlev)), maxval(v(:,:,vlev)) + enddo + endif + +end subroutine read_winds !--------------------------------------------------------------------------- ! Convert from 2-d to 3-d winds. @@ -5672,6 +6234,85 @@ subroutine convert_winds end subroutine convert_winds +!--------------------------------------------------------------------------- +! Compute grid rotation angle for non-latlon grids +!--------------------------------------------------------------------------- + +!# NG The original gridrot subroutine was specific to polar stereographic grids. +! We need to compute it for Lambert Conformal grids. So we need lat1,lat2 +! Note this follows the ncl_ncarg source code +! ncl_ncarg-6.6.2/ni/src/ncl/GetGrids.c +subroutine gridrot(lov,latin1,latin2,lon,rot) + + use model_grid, only : i_input,j_input + implicit none + + + real(esmf_kind_r4), intent(in) :: lov,latin1,latin2 + real(esmf_kind_r4), intent(inout) :: rot(i_input,j_input) + real(esmf_kind_r8), intent(in) :: lon(i_input,j_input) + + real(esmf_kind_r4) :: trot(i_input,j_input), tlon(i_input,j_input) + real(esmf_kind_r4) :: dtor = 3.14159265359/180.0_esmf_kind_r4 + real(esmf_kind_r4) :: an + !trot_tmp = real(lon,esmf_kind_r4)-lov + !trot = trot_tmp + !where(trot_tmp > 180.0) trot = trot-360.0_esmf_kind_r4 + !where(trot_tmp < -180.0) trot = trot-360.0_esmf_kind_r4 + + if ( (latin1 - latin2) .lt. 0.000001 ) then + an = sin(latin1*dtor) + else + an = log( cos(latin1*dtor) / cos(latin2*dtor) ) / & + log( tan(dtor*(90.0-latin1)/2.) / tan(dtor*(90.0-latin2)/2.)) + end if + + tlon = mod(lon - lov + 180. + 3600., 360.) - 180. + trot = an * tlon + + rot = trot * dtor + +end subroutine gridrot + +! Subroutine calcalpha_rotlatlon calculates rotation angle +! specific to rotated latlon grids, needed to convert to +! earth-relative winds +subroutine calcalpha_rotlatlon(latgrid,longrid,cenlat,cenlon,alpha) + + use model_grid, only : i_input,j_input + implicit none + + real(esmf_kind_r8), intent(in) :: latgrid(i_input,j_input), & + longrid(i_input,j_input) + real(esmf_kind_r4), intent(in) :: cenlat, cenlon + real(esmf_kind_r4), intent(out) :: alpha(i_input,j_input) + + ! Variables local to subroutine + real(esmf_kind_r8) :: D2R,lon0_r,lat0_r,sphi0,cphi0 + real(esmf_kind_r8), DIMENSION(i_input,j_input) :: tlat,tlon,tph,sinalpha + + D2R = acos(-1.0_esmf_kind_r8) / 180.0_esmf_kind_r8 + if (cenlon .lt. 0) then + lon0_r = (cenlon + 360.0)*D2R + else + lon0_r = cenlon*D2R + end if + lat0_r=cenlat*D2R + sphi0=sin(lat0_r) + cphi0=cos(lat0_r) + + ! deal with input lat/lon + tlat = latgrid * D2R + tlon = longrid * D2R + + ! Calculate alpha (rotation angle) + tlon = -tlon + lon0_r + tph = asin(cphi0*sin(tlat) - sphi0*cos(tlat)*cos(tlon)) + sinalpha = sphi0 * sin(tlon) / cos(tph) + alpha = -asin(sinalpha)/D2R + ! returns alpha in degrees +end subroutine calcalpha_rotlatlon + subroutine handle_grib_error(vname,lev,method,value,varnum, iret,var,var8,var3d) use, intrinsic :: ieee_arithmetic @@ -5748,9 +6389,13 @@ subroutine read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc) if(lsoil_input == 4) then slevs = (/character(24)::':0-0.1 m below ground:', ':0.1-0.4 m below ground:', & ':0.4-1 m below ground:', ':1-2 m below ground:'/) + elseif(lsoil_input == 9) then + slevs = (/character(26)::':0-0 m below ground',':0.01-0.01 m below ground:',':0.04-0.04 m below ground:', & + ':0.1-0.1 m below ground:',':0.3-0.3 m below ground:',':0.6-0.6 m below ground:', & + ':1-1 m below ground:',':1.6-1.6 m below ground:',':3-3 m below ground:'/) else rc = -1 - call error_handler("reading soil levels. File must have 4 soil levels.", rc) + call error_handler("reading soil levels. File must have 4 or 9 soil levels.", rc) endif call get_var_cond(vname,this_miss_var_method=method,this_miss_var_value=value, & @@ -5864,6 +6509,16 @@ subroutine cleanup_input_sfc_data call ESMF_FieldDestroy(veg_type_input_grid, rc=rc) call ESMF_FieldDestroy(z0_input_grid, rc=rc) call ESMF_FieldDestroy(terrain_input_grid, rc=rc) + if (.not. vgfrc_from_climo) then + call ESMF_FieldDestroy(veg_greenness_input_grid, rc=rc) + endif + if (.not. minmax_vgfrc_from_climo) then + call ESMF_FieldDestroy(min_veg_greenness_input_grid, rc=rc) + call ESMF_FieldDestroy(max_veg_greenness_input_grid, rc=rc) + endif + if (.not. lai_from_climo) then + call ESMF_FieldDestroy(lai_input_grid, rc=rc) + endif end subroutine cleanup_input_sfc_data diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 84e84a7a3..5f207657e 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -56,6 +56,7 @@ module model_grid !-------------------------------------------------------------------------- use esmf + use ESMF_LogPublicMod implicit none @@ -63,8 +64,12 @@ module model_grid character(len=5), allocatable, public :: tiles_target_grid(:) character(len=10), public :: inv_file = "chgres.inv" + character(len=50), public :: input_grid_type = "latlon" + !character(len=100), public :: the_file_hrrr = "./HRRR_adj_rad.grib2" - integer, parameter, public :: lsoil_target = 4 ! # soil layers + ! Made lsoil_target non-parameter to allow for RAP land surface initiation + integer, public :: lsoil_target = 4 ! # soil layers + integer, public :: i_input, j_input integer, public :: ip1_input, jp1_input integer, public :: i_target, j_target @@ -109,7 +114,7 @@ module model_grid subroutine define_input_grid(localpet, npets) - use program_setup, only : input_type + use program_setup, only : input_type, external_model implicit none @@ -120,8 +125,10 @@ subroutine define_input_grid(localpet, npets) trim(input_type) == "gfs_sigio" .or. & trim(input_type) == "gaussian_netcdf") then call define_input_grid_gaussian(localpet, npets) - elseif (trim(input_type) == "grib2") then + elseif (trim(external_model) == "GFS" .and. trim(input_type) == "grib2") then call define_input_grid_gfs_grib2(localpet,npets) + elseif (trim(input_type) == "grib2") then + call define_input_grid_grib2(localpet,npets) else call define_input_grid_mosaic(localpet, npets) endif @@ -602,10 +609,8 @@ end subroutine define_input_grid_mosaic subroutine define_input_grid_gfs_grib2(localpet, npets) - use mpi - use wgrib2api - + use mpi use program_setup, only : data_dir_input_grid, & grib2_file_input_grid @@ -615,7 +620,8 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) character(len=250) :: the_file - integer :: i, j, rc, clb(2), cub(2), ierr + integer :: i, j, rc, clb(2), cub(2) + integer :: ierr real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) @@ -633,15 +639,15 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) num_tiles_input_grid = 1 the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid - if(localpet == 0) then + if (localpet == 0) then print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) rc=grb2_mk_inv(the_file,inv_file) if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc) endif -! Wait for localpet 0 to create inventory. +! Wait for localpet 0 to create inventory call mpi_barrier(mpi_comm_world, ierr) - + rc = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & lat=lat4, lon=lon4) if (rc /= 1) call error_handler("READING GRIB2 FILE", rc) @@ -788,12 +794,311 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) end subroutine define_input_grid_gfs_grib2 + subroutine define_input_grid_grib2(localpet, npets) + + use mpi + use netcdf + use wgrib2api + use program_setup, only : grib2_file_input_grid, data_dir_input_grid, & + fix_dir_input_grid, external_model + implicit none + + character(len=500) :: the_file, temp_file + + integer, intent(in) :: localpet, npets + + integer :: error, extra, i, j, clb(2), cub(2) + + real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:), lat_corners(:,:) + real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:), lon_corners(:,:) + real(esmf_kind_r8) :: lat_target(i_target,j_target), & + lon_target(i_target,j_target) + real(esmf_kind_r8) :: deltalon, dx + integer :: ncid,id_var, id_dim + real(esmf_kind_r8), pointer :: lat_src_ptr(:,:), lon_src_ptr(:,:) + character(len=10000) :: cmdline_msg, temp_msg, temp_msg2 + character(len=10) :: temp_num = 'NA' + + num_tiles_input_grid = 1 + + !inv_file = "chgres.inv" + the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid + temp_file = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc" + + call ESMF_FieldGather(latitude_target_grid, lat_target, rootPet=0, tile=1, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + call ESMF_FieldGather(longitude_target_grid, lon_target, rootPet=0, tile=1, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + + if (localpet==0) then + print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) + error=grb2_mk_inv(the_file,inv_file) + if (error /=0) call error_handler("OPENING GRIB2 FILE",error) + error = grb2_inq(the_file, inv_file,grid_desc=temp_msg) + i = index(temp_msg, "grid_template=") + len("grid_template=") + j = index(temp_msg,":winds(") + temp_num=temp_msg(i:j-1) + endif + call MPI_BARRIER(MPI_COMM_WORLD, error) + call MPI_BCAST(temp_num,10,MPI_CHAR,0,MPI_COMM_WORLD,error) + + ! Wgrib2 can't properly read the lat/lon arrays of data on NCEP rotated lat/lon + ! grids, so read in lat/lon from fixed coordinate file + print*, 'temp num =', temp_num + if (trim(temp_num)=="3.32769" .or. trim(temp_num)=="32769") then + + input_grid_type = "rotated_latlon" + + error=nf90_open(trim(temp_file),nf90_nowrite,ncid) + call netcdf_err(error, 'opening: '//trim(temp_file)) + + error=nf90_inq_dimid(ncid, 'nx', id_dim) + call netcdf_err(error, 'reading nx id' ) + error=nf90_inquire_dimension(ncid,id_dim,len=i_input) + call netcdf_err(error, 'reading nx value' ) + + error=nf90_inq_dimid(ncid, 'ny', id_dim) + call netcdf_err(error, 'reading ny id' ) + error=nf90_inquire_dimension(ncid,id_dim,len=j_input) + call netcdf_err(error, 'reading ny value' ) + + allocate(longitude_one_tile(i_input,j_input)) + allocate(latitude_one_tile(i_input,j_input)) + + error=nf90_inq_varid(ncid, 'gridlat', id_var) + call netcdf_err(error, 'reading field id' ) + error=nf90_get_var(ncid, id_var, latitude_one_tile) + call netcdf_err(error, 'reading field' ) + + error=nf90_inq_varid(ncid, 'gridlon', id_var) + call netcdf_err(error, 'reading field id' ) + error=nf90_get_var(ncid, id_var, longitude_one_tile) + call netcdf_err(error, 'reading field' ) + + elseif (temp_num == "3.0" .or. temp_num == "3.30" .or. temp_num=="30" .or. temp_num == "0") then + + if (temp_num =="3.0" .or. temp_num == "0") input_grid_type = "latlon" + if (temp_num =="3.30" .or. temp_num=='30') input_grid_type = "lambert" + + error = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & + lat=latitude_one_tile, lon=longitude_one_tile) + if (error /= 1) call error_handler("READING FILE", error) + + + if (localpet==0) print*, "from file lon(1:10,1) = ", longitude_one_tile(1:10,1) + if (localpet==0) print*, "from file lat(1,1:10) = ", latitude_one_tile(1,1:10) + elseif (temp_num=="NA") then + error = 0 + call error_handler("Grid template number cannot be read from the input file. Please " //& + "check that the wgrib2 executable is in your path.", error) + else + error = 0 + call error_handler("Unknown input file grid template number. Must be one of: " //& + "3, 3.30, 3.32769", error) + endif + + print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input + + ip1_input = i_input + 1 + jp1_input = j_input + 1 + +!----------------------------------------------------------------------- +! Create ESMF grid object for the model grid. +!----------------------------------------------------------------------- + + extra = npets / num_tiles_input_grid + + print*,"- CALL GridCreateNoPeriDim FOR INPUT MODEL GRID" + input_grid = ESMF_GridCreateNoPeriDim(maxIndex=(/i_input,j_input/), & + indexflag=ESMF_INDEX_GLOBAL, & + rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreateNoPeriDim", error) + + +!----------------------------------------------------------------------- +! Read the mask and lat/lons. +!----------------------------------------------------------------------- + + print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." + latitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_latitude", & + rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", error) + + print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE." + longitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_longitude", & + rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", error) + + print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. " + call ESMF_FieldScatter(latitude_input_grid, real(latitude_one_tile,esmf_kind_r8), rootpet=0, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", error) + +print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." + call ESMF_FieldScatter(longitude_input_grid, real(longitude_one_tile,esmf_kind_r8), rootpet=0, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", error) + + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", error) + + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=1, & + farrayPtr=lon_src_ptr, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", error) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_src_ptr, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", error) + + do j = clb(2),cub(2) + do i = clb(1), cub(1) + lon_src_ptr(i,j)=real(longitude_one_tile(i,j),esmf_kind_r8) + lat_src_ptr(i,j)=real(latitude_one_tile(i,j),esmf_kind_r8) + enddo + enddo + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CORNER, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", error) + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=1, & + farrayPtr=lon_src_ptr, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", error) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_src_ptr, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", error) + + print*,'bounds for corners ',localpet,clb(1),cub(1),clb(2),cub(2) + + ! If we have data on a lat/lon or lambert grid, create staggered coordinates + if(trim(input_grid_type)=="latlon" .or. trim(input_grid_type) == "lambert") then + if (trim(input_grid_type) == "latlon") then + + deltalon = abs(longitude_one_tile(2,1)-longitude_one_tile(1,1)) + do j = clb(2), cub(2) + do i = clb(1), cub(1) + + if (i == ip1_input) then + lon_src_ptr(i,j) = longitude_one_tile(i_input,1) + (0.5_esmf_kind_r8*deltalon) + else + lon_src_ptr(i,j) = longitude_one_tile(i,1) - (0.5_esmf_kind_r8*deltalon) + endif + + if (j == jp1_input) then + lat_src_ptr(i,j) = latitude_one_tile(1,j_input) + (0.5_esmf_kind_r8*deltalon) + else + lat_src_ptr(i,j) = latitude_one_tile(1,j) - (0.5_esmf_kind_r8*deltalon) + endif + + enddo + enddo + else + if (localpet==0) then + !cmdline_msg = "wgrib2 "//trim(the_file)//" -d 1 -grid &> temp2.out" + !call system(cmdline_msg) + !open(4,file="temp2.out") + !do i = 1,6 + ! read(4,"(A)") temp_msg2 + !enddo + !close(4) + print*, trim(temp_msg) + i = index(temp_msg, "Dx ") + len("Dx ") + j = index(temp_msg," m Dy") + read(temp_msg(i:j-1),"(F9.6)") dx + print*, "DX = ", dx + endif + call MPI_BARRIER(MPI_COMM_WORLD,error) + call MPI_BCAST(dx,1,MPI_REAL8,0,MPI_COMM_WORLD,error) + + call get_cell_corners(real(latitude_one_tile,esmf_kind_r8), & + real(longitude_one_tile, esmf_kind_r8), & + lat_src_ptr, lon_src_ptr, dx, clb, cub) + endif + elseif (trim(input_grid_type) == "rotated_latlon") then !Read the corner coords from file + + allocate(lon_corners(ip1_input,jp1_input)) + allocate(lat_corners(ip1_input,jp1_input)) + + error=nf90_inq_varid(ncid, 'gridlon_corners', id_var) + call netcdf_err(error, 'reading field id' ) + error=nf90_get_var(ncid, id_var, lon_corners) + call netcdf_err(error, 'reading field' ) + + error=nf90_inq_varid(ncid, 'gridlat_corners', id_var) + call netcdf_err(error, 'reading field id' ) + error=nf90_get_var(ncid, id_var, lat_corners) + call netcdf_err(error, 'reading field' ) + print*, 'min, max lat corners =', minval(lat_corners), maxval(lat_corners) + print*, 'min, max lon corners =', minval(lon_corners), maxval(lon_corners) + + do j = clb(2),cub(2) + do i = clb(1), cub(1) + lon_src_ptr(i,j)=real(lon_corners(i,j),esmf_kind_r8) + lat_src_ptr(i,j)=real(lat_corners(i,j),esmf_kind_r8) + enddo + enddo + + error= nf90_close(ncid) + endif + + nullify(lon_src_ptr) + nullify(lat_src_ptr) + + deallocate(longitude_one_tile) + deallocate(latitude_one_tile) + + end subroutine define_input_grid_grib2 + subroutine define_target_grid(localpet, npets) use netcdf use program_setup, only : mosaic_file_target_grid, & orog_dir_target_grid, & - orog_files_target_grid + orog_files_target_grid, & + nsoill_out implicit none @@ -817,6 +1122,8 @@ subroutine define_target_grid(localpet, npets) real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:) real(esmf_kind_r8), allocatable :: terrain_one_tile(:,:) + lsoil_target = nsoill_out + print*,'- OPEN TARGET GRID MOSAIC FILE: ',trim(mosaic_file_target_grid) error=nf90_open(trim(mosaic_file_target_grid),nf90_nowrite,ncid) call netcdf_err(error, 'opening grid mosaic file') @@ -1200,6 +1507,84 @@ subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, & error = nf90_close(ncid) end subroutine get_model_latlons + + !---------------------------------------------------------------------------------------- +! For grids with equal cell sizes (e.g., lambert conformal), get lat and on of the grid +! cell corners +!---------------------------------------------------------------------------------------- + + subroutine get_cell_corners( latitude, longitude, latitude_sw, longitude_sw, dx,clb,cub) + implicit none + + real(esmf_kind_r8), intent(in) :: latitude(i_input,j_input) + real(esmf_kind_r8), intent(inout), pointer :: latitude_sw(:,:) + real(esmf_kind_r8), intent(in) :: longitude(i_input, j_input) + real(esmf_kind_r8), intent(inout), pointer :: longitude_sw(:,:) + real(esmf_kind_r8), intent(in) :: dx !grid cell side size (m) + + integer, intent(in) :: clb(2), cub(2) + + real(esmf_kind_r8) :: lat1, lon1, lat2, lon2, d, brng + + + real(esmf_kind_r8), parameter :: pi = 3.14159265359 + real(esmf_kind_r8), parameter :: R = 6370000.0 + real(esmf_kind_r8), parameter :: bearingInDegrees = 135.0 + + integer :: i, j + + d = sqrt((dx**2.0_esmf_kind_r8)/2.0_esmf_kind_r8) + + do j = clb(2),cub(2) + do i = clb(1), cub(1) + + if (j == jp1_input .and. i == ip1_input) then + lat1 = latitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 ) + lon1 = longitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 ) + brng = 315.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 + lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); + lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); + latitude_sw(ip1_input,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi + longitude_sw(ip1_input,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi + cycle + endif + + if (i == ip1_input) then + brng = 225.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 + lat1 = latitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 ) + lon1 = longitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 ) + lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); + lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); + latitude_sw(ip1_input,j) = lat2 * 180.0_esmf_kind_r8 / pi + longitude_sw(ip1_input,j) = lon2 * 180.0_esmf_kind_r8 / pi + cycle + endif + + if (j == jp1_input) then + brng = 45.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 + lat1 = latitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 ) + lon1 = longitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 ) + lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); + lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); + latitude_sw(i,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi + longitude_sw(i,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi + cycle + endif + + lat1 = latitude(i,j) * ( pi / 180.0_esmf_kind_r8 ) + lon1 = longitude(i,j) * ( pi / 180.0_esmf_kind_r8 ) + + brng = bearingInDegrees * ( pi / 180.0_esmf_kind_r8 ); + lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); + lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); + + latitude_sw(i,j) = lat2 * 180.0_esmf_kind_r8 / pi + longitude_sw(i,j) = lon2 * 180.0_esmf_kind_r8 / pi + + enddo + enddo + + end subroutine get_cell_corners !----------------------------------------------------------------------- ! Read the model land mask and terrain for a single tile. diff --git a/sorc/chgres_cube.fd/program_setup.f90 b/sorc/chgres_cube.fd/program_setup.f90 index ed6935457..6340e13b8 100644 --- a/sorc/chgres_cube.fd/program_setup.f90 +++ b/sorc/chgres_cube.fd/program_setup.f90 @@ -36,10 +36,6 @@ module program_setup ! target grids. ! fix_dir_target_grid Directory containing target grid ! pre-computed fixed data (ex: soil type) -! grib2_file_input_grid File name of grib2 input data. -! Assumes atmospheric and surface data -! are in a single file. 'grib2' input -! type only. ! halo_blend Number of row/cols of blending halo, ! where model tendencies and lateral ! boundary tendencies are applied. @@ -56,7 +52,7 @@ module program_setup ! nemsio files; ! (4) "gaussian_netcdf" for fv3 gaussian ! netcdf files. -! (5) "grib2" for fv3gfs grib2 files. +! (5) "grib2" for grib2 files. ! (6) "gfs_gaussian_nemsio" for spectral gfs ! gaussian nemsio files ! (7) "gfs_sigio" for spectral gfs @@ -98,16 +94,71 @@ module program_setup ! These names will be used to identify ! the tracer records in the output files. ! Follows the convention in the field table. +! FOR GRIB2 FILES: Not used. Tracers instead taken +! from the varmap file. ! tracers_input Name of each atmos tracer record in ! the input file. May be different from ! value in 'tracers'. +! FOR GRIB2 FILES: Not used. Tracers instead taken +! from the varmap file. ! use_thomp_mp_climo When true, read and process Thompson ! MP climatological tracers. False, ! when 'thomp_mp_climo_file' is NULL. ! vcoord_file_target_grid Vertical coordinate definition file ! wltsmc_input/target Wilting point soil moisture content ! input/target grids -! +! +! nsoill_out Number of soil levels desired in the output data. +! chgres_cube can interpolate from 9 input to 4 output +! levels. +! DEFAULT: 4 +! +! Variables that are relevant only for "grib2" input type: +! +! grib2_file_input_grid REQUIRED. File name of grib2 input data. +! Assumes atmospheric and surface data are in a single +! file. +! +! varmap_file REQUIRED. Full path of the relevant varmap file. +! +! external_model The model that the input data is derived from. Current +! supported options are: "GFS", "HRRR", "NAM", "RAP". +! Default: "GFS" +! +! vgtyp_from_climo If false, interpolate vegetation type from the input +! data to the target grid instead of using data from +! static data. Use with caution as vegetation categories +! can vary. +! Default: False +! +! sotyp_from_climo If false, interpolate soil type from the input +! data to the target grid instead of using data from +! static data. Use with caution as the code assumes +! input soil type use STATSGO soil categories. +! Default: False +! +! vgfrc_from_climo If false, interpolate vegetation fraction from the input +! data to the target grid instead of using data from +! static data. Use with caution as vegetation categories +! can vary. +! Default: False +! +! minmax_vgfrc_from_climo If false, interpolate min/max vegetation fraction from +! the input data to the target grid instead of using data +! from static data. Use with caution as vegetation +! categories can vary. +! Default: False +! +! lai_from_climo If false, interpolate leaf area index from the input +! data to the target grid instead of using data from +! static data. +! Default: False +! +! tg3_from_soil If false, use lowest level soil temperature for the +! base soil temperature instead of using data from +! static data. +! Default: False +! !-------------------------------------------------------------------------- implicit none @@ -124,6 +175,7 @@ module program_setup character(len=500), public :: mosaic_file_target_grid = "NULL" character(len=500), public :: nst_files_input_grid = "NULL" character(len=500), public :: grib2_file_input_grid = "NULL" + character(len=500), public :: geogrid_file_input_grid = "NULL" character(len=500), public :: orog_dir_input_grid = "NULL" character(len=500), public :: orog_files_input_grid(6) = "NULL" character(len=500), public :: orog_dir_target_grid = "NULL" @@ -134,7 +186,12 @@ module program_setup character(len=6), public :: cres_target_grid = "NULL" character(len=500), public :: atm_weight_file="NULL" character(len=25), public :: input_type="restart" - character(len=20), public :: phys_suite="GFS" !Default to gfs physics suite + character(len=20), public :: external_model="GFS" !Default assume gfs data + + + + character(len=500), public :: fix_dir_input_grid = "NULL" + integer, parameter, public :: max_tracers=100 integer, public :: num_tracers, num_tracers_input @@ -154,11 +211,20 @@ module program_setup integer, public :: regional = 0 integer, public :: halo_bndy = 0 integer, public :: halo_blend = 0 + integer, public :: nsoill_out = 4 logical, public :: convert_atm = .false. logical, public :: convert_nst = .false. logical, public :: convert_sfc = .false. - + + ! Options for replacing vegetation/soil type, veg fraction, and lai with data from the grib2 file + ! Default is to use climatology instead + logical, public :: vgtyp_from_climo = .true. + logical, public :: sotyp_from_climo = .true. + logical, public :: vgfrc_from_climo = .true. + logical, public :: minmax_vgfrc_from_climo = .true. + logical, public :: lai_from_climo = .true. + logical, public :: tg3_from_soil = .false. logical, public :: use_thomp_mp_climo=.false. real, allocatable, public :: drysmc_input(:), drysmc_target(:) @@ -199,16 +265,26 @@ subroutine read_setup_namelist atm_core_files_input_grid, & atm_tracer_files_input_grid, & grib2_file_input_grid, & + geogrid_file_input_grid, & data_dir_input_grid, & vcoord_file_target_grid, & cycle_mon, cycle_day, & cycle_hour, convert_atm, & convert_nst, convert_sfc, & + vgtyp_from_climo, & + sotyp_from_climo, & + vgfrc_from_climo, & + minmax_vgfrc_from_climo, & + lai_from_climo, tg3_from_soil, & regional, input_type, & + external_model, & atm_weight_file, tracers, & - tracers_input,phys_suite, & + tracers_input, & halo_bndy, & - halo_blend, thomp_mp_climo_file + halo_blend, & + fix_dir_input_grid, & + nsoill_out, & + thomp_mp_climo_file print*,"- READ SETUP NAMELIST" @@ -219,7 +295,7 @@ subroutine read_setup_namelist close (41) call to_lower(input_type) - call to_upper(phys_suite) +! call to_upper(phys_suite) orog_dir_target_grid = trim(orog_dir_target_grid) // '/' orog_dir_input_grid = trim(orog_dir_input_grid) // '/' @@ -229,7 +305,7 @@ subroutine read_setup_namelist !------------------------------------------------------------------------- is = index(mosaic_file_target_grid, "/", .true.) - ie = index(mosaic_file_target_grid, "_mosaic") + ie = index(mosaic_file_target_grid, "mosaic") - 1 if (is == 0 .or. ie == 0) then call error_handler("CANT DETERMINE CRES FROM MOSAIC FILE.", 1) @@ -293,6 +369,45 @@ subroutine read_setup_namelist case default call error_handler("UNRECOGNIZED INPUT DATA TYPE.", 1) end select + +!------------------------------------------------------------------------- +! Ensure proper file variable provided for grib2 input +!------------------------------------------------------------------------- + + if (trim(input_type) == "grib2") then + if (trim(grib2_file_input_grid) == "NULL" .or. trim(grib2_file_input_grid) == "") then + call error_handler("FOR GRIB2 DATA, PLEASE PROVIDE GRIB2_FILE_INPUT_GRID", 1) + endif + endif + + !------------------------------------------------------------------------- +! For grib2 input, warn about possibly unsupported external model types +!------------------------------------------------------------------------- + + if (trim(input_type) == "grib2") then + if (.not. any((/character(4)::"GFS","NAM","RAP","HRRR"/)==trim(external_model))) then + call error_handler( "KNOWN SUPPORTED external_model INPUTS ARE GFS, NAM, RAP, AND HRRR. " // & + "IF YOU WISH TO PROCESS GRIB2 DATA FROM ANOTHER MODEL, YOU MAY ATTEMPT TO DO SO AT YOUR OWN RISK. " // & + "ONE WAY TO DO THIS IS PROVIDE NAM FOR external_model AS IT IS A RELATIVELY STRAIGHT-" // & + "FORWARD REGIONAL GRIB2 FILE. YOU MAY ALSO COMMENT OUT THIS ERROR MESSAGE IN " // & + "program_setup.f90 LINE 389. NO GUARANTEE IS PROVIDED THAT THE CODE WILL WORK OR "// & + "THAT THE RESULTING DATA WILL BE CORRECT OR WORK WITH THE ATMOSPHERIC MODEL.", 1) + endif + endif + +!------------------------------------------------------------------------- +! For grib2 hrrr input without geogrid file input, warn that soil moisture interpolation +! will be less accurate +!------------------------------------------------------------------------- + + if (trim(input_type) == "grib2" .and. trim(external_model)=="HRRR") then + if (trim(geogrid_file_input_grid) == "NULL" .or. trim(grib2_file_input_grid) == "") then + print*, "HRRR DATA DOES NOT CONTAIN SOIL TYPE INFORMATION. WITHOUT & + GEOGRID_FILE_INPUT_GRID SPECIFIED, SOIL MOISTURE INTERPOLATION MAY BE LESS & + ACCURATE. " + endif + endif + return if (trim(thomp_mp_climo_file) /= "NULL") then use_thomp_mp_climo=.true. @@ -346,6 +461,12 @@ subroutine read_varmap if(trim(var_type(k))=='T') then num_tracers = num_tracers + 1 tracers_input(num_tracers)=chgres_var_names(k) + if ((trim(chgres_var_names(k)) == "ice_aero" .or. trim(chgres_var_names(k)) == "liq_aero") .and. & + trim(thomp_mp_climo_file) .ne. "NULL" .and. trim(input_type) == "grib2") then + call error_handler("VARMAP TABLE CONTAINS TRACER ENTRIES FOR THOMPSON AEROSOLS liq_aero or "// & + "ice_aero. REMOVE THESE ENTRIES OR REMOVE THE NAMELIST ENTRY FOR "// & + "thomp_mp_climo_file AND TRY AGAIN.",1) + endif endif enddo close(14) diff --git a/sorc/chgres_cube.fd/search_util.f90 b/sorc/chgres_cube.fd/search_util.f90 index 58bb3cf37..e5bec9c53 100644 --- a/sorc/chgres_cube.fd/search_util.f90 +++ b/sorc/chgres_cube.fd/search_util.f90 @@ -19,7 +19,7 @@ module search_util contains - subroutine search (field, mask, idim, jdim, tile, field_num, latitude) + subroutine search (field, mask, idim, jdim, tile, field_num, latitude, terrain_land, soilt_climo) !----------------------------------------------------------------------- ! Replace undefined values on the model grid with a valid value at @@ -43,6 +43,8 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) integer(esmf_kind_i8), intent(in) :: mask(idim,jdim) real(esmf_kind_r8), intent(in), optional :: latitude(idim,jdim) + real(esmf_kind_r8), intent(in), optional :: terrain_land(idim,jdim) + real(esmf_kind_r8), intent(in), optional :: soilt_climo(idim,jdim) real(esmf_kind_r8), intent(inout) :: field(idim,jdim) @@ -53,6 +55,7 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) real :: default_value real(esmf_kind_r8) :: field_save(idim,jdim) + integer :: repl_nearby, repl_default !----------------------------------------------------------------------- ! Set default value. @@ -90,6 +93,18 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) default_value = 0.0 case (224) ! soil type, flag value to turn off soil moisture rescaling. default_value = -99999.9 + case (225) ! vegetation type, flag value to be replaced + default_value = -99999.9 + case (226) ! vegetation fraction, flag value to be replaced + default_value = 0.5 + case (227) ! max vegetation fraction, flag value to be replaced + default_value = 0.5 + case (228) ! min vegetation fraction, flag value to be replaced + default_value = 0.5 + case (229) ! lai, flag value to be replaced + default_value = 1.0 + case (230) ! soil type on the input grid + default_value = 11.0 case default print*,'- FATAL ERROR. UNIDENTIFIED FIELD NUMBER : ', field call mpi_abort(mpi_comm_world, 77, ierr) @@ -100,9 +115,10 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) !----------------------------------------------------------------------- field_save = field - + repl_nearby = 0 + repl_default = 0 !$OMP PARALLEL DO DEFAULT(NONE), & -!$OMP SHARED(IDIM,JDIM,MASK,FIELD_SAVE,FIELD,TILE,LATITUDE,DEFAULT_VALUE,FIELD_NUM), & +!$OMP SHARED(IDIM,JDIM,MASK,FIELD_SAVE,FIELD,TILE,LATITUDE,DEFAULT_VALUE,FIELD_NUM,REPL_NEARBY,REPL_DEFAULT,SOILT_CLIMO,TERRAIN_LAND), & !$OMP PRIVATE(I,J,KRAD,ISTART,IEND,JSTART,JEND,II,JJ) J_LOOP : do j = 1, jdim @@ -132,7 +148,10 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) if (mask(ii,jj) == 1 .and. field_save(ii,jj) > -9999.0) then field(i,j) = field_save(ii,jj) - write(6,100) tile,i,j,ii,jj,field(i,j) + ! write(6,100) field_num,tile,i,j,ii,jj,field(i,j) + ! When using non-GFS data, there are a lot of these print statements even + ! when everything is working correctly. Count instead of printing each + repl_nearby = repl_nearby + 1 cycle I_LOOP endif @@ -148,22 +167,39 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) elseif (field_num == 91) then ! sea ice fract if (abs(latitude(i,j)) > 55.0) then field(i,j) = default_value + repl_default = repl_default + 1 else field(i,j) = 0.0 + repl_default = repl_default + 1 endif + elseif (field_num == 7 .and. PRESENT(terrain_land)) then + ! Terrain heights for isolated landice points never get a correct value, so replace + ! with terrain height from the input grid interpolated to the target grid + field(i,j) = terrain_land(i,j) + repl_default = repl_default + 1 + elseif (field_num == 224 .and. PRESENT(soilt_climo)) then + ! When using input soil type fields instead of climatological data on the + ! target grid, isolated land locations that exist in the target grid but + ! not the input grid don't receiving proper soil type information, so replace + ! with climatological values + field(i,j) = soilt_climo(i,j) + repl_default = repl_default + 1 else field(i,j) = default_value ! Search failed. Use default value. + repl_default = repl_default + 1 endif - write(6,101) tile,i,j,field(i,j) + !write(6,101) field_num,tile,i,j,field(i,j) endif enddo I_LOOP enddo J_LOOP !$OMP END PARALLEL DO - 100 format(1x,"- MISSING POINT TILE: ",i2," I/J: ",i5,i5," SET TO VALUE AT: ",i5,i5,". NEW VALUE IS: ",f8.3) - 101 format(1x,"- MISSING POINT TILE: ",i2," I/J: ",i5,i5," SET TO DEFAULT VALUE OF: ",f8.3) +! 100 format(1x,"- MISSING POINT FIELD ",i4," TILE: ",i2," I/J: ",i5,i5," SET TO VALUE AT: ",i5,i5,". NEW VALUE IS: ",f8.3) +! 101 format(1x,"- MISSING POINT FIELD ",i4," TILE: ",i2," I/J: ",i5,i5," SET TO DEFAULT VALUE OF: ",f8.3) + print*, "- TOTAL POINTS FOR VAR ", field_num, " REPLACED BY NEARBY VALUES: ", repl_nearby + print*, "- TOTAL POINTS FOR VAR ", field_num, " REPLACED BY DEFAULT VALUE: ", repl_default end subroutine search diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index ce1b64198..7701534aa 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -68,6 +68,8 @@ module surface ! friction velocity type(esmf_field), public :: z0_target_grid ! roughness length + type(esmf_field), public :: lai_target_grid + ! leaf area index ! nst fields type(esmf_field), public :: c_d_target_grid @@ -96,6 +98,9 @@ module surface type(esmf_field) :: terrain_from_input_grid ! terrain height interpolated ! from input grid + type(esmf_field) :: terrain_from_input_grid_land + ! terrain height interpolated + ! from input grid at all land points real, parameter, private :: blim = 5.5 ! soil 'b' parameter limit @@ -120,9 +125,10 @@ subroutine surface_driver(localpet) read_input_nst_data use program_setup, only : calc_soil_params_driver, & - convert_nst - - use static_data, only : get_static_fields, & + convert_nst, & + vgtyp_from_climo, & + sotyp_from_climo + use static_data, only : get_static_fields, & cleanup_static_fields implicit none @@ -164,26 +170,54 @@ subroutine surface_driver(localpet) !----------------------------------------------------------------------- if (convert_nst) call create_nst_esmf_fields + +!----------------------------------------------------------------------- +! Adjust soil levels of input grid !! not implemented yet +!----------------------------------------------------------------------- + + call adjust_soil_levels(localpet) !----------------------------------------------------------------------- ! Horizontally interpolate fields. !----------------------------------------------------------------------- call interp(localpet) - + + !if (.not. (vgtyp_from_climo .or. sotyp_from_climo)) then + !--------------------------------------------------------------------------------------------- + ! Check for points where smois is too high to be a land point at a land point + !--------------------------------------------------------------------------------------------- + ! + !call check_smois_water + !endif + !--------------------------------------------------------------------------------------------- -! Adjust soil/landice column temperatures for any change in elevation between the +! Adjust soil/landice column temperatures for any change in elevation between +! the ! input and target grids. !--------------------------------------------------------------------------------------------- call adjust_soilt_for_terrain - + !--------------------------------------------------------------------------------------------- ! Rescale soil moisture for changes in soil type between the input and target grids. !--------------------------------------------------------------------------------------------- call rescale_soil_moisture + + !if (.not. (vgtyp_from_climo .or. sotyp_from_climo)) then + !--------------------------------------------------------------------------------------------- + ! Check soil moisture again for mismatches after rescale_soil_moisture subroutine + !--------------------------------------------------------------------------------------------- + ! call check_smois_land + + !--------------------------------------------------------------------------------------------- + ! Replacing values various land surface parameters at points identified as mis-masked in + ! check_smois_land + !--------------------------------------------------------------------------------------------- + ! call replace_land_sfcparams(localpet) + !endif !--------------------------------------------------------------------------------------------- ! Compute liquid portion of total soil moisture. !--------------------------------------------------------------------------------------------- @@ -283,7 +317,11 @@ subroutine interp(localpet) xzts_input_grid, & z_c_input_grid, & zm_input_grid, terrain_input_grid, & - veg_type_landice_input + veg_type_landice_input, & + veg_greenness_input_grid, & + max_veg_greenness_input_grid, & + min_veg_greenness_input_grid, & + lai_input_grid use model_grid, only : input_grid, target_grid, & i_target, j_target, & @@ -293,10 +331,22 @@ subroutine interp(localpet) seamask_target_grid, & latitude_target_grid - use program_setup, only : convert_nst, input_type - + use program_setup, only : convert_nst, & + vgtyp_from_climo, & + sotyp_from_climo, & + vgfrc_from_climo, & + minmax_vgfrc_from_climo, & + lai_from_climo, & + tg3_from_soil, & + external_model, & + input_type + use static_data, only : veg_type_target_grid, & - soil_type_target_grid + soil_type_target_grid, & + veg_greenness_target_grid, & + substrate_temp_target_grid,& + min_veg_greenness_target_grid,& + max_veg_greenness_target_grid use search_util @@ -356,6 +406,10 @@ subroutine interp(localpet) real(esmf_kind_r8), pointer :: landmask_input_ptr(:,:) real(esmf_kind_r8), pointer :: veg_type_input_ptr(:,:) real(esmf_kind_r8), allocatable :: veg_type_target_one_tile(:,:) + real(esmf_kind_r8), pointer :: veg_greenness_target_ptr(:,:) + real(esmf_kind_r8), pointer :: min_veg_greenness_target_ptr(:,:) + real(esmf_kind_r8), pointer :: max_veg_greenness_target_ptr(:,:) + real(esmf_kind_r8), pointer :: lai_target_ptr(:,:) type(esmf_regridmethod_flag) :: method type(esmf_routehandle) :: regrid_bl_no_mask @@ -459,8 +513,6 @@ subroutine interp(localpet) call error_handler("IN FieldRegridRelease", rc) !----------------------------------------------------------------------- -! Next, determine the sea ice fraction on target grid. -! ! First, set the mask on the target and input grids. !----------------------------------------------------------------------- @@ -485,8 +537,12 @@ subroutine interp(localpet) farrayPtr=seamask_target_ptr, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) - - mask_target_ptr = seamask_target_ptr + + print*,"- CALL FieldGet FOR TARGET GRID LANDMASK." + call ESMF_FieldGet(landmask_target_grid, & + farrayPtr=landmask_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) print*,"- CALL GridAddItem FOR INPUT GRID SEAMASK." call ESMF_GridAddItem(input_grid, & @@ -507,25 +563,176 @@ subroutine interp(localpet) farrayPtr=mask_input_ptr, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN GridGetItem", rc) - - mask_input_ptr = 1 - where (nint(landmask_input_ptr) == 1) mask_input_ptr = 0 - -!----------------------------------------------------------------------- -! Interpolate. -!----------------------------------------------------------------------- - - if (localpet == 0) then + + if (localpet == 0) then allocate(data_one_tile(i_target,j_target)) - allocate(data_one_tile2(i_target,j_target)) allocate(data_one_tile_3d(i_target,j_target,lsoil_target)) allocate(mask_target_one_tile(i_target,j_target)) else allocate(data_one_tile(0,0)) - allocate(data_one_tile2(0,0)) allocate(data_one_tile_3d(0,0,0)) allocate(mask_target_one_tile(0,0)) endif + + !----------------------------------------------------------------------- + ! Interpolate vegetation type to target grid if chosen in namelist and terrain + ! for use in replacing isolated bad terrain values + !----------------------------------------------------------------------- + + method=ESMF_REGRIDMETHOD_NEAREST_STOD + + isrctermprocessing = 1 + + mask_input_ptr = 0 + where (nint(landmask_input_ptr) == 1) mask_input_ptr = 1 + + mask_target_ptr = 0 + where (landmask_target_ptr == 1) mask_target_ptr = 1 + + print*,"- CALL FieldCreate FOR TERRAIN FROM INPUT GRID LAND." + terrain_from_input_grid_land = ESMF_FieldCreate(target_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldRegridStore for land fields." + call ESMF_FieldRegridStore(terrain_input_grid, & + terrain_from_input_grid_land, & + srcmaskvalues=(/0/), & + dstmaskvalues=(/0/), & + polemethod=ESMF_POLEMETHOD_NONE, & + srctermprocessing=isrctermprocessing, & + unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & + normtype=ESMF_NORMTYPE_FRACAREA, & + routehandle=regrid_all_land, & + 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 TERRAIN." + call ESMF_FieldRegrid(terrain_input_grid, & + terrain_from_input_grid_land, & + routehandle=regrid_all_land, & + 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 FieldGet FOR terrain from input grid at land." + call ESMF_FieldGet(terrain_from_input_grid_land, & + farrayPtr=terrain_from_input_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + l = lbound(unmapped_ptr) + u = ubound(unmapped_ptr) + + do ij = l(1), u(1) + call ij_to_i_j(unmapped_ptr(ij), i_target, j_target, i, j) + terrain_from_input_ptr(i,j) = -9999.9 + enddo + nullify(terrain_from_input_ptr) + + do tile = 1, num_tiles_target_grid + + print*,"- CALL FieldGather FOR TARGET LANDMASK TILE: ", tile + call ESMF_FieldGather(landmask_target_grid, mask_target_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 TERRAIN FROM INPUT GRID: ", tile + call ESMF_FieldGather(terrain_from_input_grid, data_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 + 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 + call search(data_one_tile, land_target_one_tile, i_target, j_target, tile, 7) + deallocate(land_target_one_tile) + endif + + print*,"- CALL FieldScatter FOR TERRAIN FROM INPUT GRID: ", tile + call ESMF_FieldScatter(terrain_from_input_grid, data_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) + enddo + + if(.not. vgtyp_from_climo) then + + do tile = 1, num_tiles_target_grid + print*,"-CALL FieldGather VEG TYPE TARGET GRID" + call ESMF_FieldGather(veg_type_target_grid, data_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) + + data_one_tile(:,:) = 0 + enddo + + print*,"- CALL FieldRegrid VEG TYPE." + call ESMF_FieldRegrid(veg_type_input_grid, & + veg_type_target_grid, & + routehandle=regrid_all_land, & + 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 FieldGet FOR TARGET grid veg type." + call ESMF_FieldGet(veg_type_target_grid, & + farrayPtr=veg_type_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + l = lbound(unmapped_ptr) + u = ubound(unmapped_ptr) + + do ij = l(1), u(1) + call ij_to_i_j(unmapped_ptr(ij), i_target, j_target, i, j) + veg_type_target_ptr(i,j) = -9999.9 + enddo + + do tile = 1, num_tiles_target_grid + print*,"- CALL FieldGather FOR TARGET GRID VEG TYPE TILE: ", tile + call ESMF_FieldGather(veg_type_target_grid, data_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 TARGET LANDMASK TILE: ", tile + call ESMF_FieldGather(landmask_target_grid, mask_target_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 + 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 + call search(data_one_tile, land_target_one_tile, i_target, j_target, tile, 225) + deallocate(land_target_one_tile) + endif + + print*,"- CALL FieldScatter FOR TARGET GRID VEG TYPE: ", tile + call ESMF_FieldScatter(veg_type_target_grid, data_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) + enddo + nullify(veg_type_target_ptr) + endif + print*,"- CALL FieldRegridRelease." + call ESMF_FieldRegridRelease(routehandle=regrid_all_land, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldRegridRelease", rc) + +!----------------------------------------------------------------------- +! Next, determine the sea ice fraction on target grid. +! Interpolate. +!----------------------------------------------------------------------- + + mask_input_ptr = 1 + where (nint(landmask_input_ptr) == 1) mask_input_ptr = 0 + + mask_target_ptr = seamask_target_ptr method=ESMF_REGRIDMETHOD_CONSERVE @@ -582,7 +789,7 @@ subroutine interp(localpet) call ESMF_FieldGather(seaice_fract_target_grid, data_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 TARGET GRID MASK TILE: ", tile call ESMF_FieldGather(seamask_target_grid, mask_target_one_tile, rootPet=0, tile=tile, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -602,7 +809,7 @@ subroutine interp(localpet) call ESMF_FieldGather(landmask_target_grid, mask_target_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 do j = 1, j_target @@ -713,7 +920,7 @@ subroutine interp(localpet) print*,"- CALL FieldGet FOR TARGET grid snow depth." call ESMF_FieldGet(snow_depth_target_grid, & farrayPtr=snow_depth_target_ptr, 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 FieldGet", rc) print*,"- CALL Field_Regrid for snow liq equiv." @@ -769,6 +976,9 @@ subroutine interp(localpet) call error_handler("IN FieldGather", rc) if (localpet == 0) then + ! I don't think is needed anymore with the more recent fixes to fill values in input_data + !if (count(landmask_target_ptr == 2) == 0) data_one_tile(:,:) =0.0_esmf_kind_r8 + where(mask_target_one_tile == 1) mask_target_one_tile = 0 where(mask_target_one_tile == 2) mask_target_one_tile = 1 call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 92) @@ -1795,9 +2005,11 @@ subroutine interp(localpet) if (localpet == 0) then allocate (veg_type_target_one_tile(i_target,j_target)) allocate (land_target_one_tile(i_target,j_target)) + allocate (data_one_tile2(i_target,j_target)) else allocate (veg_type_target_one_tile(0,0)) allocate (land_target_one_tile(0,0)) + allocate (data_one_tile2(0,0)) endif do tile = 1, num_tiles_target_grid @@ -1827,9 +2039,14 @@ subroutine interp(localpet) call ESMF_FieldGather(terrain_from_input_grid, data_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 TERRAIN FROM INPUT GRID LAND, TILE: ", tile + call ESMF_FieldGather(terrain_from_input_grid_land, 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) if (localpet == 0) then - call search(data_one_tile, land_target_one_tile, i_target, j_target, tile, 7) + call search(data_one_tile, land_target_one_tile, i_target, j_target, tile, 7, terrain_land=data_one_tile2) endif print*,"- CALL FieldScatter FOR TERRAIN FROM INPUT GRID, TILE: ", tile @@ -1937,6 +2154,45 @@ subroutine interp(localpet) 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) + + if (.not. vgfrc_from_climo) then + print*,"- CALL Field_Regrid for veg greenness over land." + call ESMF_FieldRegrid(veg_greenness_input_grid, & + veg_greenness_target_grid, & + routehandle=regrid_land, & + 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) + endif + + if (.not. minmax_vgfrc_from_climo) then + print*,"- CALL Field_Regrid for max veg greenness over land." + call ESMF_FieldRegrid(max_veg_greenness_input_grid, & + max_veg_greenness_target_grid, & + routehandle=regrid_land, & + 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 Field_Regrid for min veg greenness over land." + call ESMF_FieldRegrid(min_veg_greenness_input_grid, & + min_veg_greenness_target_grid, & + routehandle=regrid_land, & + 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) + endif + + if (.not. lai_from_climo) then + print*,"- CALL Field_Regrid for leaf area index over land." + call ESMF_FieldRegrid(lai_input_grid, & + lai_target_grid, & + routehandle=regrid_land, & + 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) + + endif print*,"- CALL FieldGet FOR TARGET grid total soil moisture over land." call ESMF_FieldGet(soilm_tot_target_grid, & @@ -1968,9 +2224,38 @@ subroutine interp(localpet) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGet", rc) + if (.not. vgfrc_from_climo) then + print*,"- CALL FieldGet FOR TARGET veg greenness." + call ESMF_FieldGet(veg_greenness_target_grid, & + farrayPtr=veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + endif + + if (.not. minmax_vgfrc_from_climo) then + print*,"- CALL FieldGet FOR TARGET max veg greenness." + call ESMF_FieldGet(max_veg_greenness_target_grid, & + farrayPtr=max_veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET min veg greenness." + call ESMF_FieldGet(min_veg_greenness_target_grid, & + farrayPtr=min_veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldGet", rc) + endif + + if (.not. lai_from_climo) then + print*,"- CALL FieldGet FOR TARGET lai." + call ESMF_FieldGet(lai_target_grid, & + farrayPtr=lai_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__))& + call error_handler("IN FieldGet", rc) + endif + l = lbound(unmapped_ptr) u = ubound(unmapped_ptr) - do ij = l(1), u(1) call ij_to_i_j(unmapped_ptr(ij), i_target, j_target, i, j) soilm_tot_target_ptr(i,j,:) = -9999.9 @@ -1978,6 +2263,10 @@ subroutine interp(localpet) skin_temp_target_ptr(i,j) = -9999.9 terrain_from_input_ptr(i,j) = -9999.9 soil_type_from_input_ptr(i,j) = -9999.9 + veg_greenness_target_ptr(i,j) = -9999.9 + max_veg_greenness_target_ptr(i,j) = -9999.9 + min_veg_greenness_target_ptr(i,j) = -9999.9 + lai_target_ptr(i,j) = -9999.9 enddo if (localpet == 0) then @@ -2038,24 +2327,103 @@ subroutine interp(localpet) call error_handler("IN FieldGather", rc) !--------------------------------------------------------------------------------------- -! grib2 data does not have soil type. Set soil type interpolated from input +! Some grib2 data does not have soil type. Set soil type interpolated from input ! grid to the target (model) grid soil type. This turns off the soil moisture ! rescaling. !--------------------------------------------------------------------------------------- - if (localpet == 0) then - if (trim(input_type) .ne. "grib2") then + if (.not. sotyp_from_climo) then + if (localpet==0) then + call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 224,soilt_climo=data_one_tile2) + endif + else + if (localpet == 0 .and. maxval(data_one_tile) > 0 .and. (trim(external_model) .ne. "GFS" .or. trim(input_type) .ne. "grib2")) then + ! If soil type from the input grid has any non-zero points then soil type must exist for + ! use call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 224) - else + elseif (localpet == 0) then data_one_tile = data_one_tile2 endif endif + + if (.not. sotyp_from_climo) then + print*,"- CALL FieldScatter FOR SOIL TYPE TARGET GRID, TILE: ", tile + call ESMF_FieldScatter(soil_type_target_grid, data_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 print*,"- CALL FieldScatter FOR SOIL TYPE FROM INPUT GRID, TILE: ", tile call ESMF_FieldScatter(soil_type_from_input_grid, data_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) + + + + if (.not. vgfrc_from_climo) then + print*,"- CALL FieldGather FOR TARGET GRID VEG GREENNESS, TILE: ", tile + call ESMF_FieldGather(veg_greenness_target_grid, data_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 .and. maxval(data_one_tile) > 0.0) then + call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 226) + endif + + print*,"- CALL FieldScatter FOR VEG GREENNESS TARGET GRID, TILE: ", tile + call ESMF_FieldScatter(veg_greenness_target_grid, data_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 + + if (.not. minmax_vgfrc_from_climo) then + print*,"- CALL FieldGather FOR TARGET GRID MAX VEG GREENNESS, TILE: ", tile + call ESMF_FieldGather(max_veg_greenness_target_grid, data_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 .and. maxval(data_one_tile) > 0.0) then + call search(data_one_tile, mask_target_one_tile, i_target, j_target,tile, 227) + endif + + print*,"- CALL FieldScatter FOR MAX VEG GREENNESS TARGET GRID, TILE: ", tile + call ESMF_FieldScatter(max_veg_greenness_target_grid, data_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 FieldGather FOR TARGET GRID MIN VEG GREENNESS, TILE: ", tile + call ESMF_FieldGather(min_veg_greenness_target_grid, data_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 .and. maxval(data_one_tile) > 0.0) then + call search(data_one_tile, mask_target_one_tile, i_target, j_target,tile,228) + endif + + + print*,"- CALL FieldScatter FOR MIN VEG GREENNESS TARGET GRID, TILE: ",tile + call ESMF_FieldScatter(min_veg_greenness_target_grid, data_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 + + if (.not. lai_from_climo) then + print*,"- CALL FieldGather FOR TARGET GRID LEAF AREA INDEX, TILE: ", tile + call ESMF_FieldGather(lai_target_grid, data_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 .and. maxval(data_one_tile) > 0.0) then + call search(data_one_tile, mask_target_one_tile, i_target, j_target,tile, 229) + endif + + print*,"- CALL FieldScatter FOR LEAF AREA INDEX TARGET GRID, TILE: ", tile + call ESMF_FieldScatter(lai_target_grid, data_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 + 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) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -2091,6 +2459,13 @@ subroutine interp(localpet) call ESMF_FieldScatter(soil_temp_target_grid, data_one_tile_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) + + if (tg3_from_soil) then + print*,"- CALL FieldScatter FOR TARGET GRID SUBSTRATE TEMPERATURE, TILE: ", tile + call ESMF_FieldScatter(substrate_temp_target_grid, data_one_tile_3d(:,:,lsoil_target), 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 enddo @@ -2233,6 +2608,518 @@ subroutine calc_liq_soil_moisture end subroutine calc_liq_soil_moisture +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! This subroutine was previously used to correct for points with soil moisture that was +!! too high for land points. This happened when there was a mismatch between input and +!! target landmasks, specifically when input vegetation type was used to replace that +!! in the geogrid file. The functions performed by this subroutine are now performed +!! when the data is read in. Note that the target grid landmask is no longer modified +!! anywhere in the code, and the input data is modified instead. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!subroutine check_smois_water +! +!use model_grid, only : landmask_target_grid +! +!use static_data, only : veg_type_target_grid, veg_greenness_target_grid, & +! soil_type_target_grid, max_veg_greenness_target_grid,& +! min_veg_greenness_target_grid, mxsno_albedo_target_grid, & +! alvsf_target_grid,alvwf_target_grid,& +! alnsf_target_grid,alnwf_target_grid +! +! implicit none +! +! integer :: clb(3), cub(3), i, j, rc +! +! integer(esmf_kind_r8), pointer :: landmask_ptr(:,:) +! +! real(esmf_kind_r8), pointer :: soilm_target_ptr(:,:,:), & +! alvsf_target_ptr(:,:), & +! alnsf_target_ptr(:,:), & +! alvwf_target_ptr(:,:), & +! alnwf_target_ptr(:,:), & +! veg_greenness_target_ptr(:,:), & +! min_veg_greenness_target_ptr(:,:), & +! max_veg_greenness_target_ptr(:,:), & +! canopy_mc_target_ptr(:,:), & +! mxsno_albedo_target_ptr(:,:), & +! soil_type_target_ptr(:,:), & +! veg_type_target_ptr(:,:) +! +! +! print*,"- CALL FieldGet FOR TARGET GRID LAND-SEA MASK." +! call ESMF_FieldGet(landmask_target_grid, & +! farrayPtr=landmask_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR SOIL MOIS TARGET GRID." +! call ESMF_FieldGet(soilm_tot_target_grid, & +! computationalLBound=clb, & +! computationalUBound=cub, & +! farrayPtr=soilm_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID SOIL TYPE." +! call ESMF_FieldGet(soil_type_target_grid, & +! farrayPtr=soil_type_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID VEG TYPE." +! call ESMF_FieldGet(veg_type_target_grid, & +! farrayPtr=veg_type_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID ALVSF." +! call ESMF_FieldGet(alvsf_target_grid, & +! farrayPtr=alvsf_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID ALNSF." +! call ESMF_FieldGet(alnsf_target_grid, & +! farrayPtr=alnsf_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID ALVWF." +! call ESMF_FieldGet(alvwf_target_grid, & +! farrayPtr=alvwf_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID ALNWF." +! call ESMF_FieldGet(alnwf_target_grid, & +! farrayPtr=alnwf_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID VEG FRAC." +! call ESMF_FieldGet(veg_greenness_target_grid, & +! farrayPtr=veg_greenness_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID MAX VEG FRAC." +! call ESMF_FieldGet(max_veg_greenness_target_grid, & +! farrayPtr=max_veg_greenness_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID MIN VEG FRAC." +! call ESMF_FieldGet(min_veg_greenness_target_grid, & +! farrayPtr=min_veg_greenness_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID CANOPY MC." +! call ESMF_FieldGet(canopy_mc_target_grid, & +! farrayPtr=canopy_mc_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! print*,"- CALL FieldGet FOR TARGET GRID SNOW ALBEDO." +! call ESMF_FieldGet(mxsno_albedo_target_grid, & +! farrayPtr=mxsno_albedo_target_ptr, rc=rc) +! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & +! call error_handler("IN FieldGet", rc) +! +! do i =clb(1),cub(1) +! do j = clb(2),cub(2) +! if (landmask_ptr(i,j)==1 .and. soilm_target_ptr(i,j,1) > 0.75 .and. nint(veg_type_target_ptr(i,j)) /= veg_type_landice_target) then +! !write(*,'(a,2i5,a,i3,f5.2,2i3)') "CORRECTING G.P.",i,j," FROM LAND TO SEA VALUES; & +! ! curr landmask, soilm, stype, vtype = ",landmask_ptr(i,j),& +! ! soilm_target_ptr(i,j,1),nint(soil_type_target_ptr(i,j)), & +! ! nint(veg_type_target_ptr(i,j)) +! soil_type_target_ptr(i,j) = 0.0 +! veg_type_target_ptr(i,J) = 0.0 +! landmask_ptr(i,j) = 0 +! alvsf_target_ptr(i,j) = 0.06 +! alvwf_target_ptr(i,j) = 0.06 +! alnsf_target_ptr(i,j) = 0.06 +! alnwf_target_ptr(i,j) = 0.06 +! min_veg_greenness_target_ptr(i,j) = 0.0 +! max_veg_greenness_target_ptr(i,j) = 0.0 +! veg_greenness_target_ptr(i,j) = 0.0 +! mxsno_albedo_target_ptr(i,j) = 0.0 +! canopy_mc_target_ptr(i,j) = 0.0 +! endif +! enddo +! enddo +! +!end subroutine check_smois_water + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! When using vegetation type from the input data instead of the orography file, there +!! are frequently points with ~0 soil moisture at land points. For these points, set +!! values in all relevant target grid surface arrays to fill values (done in +!! check_smois_land) then run the search routine again to fill with appropriate values +!! from nearby points (done in replace_land_sfcparams). +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine check_smois_land + +use model_grid, only : landmask_target_grid + +use static_data, only : veg_type_target_grid, veg_greenness_target_grid, & + soil_type_target_grid, max_veg_greenness_target_grid,& + min_veg_greenness_target_grid, mxsno_albedo_target_grid, & + alvsf_target_grid,alvwf_target_grid,& + alnsf_target_grid,alnwf_target_grid,& + facsf_target_grid,facwf_target_grid, & + slope_type_target_grid + + implicit none + + integer :: clb(3), cub(3), i, j, rc + + integer(esmf_kind_r8), pointer :: landmask_ptr(:,:) + + real(esmf_kind_r8), pointer :: soilm_target_ptr(:,:,:), & + soilt_target_ptr(:,:,:), & + alvsf_target_ptr(:,:), & + alnsf_target_ptr(:,:), & + alvwf_target_ptr(:,:), & + alnwf_target_ptr(:,:), & + veg_greenness_target_ptr(:,:), & + min_veg_greenness_target_ptr(:,:), & + max_veg_greenness_target_ptr(:,:), & + canopy_mc_target_ptr(:,:), & + mxsno_albedo_target_ptr(:,:), & + soil_type_target_ptr(:,:), & + veg_type_target_ptr(:,:), & + facsf_target_ptr(:,:), & + facwf_target_ptr(:,:), & + slope_type_target_ptr(:,:) + + + print*,"- CALL FieldGet FOR TARGET GRID LAND-SEA MASK." + call ESMF_FieldGet(landmask_target_grid, & + farrayPtr=landmask_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR SOIL MOIS TARGET GRID." + call ESMF_FieldGet(soilm_tot_target_grid, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=soilm_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR SOIL TEMP TARGET GRID." + call ESMF_FieldGet(soil_temp_target_grid, & + farrayPtr=soilt_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID SOIL TYPE." + call ESMF_FieldGet(soil_type_target_grid, & + farrayPtr=soil_type_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID SOIL TYPE." + call ESMF_FieldGet(slope_type_target_grid, & + farrayPtr=slope_type_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID VEG TYPE." + call ESMF_FieldGet(veg_type_target_grid, & + farrayPtr=veg_type_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID ALVSF." + call ESMF_FieldGet(alvsf_target_grid, & + farrayPtr=alvsf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID ALNSF." + call ESMF_FieldGet(alnsf_target_grid, & + farrayPtr=alnsf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID ALVWF." + call ESMF_FieldGet(alvwf_target_grid, & + farrayPtr=alvwf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID ALNWF." + call ESMF_FieldGet(alnwf_target_grid, & + farrayPtr=alnwf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + +print*,"- CALL FieldGet FOR TARGET GRID FACSF." + call ESMF_FieldGet(facsf_target_grid, & + farrayPtr=facsf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID FACWF." + call ESMF_FieldGet(facwf_target_grid, & + farrayPtr=facwf_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID VEG FRAC." + call ESMF_FieldGet(veg_greenness_target_grid, & + farrayPtr=veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID MAX VEG FRAC." + call ESMF_FieldGet(max_veg_greenness_target_grid, & + farrayPtr=max_veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID MIN VEG FRAC." + call ESMF_FieldGet(min_veg_greenness_target_grid, & + farrayPtr=min_veg_greenness_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID CANOPY MC." + call ESMF_FieldGet(canopy_mc_target_grid, & + farrayPtr=canopy_mc_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + print*,"- CALL FieldGet FOR TARGET GRID SNOW ALBEDO." + call ESMF_FieldGet(mxsno_albedo_target_grid, & + farrayPtr=mxsno_albedo_target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + do i =clb(1),cub(1) + do j = clb(2),cub(2) + if (landmask_ptr(i,j)==1 .and. soilm_target_ptr(i,j,1) < 0.001 .and. nint(veg_type_target_ptr(i,j)) /= veg_type_landice_target) then !.and. & + WRITE(*,'(a,2i5,a,2i3)'), " CORRECTING G.P. ",i,j," PARAMS FROM SEA TO LAND & + VALUES; curr stype,vtype=", nint(soil_type_target_ptr(i,j)),nint(veg_type_target_ptr(i,j)) + ! Set values to missing so that search function can then replace + ! them with nearby point values (see replace_land_sfcparams) + ! subroutine) + soilm_target_ptr(i,j,:) = -99999.9 + soilt_target_ptr(i,j,:) = -99999.9 + + soil_type_target_ptr(i,j) = -99999.9 + veg_type_target_ptr(i,J) = -99999.9 + + alvsf_target_ptr(i,j) = -99999.9 + alvwf_target_ptr(i,j) = -99999.9 + alnsf_target_ptr(i,j) = -99999.9 + alnwf_target_ptr(i,j) = -99999.9 + facsf_target_ptr(i,j) = -99999.9 + facwf_target_ptr(i,j) = -99999.9 + min_veg_greenness_target_ptr(i,j) = -99999.9 + max_veg_greenness_target_ptr(i,j) = -99999.9 + veg_greenness_target_ptr(i,j) = -99999.9 + mxsno_albedo_target_ptr(i,j) = -99999.9 + canopy_mc_target_ptr(i,j) = -99999.9 + slope_type_target_ptr(i,j) = -99999.9 + end if + enddo + enddo + !search (field, mask, idim, jdim, tile, field_num, + !call search(soilm_target_ptr(clb(1):cub(1),clb(2):cub(2), +end subroutine check_smois_land + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! When using vegetation type from the input data instead of the orography file, there +!! are frequently points with ~0 soil moisture at land points. For these points, set +!! values in all relevant target grid surface arrays to fill values (done in +!! check_smois_land) then run the search routine again to fill with appropriate values +!! from nearby points (done in replace_land_sfcparams). +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine replace_land_sfcparams(localpet) + + use search_util + use model_grid, only : lsoil_target, i_target, j_target, & + landmask_target_grid + use static_data, only : soil_type_target_grid, & + alvsf_target_grid, & + alvwf_target_grid, & + alnsf_target_grid, & + alnwf_target_grid, & + facsf_target_grid, & + facwf_target_grid, & + mxsno_albedo_target_grid, & + max_veg_greenness_target_grid, & + min_veg_greenness_target_grid, & + slope_type_target_grid, & + veg_greenness_target_grid, & + veg_type_target_grid + +!i_input, j_input, input_grid + implicit none + integer, intent(in) :: localpet + !character(len=1000) :: msg + integer :: rc, k_soil + integer(esmf_kind_i8) :: maskdata_one_tile(i_target,j_target) + real(esmf_kind_r8) :: soiltdata_one_tile_3d(i_target,j_target,lsoil_target), & + soilmdata_one_tile_3d(i_target,j_target,lsoil_target), & + soiltype_one_tile(i_target,j_target), & + alvsf_one_tile(i_target,j_target), & + alvwf_one_tile(i_target,j_target), & + alnsf_one_tile(i_target,j_target), & + alnwf_one_tile(i_target,j_target), & + facsf_one_tile(i_target,j_target), & + facwf_one_tile(i_target,j_target), & + mxsno_albedo_one_tile(i_target,j_target), & + max_veg_greenness_one_tile(i_target,j_target), & + min_veg_greenness_one_tile(i_target,j_target), & + slope_type_one_tile(i_target,j_target), & + veg_greenness_one_tile(i_target,j_target), & + veg_type_one_tile(i_target,j_target), & + canopy_mc_one_tile(i_target,j_target), & + tmp(i_target,j_target) + + !if (localpet==0) PRINT *, "STARTING SUBROUTINE replace_land_sfcparams" + call ESMF_FieldGather(landmask_target_grid, maskdata_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + + +! Now get all variables to call search routine for + call ESMF_FieldGather(soil_temp_target_grid, soiltdata_one_tile_3d, rootPet=0, tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(soilm_tot_target_grid, soilmdata_one_tile_3d, rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(soil_type_target_grid, soiltype_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(alvsf_target_grid, alvsf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(alvwf_target_grid, alvwf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(alnsf_target_grid, alnsf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(alnwf_target_grid, alnwf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(facsf_target_grid, facsf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(facwf_target_grid, facwf_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(mxsno_albedo_target_grid,mxsno_albedo_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(max_veg_greenness_target_grid, max_veg_greenness_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(min_veg_greenness_target_grid, min_veg_greenness_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(slope_type_target_grid, slope_type_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(veg_greenness_target_grid, veg_greenness_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(veg_type_target_grid, veg_type_one_tile,rootPet=0,tile=1,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + call ESMF_FieldGather(canopy_mc_target_grid,canopy_mc_one_tile,rootPet=0,tile=1,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 + ! 2d vars + call search(soiltype_one_tile,maskdata_one_tile,i_target,j_target,1,224) + call search(alvsf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(alvwf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(alnsf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(alnwf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(facsf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(facwf_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(mxsno_albedo_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(max_veg_greenness_one_tile,maskdata_one_tile,i_target,j_target,1,226) + call search(min_veg_greenness_one_tile,maskdata_one_tile,i_target,j_target,1,226) + call search(slope_type_one_tile,maskdata_one_tile,i_target,j_target,1,0) + call search(veg_greenness_one_tile,maskdata_one_tile,i_target,j_target,1,226) + call search(veg_type_one_tile,maskdata_one_tile,i_target,j_target,1,225) + call search(canopy_mc_one_tile,maskdata_one_tile,i_target,j_target,1,223) + ! 3d vars + do k_soil = 1, lsoil_target + tmp = soiltdata_one_tile_3d(:,:,k_soil) + call search(tmp,maskdata_one_tile,i_target,j_target,1,85) + soiltdata_one_tile_3d(:,:,k_soil) = tmp + + tmp = soilmdata_one_tile_3d(:,:,k_soil) + call search(tmp,maskdata_one_tile,i_target,j_target,1,86) + soilmdata_one_tile_3d(:,:,k_soil) = tmp + end do + end if ! localpet + +! scatter data back to procs + call ESMF_FieldScatter(soilm_tot_target_grid, soilmdata_one_tile_3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(soil_temp_target_grid, soiltdata_one_tile_3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(soil_type_target_grid, soiltype_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(alvsf_target_grid, alvsf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(alvwf_target_grid, alvwf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(alnsf_target_grid, alnsf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(alnwf_target_grid, alnwf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(facsf_target_grid, facsf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(facwf_target_grid, facwf_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(mxsno_albedo_target_grid, mxsno_albedo_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(max_veg_greenness_target_grid, max_veg_greenness_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(min_veg_greenness_target_grid, min_veg_greenness_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(slope_type_target_grid, slope_type_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(veg_greenness_target_grid, veg_greenness_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(veg_type_target_grid, veg_type_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + call ESMF_FieldScatter(canopy_mc_target_grid, canopy_mc_one_tile, rootpet=0,rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + +end subroutine replace_land_sfcparams + + FUNCTION FRH2O (TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS) !$$$ function documentation block ! @@ -2500,7 +3387,6 @@ subroutine rescale_soil_moisture !--------------------------------------------------------------------------------------------- if (soilt_target /= soilt_input) then - !--------------------------------------------------------------------------------------------- ! Rescale top layer. First, determine direct evaporation part: !--------------------------------------------------------------------------------------------- @@ -2653,6 +3539,122 @@ subroutine adjust_soilt_for_terrain end subroutine adjust_soilt_for_terrain +!--------------------------------------------------------------------------------------------- +! Adjust soil levels of the input grid if there's a mismatch between input and +! target grids. Presently can only convert from 9 to 4 levels. +!--------------------------------------------------------------------------------------------- + + subroutine adjust_soil_levels(localpet) + use model_grid, only : lsoil_target, i_input, j_input, input_grid + use input_data, only : lsoil_input, soil_temp_input_grid, & + soilm_liq_input_grid, soilm_tot_input_grid + implicit none + integer, intent(in) :: localpet + character(len=1000) :: msg + integer :: rc + real(esmf_kind_r8) :: tmp(i_input,j_input), & + data_one_tile(i_input,j_input,lsoil_input), & + tmp3d(i_input,j_input,lsoil_target) + if (lsoil_input == 9 .and. lsoil_target == 4) then + print*, "CONVERTING FROM 9 INPUT SOIL LEVELS TO 4 TARGET SOIL LEVELS" + call ESMF_FieldGather(soil_temp_input_grid, data_one_tile, rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + + call ESMF_FieldDestroy(soil_temp_input_grid,rc=rc) + soil_temp_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_target/), rc=rc) + + if(localpet==0)then + tmp3d(:,:,1)= (data_one_tile(:,:,1) + data_one_tile(:,:,2))/2.0 * 0.1 + & + (data_one_tile(:,:,2) + data_one_tile(:,:,3))/2.0 * 0.3 + & + (data_one_tile(:,:,3) + data_one_tile(:,:,4))/2.0 * 0.6 + tmp = (data_one_tile(:,:,6) - data_one_tile(:,:,5)) / 30.0 * 10.0 + data_one_tile(:,:,5) !Linear approx. of 40 cm obs + tmp3d(:,:,2)= (data_one_tile(:,:,4) + data_one_tile(:,:,5)) / 2.0 * 0.75 + & + (data_one_tile(:,:,5) + tmp) / 2.0 * 0.25 + tmp3d(:,:,3)= (tmp + data_one_tile(:,:,6)) /2.0 * (1.0/3.0) + & + (data_one_tile(:,:,6) + data_one_tile(:,:,7)) / 2.0 * (2.0/3.0) + tmp = (data_one_tile(:,:,9) - data_one_tile(:,:,9)) / 140.0 * 40.0 + data_one_tile(:,:,8) !Linear approx of 200 cm obs + tmp3d(:,:,4)= (data_one_tile(:,:,7) + data_one_tile(:,:,8)) / 2.0 * 0.6 + & + (data_one_tile(:,:,8) + tmp) / 2.0 * 0.4 + endif + + call ESMF_FieldScatter(soil_temp_input_grid, tmp3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + call ESMF_FieldGather(soilm_tot_input_grid, data_one_tile, rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + + call ESMF_FieldDestroy(soilm_tot_input_grid,rc=rc) + soilm_tot_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_target/), rc=rc) + + if(localpet==0) then + tmp3d(:,:,1)= (data_one_tile(:,:,1) + data_one_tile(:,:,2))/2.0 * 0.1 + & + (data_one_tile(:,:,2) + data_one_tile(:,:,3))/2.0 * 0.3 + & + (data_one_tile(:,:,3) + data_one_tile(:,:,4))/2.0 * 0.6 + tmp = (data_one_tile(:,:,6) - data_one_tile(:,:,5)) / 30.0 * 10.0 + data_one_tile(:,:,5) !Linear approx. of 40 cm obs + tmp3d(:,:,2)= (data_one_tile(:,:,4) + data_one_tile(:,:,5)) / 2.0 * 0.75 + & + (data_one_tile(:,:,5) + tmp) / 2.0 * 0.25 + tmp3d(:,:,3)= (tmp + data_one_tile(:,:,6)) /2.0 * (1.0/3.0) + & + (data_one_tile(:,:,6) + data_one_tile(:,:,7)) / 2.0 * (2.0/3.0) + tmp = (data_one_tile(:,:,9) - data_one_tile(:,:,9)) / 140.0 * 40.0 + data_one_tile(:,:,8) !Linear approx of 200 cm obs + tmp3d(:,:,4)= (data_one_tile(:,:,7) + data_one_tile(:,:,8)) / 2.0 * 0.6 + & + (data_one_tile(:,:,8) + tmp) / 2.0 * 0.4 + endif + + call ESMF_FieldScatter(soilm_tot_input_grid, tmp3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + call ESMF_FieldGather(soilm_liq_input_grid, data_one_tile, rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", rc) + + call ESMF_FieldDestroy(soilm_liq_input_grid,rc=rc) + soilm_liq_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lsoil_target/), rc=rc) + if(localpet==0) then + tmp3d(:,:,1)= (data_one_tile(:,:,1) + data_one_tile(:,:,2))/2.0 * 0.1 + & + (data_one_tile(:,:,2) + data_one_tile(:,:,3))/2.0 * 0.3 + & + (data_one_tile(:,:,3) + data_one_tile(:,:,4))/2.0 * 0.6 + tmp = (data_one_tile(:,:,6) - data_one_tile(:,:,5)) / 30.0 * 10.0 + data_one_tile(:,:,5) !Linear approx. of 40 cm obs + tmp3d(:,:,2)= (data_one_tile(:,:,4) + data_one_tile(:,:,5)) / 2.0 * 0.75 + & + (data_one_tile(:,:,5) + tmp) / 2.0 * 0.25 + tmp3d(:,:,3)= (tmp + data_one_tile(:,:,6)) /2.0 * (1.0/3.0) + & + (data_one_tile(:,:,6) + data_one_tile(:,:,7)) / 2.0 * (2.0/3.0) + tmp = (data_one_tile(:,:,9) - data_one_tile(:,:,9)) / 140.0 * 40.0 + data_one_tile(:,:,8) !Linear approx of 200 cm obs + tmp3d(:,:,4)= (data_one_tile(:,:,7) + data_one_tile(:,:,8)) / 2.0 * 0.6 + & + (data_one_tile(:,:,8) + tmp) / 2.0 * 0.4 + endif + + call ESMF_FieldScatter(soilm_liq_input_grid, tmp3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + elseif (lsoil_input /= lsoil_target) then + rc = -1 + + write(msg,'("NUMBER OF SOIL LEVELS IN INPUT (",I2,") and OUPUT & + (",I2,") MUST EITHER BE EQUAL OR 9 AND 4, RESPECTIVELY")') & + lsoil_input, lsoil_target + + call error_handler(trim(msg), rc) + endif + + end subroutine adjust_soil_levels + !--------------------------------------------------------------------------------------------- ! Set roughness at land and sea ice. !--------------------------------------------------------------------------------------------- @@ -3474,6 +4476,21 @@ subroutine create_surface_esmf_fields call error_handler("IN FieldGet", rc) target_ptr = init_val + + print*,"- CALL FieldCreate FOR TARGET GRID LEAF AREA INDEX." + lai_target_grid = ESMF_FieldCreate(target_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- INITIALIZE TARGET leaf area index." + call ESMF_FieldGet(lai_target_grid, & + farrayPtr=target_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + target_ptr = init_val print*,"- CALL FieldCreate FOR TARGET GRID Z0." z0_target_grid = ESMF_FieldCreate(target_grid, & @@ -3758,6 +4775,7 @@ subroutine cleanup_target_sfc_data call ESMF_FieldDestroy(canopy_mc_target_grid, rc=rc) call ESMF_FieldDestroy(z0_target_grid, rc=rc) call ESMF_FieldDestroy(terrain_from_input_grid, rc=rc) + call ESMF_FieldDestroy(terrain_from_input_grid_land, rc=rc) call ESMF_FieldDestroy(soil_type_from_input_grid, rc=rc) call ESMF_FieldDestroy(soil_temp_target_grid, rc=rc) call ESMF_FieldDestroy(soilm_tot_target_grid, rc=rc) diff --git a/sorc/chgres_cube.fd/write_data.F90 b/sorc/chgres_cube.fd/write_data.F90 index 0e04e7330..a05c824d8 100644 --- a/sorc/chgres_cube.fd/write_data.F90 +++ b/sorc/chgres_cube.fd/write_data.F90 @@ -1198,7 +1198,8 @@ subroutine write_fv3_atm_data_netcdf(localpet) use program_setup, only : halo=>halo_bndy, & input_type, tracers, num_tracers, & - use_thomp_mp_climo + use_thomp_mp_climo, & + regional use atmosphere, only : lev_target, & levp1_target, & @@ -1283,7 +1284,11 @@ subroutine write_fv3_atm_data_netcdf(localpet) HEADER : if (localpet < num_tiles_target_grid) then tile = localpet + 1 - WRITE(OUTFILE, '(A, I1, A)') 'out.atm.tile', tile, '.nc' + if (regional > 0) then + outfile = "out.atm.tile7.nc" + else + WRITE(OUTFILE, '(A, I1, A)') 'out.atm.tile', tile, '.nc' + endif !--- open the file error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & @@ -1530,6 +1535,7 @@ subroutine write_fv3_atm_data_netcdf(localpet) if (localpet < num_tiles_target_grid) then dum3d(:,:,:) = data_one_tile_3d(i_start:i_end,j_start:j_end,:) dum3d(:,:,1:lev_target) = dum3d(:,:,lev_target:1:-1) + print*,"MIN MAX W AT WRITE = ", minval(dum3d(:,:,:)), maxval(dum3d(:,:,:)) error = nf90_put_var( ncid, id_w, dum3d) call netcdf_err(error, 'WRITING VERTICAL VELOCITY RECORD' ) endif @@ -1680,6 +1686,7 @@ subroutine write_fv3_atm_data_netcdf(localpet) if (localpet < num_tiles_target_grid) then dum3d(:,:,:) = data_one_tile_3d(i_start:i_end,j_start:jp1_end,:) dum3d(:,:,1:lev_target) = dum3d(:,:,lev_target:1:-1) + print*,"MIN MAX US AT WRITE = ", minval(dum3d(:,:,:)), maxval(dum3d(:,:,:)) error = nf90_put_var( ncid, id_u_s, dum3d) call netcdf_err(error, 'WRITING U_S RECORD' ) endif @@ -1696,6 +1703,7 @@ subroutine write_fv3_atm_data_netcdf(localpet) if (localpet < num_tiles_target_grid) then dum3d(:,:,:) = data_one_tile_3d(i_start:i_end,j_start:jp1_end,:) dum3d(:,:,1:lev_target) = dum3d(:,:,lev_target:1:-1) + print*,"MIN MAX VS AT WRITE = ", minval(dum3d(:,:,:)), maxval(dum3d(:,:,:)) error = nf90_put_var( ncid, id_v_s, dum3d) call netcdf_err(error, 'WRITING V_S RECORD' ) endif @@ -1760,6 +1768,7 @@ subroutine write_fv3_atm_data_netcdf(localpet) if (localpet < num_tiles_target_grid) then dum3d(:,:,:) = data_one_tile_3d(i_start:ip1_end,j_start:j_end,:) dum3d(:,:,1:lev_target) = dum3d(:,:,lev_target:1:-1) + print*,"MIN MAX UW AT WRITE = ", minval(dum3d(:,:,:)), maxval(dum3d(:,:,:)) error = nf90_put_var( ncid, id_u_w, dum3d) call netcdf_err(error, 'WRITING U_W RECORD' ) endif @@ -1776,6 +1785,7 @@ subroutine write_fv3_atm_data_netcdf(localpet) if (localpet < num_tiles_target_grid) then dum3d(:,:,:) = data_one_tile_3d(i_start:ip1_end,j_start:j_end,:) dum3d(:,:,1:lev_target) = dum3d(:,:,lev_target:1:-1) + print*,"MIN MAX VW AT WRITE = ", minval(dum3d(:,:,:)), maxval(dum3d(:,:,:)) error = nf90_put_var( ncid, id_v_w, dum3d) call netcdf_err(error, 'WRITING V_W RECORD' ) endif @@ -1804,7 +1814,8 @@ subroutine write_fv3_sfc_data_netcdf(localpet) longitude_target_grid, & i_target, j_target, lsoil_target - use program_setup, only : convert_nst, halo=>halo_bndy + use program_setup, only : convert_nst, halo=>halo_bndy, & + regional, lai_from_climo use surface, only : canopy_mc_target_grid, & f10m_target_grid, & @@ -1824,6 +1835,7 @@ subroutine write_fv3_sfc_data_netcdf(localpet) tprcp_target_grid, & ustar_target_grid, & z0_target_grid, & + lai_target_grid, & c_d_target_grid, & c_0_target_grid, & d_conv_target_grid, & @@ -1880,6 +1892,7 @@ subroutine write_fv3_sfc_data_netcdf(localpet) integer :: id_fice, id_tisfc, id_tprcp integer :: id_srflag, id_snwdph, id_shdmin integer :: id_shdmax, id_slope, id_snoalb + integer :: id_lai integer :: id_stc, id_smc, id_slc integer :: id_tref, id_z_c, id_c_0 integer :: id_c_d, id_w_0, id_w_d @@ -1947,7 +1960,11 @@ subroutine write_fv3_sfc_data_netcdf(localpet) LOCAL_PET : if (localpet == 0) then - WRITE(OUTFILE, '(A, I1, A)') 'out.sfc.tile', tile, '.nc' + if (regional > 0) then + outfile = "out.sfc.tile7.nc" + else + WRITE(OUTFILE, '(A, I1, A)') 'out.sfc.tile', tile, '.nc' + endif !--- open the file error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & @@ -2293,6 +2310,17 @@ subroutine write_fv3_sfc_data_netcdf(localpet) call netcdf_err(error, 'DEFINING SNOALB UNITS' ) error = nf90_put_att(ncid, id_snoalb, "coordinates", "geolon geolat") call netcdf_err(error, 'DEFINING SNOALB COORD' ) + + if (.not. lai_from_climo) then + error = nf90_def_var(ncid, 'lai', NF90_DOUBLE, (/dim_x,dim_y,dim_time/), id_lai) + call netcdf_err(error, 'DEFINING LAI' ) + error = nf90_put_att(ncid, id_lai, "long_name", "lai") + call netcdf_err(error, 'DEFINING LAI LONG NAME' ) + error = nf90_put_att(ncid, id_lai, "units", "none") + call netcdf_err(error, 'DEFINING LAI UNITS' ) + error = nf90_put_att(ncid, id_lai, "coordinates", "geolon geolat") + call netcdf_err(error, 'DEFINING LAI COORD' ) + endif error = nf90_def_var(ncid, 'stc', NF90_DOUBLE, (/dim_x,dim_y,dim_lsoil,dim_time/), id_stc) call netcdf_err(error, 'DEFINING STC' ) @@ -2580,7 +2608,20 @@ subroutine write_fv3_sfc_data_netcdf(localpet) error = nf90_put_var( ncid, id_snoalb, dum2d, start=(/1,1,1/), count=(/i_target_out,j_target_out,1/)) call netcdf_err(error, 'WRITING MAX SNOW ALBEDO RECORD' ) endif + + if (.not. lai_from_climo) then + print*,"- CALL FieldGather FOR TARGET GRID LEAF AREA INDEX FOR TILE: ", tile + call ESMF_FieldGather(lai_target_grid, data_one_tile, rootPet=0, tile=tile, rc=error) + if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGather", error) + if (localpet == 0) then + dum2d(:,:) = data_one_tile(istart:iend, jstart:jend) + error = nf90_put_var( ncid, id_lai, dum2d, start=(/1,1,1/), count=(/i_target_out,j_target_out,1/)) + call netcdf_err(error, 'WRITING LEAF AREA INDEX RECORD' ) + endif + endif + print*,"- CALL FieldGather FOR TARGET GRID SOIL TYPE FOR TILE: ", tile call ESMF_FieldGather(soil_type_target_grid, data_one_tile, rootPet=0, tile=tile, rc=error) if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & diff --git a/sorc/emcsfc_snow2mdl.fd/.gitignore b/sorc/emcsfc_snow2mdl.fd/.gitignore deleted file mode 100644 index 80df0fa54..000000000 --- a/sorc/emcsfc_snow2mdl.fd/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.o -*.mod -emcsfc_snow2mdl diff --git a/sorc/global_chgres.fd/.gitignore b/sorc/global_chgres.fd/.gitignore deleted file mode 100644 index 013de5d46..000000000 --- a/sorc/global_chgres.fd/.gitignore +++ /dev/null @@ -1 +0,0 @@ -global_chgres diff --git a/sorc/machine-setup.sh b/sorc/machine-setup.sh index 08e500537..6664446fe 100644 --- a/sorc/machine-setup.sh +++ b/sorc/machine-setup.sh @@ -78,12 +78,12 @@ elif [[ -L /usrx && "$( readlink /usrx 2> /dev/null )" =~ dell ]] ; then target=wcoss_dell_p3 module purge elif [[ -d /glade ]] ; then - # We are on NCAR Yellowstone + # We are on NCAR Cheyenne if ( ! eval module help > /dev/null 2>&1 ) ; then echo load the module command 1>&2 - . /usr/share/Modules/init/$__ms_shell + . /glade/u/apps/ch/opt/lmod/8.1.7/lmod/8.1.7/init/sh fi - target=yellowstone + target=cheyenne module purge elif [[ -d /lustre && -d /ncrc ]] ; then # We are on GAEA. @@ -103,6 +103,9 @@ elif [[ "$(hostname)" =~ "Orion" ]]; then module purge elif [[ "$(hostname)" =~ "odin" ]]; then target="odin" +elif [[ -d /work/00315 && -d /scratch/00315 ]] ; then + target=stampede + module purge else echo WARNING: UNKNOWN PLATFORM 1>&2 fi diff --git a/sorc/nst_tf_chg.fd/.gitignore b/sorc/nst_tf_chg.fd/.gitignore deleted file mode 100644 index 3eff26396..000000000 --- a/sorc/nst_tf_chg.fd/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.o -*.mod -nst_tf_chg.x diff --git a/sorc/sfc_climo_gen.fd/.gitignore b/sorc/sfc_climo_gen.fd/.gitignore deleted file mode 100644 index 3c4362860..000000000 --- a/sorc/sfc_climo_gen.fd/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.o -*.mod -gridgen_sfc diff --git a/util/gdas_init/driver.cray.sh b/util/gdas_init/driver.cray.sh index 93cae342b..28d92a63a 100755 --- a/util/gdas_init/driver.cray.sh +++ b/util/gdas_init/driver.cray.sh @@ -9,7 +9,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list PROJECT_CODE=GFS-DEV diff --git a/util/gdas_init/driver.dell.sh b/util/gdas_init/driver.dell.sh index ec18e09a9..f6a0ccb49 100755 --- a/util/gdas_init/driver.dell.sh +++ b/util/gdas_init/driver.dell.sh @@ -9,7 +9,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list PROJECT_CODE=GFS-DEV diff --git a/util/gdas_init/driver.hera.sh b/util/gdas_init/driver.hera.sh index 3fee1147e..9bf57791b 100755 --- a/util/gdas_init/driver.hera.sh +++ b/util/gdas_init/driver.hera.sh @@ -9,7 +9,9 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target +module list # Needed for NDATE utility module use -a /scratch2/NCEPDEV/nwprod/NCEPLIBS/modulefiles diff --git a/util/vcoord_gen/run.cray.sh b/util/vcoord_gen/run.cray.sh index b6f9c7e13..451fc7c3c 100755 --- a/util/vcoord_gen/run.cray.sh +++ b/util/vcoord_gen/run.cray.sh @@ -21,7 +21,8 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 -source ../../modulefiles/build.$target +module use ../../modulefiles +module load build.$target module list outfile="./global_hyblev.txt" From bf70549328bc3151955eeab62f09cb922081478f Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Tue, 20 Oct 2020 10:23:21 -0400 Subject: [PATCH 2/2] Remove landsfcutil library dependency from snow2mdl program Incorporate landsfcutil library routines 'uninterpred' and 'intlon' into the snow2mdl program. The goal is to eventually retire the landsfcutil library. For details, see issue #169. --- sorc/emcsfc_snow2mdl.fd/CMakeLists.txt | 1 - sorc/emcsfc_snow2mdl.fd/snow2mdl.f | 72 +++++++++++++++++++++++++- 2 files changed, 70 insertions(+), 3 deletions(-) diff --git a/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt b/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt index e884e1d7a..c3cba21da 100644 --- a/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt +++ b/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt @@ -19,7 +19,6 @@ target_link_libraries( g2::g2_d ip::ip_d sp::sp_d - landsfcutil::landsfcutil_d bacio::bacio_4 w3nco::w3nco_d) if(OpenMP_Fortran_FOUND) diff --git a/sorc/emcsfc_snow2mdl.fd/snow2mdl.f b/sorc/emcsfc_snow2mdl.fd/snow2mdl.f index 4102dd045..200fd45a8 100755 --- a/sorc/emcsfc_snow2mdl.fd/snow2mdl.f +++ b/sorc/emcsfc_snow2mdl.fd/snow2mdl.f @@ -87,8 +87,6 @@ module snow2mdl kgds_autosnow, & bad_afwa_nh, bad_afwa_sh - use read_write_utils, only : uninterpred - private real, allocatable :: snow_cvr_mdl(:,:) ! cover in % on mdl grid @@ -1145,4 +1143,74 @@ subroutine write_grib1 end subroutine write_grib1 +!----------------------------------------------------------------------- +! fills out full grid using thinned grid data. use an iord of +! "1" to use a nearest neighbor approach. +!----------------------------------------------------------------------- + + subroutine uninterpred(iord,kmsk,fi,f,lonl,latd,len,lonsperlat) + + implicit none + + integer, intent(in) :: len + integer, intent(in) :: iord + integer, intent(in) :: lonl + integer, intent(in) :: latd + integer, intent(in) :: lonsperlat(latd/2) + integer, intent(in) :: kmsk(lonl*latd) + integer :: j,lons,jj,latd2,ii,i + + real, intent(in) :: fi(len) + real, intent(out) :: f(lonl,latd) + + latd2 = latd / 2 + ii = 1 + + do j=1,latd + + jj = j + if (j .gt. latd2) jj = latd - j + 1 + lons=lonsperlat(jj) + + if(lons.ne.lonl) then + call intlon(iord,1,1,lons,lonl,kmsk(ii),fi(ii),f(1,j)) + else + do i=1,lonl + f(i,j) = fi(ii+i-1) + enddo + endif + + ii = ii + lons + + enddo + + end subroutine uninterpred + + subroutine intlon(iord,imon,imsk,m1,m2,k1,f1,f2) + + implicit none + + integer,intent(in) :: iord,imon,imsk,m1,m2 + integer,intent(in) :: k1(m1) + integer :: i2,in,il,ir + + real,intent(in) :: f1(m1) + real,intent(out) :: f2(m2) + real :: r,x1 + + r=real(m1)/real(m2) + do i2=1,m2 + x1=(i2-1)*r + il=int(x1)+1 + ir=mod(il,m1)+1 + if(iord.eq.2.and.(imsk.eq.0.or.k1(il).eq.k1(ir))) then + f2(i2)=f1(il)*(il-x1)+f1(ir)*(x1-il+1) + else + in=mod(nint(x1),m1)+1 + f2(i2)=f1(in) + endif + enddo + + end subroutine intlon + end module snow2mdl