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/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 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"