From 15bab55477b40503ec041a3b3a1445b8da1c9537 Mon Sep 17 00:00:00 2001 From: mdtoy Date: Fri, 9 Apr 2021 15:13:55 -0600 Subject: [PATCH] GSL orographic drag parameter computation (#273) Add new program to compute parameters needed for the GSL subgrid-scale orographic drag suite. Update the grid generation scripts to make the creation of these parameters an option. Add a regression test for this new option. Fixes #253 --- driver_scripts/driver_grid.cray.sh | 2 + driver_scripts/driver_grid.dell.sh | 2 + driver_scripts/driver_grid.hera.sh | 8 +- driver_scripts/driver_grid.jet.sh | 2 + driver_scripts/driver_grid.orion.sh | 2 + reg_tests/grid_gen/driver.cray.sh | 9 +- reg_tests/grid_gen/driver.dell.sh | 9 +- reg_tests/grid_gen/driver.hera.sh | 9 +- reg_tests/grid_gen/driver.jet.sh | 9 +- reg_tests/grid_gen/driver.orion.sh | 9 +- reg_tests/grid_gen/regional.gsl.gwd.sh | 69 + sorc/orog_mask_tools.fd/CMakeLists.txt | 1 + .../orog_gsl.fd/CMakeLists.txt | 21 + sorc/orog_mask_tools.fd/orog_gsl.fd/README | 18 + .../orog_gsl.fd/gsl_oro_data.f90 | 86 ++ .../module_gsl_oro_data_lg_scale.f90 | 1344 +++++++++++++++++ .../module_gsl_oro_data_sm_scale.f90 | 1270 ++++++++++++++++ ush/fv3gfs_driver_grid.sh | 76 +- ush/fv3gfs_make_orog_gsl.sh | 70 + 19 files changed, 3004 insertions(+), 12 deletions(-) create mode 100755 reg_tests/grid_gen/regional.gsl.gwd.sh create mode 100644 sorc/orog_mask_tools.fd/orog_gsl.fd/CMakeLists.txt create mode 100644 sorc/orog_mask_tools.fd/orog_gsl.fd/README create mode 100644 sorc/orog_mask_tools.fd/orog_gsl.fd/gsl_oro_data.f90 create mode 100644 sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_lg_scale.f90 create mode 100644 sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_sm_scale.f90 create mode 100755 ush/fv3gfs_make_orog_gsl.sh diff --git a/driver_scripts/driver_grid.cray.sh b/driver_scripts/driver_grid.cray.sh index 1631bba73..6e39f266a 100755 --- a/driver_scripts/driver_grid.cray.sh +++ b/driver_scripts/driver_grid.cray.sh @@ -68,6 +68,8 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg' +export make_gsl_orog=false # 'true' if user needs 'oro' files for GSL + # orographic drag suite if [ $gtype = uniform ]; then export res=96 export add_lake=false # Add lake frac and depth to orography data. diff --git a/driver_scripts/driver_grid.dell.sh b/driver_scripts/driver_grid.dell.sh index 308a9ec06..3f04da922 100755 --- a/driver_scripts/driver_grid.dell.sh +++ b/driver_scripts/driver_grid.dell.sh @@ -70,6 +70,8 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg' +export make_gsl_orog=false # 'true' if user needs 'oro' files for GSL + # orographic drag suite if [ $gtype = uniform ]; then export res=96 export add_lake=false # Add lake frac and depth to orography data. diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh index 529f14b07..9608b2220 100755 --- a/driver_scripts/driver_grid.hera.sh +++ b/driver_scripts/driver_grid.hera.sh @@ -1,7 +1,7 @@ #!/bin/bash #SBATCH -J fv3_grid_driver -#SBATCH -A fv3-cpu +#SBATCH -A wrfruc #SBATCH --open-mode=truncate #SBATCH -o log.fv3_grid_driver #SBATCH -e log.fv3_grid_driver @@ -71,6 +71,8 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg' +export make_gsl_orog=false # 'true' if user needs 'oro' files for GSL + # orographic drag suite if [ $gtype = uniform ]; then export res=96 export add_lake=false # Add lake frac and depth to orography data. @@ -92,7 +94,7 @@ elif [ $gtype = nest ] || [ $gtype = regional_gfdl ]; then export jstart_nest=331 # Starting j-direction index of nest grid in parent tile supergrid export iend_nest=1402 # Ending i-direction index of nest grid in parent tile supergrid export jend_nest=1194 # Ending j-direction index of nest grid in parent tile supergrid - export halo=3 # Lateral boundary halo + export halo=4 # Lateral boundary halo elif [ $gtype = regional_esg ] ; then export res=-999 # equivalent resolution is computed export target_lon=-97.5 # Center longitude of grid @@ -105,7 +107,7 @@ elif [ $gtype = regional_esg ] ; then # direction is related to delx as follows: # distance = 2*delx*(circumf_Earth/360 deg) export dely=0.0585 # Grid spacing (in degrees) in the 'j' direction. - export halo=3 # number of row/cols for halo + export halo=4 # number of row/cols for halo fi #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh index a250b0cd6..853c15414 100755 --- a/driver_scripts/driver_grid.jet.sh +++ b/driver_scripts/driver_grid.jet.sh @@ -71,6 +71,8 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg' +export make_gsl_orog=false # 'true' if user needs 'oro' files for GSL + # orographic drag suite if [ $gtype = uniform ]; then export res=96 export add_lake=false # Add lake frac and depth to orography data. diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh index f4ac50a39..a4d9812ef 100755 --- a/driver_scripts/driver_grid.orion.sh +++ b/driver_scripts/driver_grid.orion.sh @@ -70,6 +70,8 @@ module list export gtype=uniform # 'uniform', 'stretch', 'nest', # 'regional_gfdl', 'regional_esg' +export make_gsl_orog=false # 'true' if user needs 'oro' files for GSL + # orographic drag suite if [ $gtype = uniform ]; then export res=96 export add_lake=false # Add lake frac and depth to orography data. diff --git a/reg_tests/grid_gen/driver.cray.sh b/reg_tests/grid_gen/driver.cray.sh index d19500ac4..f3d55b731 100755 --- a/reg_tests/grid_gen/driver.cray.sh +++ b/reg_tests/grid_gen/driver.cray.sh @@ -73,10 +73,17 @@ bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J gfdl.regional -W 0: bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J esg.regional -W 0:10 -M 2400 \ -w 'ended(gfdl.regional)' -extsched 'CRAYLINUX[]' "export NODES=1; $PWD/esg.regional.sh" +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag. +#----------------------------------------------------------------------------- + +bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J reg.gsl.gwd -W 0:08 -M 2400 \ + -w 'ended(esg.regional)' -extsched 'CRAYLINUX[]' "export NODES=1; $PWD/regional.gsl.gwd.sh" + #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- -bsub -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J summary -R "rusage[mem=100]" -W 0:01 -w 'ended(esg.regional)' "grep -a '<<<' $LOG_FILE >> $SUM_FILE" +bsub -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J summary -R "rusage[mem=100]" -W 0:01 -w 'ended(reg.gsl.gwd)' "grep -a '<<<' $LOG_FILE >> $SUM_FILE" exit diff --git a/reg_tests/grid_gen/driver.dell.sh b/reg_tests/grid_gen/driver.dell.sh index 45880c8a7..3495e9e01 100755 --- a/reg_tests/grid_gen/driver.dell.sh +++ b/reg_tests/grid_gen/driver.dell.sh @@ -70,9 +70,16 @@ bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J gfdl.regional -W 0: bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J esg.regional -W 0:10 -x -n 24 -w 'ended(gfdl.regional)' \ -R "span[ptile=24]" -R "affinity[core(1):distribute=balance]" "$PWD/esg.regional.sh" +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag. +#----------------------------------------------------------------------------- + +bsub -e $LOG_FILE -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J reg.gsl.gwd -W 0:08 -x -n 24 -w 'ended(esg.regional)' \ + -R "span[ptile=24]" -R "affinity[core(1):distribute=balance]" "$PWD/regional.gsl.gwd.sh" + #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- bsub -o $LOG_FILE -q $QUEUE -P $PROJECT_CODE -J summary -R "affinity[core(1)]" -R "rusage[mem=100]" -W 0:01 \ - -w 'ended(esg.regional)' "grep -a '<<<' $LOG_FILE >> $SUM_FILE" + -w 'ended(reg.gsl.gwd)' "grep -a '<<<' $LOG_FILE >> $SUM_FILE" diff --git a/reg_tests/grid_gen/driver.hera.sh b/reg_tests/grid_gen/driver.hera.sh index 1573a6c5a..31a17c0fb 100755 --- a/reg_tests/grid_gen/driver.hera.sh +++ b/reg_tests/grid_gen/driver.hera.sh @@ -75,12 +75,19 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional \ -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST2 ./esg.regional.sh) +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag test. +#----------------------------------------------------------------------------- + +TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST3 ./regional.gsl.gwd.sh) + #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - --open-mode=append -q $QUEUE -d afterok:$TEST3 << EOF + --open-mode=append -q $QUEUE -d afterok:$TEST4 << EOF #!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.jet.sh b/reg_tests/grid_gen/driver.jet.sh index 6cd9c53d1..ed2cb3492 100755 --- a/reg_tests/grid_gen/driver.jet.sh +++ b/reg_tests/grid_gen/driver.jet.sh @@ -73,12 +73,19 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional \ --partition=xjet -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST2 ./esg.regional.sh) +#----------------------------------------------------------------------------- +# Regional GSL gravity wave drag. +#----------------------------------------------------------------------------- + +TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + --partition=xjet -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST3 ./regional.gsl.gwd.sh) + #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --partition=xjet --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - --open-mode=append -q $QUEUE -d afterok:$TEST3 << EOF + --open-mode=append -q $QUEUE -d afterok:$TEST4 << EOF #!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.orion.sh b/reg_tests/grid_gen/driver.orion.sh index 01c812af4..a18bac7ff 100755 --- a/reg_tests/grid_gen/driver.orion.sh +++ b/reg_tests/grid_gen/driver.orion.sh @@ -69,12 +69,19 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ TEST3=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J esg.regional \ --open-mode=append -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST2 ./esg.regional.sh) +#----------------------------------------------------------------------------- +# Regional grid with GSL gravity wave drag fields. +#----------------------------------------------------------------------------- + +TEST4=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J reg.gsl.gwd \ + --open-mode=append -o $LOG_FILE -e $LOG_FILE -d afterok:$TEST3 ./regional.gsl.gwd.sh) + #----------------------------------------------------------------------------- # Create summary log. #----------------------------------------------------------------------------- sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ - --open-mode=append -q $QUEUE -d afterok:$TEST3 << EOF + --open-mode=append -q $QUEUE -d afterok:$TEST4 << EOF #!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/regional.gsl.gwd.sh b/reg_tests/grid_gen/regional.gsl.gwd.sh new file mode 100755 index 000000000..0a66daeeb --- /dev/null +++ b/reg_tests/grid_gen/regional.gsl.gwd.sh @@ -0,0 +1,69 @@ +#!/bin/bash + +#----------------------------------------------------------------------- +# Create a regional esg grid with GSL gravity wave drag files. +# Compare output to a set of baseline files using the 'nccmp' utility. +# This script is run by the machine specific driver script. +#----------------------------------------------------------------------- + +set -x + +export TEMP_DIR=${WORK_DIR}/regional.gsl.gwd.work +export out_dir=${WORK_DIR}/regional.gsl.gwd + +export gtype=regional_esg +export make_gsl_orog=true # Create GSL gravity wave drag fields +export target_lon=-97.5 # Center longitude of the highest resolution tile +export target_lat=35.5 # Center latitude of the highest resolution tile +export idim=301 # Dimension of grid in 'i' direction +export jdim=200 # Dimension of grid in 'j' direction +export delx=0.0585 # Grid spacing in degrees in 'i' direction +export dely=0.0585 # Grid spacing in degrees in 'j' direction +export halo=4 + +NCCMP=${NCCMP:-$(which nccmp)} + +#----------------------------------------------------------------------- +# Start script. +#----------------------------------------------------------------------- + +echo "Starting at: " `date` + +$home_dir/ush/fv3gfs_driver_grid.sh + +iret=$? +if [ $iret -ne 0 ]; then + set +x + echo "<<< REGIONAL GSL GWD TEST FAILED. <<<" + exit $iret +fi + +echo "Ending at: " `date` + +#----------------------------------------------------------------------------- +# Compare output to baseline set of data. +#----------------------------------------------------------------------------- + +cd $out_dir/C772 + +test_failed=0 +for files in *tile*.nc ./fix_sfc/*tile*.nc +do + if [ -f $files ]; then + echo CHECK $files + $NCCMP -dmfqS $files $HOMEreg/regional.gsl.gwd/$files + iret=$? + if [ $iret -ne 0 ]; then + test_failed=1 + fi + fi +done + +set +x +if [ $test_failed -ne 0 ]; then + echo "<<< REGIONAL GSL GWD TEST FAILED. >>>" +else + echo "<<< REGIONAL GSL GWD TEST PASSED. >>>" +fi + +exit 0 diff --git a/sorc/orog_mask_tools.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/CMakeLists.txt index 39d57220d..59f9635dd 100644 --- a/sorc/orog_mask_tools.fd/CMakeLists.txt +++ b/sorc/orog_mask_tools.fd/CMakeLists.txt @@ -4,6 +4,7 @@ # George Gayno, Mark Potts add_subdirectory(orog.fd) +add_subdirectory(orog_gsl.fd) add_subdirectory(lake.fd) add_subdirectory(inland.fd) diff --git a/sorc/orog_mask_tools.fd/orog_gsl.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/orog_gsl.fd/CMakeLists.txt new file mode 100644 index 000000000..940c86c09 --- /dev/null +++ b/sorc/orog_mask_tools.fd/orog_gsl.fd/CMakeLists.txt @@ -0,0 +1,21 @@ +set(fortran_src + module_gsl_oro_data_lg_scale.f90 + module_gsl_oro_data_sm_scale.f90 + gsl_oro_data.f90) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -convert big_endian -assume byterecl") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8 -fconvert=big-endian -fno-range-check") + if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch -fallow-invalid-boz") + endif() +endif() + +set(exe_name orog_gsl) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries( + ${exe_name} + NetCDF::NetCDF_Fortran) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/orog_mask_tools.fd/orog_gsl.fd/README b/sorc/orog_mask_tools.fd/orog_gsl.fd/README new file mode 100644 index 000000000..6b8e0397c --- /dev/null +++ b/sorc/orog_mask_tools.fd/orog_gsl.fd/README @@ -0,0 +1,18 @@ +This code creates orographic statistics fields required for the orographic drag suite developed by NOAA's Global Systems Laboratory (GSL). The fields are a subset of the ones calculated by "orog.fd" except that they are calculated in a different manner. The fields are: +stddev -- standard deviation of subgrid topography +convexity -- convexity of subgrid topography +oa1,oa2,oa3,oa4 -- orographic asymmetry of subgrid topography (for 4 orientations: + 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly) +ol1,ol2,ol3,ol4 -- orographic effective length of subgrid topography (for 4 orientations: + 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly) + +Two Cxxx_oro_data files are created: + - an "ls" (large-scale) file for the gravity wave drag and blocking schemes of Kim and Doyle (QJRMS, 2005) + - a "ss" (small-scale) file for the small-scale gravity wave drag scheme of Tsiringakis et al. (QJRMS, 2017) and the turbulent orographic form drag (TOFD) scheme of Beljaars et al. (QJRMS, 2004) + +The inputs to be passed to the executable are tile number (1-6 for global, 7 for stand-alone regional) and grid resolution, e.g., 768 for C768. + +The source data are the following two files to be located in the "fix" directory: + - geo_em.d01.lat-lon.2.5m.HGT_M.nc -- global topographic data on 2.5-minute lat-lon grid (interpolated from GMTED2010 30-second topographic data) + - HGT.Beljaars_filtered.lat-lon.30s_res.nc -- global topographic data on 30-second lat-lon grid (GMTED2010 data smoothed according to Beljaars et al. (QJRMS, 2004)) + diff --git a/sorc/orog_mask_tools.fd/orog_gsl.fd/gsl_oro_data.f90 b/sorc/orog_mask_tools.fd/orog_gsl.fd/gsl_oro_data.f90 new file mode 100644 index 000000000..e629dd2de --- /dev/null +++ b/sorc/orog_mask_tools.fd/orog_gsl.fd/gsl_oro_data.f90 @@ -0,0 +1,86 @@ +!> @file +!! @brief Create orographic (oro_data) files for use by GSL drag suite +!! @author Michael Toy, NOAA/GSL +!! @date 2021-03-12 +!! +!! Program GSL_ORO_DATA +!! +!! This program calls subroutines which calculate the parameters +!! required for the GSL subgrid-scale orographic gravity-wave drag (GWDO) +!! suite on the FV3 grid. These parameters are for the small-scale +!! GWD (Tsiringakis et al., 2017) and turbulent orographic form drag (TOFD) +!! (Beljaars et al., 2004) schemes of the GSL drag suite. +!! The output fields are: +!! - stddev standard deviation of subgrid-scale topograpy +!! - convexity convexity (kurtosis) of subgrid-scale topography +!! - ol{1,2,3,4} orographic effective lengths of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! - oa{1,2,3,4} orographic asymmetries of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! +!! Note: This program works for both the global FV3GFS cubed +!! sphere, i.e., for tiles 1 through 6, (and 7 if nested +!! grid) (halo.eq.-999 for no halo), and for the stand-alone +!! regional lam (tile 7 and halo.ne.-999) +!! If a halo number is given, this is only to specify the +!! Cxxx_grid.halox data used for input. The oro_data files +!! are always "halo0" output. +!! +!! Based on code by Michael Duda provided by NCAR/MMM + +!> Brief description of program: Creates orographic (oro_data) files +!! needed by the GSL drag suite physics parameterization +!! +!! @author Michaei Toy, NOAA/GSL +!! @return 0 for success, error code otherwise. +program gsl_oro_data + +use gsl_oro_data_sm_scale, only: calc_gsl_oro_data_sm_scale +use gsl_oro_data_lg_scale, only: calc_gsl_oro_data_lg_scale + +implicit none + + +character(len=2) :: tile_num ! tile number entered by user +character(len=7) :: res_indx ! grid-resolution index, e.g., 96, 192, 384, 768, + ! etc. entered by user +character(len=4) :: halo ! halo value entered by user (for input grid data) + +logical :: duplicate_oro_data_file ! flag for whether oro_data_ls file is a duplicate + ! of oro_data_ss due to minimum grid size being less than 7.5km + + + +! Read in FV3GFS grid info +print * +print *, "Enter tile number:" +read (5,*) tile_num +print * +print *, "Enter grid-resolution index:" +read (5,*) res_indx +print * +print *, "Enter halo number (-999 for no halo):" +read (5,*) halo +print * +print *, "Creating tile oro_data for tile number: ", tile_num +print *, "Grid resolution = ", res_indx +print *, "Halo = ", halo +print * + + +call calc_gsl_oro_data_sm_scale(tile_num,res_indx,halo,duplicate_oro_data_file) + +print *, "duplicate_oro_data_file =", duplicate_oro_data_file +print * + +if ( .not.duplicate_oro_data_file ) then + call calc_gsl_oro_data_lg_scale(tile_num,res_indx,halo) +end if + + +print * +print *, "End program gsl_oro_data" +print * + + +end program gsl_oro_data diff --git a/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_lg_scale.f90 b/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_lg_scale.f90 new file mode 100644 index 000000000..eb3c2a56d --- /dev/null +++ b/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_lg_scale.f90 @@ -0,0 +1,1344 @@ +!> @file +!! @brief Calculates large-scale GWD orographic stats for FV3GFS grids +!! @author Michael Toy, NOAA/GSL +!! @date 2021-03-12 +!! +!! This module calculates the parameters required for the subgrid- +!! scale orographic gravity-wave drag (GWDO) scheme on the FV3 +!! grid. These parameters are for the large-scale GWDO and blocking +!! schemes of the GSL drag suite. 2.5minute (~5km) global topography +!! is used. The topographic data comes from the 'fix' file +!! geo_em.d01.lat-lon.2.5m.HGT_M.nc. +!! The output fields are: +!! - stddev standard deviation of subgrid-scale topograpy +!! - convexity convexity (kurtosis) of subgrid-scale topography +!! - ol{1,2,3,4} orographic effective lengths of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! - oa{1,2,3,4} orographic asymmetries of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! +!! Based on code by Michael Duda provided by NCAR/MMM +!! +module gsl_oro_data_lg_scale + +implicit none + +integer, parameter :: real_kind = selected_real_kind(6) !< single precision +integer, parameter :: dbl_kind = selected_real_kind(13) !< double precision + +real, parameter :: pi = 3.1415926535897_real_kind !< pi +integer :: dimX_fine !< x-dimension of fine grid +integer :: dimY_fine !< y-dimension of fine grid + +real (kind = real_kind), allocatable :: lat1d_fine(:) !< latitude of fine grid pts +real (kind = real_kind), allocatable :: lon1d_fine(:) !< longitude of fine grid pts + +real (kind = real_kind), parameter :: p5 = 0.5_real_kind !< one half + +real (kind = real_kind), allocatable :: HGT_M_fine(:,:) !< Height of fine grid pts (m) +real (kind = real_kind), parameter :: HGT_missing = 1.E+10 !< Flag for missing data + +contains + +!> Subroutine to compute orographic statistics needed for large-scale +!! orograhic drag (gravity wave and blocking) schemes +!! +!! @param[in] tile_num (tile number) +!! @param[in] res_indx (resolution) +!! @param[in] halo (halo number) +!! @author Michael Toy, NOAA/GSL +subroutine calc_gsl_oro_data_lg_scale(tile_num,res_indx,halo) + +use netcdf + +implicit none + +character(len=2), intent(in) :: tile_num ! tile number +character(len=7), intent(in) :: res_indx ! grid-resolution index, e.g., 96, 192, 384, etc +character(len=4), intent(in) :: halo ! halo value (for input grid data) + +integer :: i,j,ii,jj + +integer :: ncid_in,ncid_out,err +integer :: varid +integer :: dimid,latid,lonid +integer, dimension (2) :: dimids + +integer :: nfinepoints ! number of fine grid points in each coarse grid cell + +real (kind = real_kind) :: sum2, sum4, var + + +real (kind = real_kind), allocatable :: & + zs(:,:) + +logical :: zs_accum + +real (kind = real_kind) :: zs_mean +real (kind = real_kind), allocatable :: & + std_dev(:),convexity(:), & + OA1(:),OA2(:),OA3(:),OA4(:), & + OL1(:),OL2(:),OL3(:),OL4(:) + +real (kind = real_kind), parameter :: max_convexity = 10._real_kind ! max value for convexity + +integer :: nu, nd, nw, nt +real (kind = real_kind) :: ratio + + +real, parameter :: ae = 6371220._real_kind ! Earth radius in meters + +character(len=35) :: FV3_grid_input_file_name +character(len=150) :: fine_topo_source_file_name +character(len=35) :: oro_data_output_file_name + +integer :: temp_int, dimX_FV3, dimY_FV3 +real (kind = dbl_kind), allocatable :: lat_FV3_raw(:,:), lon_FV3_raw(:,:) +real (kind = real_kind), allocatable :: lat_FV3(:,:), lon_FV3(:,:) +real (kind = real_kind), allocatable :: lat_FV3_deg(:,:) ! saved version of latitude for output +real (kind = real_kind), allocatable :: lon_FV3_deg(:,:) ! saved version of longitude for output +real (kind = dbl_kind), allocatable :: area_FV3_qtr(:,:) ! meters squared +real (kind = real_kind), allocatable :: area_FV3(:,:) ! meters squared +real (kind = real_kind) :: dlta_lat, dlta_lon + +integer :: i_blk, j_blk +integer :: ii_loc, jj_loc, ii_m, jj_m +integer, dimension(3) :: s_ii, e_ii, s_jj, e_jj +real (kind = real_kind), dimension(3) :: lat_blk, lon_blk +real (kind = real_kind), dimension(3,3) :: HGT_M_coarse +real (kind = real_kind), allocatable :: HGT_M_coarse_on_fine(:,:) +integer :: cell_count ! allows for use of 1D arrays for GWD statistics fields +integer :: halo_int ! integer form of halo + +logical :: fexist + +logical, parameter :: detrend_topography = .true. ! switch for using detrended + ! or actual fine-grid topography + ! to represent subgrid terrain + + +print *, "Creating oro_data_ls file" +print * + +! File name for file that contains grid information +if ( halo.eq."-999" ) then ! global or nested tile + FV3_grid_input_file_name = "C" // trim(res_indx) // "_grid.tile" // & + trim(tile_num) // ".nc" +else ! stand-alone regional tile + FV3_grid_input_file_name = "C" // trim(res_indx) // "_grid.tile" // & + trim(tile_num) // ".halo" // trim(halo) // ".nc" +end if + +print *, "Reading from file: ", FV3_grid_input_file_name + +! Check that input file exists -- exit with error message if not +inquire(file=FV3_grid_input_file_name,exist=fexist) +if (.not.fexist) then + print * + print *, "Fatal error: Input file "//trim(FV3_grid_input_file_name)// & + " does not exist" + print *, "Exiting program gsl_oro_data.f90" + print * + call exit(4) +end if + + +! In preparation for reading in grid data, account for existence +! of halo points +if ( halo.eq."-999" ) then ! global or nested tile + halo_int = 0 +else + read(halo,*) halo_int ! integer form of halo +end if + + +! Open Cxxx_grid netCDF file for input and get dimensions +err = nf90_open(FV3_grid_input_file_name,nf90_nowrite,ncid_in) +call netcdf_err(err, 'opening: '//FV3_grid_input_file_name) + +err = nf90_inq_dimid(ncid_in,'nx',dimid) +call netcdf_err(err, 'reading nx id') +err = nf90_inquire_dimension(ncid_in,dimid,len=temp_int) +call netcdf_err(err, 'reading nx value') +dimX_FV3 = temp_int/2 - 2*halo_int ! shaving off any halo points from edges + +err = nf90_inq_dimid(ncid_in,'ny',dimid) +call netcdf_err(err, 'reading ny id') +err = nf90_inquire_dimension(ncid_in,dimid,len=temp_int) +call netcdf_err(err, 'reading ny value') +dimY_FV3 = temp_int/2 - 2*halo_int ! shaving off any halo points from edges + +print *, "dimX_FV3 =", dimX_FV3 ! number of model cells in x-direction +print *, "dimY_FV3 =", dimY_FV3 ! number of model cells in y-direction +print * + +! Read in lat/lon (in degrees) +allocate (lat_FV3_raw((2*dimX_FV3+1),(2*dimY_FV3+1))) +err = nf90_inq_varid(ncid_in,'y',varid) +call netcdf_err(err, 'reading y id') +err = nf90_get_var(ncid_in,varid,lat_FV3_raw, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3+1,2*dimY_FV3+1/)) +call netcdf_err(err, 'reading y') + +allocate (lon_FV3_raw((2*dimX_FV3+1),(2*dimY_FV3+1))) +err = nf90_inq_varid(ncid_in,'x',varid) +call netcdf_err(err, 'reading x id') +err = nf90_get_var(ncid_in,varid,lon_FV3_raw, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3+1,2*dimY_FV3+1/)) +call netcdf_err(err, 'reading x') + +! Read in quarter grid-cell areas +allocate (area_FV3_qtr((2*dimX_FV3),(2*dimY_FV3))) +err = nf90_inq_varid(ncid_in,'area',varid) +call netcdf_err(err, 'reading area id') +err = nf90_get_var(ncid_in,varid,area_FV3_qtr, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3,2*dimY_FV3/)) +call netcdf_err(err, 'reading area') + +! Calculate lat/lon at mass points (cell-centers) +! Stride by 2 starting with 2nd point +! NOTE: "Converting" from dbl_kind to real_kind +allocate (lat_FV3(dimX_FV3,dimY_FV3)) +allocate (lon_FV3(dimX_FV3,dimY_FV3)) +allocate (lat_FV3_deg(dimX_FV3,dimY_FV3)) +allocate (lon_FV3_deg(dimX_FV3,dimY_FV3)) +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + lat_FV3(i,j) = lat_FV3_raw(2*i,2*j) + lon_FV3(i,j) = lon_FV3_raw(2*i,2*j) + end do +end do +lat_FV3_deg(:,:) = lat_FV3(:,:) +lon_FV3_deg(:,:) = lon_FV3(:,:) +deallocate(lat_FV3_raw) +deallocate(lon_FV3_raw) + +! Convert lat/lon to radians +lat_FV3 = lat_FV3*pi/180._real_kind +lon_FV3 = lon_FV3*pi/180._real_kind + +! Create full grid-cell areas (4 raw areas per grid cell area) +! NOTE: "Converting" from dbl_kind to real_kind +allocate(area_FV3(dimX_FV3,dimY_FV3)) +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + area_FV3(i,j) = area_FV3_qtr(2*i-1,2*j-1) + area_FV3_qtr(2*i-1,2*j) + & + area_FV3_qtr(2*i,2*j-1) + area_FV3_qtr(2*i,2*j) + end do +end do +deallocate(area_FV3_qtr) + +err = nf90_close(ncid_in) + + + +! Open file containing 2.5min topo data (fine grid) +fine_topo_source_file_name = "geo_em.d01.lat-lon.2.5m.HGT_M.nc" +! Check that input file exists -- exit with error message if not +inquire(file=fine_topo_source_file_name,exist=fexist) +if (.not.fexist) then + print * + print *, "Fatal error: Topo source file "// & + trim(fine_topo_source_file_name)//" does not exist" + print *, "Exiting program gsl_oro_data.f90" + print * + call exit(4) +end if +err = nf90_open(trim(fine_topo_source_file_name),nf90_nowrite,ncid_in) +call netcdf_err(err, 'opening: '//trim(fine_topo_source_file_name)) + +! Get dimensions +err = nf90_inq_dimid(ncid_in,'west_east',dimid) +call netcdf_err(err, 'reading west_east id') +err = nf90_inquire_dimension(ncid_in,dimid,len=dimX_fine) +call netcdf_err(err, 'reading west_east value') + +err = nf90_inq_dimid(ncid_in,'south_north',dimid) +call netcdf_err(err, 'reading south_north id') +err = nf90_inquire_dimension(ncid_in,dimid,len=dimY_fine) +call netcdf_err(err, 'reading south_north value') + +print *, "Source file for high-resolution topography: ", & + trim(fine_topo_source_file_name) +print *, "dimX_fine =", dimX_fine +print *, "dimY_fine =", dimY_fine +print * + + +! Read in lat/lon of fine grid +allocate(lat1d_fine(dimY_fine)) +allocate(lon1d_fine(dimX_fine)) +err = nf90_inq_varid(ncid_in,'CLAT',varid) +call netcdf_err(err, 'reading CLAT id') +err = nf90_get_var(ncid_in,varid,lat1d_fine,start=(/1,1/), & + count=(/1,dimY_fine/)) +call netcdf_err(err, 'reading CLAT') + +err = nf90_inq_varid(ncid_in,'CLONG',varid) +call netcdf_err(err, 'reading CLONG id') +err = nf90_get_var(ncid_in,varid,lon1d_fine,start=(/1,1/), & + count=(/dimX_fine,1/)) +call netcdf_err(err, 'reading CLONG') + +! Convert lat/lon to radians +lat1d_fine = lat1d_fine*pi/180._real_kind +lon1d_fine = lon1d_fine*pi/180._real_kind + + +! Reassign FV3 longitude to vary from -pi to pi to match lon1d_fine range +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + if ( lon_FV3(i,j).gt.pi ) then + lon_FV3(i,j) = lon_FV3(i,j) - 2*pi + end if + end do +end do + + +! Read in fine-scale topography +allocate(HGT_M_fine(dimX_fine,dimY_fine)) +err = nf90_inq_varid(ncid_in,'HGT_M',varid) +call netcdf_err(err, 'reading HGT_M id') +err = nf90_get_var(ncid_in,varid,HGT_M_fine,start=(/1,1/), & + count=(/dimX_fine,dimY_fine/)) +call netcdf_err(err, 'reading HGT_M') + + +err = nf90_close(ncid_in) + + +! Allocate GWD statistics fields +allocate (std_dev(dimX_FV3*dimY_FV3)) +allocate (convexity(dimX_FV3*dimY_FV3)) +allocate (OA1(dimX_FV3*dimY_FV3)) +allocate (OA2(dimX_FV3*dimY_FV3)) +allocate (OA3(dimX_FV3*dimY_FV3)) +allocate (OA4(dimX_FV3*dimY_FV3)) +allocate (OL1(dimX_FV3*dimY_FV3)) +allocate (OL2(dimX_FV3*dimY_FV3)) +allocate (OL3(dimX_FV3*dimY_FV3)) +allocate (OL4(dimX_FV3*dimY_FV3)) + +! Initialize GWD statistics fields +std_dev(:) = 0._real_kind +convexity(:) = 0._real_kind +OA1(:) = 0._real_kind +OA2(:) = 0._real_kind +OA3(:) = 0._real_kind +OA4(:) = 0._real_kind +OL1(:) = 0._real_kind +OL2(:) = 0._real_kind +OL3(:) = 0._real_kind +OL4(:) = 0._real_kind + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! This is a loop over all the FV3 (coarse) grid cells +! The subgrid-scale topographic variables needed for the large-scale +! orographic gravity wave drag schemes are calculated by the following steps: +! 1) Sample the fine-scale (2.5min) topography contained within each +! coarse grid cell. +! 2) Detrend the topography by subtracting a bilinear-interpolated height field +! from the fine-scale height field (if detrend_topography = .true.), +! otherwise actual fine-scale height field is used to calculate statistics +! 3) Calculate the orographic statistics: stddev,convexity,oa1,...oa4, +! ol1,...,ol4 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +cell_count = 1 + +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + + ! Calculate approximate side-lengths of square lat-long "coarse" grid + ! cell centered on FV3 cell (units = radians) + dlta_lat = sqrt(area_FV3(i,j))/ae + dlta_lon = sqrt(area_FV3(i,j))/(ae*COS(lat_FV3(i,j))) + + ! Determine lat/lon of 9 lat-lon block centers + ! Note: lat_blk(2)/lon_blk(2) = lat_FV3(i,j)/lon_FV3(i,j) + ! Note: abs(lon_blk) may exceed pi + do i_blk = 1,3 + lon_blk(i_blk) = lon_FV3(i,j) + (i_blk-2)*dlta_lon + end do + ! Note: abs(lat_blk) may exceed pi/2 (90 degrees) + do j_blk = 1,3 + lat_blk(j_blk) = lat_FV3(i,j) + (j_blk-2)*dlta_lat + end do + + ! Find starting and ending fine-grid i,j indices for each + ! of the 9 "coarse-grid" blocks + ! Note: Index value of -999 is returned if latitude of grid points + ! exceed 90 degrees north or south + do i_blk = 1,3 + s_ii(i_blk) = nearest_i_east(lon_blk(i_blk)-p5*dlta_lon) + e_ii(i_blk) = nearest_i_west(lon_blk(i_blk)+p5*dlta_lon) + end do + do j_blk = 1,3 + s_jj(j_blk) = nearest_j_north(lat_blk(j_blk)-p5*dlta_lat) + e_jj(j_blk) = nearest_j_south(lat_blk(j_blk)+p5*dlta_lat) + end do + + + ! Calculate mean topographic height in each "coarse grid" block + ! Note: We only do the mean-height calculation if we are detrending + ! the subgrid topography, otherwise, we still need the + ! fine-grid indices for the block limits -- s_ii, etc. + do i_blk = 1,3 + + ! "Shave" blocks on north or south due to proximity to poles + ! if necessary + j_blk = 1 ! southern row + ! Check for "shaved" block due to proximity to south pole + if ( (s_jj(j_blk).eq.-999).and.(e_jj(j_blk).ne.-999) ) then + s_jj(j_blk) = 1 ! southern boundary of shaved block + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(1)+lat1d_fine(e_jj(j_blk))) + end if + + j_blk = 2 ! center row + ! Check for "shaved" block due to proximity to south or north pole + ! Note: We're assuming e_jj(2) and s_jj(2) can't both be -999 + if ( s_jj(j_blk).eq.-999 ) then + s_jj(j_blk) = 1 ! block shaved on the south + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(1)+lat1d_fine(e_jj(j_blk))) + end if + if ( e_jj(j_blk).eq.-999 ) then + e_jj(j_blk) = dimY_fine ! block shaved on the north + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(s_jj(j_blk))+lat1d_fine(dimY_fine)) + end if + + j_blk = 3 ! northern row + ! Check for "shaved" block due to proximity to north pole + if ( (e_jj(j_blk).eq.-999).and.(s_jj(j_blk).ne.-999) ) then + e_jj(j_blk) = dimY_fine ! northern boundary of shaved block + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(s_jj(j_blk))+lat1d_fine(dimY_fine)) + end if + + if ( detrend_topography ) then + do j_blk = 1,3 + call calc_mean_HGT(s_ii(i_blk),e_ii(i_blk), & + s_jj(j_blk),e_jj(j_blk),HGT_M_coarse(i_blk,j_blk)) + ! Note: If there is no block because s_jj and e_jj are + ! both -999 HGT_M_coarse returns with a "missing" + ! value of HGT_missing = 1.E+10 + end do + end if + + end do + + + ! Calculate number of fine-grid points within center coarse block (2,2) + ! Check if center block straddles date line + if ( s_ii(2).gt.e_ii(2) ) then + ii_m = dimX_fine - s_ii(2) + 1 + e_ii(2) + else + ii_m = e_ii(2) - s_ii(2) + 1 + end if + jj_m = e_jj(2) - s_jj(2) + 1 + + + ! Bilinearly interpolate coarse-grid topography of the 9 blocks to + ! fine grid for the purpose of detrending the fine-grid topography + ! to represent the sub-grid topography + ! Note: The detrending only occurs within the center coarse block (2,2) + if ( detrend_topography ) then + + ! i,j indices of HGT_M_coarse_on_fine range from 1,ii_m and 1,jj_m + ! i.e., a "local" index system + allocate (HGT_M_coarse_on_fine(ii_m,jj_m)) + + do jj = s_jj(2), e_jj(2) + jj_loc = jj - s_jj(2) + 1 ! local j-index (1 ... jj_m) + ! Check if block straddles the date line + if ( s_ii(2).gt.e_ii(2) ) then + do ii = s_ii(2), dimX_fine ! west of the date line + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + call HGT_interpolate(lat1d_fine(jj),lon1d_fine(ii), & + lat_blk(:),lon_blk(:),HGT_M_coarse(:,:), & + HGT_M_coarse_on_fine(ii_loc,jj_loc)) + end do + do ii = 1, e_ii(2) ! east of the date line + ii_loc = ii_loc + 1 ! local i-index ( 1 ... ii_m ) + call HGT_interpolate(lat1d_fine(jj),lon1d_fine(ii), & + lat_blk(:),lon_blk(:),HGT_M_coarse(:,:), & + HGT_M_coarse_on_fine(ii_loc,jj_loc)) + end do + else ! no crossing of the date line + do ii = s_ii(2), e_ii(2) + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + call HGT_interpolate(lat1d_fine(jj),lon1d_fine(ii), & + lat_blk(:),lon_blk(:),HGT_M_coarse(:,:), & + HGT_M_coarse_on_fine(ii_loc,jj_loc)) + end do + end if + end do + + end if + + + ! Assign values to "zs", which is the fine-grid surface topography field + ! that we will calculate statistics on, i.e, stddev, convexity, etc. + ! This will either be the detrended values (detrend_topography = .true.) + ! or the actual values (detrend_topography = .false.) + allocate (zs(ii_m,jj_m)) + + do jj = s_jj(2), e_jj(2) + jj_loc = jj - s_jj(2) + 1 ! local j-index (1 ... jj_m) + ! Check if block straddles the date line + if ( s_ii(2).gt.e_ii(2) ) then + do ii = s_ii(2), dimX_fine ! west of the date line + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + if ( detrend_topography ) then + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) - & + HGT_M_coarse_on_fine(ii_loc,jj_loc) + else + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end if + end do + do ii = 1, e_ii(2) ! east of the date line + ii_loc = ii_loc + 1 ! local i-index ( 1 ... ii_m ) + if ( detrend_topography ) then + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) - & + HGT_M_coarse_on_fine(ii_loc,jj_loc) + else + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end if + end do + else ! no crossing of the date line + do ii = s_ii(2), e_ii(2) + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + if ( detrend_topography ) then + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) - & + HGT_M_coarse_on_fine(ii_loc,jj_loc) + else + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end if + end do + end if + end do + + + + ! + ! Finally, we can now calculate the topographic statistics fields needed + ! for the gravity wave drag scheme + ! + + ! Make sure statistics are zero if there is no terrain in the grid cell + zs_accum = .false. + do jj = 1,jj_m + do ii = 1,ii_m + if ( abs(zs(ii,jj)).gt.1.E-3 ) zs_accum = .true. + end do + end do + if ( .not.zs_accum ) then ! no terrain in the grid cell + std_dev(cell_count) = 0._real_kind + convexity(cell_count) = 0._real_kind + OA1(cell_count) = 0._real_kind + OA2(cell_count) = 0._real_kind + OA3(cell_count) = 0._real_kind + OA4(cell_count) = 0._real_kind + OL1(cell_count) = 0._real_kind + OL2(cell_count) = 0._real_kind + OL3(cell_count) = 0._real_kind + OL4(cell_count) = 0._real_kind + if ( detrend_topography ) deallocate (HGT_M_coarse_on_fine) + deallocate(zs) + cell_count = cell_count + 1 + cycle ! move on to next (coarse) grid cell + end if + + + ! + ! Calculate standard deviation of subgrid-scale terrain height + ! + + ! Calculate mean height + sum2 = 0._real_kind + nfinepoints = ii_m*jj_m + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + zs(ii,jj) + end do + end do + zs_mean = sum2 / real(nfinepoints,real_kind) + + ! Calculate standard deviation + sum2 = 0._real_kind + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + ( zs(ii,jj) - zs_mean )**2 + end do + end do + std_dev(cell_count) = sqrt( sum2/real(nfinepoints,real_kind) ) + + + ! + ! Calculate convexity of sub-grid-scale terrain + ! + + sum2 = 0._real_kind + sum4 = 0._real_kind + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + ( zs(ii,jj) - zs_mean )**2 + sum4 = sum4 + ( zs(ii,jj) - zs_mean )**4 + end do + end do + + var = sum2 / real(nfinepoints,real_kind) + if ( abs(var) < 1.0E-05_real_kind ) then + convexity(cell_count) = 0._real_kind + else + convexity(cell_count) = min( sum4 / ( var**2 * & + real(nfinepoints,real_kind) ), max_convexity ) + end if + + + ! + ! Calculate orographic asymmetries + ! + + ! OA1 -- orographic asymmetry in West direction + nu = 0 + nd = 0 + do jj = 1,jj_m + do ii = 1,ii_m/2 ! left half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + end do + do ii = ii_m/2 + 1, ii_m ! right half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end do + end do + if ( nu + nd > 0 ) then + OA1(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA1(cell_count) = 0._real_kind + end if + + ! OA2 -- orographic asymmetry in South direction + nu = 0 + nd = 0 + do jj = 1,jj_m/2 ! bottom half of box + do ii = 1,ii_m + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + end do + end do + do jj = jj_m/2 + 1,jj_m ! top half of box + do ii = 1, ii_m + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end do + end do + if ( nu + nd > 0 ) then + OA2(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA2(cell_count) = 0._real_kind + end if + + ! OA3 -- orographic asymmetry in South-West direction + nu = 0 + nd = 0 + ratio = real(jj_m,real_kind)/real(ii_m,real_kind) + do jj = 1,jj_m + do ii = 1,ii_m + if ( nint(real(ii,real_kind)*ratio) < (jj_m - jj) ) then + ! south-west half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + else ! north-east half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end if + end do + end do + if ( nu + nd > 0 ) then + OA3(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA3(cell_count) = 0._real_kind + end if + + ! OA4 -- orographic asymmetry in North-West direction + nu = 0 + nd = 0 + ratio = real(jj_m,real_kind)/real(ii_m,real_kind) + do jj = 1,jj_m + do ii = 1,ii_m + if ( nint(real(ii,real_kind)*ratio) < jj ) then + ! north-west half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + else ! south-east half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end if + end do + end do + if ( nu + nd > 0 ) then + OA4(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA4(cell_count) = 0._real_kind + end if + + + ! + ! Calculate orographic effective lengths + ! + + ! OL1 -- orographic effective length for Westerly flow + nw = 0 + nt = 0 + do jj = max(jj_m/4,1), 3*jj_m/4 + ! within central east-west band of box + do ii = 1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL1(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL1(cell_count) = 0._real_kind + end if + + ! OL2 -- orographic effective length for Southerly flow + nw = 0 + nt = 0 + do jj = 1, jj_m + do ii = max(ii_m/4,1), 3*ii_m/4 + ! within central north-south band of box + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL2(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL2(cell_count) = 0._real_kind + end if + + ! OL3 -- orographic effective length for South-Westerly flow + nw = 0 + nt = 0 + do jj = 1, jj_m/2 + do ii = 1, ii_m/2 + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + do jj = jj_m/2+1, jj_m + do ii = ii_m/2+1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL3(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL3(cell_count) = 0._real_kind + end if + + ! OL4 -- orographic effective length for North-Westerly flow + nw = 0 + nt = 0 + do jj = jj_m/2+1, jj_m + do ii = 1, ii_m/2 + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + do jj = 1, jj_m/2 + do ii = ii_m/2+1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL4(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL4(cell_count) = 0._real_kind + end if + + + + if ( detrend_topography ) deallocate (HGT_M_coarse_on_fine) + deallocate (zs) + + cell_count = cell_count + 1 + + end do ! j = 1,dimY_FV3 +end do ! i = 1,dimX_FV3 + + + +! +! Output GWD statistics fields to netCDF file +! + + +if ( halo.eq."-999" ) then ! global or nested tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ls.tile" & + // trim(tile_num) // ".nc" +else ! stand-alone regional tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ls.tile" & + // trim(tile_num) // ".halo0.nc" +end if + +! Open netCDF file for output +err = nf90_create(oro_data_output_file_name, NF90_CLOBBER, ncid_out) +call netcdf_err(err, 'creating: '//oro_data_output_file_name) + +err = nf90_redef(ncid_out) + +! Define dimensions +err = nf90_def_dim(ncid_out,'lon',dimX_FV3,lonid) +call netcdf_err(err, 'defining lon dimension') +err = nf90_def_dim(ncid_out,'lat',dimY_FV3,latid) +call netcdf_err(err, 'defining lat dimension') + +! Define the 'dimensions vector' dimids to be used for writing +! the 2-dimensional variables to the netCDF file +dimids(1) = lonid +dimids(2) = latid + +! Define variables and attributes to put in the netCDF file +err = nf90_def_var(ncid_out,'geolon',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining geolon') +err = nf90_put_att(ncid_out,varid,'units','degrees') +err = nf90_put_att(ncid_out,varid,'description','longitude') +err = nf90_def_var(ncid_out,'geolat',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining geolat') +err = nf90_put_att(ncid_out,varid,'units','degrees') +err = nf90_put_att(ncid_out,varid,'description','latitude') +err = nf90_def_var(ncid_out,'stddev',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'stddev') +err = nf90_put_att(ncid_out,varid,'units','meters') +err = nf90_put_att(ncid_out,varid,'description', & + 'standard deviation of subgrid topography') +err = nf90_def_var(ncid_out,'convexity',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining convexity') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'convexity of subgrid topography') +err = nf90_def_var(ncid_out,'oa1',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa1') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in west direction') +err = nf90_def_var(ncid_out,'oa2',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa2') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south direction') +err = nf90_def_var(ncid_out,'oa3',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa3') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south-west direction') +err = nf90_def_var(ncid_out,'oa4',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa4') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in north-west direction') +err = nf90_def_var(ncid_out,'ol1',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol1') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for westerly flow') +err = nf90_def_var(ncid_out,'ol2',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol2') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for southerly flow') +err = nf90_def_var(ncid_out,'ol3',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol3') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for south-westerly flow') +err = nf90_def_var(ncid_out,'ol4',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol4') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for north-westerly flow') + +! Add global attributes +err = nf90_put_att(ncid_out,nf90_global, & + 'source_file_for_high-resolution_topography', & + trim(fine_topo_source_file_name)) +if ( detrend_topography ) then + err = nf90_put_att(ncid_out,nf90_global, & + 'high-res_topography_detrended','.TRUE.') +else + err = nf90_put_att(ncid_out,nf90_global, & + 'high-res_topography_detrended','.FALSE.') +end if + +err = nf90_enddef(ncid_out) + + +! Write data to output netCDF file +err = nf90_inq_varid(ncid_out,'geolon',varid) +call netcdf_err(err, 'reading geolon id') +err = nf90_put_var(ncid_out,varid,lon_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing geolon') +err = nf90_inq_varid(ncid_out,'geolat',varid) +call netcdf_err(err, 'reading geolat id') +err = nf90_put_var(ncid_out,varid,lat_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing geolat') +err = nf90_inq_varid(ncid_out,'stddev',varid) +call netcdf_err(err, 'reading stddev id') +err = nf90_put_var(ncid_out,varid,std_dev,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing stddev') +err = nf90_inq_varid(ncid_out,'convexity',varid) +call netcdf_err(err, 'reading convexity id') +err = nf90_put_var(ncid_out,varid,convexity,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing convexity') +err = nf90_inq_varid(ncid_out,'oa1',varid) +call netcdf_err(err, 'reading oa1 id') +err = nf90_put_var(ncid_out,varid,OA1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa1') +err = nf90_inq_varid(ncid_out,'oa2',varid) +call netcdf_err(err, 'reading oa2 id') +err = nf90_put_var(ncid_out,varid,OA2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa2') +err = nf90_inq_varid(ncid_out,'oa3',varid) +call netcdf_err(err, 'reading oa3 id') +err = nf90_put_var(ncid_out,varid,OA3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa3') +err = nf90_inq_varid(ncid_out,'oa4',varid) +call netcdf_err(err, 'reading oa4 id') +err = nf90_put_var(ncid_out,varid,OA4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa4') +err = nf90_inq_varid(ncid_out,'ol1',varid) +call netcdf_err(err, 'reading ol1 id') +err = nf90_put_var(ncid_out,varid,OL1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol1') +err = nf90_inq_varid(ncid_out,'ol2',varid) +call netcdf_err(err, 'reading ol2 id') +err = nf90_put_var(ncid_out,varid,OL2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol2') +err = nf90_inq_varid(ncid_out,'ol3',varid) +call netcdf_err(err, 'reading ol3 id') +err = nf90_put_var(ncid_out,varid,OL3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol3') +err = nf90_inq_varid(ncid_out,'ol4',varid) +call netcdf_err(err, 'reading ol4 id') +err = nf90_put_var(ncid_out,varid,OL4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol4') + +err = nf90_close(ncid_out) + + + +! Deallocate arrays +deallocate(lat_FV3) +deallocate(lon_FV3) +deallocate(lat_FV3_deg) +deallocate(lon_FV3_deg) +deallocate(area_FV3) +deallocate(lat1d_fine) +deallocate(lon1d_fine) +deallocate(HGT_M_fine) +deallocate(std_dev) +deallocate(convexity) +deallocate(OA1) +deallocate(OA2) +deallocate(OA3) +deallocate(OA4) +deallocate(OL1) +deallocate(OL2) +deallocate(OL3) +deallocate(OL4) + +end subroutine calc_gsl_oro_data_lg_scale + +!> Calculates average terrain height within coarse grid cell ("block") +!! +!! @param[in] s_ii Fine grid starting i-index +!! @param[in] e_ii Fine grid ending i-index +!! @param[in] s_jj Fine grid starting j-index +!! @param[in] e_jj Fine grid ending j-index +!! @param[out] hgt Fine grid height (m) +!! @author Michael Toy, NOAA/GSL +subroutine calc_mean_HGT(s_ii,e_ii,s_jj,e_jj,HGT) + +! This subroutine calculates the average terrain height within +! coarse grid cell ("block") + +implicit none + +integer :: s_ii, & ! starting fine-grid i-index + e_ii, & ! ending fine-grid i-index + s_jj, & ! starting fine-grid j-index + e_jj ! ending fine-grid j-index +real (kind=real_kind), intent(out) :: HGT + +! Local variables +integer :: i,j,grid_pt_count +real (kind=real_kind) :: HGT_sum + + +! Return a value of 0 if s_jj and e_jj are both -999, +! i.e., if there is no block adjoining the center row +! due to proximity to one of the poles +! Note: The HGT value of the block will be ignored +if ( (s_jj.eq.-999).and.(e_jj.eq.-999) ) then + HGT = HGT_missing + return +end if + +grid_pt_count = 0 +HGT_sum = 0._real_kind +do j = s_jj, e_jj + ! Note: If the grid block straddles the date line, then s_ii > e_ii + ! We need to correct for this + if ( s_ii.gt.e_ii ) then ! straddling the date line + do i = s_ii, dimX_fine ! west of the date line + HGT_sum = HGT_sum + HGT_M_fine(i,j) + grid_pt_count = grid_pt_count + 1 + end do + do i = 1, e_ii ! east of the date line + HGT_sum = HGT_sum + HGT_M_fine(i,j) + grid_pt_count = grid_pt_count + 1 + end do + else ! no crossing of the date line + do i = s_ii, e_ii + HGT_sum = HGT_sum + HGT_M_fine(i,j) + grid_pt_count = grid_pt_count + 1 + end do + end if +end do +HGT = HGT_sum/grid_pt_count + +end subroutine calc_mean_HGT + +!> Interpolates height from coarse grid on to fine grid points +!! +!! @param[in] lat Latitude of fine grid point. +!! @param[in] lon_in Longitude of fine grid point. +!! @param[in] lat_blk Latitudes of neighboring coarse grid points. +!! @param[in] lon_blk Longitudes of neighboring coarse grid points. +!! @param[in] hgt_coarse Topographic heights on coarse grid +!! @param[out] hgt_coarse_on_fine Coarse grid heights interpolated on to fine grid +!! @author Michael Toy, NOAA/GSL +subroutine HGT_interpolate(lat,lon_in,lat_blk,lon_blk,HGT_coarse, & + HGT_coarse_on_fine) + +! This subroutine bilinearly interpolates neighboring coarse-grid terrain +! heights (HGT_coarse) to fine-grid points (HGT_coarse_on_fine) +! (extrapolates in the case near poles) +! Note: Bilinear interpolation is done by calling a 1D interpolation +! function of a 1D interpolation function + +implicit none + +real (kind = real_kind), intent(in) :: & + lat, & ! latitude of fine grid point + lon_in ! longitude of fine grid point +real (kind = real_kind), dimension(3),intent(in) :: & + lat_blk, & ! latitudes of neighboring coarse grid points + lon_blk ! longitudes of neighboring coarse grid points +real (kind = real_kind), dimension(3,3), intent(in) :: HGT_coarse +real (kind = real_kind), intent(out) :: HGT_coarse_on_fine +real (kind = real_kind) :: lon + + +lon = lon_in +! We need to make sure that if we're straddling the date line, that +! we remove the possible 2*pi discontinuity between lon and +! {lon_blk(1),lon_blk(2),lon_blk(3)) for interpolation purposes +! This will line the 4 longitudes up monotonically +if ( abs(lon_in-lon_blk(2)).gt.pi ) then ! discontinuity exists + if ( lon_in.gt.0. ) lon = lon - 2*pi ! lon_in lies west of date line + if ( lon_in.lt.0. ) lon = lon + 2*pi ! lon_in lies east of date line +end if + + +! Check for need to extrapolate if top or bottom block rows +! have height = HGT_missing + +! Check for missing north row +if ( (HGT_coarse(1,3).eq.HGT_missing).or.(HGT_coarse(2,3).eq.HGT_missing) & + .or.(HGT_coarse(3,3).eq.HGT_missing) ) then + + ! Determine which quadrant of the coarse grid cell we are in + if ( (lat.ge.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant I + ! Extrapolate from lat_blk(1) and lat_blk(2) + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(3,1),HGT_coarse(3,2)) ) + elseif ( (lat.ge.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant II + ! Extrapolate from lat_blk(1) and lat_blk(2) + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(1,1),HGT_coarse(1,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)) ) + elseif ( (lat.lt.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant III + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(1,1),HGT_coarse(1,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)) ) + elseif ( (lat.lt.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant IV + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(3,1),HGT_coarse(3,2)) ) + end if + + return +end if + +! Check for missing south row +if ( (HGT_coarse(1,1).eq.HGT_missing).or.(HGT_coarse(2,1).eq.HGT_missing) & + .or.(HGT_coarse(3,1).eq.HGT_missing) ) then + + ! Determine which quadrant of the coarse grid cell we are in + if ( (lat.ge.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant I + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(3,2),HGT_coarse(3,3)) ) + elseif ( (lat.ge.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant II + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(1,2),HGT_coarse(1,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)) ) + elseif ( (lat.lt.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant III + ! Extrapolate from lat_blk(2) and lat_blk(3) + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(1,2),HGT_coarse(1,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)) ) + elseif ( (lat.lt.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant IV + ! Extrapolate from lat_blk(2) and lat_blk(3) + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(3,2),HGT_coarse(3,3)) ) + end if + + return +end if + +! Interpolation only +! Determine which quadrant of the coarse grid cell we are in +if ( (lat.ge.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant I + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(3,2),HGT_coarse(3,3)) ) +elseif ( (lat.ge.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant II + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(1,2),HGT_coarse(1,3)), & + interp_1d(lat,lat_blk(2),lat_blk(3),HGT_coarse(2,2),HGT_coarse(2,3)) ) +elseif ( (lat.lt.lat_blk(2)).and.(lon.lt.lon_blk(2)) ) then ! Quadrant III + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(1),lon_blk(2), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(1,1),HGT_coarse(1,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)) ) +elseif ( (lat.lt.lat_blk(2)).and.(lon.ge.lon_blk(2)) ) then ! Quadrant IV + HGT_coarse_on_fine = interp_1d( & + lon,lon_blk(2),lon_blk(3), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(2,1),HGT_coarse(2,2)), & + interp_1d(lat,lat_blk(1),lat_blk(2),HGT_coarse(3,1),HGT_coarse(3,2)) ) +end if + +end subroutine HGT_interpolate + +!> Finds nearest fine-grid i index to the east of a given longitude +!! +!! @param[in] lon_in longitude (radians) +!! @return nearest_i_east Nearest grid point i-index east of selected point +!! @author Michael Toy, NOAA/GSL +function nearest_i_east(lon_in) +! Calculates nearest fine-grid i index to the east of (or on) a given longitude +implicit none + +integer :: nearest_i_east +real (kind=real_kind), intent(in) :: lon_in +real (kind=real_kind) :: lon +integer :: i + +lon = lon_in +! Make sure longitude is between -pi and pi +do while ( (lon.lt.(-pi)).or.(lon.gt.pi) ) + if ( lon.lt.(-pi) ) lon = lon + 2*pi + if ( lon.gt.pi ) lon = lon - 2*pi +end do + +if ( lon.gt.lon1d_fine(dimX_fine) ) then + nearest_i_east = 1 +else + i = 1 + do while ( lon1d_fine(i).lt.lon ) + i = i + 1 + end do + nearest_i_east = i +end if + +end function nearest_i_east + +!> Finds nearest fine-grid i index to the west of a given longitude +!! +!! @param[in] lon_in longitude (radians) +!! @return nearest_i_west Nearest grid point i-index west of selected point +!! @author Michael Toy, NOAA/GSL +function nearest_i_west(lon_in) +! Calculates nearest fine-grid i index to the west of a given longitude +implicit none + +integer :: nearest_i_west +real (kind=real_kind), intent(in) :: lon_in +real (kind=real_kind) :: lon +integer :: i + +lon = lon_in +! Make sure longitude is between -pi and pi +do while ( (lon.lt.(-pi)).or.(lon.gt.pi) ) + if ( lon.lt.(-pi) ) lon = lon + 2*pi + if ( lon.gt.pi ) lon = lon - 2*pi +end do + +if ( (lon.lt.lon1d_fine(1)).or.(lon.ge.lon1d_fine(dimX_fine)) ) then + nearest_i_west = dimX_fine +else + i = 1 + do while ( lon1d_fine(i).le.lon ) + i = i + 1 + end do + nearest_i_west = i - 1 +end if + +end function nearest_i_west + +!> Calculates nearest fine-grid j index to the north of a given latitude +!! +!! @param[in] lat_in Latitude (radians) +!! @return nearest_j_north Nearest fine-grid j index to the north of a given latitude +!! @author Michael Toy, NOAA/GSL +function nearest_j_north(lat_in) +! Calculates nearest fine-grid j index to the north of a given latitude +! Note: If the abs(latitude) is greater than pi/2 (90 degrees) then +! the value -999 is returned +implicit none + +integer :: nearest_j_north +real (kind=real_kind), intent(in) :: lat_in +real (kind=real_kind) :: lat +integer :: j + +lat = lat_in +if ( abs(lat_in).gt.p5*pi ) then + nearest_j_north = -999 +else + j = 1 + do while ( (lat1d_fine(j).lt.lat).and.(j.lt.dimY_fine) ) + j = j + 1 + end do + nearest_j_north = j +end if + +end function nearest_j_north + +!> Calculates nearest fine-grid j index to the south of a given latitude +!! +!! @param[in] lat_in Latitude (radians) +!! @return nearest_j_south Nearest fine-grid j index to the south of a given latitude +!! @author Michael Toy, NOAA/GSL +function nearest_j_south(lat_in) +! Calculates nearest fine-grid j index to the south of a given latitude +! Note: If the abs(latitude) is greater than pi/2 (90 degrees) then +! the value -999 is returned +implicit none + +integer :: nearest_j_south +real (kind=real_kind), intent(in) :: lat_in +real (kind=real_kind) :: lat +integer :: j + +lat = lat_in +if ( abs(lat_in).gt.p5*pi ) then + nearest_j_south = -999 +elseif ( lat_in.le.lat1d_fine(1) ) then + nearest_j_south = 1 +else + j = 2 + do while ( (lat1d_fine(j).le.lat).and.(j.le.dimY_fine) ) + j = j + 1 + end do + nearest_j_south = j - 1 +end if + +end function nearest_j_south + +!> Interpolates (or extrapolates) linear function y = y(x) +!! +!! @param[in] x Input "x" value +!! @param[in] x1 Known point 1 +!! @param[in] x2 Known point 2 +!! @param[in] y1 Known y(x1) +!! @param[in] y2 Known y(x2) +!! @return interp_1d Interpolated y value at x +!! @author Michael Toy, NOAA/GSL +function interp_1d(x,x1,x2,y1,y2) +! Interpolates (or extrapolates) linear function y = y(x) +! to x given y1 = y(x1) and y2 = y(x2) +implicit none + +real (kind=real_kind) :: interp_1d +real (kind=real_kind), intent(in) :: x,x1,x2,y1,y2 +real (kind=real_kind) :: slope + +! Formula for a line: y = y1 + slope*(x - x1) +slope = (y2-y1)/(x2-x1) +interp_1d = y1 + slope*(x-x1) + +end function interp_1d + +!> Returns netCDF error given input err code +!! +!! @param[in] err Error code from netCDF routine +!! @param[in] string Portion of error message +!! @author Michael Toy, NOAA/GSL +subroutine netcdf_err(err,string) + +use netcdf + +implicit none + +integer, intent(in) :: err +character(len=*), intent(in) :: string +character(len=256) :: errmsg + +if (err.eq.NF90_NOERR ) return +errmsg = NF90_STRERROR(err) +print *, "" +print *, "FATAL ERROR: ", trim(string), ": ", trim(errmsg) +print *, "STOP." +call exit(4) + +return +end subroutine netcdf_err + +end module gsl_oro_data_lg_scale diff --git a/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_sm_scale.f90 b/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_sm_scale.f90 new file mode 100644 index 000000000..ce059df5d --- /dev/null +++ b/sorc/orog_mask_tools.fd/orog_gsl.fd/module_gsl_oro_data_sm_scale.f90 @@ -0,0 +1,1270 @@ +!> @file +!! @brief Calculates small-scale GWD orographic stats for FV3GFS grids +!! @author Michael Toy, NOAA/GSL +!! @date 2021-03-12 +!! +!! This module calculates the parameters required for the subgrid- +!! scale orographic gravity-wave drag (GWDO) scheme on the FV3 +!! grid. These parameters are for the small-scale GWD (Tsiringakis et al., +!! 2017) and the turbulent orographic form drag (TOFD) (Beljaars, 2004) +!! schemes of the GSL drag suite. 30 second (~1km) global topography +!! is used. The topographic data comes from the 'fix' file +!! HGT.Beljaars_filtered.lat-lon.30s_res.nc. +!! The output fields are: +!! - stddev standard deviation of subgrid-scale topograpy +!! - convexity convexity (kurtosis) of subgrid-scale topography +!! - ol{1,2,3,4} orographic effective lengths of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! - oa{1,2,3,4} orographic asymmetries of subgrid-scale topography +!! for 4 orientations: 1-westerly, 2-southerly, 3-southwesterly, 4-northwesterly +!! +!! Based on code by Michael Duda provided by NCAR/MMM +!! +module gsl_oro_data_sm_scale + +implicit none + +integer, parameter :: real_kind = selected_real_kind(6) !< single precision +integer, parameter :: dbl_kind = selected_real_kind(13) !< double precision + +real, parameter :: pi = 3.1415926535897_real_kind !< pi +integer :: dimX_fine !< x-dimension of fine grid +integer :: dimY_fine !< y-dimension of fine grid + +real (kind = real_kind), allocatable :: lat1d_fine(:) !< latitude of fine grid pts +real (kind = real_kind), allocatable :: lon1d_fine(:) !< longitude of fine grid pts + +real (kind = real_kind), parameter :: p5 = 0.5_real_kind !< one half + +contains + +!> Subroutine to compute orographic statistics needed for small-scale +!! orograhic drag (gravity wave and form drag) schemes +!! +!! @param[in] tile_num (tile number) +!! @param[in] res_indx (resolution) +!! @param[in] halo (halo number) +!! @param[out] duplicate_oro_data_file (equals .true. if min grid size <= 7.5km) +!! @author Michael Toy, NOAA/GSL +subroutine calc_gsl_oro_data_sm_scale(tile_num,res_indx,halo, & + duplicate_oro_data_file) + +use netcdf + +implicit none + +character(len=2), intent(in) :: tile_num ! tile number +character(len=7), intent(in) :: res_indx ! grid-resolution index, e.g., 96, 192, 384, etc +character(len=4), intent(in) :: halo ! halo value (for input grid data) + +logical, intent(out) :: duplicate_oro_data_file ! flag to let main program know that + ! oro_data_ls file was created by this subroutine + ! as dupliate of oro_data_ss file due to grid size + ! being below 7.5km + +real (kind = real_kind) :: min_area_FV3 ! minimum grid area in m^2 +real (kind = real_kind) :: min_DX ! minimum grid size in km + +integer :: i,j,ii,jj + +integer :: ncid_in,ncid_out,err +integer :: varid +integer :: dimid,latid,lonid +integer, dimension (2) :: dimids + +integer :: nfinepoints ! number of fine grid points in each coarse grid cell + +real (kind = real_kind) :: sum2, sum4, var + + +real (kind = real_kind), allocatable :: & + zs(:,:) + +logical :: zs_accum + +real (kind = real_kind) :: zs_mean +real (kind = real_kind), allocatable :: & + std_dev(:),convexity(:), & + OA1(:),OA2(:),OA3(:),OA4(:), & + OL1(:),OL2(:),OL3(:),OL4(:) + +real (kind = real_kind), parameter :: max_convexity = 10._real_kind ! max value for convexity + +integer :: nu, nd, nw, nt +real (kind = real_kind) :: ratio + + +real, parameter :: ae = 6371220._real_kind ! Earth radius in meters + +character(len=35) :: FV3_grid_input_file_name +character(len=150) :: fine_topo_source_file_name +character(len=35) :: oro_data_output_file_name + +integer :: temp_int, dimX_FV3, dimY_FV3 +real (kind = dbl_kind), allocatable :: lat_FV3_raw(:,:), lon_FV3_raw(:,:) +real (kind = real_kind), allocatable :: lat_FV3(:,:), lon_FV3(:,:) +real (kind = real_kind), allocatable :: lat_FV3_deg(:,:) ! saved version of latitude for output +real (kind = real_kind), allocatable :: lon_FV3_deg(:,:) ! saved version of longitude for output +real (kind = dbl_kind), allocatable :: area_FV3_qtr(:,:) ! meters squared +real (kind = real_kind), allocatable :: area_FV3(:,:) ! meters squared +real (kind = real_kind), allocatable :: HGT_M_fine(:,:) +real (kind = real_kind) :: dlta_lat, dlta_lon + +integer :: i_blk, j_blk +integer :: ii_loc, jj_loc, ii_m, jj_m +integer, dimension(3) :: s_ii, e_ii, s_jj, e_jj +real (kind = real_kind), dimension(3) :: lat_blk, lon_blk +integer :: cell_count ! allows for use of 1D arrays for GWD statistics fields +integer :: halo_int ! integer form of halo + +logical :: fexist + + +print *, "Creating oro_data_ss file" +print * + +! File name for file that contains grid information +if ( halo.eq."-999" ) then ! global or nested tile + FV3_grid_input_file_name = "C" // trim(res_indx) // "_grid.tile" // & + trim(tile_num) // ".nc" +else ! stand-alone regional tile + FV3_grid_input_file_name = "C" // trim(res_indx) // "_grid.tile" // & + trim(tile_num) // ".halo" // trim(halo) // ".nc" +end if + +print *, "Reading from file: ", FV3_grid_input_file_name + + +! Check that input file exists -- exit with error message if not +inquire(file=FV3_grid_input_file_name,exist=fexist) +if (.not.fexist) then + print * + print *, "Fatal error: Input file "//trim(FV3_grid_input_file_name)// & + " does not exist" + print *, "Exiting program gsl_oro_data.f90" + print * + call exit(4) +end if + + +! In preparation for reading in grid data, account for existence +! of halo points +if ( halo.eq."-999" ) then ! global or nested tile + halo_int = 0 +else + read(halo,*) halo_int ! integer form of halo +end if + + +! Open Cxxx_grid netCDF file for input and get dimensions +err = nf90_open(FV3_grid_input_file_name,nf90_nowrite,ncid_in) +call netcdf_err(err, 'opening: '//FV3_grid_input_file_name) + +err = nf90_inq_dimid(ncid_in,'nx',dimid) +call netcdf_err(err, 'reading nx id') +err = nf90_inquire_dimension(ncid_in,dimid,len=temp_int) +call netcdf_err(err, 'reading nx value') +dimX_FV3 = temp_int/2 - 2*halo_int ! shaving off any halo points from edges + +err = nf90_inq_dimid(ncid_in,'ny',dimid) +call netcdf_err(err, 'reading ny id') +err = nf90_inquire_dimension(ncid_in,dimid,len=temp_int) +call netcdf_err(err, 'reading ny value') +dimY_FV3 = temp_int/2 - 2*halo_int ! shaving off any halo points from edges + +print *, "dimX_FV3 =", dimX_FV3 ! number of model cells in x-direction (halo0) +print *, "dimY_FV3 =", dimY_FV3 ! number of model cells in y-direction (halo0) +print * + +! Read in lat/lon (in degrees) +allocate (lat_FV3_raw((2*dimX_FV3+1),(2*dimY_FV3+1))) +err = nf90_inq_varid(ncid_in,'y',varid) +call netcdf_err(err, 'reading y id') +err = nf90_get_var(ncid_in,varid,lat_FV3_raw, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3+1,2*dimY_FV3+1/)) +call netcdf_err(err, 'reading y') + +allocate (lon_FV3_raw((2*dimX_FV3+1),(2*dimY_FV3+1))) +err = nf90_inq_varid(ncid_in,'x',varid) +call netcdf_err(err, 'reading x id') +err = nf90_get_var(ncid_in,varid,lon_FV3_raw, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3+1,2*dimY_FV3+1/)) +call netcdf_err(err, 'reading x') + +! Read in quarter grid-cell areas +allocate (area_FV3_qtr((2*dimX_FV3),(2*dimY_FV3))) +err = nf90_inq_varid(ncid_in,'area',varid) +call netcdf_err(err, 'reading area id') +err = nf90_get_var(ncid_in,varid,area_FV3_qtr, & + start=(/1+2*halo_int,1+2*halo_int/), & + count=(/2*dimX_FV3,2*dimY_FV3/)) +call netcdf_err(err, 'reading area') + +! Calculate lat/lon at mass points (cell-centers) +! Stride by 2 starting with 2nd point +! NOTE: "Converting" from dbl_kind to real_kind +allocate (lat_FV3(dimX_FV3,dimY_FV3)) +allocate (lon_FV3(dimX_FV3,dimY_FV3)) +allocate (lat_FV3_deg(dimX_FV3,dimY_FV3)) +allocate (lon_FV3_deg(dimX_FV3,dimY_FV3)) +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + lat_FV3(i,j) = lat_FV3_raw(2*i,2*j) + lon_FV3(i,j) = lon_FV3_raw(2*i,2*j) + end do +end do +lat_FV3_deg(:,:) = lat_FV3(:,:) +lon_FV3_deg(:,:) = lon_FV3(:,:) +deallocate(lat_FV3_raw) +deallocate(lon_FV3_raw) + +! Convert lat/lon to radians +lat_FV3 = lat_FV3*pi/180._real_kind +lon_FV3 = lon_FV3*pi/180._real_kind + +! Create full grid-cell areas (4 raw areas per grid cell area) +! NOTE: "Converting" from dbl_kind to real_kind +allocate(area_FV3(dimX_FV3,dimY_FV3)) +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + area_FV3(i,j) = area_FV3_qtr(2*i-1,2*j-1) + area_FV3_qtr(2*i-1,2*j) + & + area_FV3_qtr(2*i,2*j-1) + area_FV3_qtr(2*i,2*j) + end do +end do +deallocate(area_FV3_qtr) + +err = nf90_close(ncid_in) + + + +! Open file containing 30sec topo data (fine grid) +fine_topo_source_file_name = "HGT.Beljaars_filtered.lat-lon.30s_res.nc" +! Check that input file exists -- exit with error message if not +inquire(file=fine_topo_source_file_name,exist=fexist) +if (.not.fexist) then + print * + print *, "Fatal error: Topo source file "// & + trim(fine_topo_source_file_name)//" does not exist" + print *, "Exiting program gsl_oro_data.f90" + print * + call exit(4) +end if +err = nf90_open(trim(fine_topo_source_file_name),nf90_nowrite,ncid_in) +call netcdf_err(err, 'opening: '//trim(fine_topo_source_file_name)) + +! Get dimensions +err = nf90_inq_dimid(ncid_in,'west_east',dimid) +call netcdf_err(err, 'reading west_east id') +err = nf90_inquire_dimension(ncid_in,dimid,len=dimX_fine) +call netcdf_err(err, 'reading west_east value') + +err = nf90_inq_dimid(ncid_in,'south_north',dimid) +call netcdf_err(err, 'reading south_north id') +err = nf90_inquire_dimension(ncid_in,dimid,len=dimY_fine) +call netcdf_err(err, 'reading south_north value') + +print *, "Source file for high-resolution topography: ", & + trim(fine_topo_source_file_name) +print *, "dimX_fine =", dimX_fine +print *, "dimY_fine =", dimY_fine +print * + + +! Read in lat/lon of fine grid +allocate(lat1d_fine(dimY_fine)) +allocate(lon1d_fine(dimX_fine)) +err = nf90_inq_varid(ncid_in,'CLAT',varid) +call netcdf_err(err, 'reading CLAT id') +err = nf90_get_var(ncid_in,varid,lat1d_fine,start=(/1,1/), & + count=(/1,dimY_fine/)) +call netcdf_err(err, 'reading CLAT') + +err = nf90_inq_varid(ncid_in,'CLONG',varid) +call netcdf_err(err, 'reading CLONG id') +err = nf90_get_var(ncid_in,varid,lon1d_fine,start=(/1,1/), & + count=(/dimX_fine,1/)) +call netcdf_err(err, 'reading CLONG') + +! Convert lat/lon to radians +lat1d_fine = lat1d_fine*pi/180._real_kind +lon1d_fine = lon1d_fine*pi/180._real_kind + + +! Reassign FV3 longitude to vary from -pi to pi to match lon1d_fine range +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + if ( lon_FV3(i,j).gt.pi ) then + lon_FV3(i,j) = lon_FV3(i,j) - 2*pi + end if + end do +end do + + +! Read in fine-scale topography +allocate(HGT_M_fine(dimX_fine,dimY_fine)) +err = nf90_inq_varid(ncid_in,'HGT_M',varid) +call netcdf_err(err, 'reading HGT_M id') +err = nf90_get_var(ncid_in,varid,HGT_M_fine,start=(/1,1/), & + count=(/dimX_fine,dimY_fine/)) +call netcdf_err(err, 'reading HGT_M') + + +err = nf90_close(ncid_in) + + +! Allocate GWD statistics fields +allocate (std_dev(dimX_FV3*dimY_FV3)) +allocate (convexity(dimX_FV3*dimY_FV3)) +allocate (OA1(dimX_FV3*dimY_FV3)) +allocate (OA2(dimX_FV3*dimY_FV3)) +allocate (OA3(dimX_FV3*dimY_FV3)) +allocate (OA4(dimX_FV3*dimY_FV3)) +allocate (OL1(dimX_FV3*dimY_FV3)) +allocate (OL2(dimX_FV3*dimY_FV3)) +allocate (OL3(dimX_FV3*dimY_FV3)) +allocate (OL4(dimX_FV3*dimY_FV3)) + +! Initialize GWD statistics fields +std_dev(:) = 0._real_kind +convexity(:) = 0._real_kind +OA1(:) = 0._real_kind +OA2(:) = 0._real_kind +OA3(:) = 0._real_kind +OA4(:) = 0._real_kind +OL1(:) = 0._real_kind +OL2(:) = 0._real_kind +OL3(:) = 0._real_kind +OL4(:) = 0._real_kind + + + +! Calculate the minimum coarse grid cell size as implied by the cell area +min_area_FV3 = 1.E+14 +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + min_area_FV3 = min(min_area_FV3,area_FV3(i,j)) + end do +end do +! The square root of min_area_FV3 will count as the minimum cell size +min_DX = sqrt(min_area_FV3)/1000._real_kind ! grid size in km +! NOTE: min_DX will be used after the big loop below to determine whether +! to copy topographic statistics to "large-scale" file + + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! This is a loop over all the FV3 (coarse) grid cells +! The subgrid-scale topographic variables needed for the large-scale +! orographic gravity wave drag schemes are calculated by the following steps: +! 1) Sample the fine-scale (30sec) topography contained within each +! coarse grid cell. +! 2) Calculate the orographic statistics: stddev,convexity,oa1,...oa4, +! ol1,...,ol4 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +cell_count = 1 + +do j = 1,dimY_FV3 + do i = 1,dimX_FV3 + + ! Calculate approximate side-lengths of square lat-long "coarse" grid + ! cell centered on FV3 cell (units = radians) + dlta_lat = sqrt(area_FV3(i,j))/ae + dlta_lon = sqrt(area_FV3(i,j))/(ae*COS(lat_FV3(i,j))) + + ! Determine lat/lon of 9 lat-lon block centers + ! Note: lat_blk(2)/lon_blk(2) = lat_FV3(i,j)/lon_FV3(i,j) + ! Note: abs(lon_blk) may exceed pi + do i_blk = 1,3 + lon_blk(i_blk) = lon_FV3(i,j) + (i_blk-2)*dlta_lon + end do + ! Note: abs(lat_blk) may exceed pi/2 (90 degrees) + do j_blk = 1,3 + lat_blk(j_blk) = lat_FV3(i,j) + (j_blk-2)*dlta_lat + end do + + ! Find starting and ending fine-grid i,j indices for each + ! of the 9 "coarse-grid" blocks + ! Note: Index value of -999 is returned if latitude of grid points + ! exceed 90 degrees north or south + do i_blk = 1,3 + s_ii(i_blk) = nearest_i_east(lon_blk(i_blk)-p5*dlta_lon) + e_ii(i_blk) = nearest_i_west(lon_blk(i_blk)+p5*dlta_lon) + end do + do j_blk = 1,3 + s_jj(j_blk) = nearest_j_north(lat_blk(j_blk)-p5*dlta_lat) + e_jj(j_blk) = nearest_j_south(lat_blk(j_blk)+p5*dlta_lat) + end do + + + ! Calculate lat/lon relevant to each "coarse grid" block + do i_blk = 1,3 + + ! "Shave" blocks on north or south due to proximity to poles + ! if necessary + j_blk = 1 ! southern row + ! Check for "shaved" block due to proximity to south pole + if ( (s_jj(j_blk).eq.-999).and.(e_jj(j_blk).ne.-999) ) then + s_jj(j_blk) = 1 ! southern boundary of shaved block + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(1)+lat1d_fine(e_jj(j_blk))) + end if + + j_blk = 2 ! center row + ! Check for "shaved" block due to proximity to south or north pole + ! Note: We're assuming e_jj(2) and s_jj(2) can't both be -999 + if ( s_jj(j_blk).eq.-999 ) then + s_jj(j_blk) = 1 ! block shaved on the south + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(1)+lat1d_fine(e_jj(j_blk))) + end if + if ( e_jj(j_blk).eq.-999 ) then + e_jj(j_blk) = dimY_fine ! block shaved on the north + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(s_jj(j_blk))+lat1d_fine(dimY_fine)) + end if + + j_blk = 3 ! northern row + ! Check for "shaved" block due to proximity to north pole + if ( (e_jj(j_blk).eq.-999).and.(s_jj(j_blk).ne.-999) ) then + e_jj(j_blk) = dimY_fine ! northern boundary of shaved block + ! Reassign latitude of block center + lat_blk(j_blk) = p5*(lat1d_fine(s_jj(j_blk))+lat1d_fine(dimY_fine)) + end if + + end do + + + ! Calculate number of fine-grid points within center coarse block (2,2) + ! Check if center block straddles date line + if ( s_ii(2).gt.e_ii(2) ) then + ii_m = dimX_fine - s_ii(2) + 1 + e_ii(2) + else + ii_m = e_ii(2) - s_ii(2) + 1 + end if + jj_m = e_jj(2) - s_jj(2) + 1 + + + ! Assign values to "zs", which is the fine-grid surface topography field + ! that we will calculate statistics on, i.e, stddev, convexity, etc. + allocate (zs(ii_m,jj_m)) + + do jj = s_jj(2), e_jj(2) + jj_loc = jj - s_jj(2) + 1 ! local j-index (1 ... jj_m) + ! Check if block straddles the date line + if ( s_ii(2).gt.e_ii(2) ) then + do ii = s_ii(2), dimX_fine ! west of the date line + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end do + do ii = 1, e_ii(2) ! east of the date line + ii_loc = ii_loc + 1 ! local i-index ( 1 ... ii_m ) + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end do + else ! no crossing of the date line + do ii = s_ii(2), e_ii(2) + ii_loc = ii - s_ii(2) + 1 ! local i-index ( 1 ... ii_m) + zs(ii_loc,jj_loc) = HGT_M_fine(ii,jj) + end do + end if + end do + + + + ! + ! Finally, we can now calculate the topographic statistics fields needed + ! for the gravity wave drag scheme + ! + + ! Make sure statistics are zero if there is no terrain in the grid cell + ! Note: This is a proxy for a landmask + zs_accum = .false. + do jj = 1,jj_m + do ii = 1,ii_m + if ( abs(zs(ii,jj)).gt.1.E-1 ) zs_accum = .true. + end do + end do + if ( .not.zs_accum ) then ! no terrain in the grid cell + std_dev(cell_count) = 0._real_kind + convexity(cell_count) = 0._real_kind + OA1(cell_count) = 0._real_kind + OA2(cell_count) = 0._real_kind + OA3(cell_count) = 0._real_kind + OA4(cell_count) = 0._real_kind + OL1(cell_count) = 0._real_kind + OL2(cell_count) = 0._real_kind + OL3(cell_count) = 0._real_kind + OL4(cell_count) = 0._real_kind + deallocate(zs) + cell_count = cell_count + 1 + cycle ! move on to next (coarse) grid cell + end if + + + ! + ! Calculate standard deviation of subgrid-scale terrain height + ! + + ! Calculate mean height + sum2 = 0._real_kind + nfinepoints = ii_m*jj_m + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + zs(ii,jj) + end do + end do + zs_mean = sum2 / real(nfinepoints,real_kind) + + ! Calculate standard deviation + sum2 = 0._real_kind + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + ( zs(ii,jj) - zs_mean )**2 + end do + end do + std_dev(cell_count) = sqrt( sum2/real(nfinepoints,real_kind) ) + + + ! + ! Calculate convexity of sub-grid-scale terrain + ! + + sum2 = 0._real_kind + sum4 = 0._real_kind + do jj = 1,jj_m + do ii = 1,ii_m + sum2 = sum2 + ( zs(ii,jj) - zs_mean )**2 + sum4 = sum4 + ( zs(ii,jj) - zs_mean )**4 + end do + end do + + var = sum2 / real(nfinepoints,real_kind) + if ( abs(var) < 1.0E-05_real_kind ) then + convexity(cell_count) = 0._real_kind + else + convexity(cell_count) = min( sum4 / ( var**2 * & + real(nfinepoints,real_kind) ), max_convexity ) + end if + + + ! + ! Calculate orographic asymmetries + ! + + ! OA1 -- orographic asymmetry in West direction + nu = 0 + nd = 0 + do jj = 1,jj_m + do ii = 1,ii_m/2 ! left half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + end do + do ii = ii_m/2 + 1, ii_m ! right half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end do + end do + if ( nu + nd > 0 ) then + OA1(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA1(cell_count) = 0._real_kind + end if + + ! OA2 -- orographic asymmetry in South direction + nu = 0 + nd = 0 + do jj = 1,jj_m/2 ! bottom half of box + do ii = 1,ii_m + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + end do + end do + do jj = jj_m/2 + 1,jj_m ! top half of box + do ii = 1, ii_m + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end do + end do + if ( nu + nd > 0 ) then + OA2(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA2(cell_count) = 0._real_kind + end if + + ! OA3 -- orographic asymmetry in South-West direction + nu = 0 + nd = 0 + ratio = real(jj_m,real_kind)/real(ii_m,real_kind) + do jj = 1,jj_m + do ii = 1,ii_m + if ( nint(real(ii,real_kind)*ratio) < (jj_m - jj) ) then + ! south-west half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + else ! north-east half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end if + end do + end do + if ( nu + nd > 0 ) then + OA3(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA3(cell_count) = 0._real_kind + end if + + ! OA4 -- orographic asymmetry in North-West direction + nu = 0 + nd = 0 + ratio = real(jj_m,real_kind)/real(ii_m,real_kind) + do jj = 1,jj_m + do ii = 1,ii_m + if ( nint(real(ii,real_kind)*ratio) < jj ) then + ! north-west half of box + if ( zs(ii,jj) > zs_mean ) nu = nu + 1 + else ! south-east half of box + if ( zs(ii,jj) > zs_mean ) nd = nd + 1 + end if + end do + end do + if ( nu + nd > 0 ) then + OA4(cell_count) = real((nu - nd),real_kind) / & + real((nu + nd),real_kind) + else + OA4(cell_count) = 0._real_kind + end if + + + ! + ! Calculate orographic effective lengths + ! + + ! OL1 -- orographic effective length for Westerly flow + nw = 0 + nt = 0 + do jj = max(jj_m/4,1), 3*jj_m/4 + ! within central east-west band of box + do ii = 1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL1(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL1(cell_count) = 0._real_kind + end if + + ! OL2 -- orographic effective length for Southerly flow + nw = 0 + nt = 0 + do jj = 1, jj_m + do ii = max(ii_m/4,1), 3*ii_m/4 + ! within central north-south band of box + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL2(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL2(cell_count) = 0._real_kind + end if + + ! OL3 -- orographic effective length for South-Westerly flow + nw = 0 + nt = 0 + do jj = 1, jj_m/2 + do ii = 1, ii_m/2 + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + do jj = jj_m/2+1, jj_m + do ii = ii_m/2+1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL3(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL3(cell_count) = 0._real_kind + end if + + ! OL4 -- orographic effective length for North-Westerly flow + nw = 0 + nt = 0 + do jj = jj_m/2+1, jj_m + do ii = 1, ii_m/2 + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + do jj = 1, jj_m/2 + do ii = ii_m/2+1, ii_m + if ( zs(ii,jj) > zs_mean ) nw = nw + 1 + nt = nt + 1 + end do + end do + if ( nt /= 0 ) then + OL4(cell_count) = real(nw,real_kind) / real(nt,real_kind) + else + OL4(cell_count) = 0._real_kind + end if + + + + deallocate (zs) + + cell_count = cell_count + 1 + + end do ! j = 1,dimY_FV3 +end do ! i = 1,dimX_FV3 + + + +! +! Output GWD statistics fields to netCDF file +! + + +if ( halo.eq."-999" ) then ! global or nested tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ss.tile" & + // trim(tile_num) // ".nc" +else ! stand-alone regional tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ss.tile" & + // trim(tile_num) // ".halo0.nc" +end if + +! Open netCDF file for output +err = nf90_create(oro_data_output_file_name, NF90_CLOBBER, ncid_out) +call netcdf_err(err, 'creating: '//oro_data_output_file_name) + +err = nf90_redef(ncid_out) + +! Define dimensions +err = nf90_def_dim(ncid_out,'lon',dimX_FV3,lonid) +call netcdf_err(err, 'defining lon dimension') +err = nf90_def_dim(ncid_out,'lat',dimY_FV3,latid) +call netcdf_err(err, 'defining lat dimension') + +! Define the 'dimensions vector' dimids to be used for writing +! the 2-dimensional variables to the netCDF file +dimids(1) = lonid +dimids(2) = latid + +! Define variables and attributes to put in the netCDF file +err = nf90_def_var(ncid_out,'geolon',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining geolon') +err = nf90_put_att(ncid_out,varid,'units','degrees') +err = nf90_put_att(ncid_out,varid,'description','longitude') +err = nf90_def_var(ncid_out,'geolat',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining geolat') +err = nf90_put_att(ncid_out,varid,'units','degrees') +err = nf90_put_att(ncid_out,varid,'description','latitude') +err = nf90_def_var(ncid_out,'stddev',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'stddev') +err = nf90_put_att(ncid_out,varid,'units','meters') +err = nf90_put_att(ncid_out,varid,'description', & + 'standard deviation of subgrid topography') +err = nf90_def_var(ncid_out,'convexity',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining convexity') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'convexity of subgrid topography') +err = nf90_def_var(ncid_out,'oa1',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa1') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in west direction') +err = nf90_def_var(ncid_out,'oa2',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa2') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south direction') +err = nf90_def_var(ncid_out,'oa3',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa3') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south-west direction') +err = nf90_def_var(ncid_out,'oa4',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining oa4') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in north-west direction') +err = nf90_def_var(ncid_out,'ol1',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol1') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for westerly flow') +err = nf90_def_var(ncid_out,'ol2',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol2') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for southerly flow') +err = nf90_def_var(ncid_out,'ol3',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol3') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for south-westerly flow') +err = nf90_def_var(ncid_out,'ol4',NF90_FLOAT,dimids,varid) +call netcdf_err(err, 'defining ol4') +err = nf90_put_att(ncid_out,varid,'units','-') +err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for north-westerly flow') + +! Add global attributes +err = nf90_put_att(ncid_out,nf90_global, & + 'source_file_for_high-resolution_topography', & + trim(fine_topo_source_file_name)) + +err = nf90_enddef(ncid_out) + + +! Write data to output netCDF file +err = nf90_inq_varid(ncid_out,'geolon',varid) +call netcdf_err(err, 'reading geolon id') +err = nf90_put_var(ncid_out,varid,lon_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing geolon') +err = nf90_inq_varid(ncid_out,'geolat',varid) +call netcdf_err(err, 'reading geolat id') +err = nf90_put_var(ncid_out,varid,lat_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing geolat') +err = nf90_inq_varid(ncid_out,'stddev',varid) +call netcdf_err(err, 'reading stddev id') +err = nf90_put_var(ncid_out,varid,std_dev,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing stddev') +err = nf90_inq_varid(ncid_out,'convexity',varid) +call netcdf_err(err, 'reading convexity id') +err = nf90_put_var(ncid_out,varid,convexity,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing convexity') +err = nf90_inq_varid(ncid_out,'oa1',varid) +call netcdf_err(err, 'reading oa1 id') +err = nf90_put_var(ncid_out,varid,OA1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa1') +err = nf90_inq_varid(ncid_out,'oa2',varid) +call netcdf_err(err, 'reading oa2 id') +err = nf90_put_var(ncid_out,varid,OA2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa2') +err = nf90_inq_varid(ncid_out,'oa3',varid) +call netcdf_err(err, 'reading oa3 id') +err = nf90_put_var(ncid_out,varid,OA3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa3') +err = nf90_inq_varid(ncid_out,'oa4',varid) +call netcdf_err(err, 'reading oa4 id') +err = nf90_put_var(ncid_out,varid,OA4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing oa4') +err = nf90_inq_varid(ncid_out,'ol1',varid) +call netcdf_err(err, 'reading ol1 id') +err = nf90_put_var(ncid_out,varid,OL1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol1') +err = nf90_inq_varid(ncid_out,'ol2',varid) +call netcdf_err(err, 'reading ol2 id') +err = nf90_put_var(ncid_out,varid,OL2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol2') +err = nf90_inq_varid(ncid_out,'ol3',varid) +call netcdf_err(err, 'reading ol3 id') +err = nf90_put_var(ncid_out,varid,OL3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol3') +err = nf90_inq_varid(ncid_out,'ol4',varid) +call netcdf_err(err, 'reading ol4 id') +err = nf90_put_var(ncid_out,varid,OL4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) +call netcdf_err(err, 'writing ol4') + +err = nf90_close(ncid_out) + + + +! Determine whether grid size falls below threshold for use by large-scale +! orographic gravity wave drag and blocking scheme. If it is, then +! create "dummy" oro_data_ls file containing the same data as oro_data_ss file + +duplicate_oro_data_file = .false. + +if ( min_DX.le.7.5 ) then + + duplicate_oro_data_file = .true. + print *, "Creating oro_data_ls file as duplicate of oro_data_ss" + print *, "Minimum grid cell size = ", min_DX, " km" + print * + + + if ( halo.eq."-999" ) then ! global or nested tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ls.tile" & + // trim(tile_num) // ".nc" + else ! stand-alone regional tile + oro_data_output_file_name = "C" // trim(res_indx) // "_oro_data_ls.tile" & + // trim(tile_num) // ".halo0.nc" + end if + + ! Open netCDF file for output + err = nf90_create(oro_data_output_file_name, NF90_CLOBBER, ncid_out) + call netcdf_err(err, 'creating: '//oro_data_output_file_name) + + err = nf90_redef(ncid_out) + + ! Define dimensions + err = nf90_def_dim(ncid_out,'lon',dimX_FV3,lonid) + call netcdf_err(err, 'defining lon dimension') + err = nf90_def_dim(ncid_out,'lat',dimY_FV3,latid) + call netcdf_err(err, 'defining lat dimension') + + ! Define the 'dimensions vector' dimids to be used for writing + ! the 2-dimensional variables to the netCDF file + dimids(1) = lonid + dimids(2) = latid + + ! Define variables and attributes to put in the netCDF file + err = nf90_def_var(ncid_out,'geolon',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining geolon') + err = nf90_put_att(ncid_out,varid,'units','degrees') + err = nf90_put_att(ncid_out,varid,'description','longitude') + err = nf90_def_var(ncid_out,'geolat',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining geolat') + err = nf90_put_att(ncid_out,varid,'units','degrees') + err = nf90_put_att(ncid_out,varid,'description','latitude') + err = nf90_def_var(ncid_out,'stddev',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining stddev') + err = nf90_put_att(ncid_out,varid,'units','meters') + err = nf90_put_att(ncid_out,varid,'description', & + 'standard deviation of subgrid topography') + err = nf90_def_var(ncid_out,'convexity',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining convexity') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'convexity of subgrid topography') + err = nf90_def_var(ncid_out,'oa1',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining oa1') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in west direction') + err = nf90_def_var(ncid_out,'oa2',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining oa2') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south direction') + err = nf90_def_var(ncid_out,'oa3',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining oa3') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in south-west direction') + err = nf90_def_var(ncid_out,'oa4',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining oa4') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic asymmetry in north-west direction') + err = nf90_def_var(ncid_out,'ol1',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining ol1') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for westerly flow') + err = nf90_def_var(ncid_out,'ol2',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining ol2') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for southerly flow') + err = nf90_def_var(ncid_out,'ol3',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining ol3') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for south-westerly flow') + err = nf90_def_var(ncid_out,'ol4',NF90_FLOAT,dimids,varid) + call netcdf_err(err, 'defining ol4') + err = nf90_put_att(ncid_out,varid,'units','-') + err = nf90_put_att(ncid_out,varid,'description', & + 'orographic effective length for north-westerly flow') + + ! Add global attributes + err = nf90_put_att(ncid_out,nf90_global, & + 'NOTE','This is a duplicate of the oro_data_ss file') + + err = nf90_enddef(ncid_out) + + + ! Write data to output netCDF file + err = nf90_inq_varid(ncid_out,'geolon',varid) + call netcdf_err(err, 'reading geolon id') + err = nf90_put_var(ncid_out,varid,lon_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing geolon') + err = nf90_inq_varid(ncid_out,'geolat',varid) + call netcdf_err(err, 'reading geolat id') + err = nf90_put_var(ncid_out,varid,lat_FV3_deg,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing geolat') + err = nf90_inq_varid(ncid_out,'stddev',varid) + call netcdf_err(err, 'reading stddev id') + err = nf90_put_var(ncid_out,varid,std_dev,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing stddev') + err = nf90_inq_varid(ncid_out,'convexity',varid) + call netcdf_err(err, 'reading convexity id') + err = nf90_put_var(ncid_out,varid,convexity,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing convexity') + err = nf90_inq_varid(ncid_out,'oa1',varid) + call netcdf_err(err, 'reading oa1 id') + err = nf90_put_var(ncid_out,varid,OA1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing oa1') + err = nf90_inq_varid(ncid_out,'oa2',varid) + call netcdf_err(err, 'reading oa2 id') + err = nf90_put_var(ncid_out,varid,OA2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing oa2') + err = nf90_inq_varid(ncid_out,'oa3',varid) + call netcdf_err(err, 'reading oa3 id') + err = nf90_put_var(ncid_out,varid,OA3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing oa3') + err = nf90_inq_varid(ncid_out,'oa4',varid) + call netcdf_err(err, 'reading oa4 id') + err = nf90_put_var(ncid_out,varid,OA4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing oa4') + err = nf90_inq_varid(ncid_out,'ol1',varid) + call netcdf_err(err, 'reading ol1 id') + err = nf90_put_var(ncid_out,varid,OL1,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing ol1') + err = nf90_inq_varid(ncid_out,'ol2',varid) + call netcdf_err(err, 'reading ol2 id') + err = nf90_put_var(ncid_out,varid,OL2,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing ol2') + err = nf90_inq_varid(ncid_out,'ol3',varid) + call netcdf_err(err, 'reading ol3 id') + err = nf90_put_var(ncid_out,varid,OL3,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing ol3') + err = nf90_inq_varid(ncid_out,'ol4',varid) + call netcdf_err(err, 'reading ol4 id') + err = nf90_put_var(ncid_out,varid,OL4,start=(/1,1/), & + count=(/dimX_FV3,dimY_FV3/)) + call netcdf_err(err, 'writing ol4') + + err = nf90_close(ncid_out) + +else + + print *, "Minimum grid cell size = ", min_DX, " km" + print * + +end if + + + +! Deallocate arrays +deallocate(lat_FV3) +deallocate(lon_FV3) +deallocate(lat_FV3_deg) +deallocate(lon_FV3_deg) +deallocate(area_FV3) +deallocate(lat1d_fine) +deallocate(lon1d_fine) +deallocate(HGT_M_fine) +deallocate(std_dev) +deallocate(convexity) +deallocate(OA1) +deallocate(OA2) +deallocate(OA3) +deallocate(OA4) +deallocate(OL1) +deallocate(OL2) +deallocate(OL3) +deallocate(OL4) + + +end subroutine calc_gsl_oro_data_sm_scale + +!> Finds nearest fine-grid i index to the east of a given longitude +!! +!! @param[in] lon_in longitude (radians) +!! @return nearest_i_east Nearest grid point i-index east of selected point +!! @author Michael Toy, NOAA/GSL +function nearest_i_east(lon_in) +! Calculates nearest fine-grid i index to the east of (or on) a given longitude +implicit none + +integer :: nearest_i_east +real (kind=real_kind), intent(in) :: lon_in +real (kind=real_kind) :: lon +integer :: i + +lon = lon_in +! Make sure longitude is between -pi and pi +do while ( (lon.lt.(-pi)).or.(lon.gt.pi) ) + if ( lon.lt.(-pi) ) lon = lon + 2*pi + if ( lon.gt.pi ) lon = lon - 2*pi +end do + +if ( lon.gt.lon1d_fine(dimX_fine) ) then + nearest_i_east = 1 +else + i = 1 + do while ( lon1d_fine(i).lt.lon ) + i = i + 1 + end do + nearest_i_east = i +end if + +end function nearest_i_east + +!> Finds nearest fine-grid i index to the west of a given longitude +!! +!! @param[in] lon_in longitude (radians) +!! @return nearest_i_west Nearest grid point i-index west of selected point +!! @author Michael Toy, NOAA/GSL +function nearest_i_west(lon_in) +! Calculates nearest fine-grid i index to the west of a given longitude +implicit none + +integer :: nearest_i_west +real (kind=real_kind), intent(in) :: lon_in +real (kind=real_kind) :: lon +integer :: i + +lon = lon_in +! Make sure longitude is between -pi and pi +do while ( (lon.lt.(-pi)).or.(lon.gt.pi) ) + if ( lon.lt.(-pi) ) lon = lon + 2*pi + if ( lon.gt.pi ) lon = lon - 2*pi +end do + +if ( (lon.lt.lon1d_fine(1)).or.(lon.ge.lon1d_fine(dimX_fine)) ) then + nearest_i_west = dimX_fine +else + i = 1 + do while ( lon1d_fine(i).le.lon ) + i = i + 1 + end do + nearest_i_west = i - 1 +end if + +end function nearest_i_west + +!> Calculates nearest fine-grid j index to the north of a given latitude +!! +!! @param[in] lat_in Latitude (radians) +!! @return nearest_j_north Nearest fine-grid j index to the north of a given latitude +!! @author Michael Toy, NOAA/GSL +function nearest_j_north(lat_in) +! Calculates nearest fine-grid j index to the north of a given latitude +! Note: If the abs(latitude) is greater than pi/2 (90 degrees) then +! the value -999 is returned +implicit none + +integer :: nearest_j_north +real (kind=real_kind), intent(in) :: lat_in +real (kind=real_kind) :: lat +integer :: j + +lat = lat_in +if ( abs(lat_in).gt.p5*pi ) then + nearest_j_north = -999 +else + j = 1 + do while ( (lat1d_fine(j).lt.lat).and.(j.lt.dimY_fine) ) + j = j + 1 + end do + nearest_j_north = j +end if + +end function nearest_j_north + +!> Calculates nearest fine-grid j index to the south of a given latitude +!! +!! @param[in] lat_in Latitude (radians) +!! @return nearest_j_south Nearest fine-grid j index to the south of a given latitude +!! @author Michael Toy, NOAA/GSL +function nearest_j_south(lat_in) +! Calculates nearest fine-grid j index to the south of a given latitude +! Note: If the abs(latitude) is greater than pi/2 (90 degrees) then +! the value -999 is returned +implicit none + +integer :: nearest_j_south +real (kind=real_kind), intent(in) :: lat_in +real (kind=real_kind) :: lat +integer :: j + +lat = lat_in +if ( abs(lat_in).gt.p5*pi ) then + nearest_j_south = -999 +elseif ( lat_in.le.lat1d_fine(1) ) then + nearest_j_south = 1 +else + j = 2 + do while ( (lat1d_fine(j).le.lat).and.(j.le.dimY_fine) ) + j = j + 1 + end do + nearest_j_south = j - 1 +end if + +end function nearest_j_south + +!> Interpolates (or extrapolates) linear function y = y(x) +!! +!! @param[in] x Input "x" value +!! @param[in] x1 Known point 1 +!! @param[in] x2 Known point 2 +!! @param[in] y1 Known y(x1) +!! @param[in] y2 Known y(x2) +!! @return interp_1d Interpolated y value at x +!! @author Michael Toy, NOAA/GSL +function interp_1d(x,x1,x2,y1,y2) +! Interpolates (or extrapolates) linear function y = y(x) +! to x given y1 = y(x1) and y2 = y(x2) +implicit none + +real (kind=real_kind) :: interp_1d +real (kind=real_kind), intent(in) :: x,x1,x2,y1,y2 +real (kind=real_kind) :: slope + +! Formula for a line: y = y1 + slope*(x - x1) +slope = (y2-y1)/(x2-x1) +interp_1d = y1 + slope*(x-x1) + +end function interp_1d + +!> Returns netCDF error given input err code +!! +!! @param[in] err Error code from netCDF routine +!! @param[in] string Portion of error message +!! @author Michael Toy, NOAA/GSL +subroutine netcdf_err(err,string) + +use netcdf + +implicit none + +integer, intent(in) :: err +character(len=*), intent(in) :: string +character(len=256) :: errmsg + +if (err.eq.NF90_NOERR ) return +errmsg = NF90_STRERROR(err) +print *, "" +print *, "FATAL ERROR: ", trim(string), ": ", trim(errmsg) +print *, "STOP." +call exit(4) + +return +end subroutine netcdf_err + + + +end module gsl_oro_data_sm_scale diff --git a/ush/fv3gfs_driver_grid.sh b/ush/fv3gfs_driver_grid.sh index 209f8de05..e5c898e7e 100755 --- a/ush/fv3gfs_driver_grid.sh +++ b/ush/fv3gfs_driver_grid.sh @@ -15,14 +15,18 @@ # records that describe the model grid. # 2) 'oro' files containing land mask, terrain and gravity # wave drag fields. -# 3) surface climo fields, such as soil type, vegetation +# 3) 'oro' files ('oro_data_ls' and 'oro_data_ss') specific to +# the GSL drag suite physics parameterization (only if +# flag make_gsl_orog = true) +# 4) surface climo fields, such as soil type, vegetation # greenness and albedo. # # Calls the following scripts: # 1) fv3gfs_make_grid.sh (make 'grid' files) -# 2) fv3gfs_mask_orog.sh (make 'oro' files) -# 3) fv3gfs_filter_topo.sh (filter topography) -# 4) sfc_climo_gen.sh (create surface climo fields) +# 2) fv3gfs_make_orog.sh (make 'oro' files) +# 3) fv3gfs_make_orog_gsl.sh (make gsl drag 'oro' files) +# 4) fv3gfs_filter_topo.sh (filter topography) +# 5) sfc_climo_gen.sh (create surface climo fields) # # Note: The sfc_climo_gen program only runs with an # mpi task count that is a multiple of six. This is @@ -49,6 +53,8 @@ export add_lake=${add_lake:-false} # add lake fraction and depth. uniform export lake_cutoff=${lake_cutoff:-0.20} # lake fractions less than lake_cutoff # are ignored. +export make_gsl_orog=${make_gsl_orog:-false} # when true, create GSL drag suite orog files. + if [ $gtype = uniform ]; then echo "Creating global uniform grid" elif [ $gtype = stretch ]; then @@ -99,6 +105,7 @@ export home_dir=${home_dir:-"$PWD/../"} export script_dir=$home_dir/ush export exec_dir=${exec_dir:-"$home_dir/exec"} export topo=$home_dir/fix/fix_orog +export topo_am=$home_dir/fix/fix_am export NCDUMP=${NCDUMP:-ncdump} @@ -170,9 +177,16 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then if [ $machine = WCOSS_C ]; then touch $TEMP_DIR/orog.file1 + if [ $make_gsl_orog = true ]; then + touch $TEMP_DIR/orog_gsl.file1 + fi tile=1 while [ $tile -le $ntiles ]; do echo "$script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo " >>$TEMP_DIR/orog.file1 + if [ $make_gsl_orog = true ]; then + export halo_tmp="-999" # no halo + echo $script_dir/fv3gfs_make_orog_gsl.sh $res $tile $halo_tmp $grid_dir $orog_dir $topo_am >>$TEMP_DIR/orog_gsl.file1 + fi tile=$(( $tile + 1 )) done aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TEMP_DIR/orog.file1 @@ -181,6 +195,14 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then exit $err fi rm $TEMP_DIR/orog.file1 + if [ $make_gsl_orog = true ]; then + aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TEMP_DIR/orog_gsl.file1 + err=$? + if [ $err != 0 ]; then + exit $err + fi + rm $TEMP_DIR/orog_gsl.file1 + fi else tile=1 while [ $tile -le $ntiles ]; do @@ -194,6 +216,19 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then if [ $err != 0 ]; then exit $err fi + if [ $make_gsl_orog = true ]; then + set +x + echo + echo "............ Execute fv3gfs_make_orog_gsl.sh for tile $tile .................." + echo + set -x + export halo_tmp="-999" # no halo + $script_dir/fv3gfs_make_orog_gsl.sh $res $tile $halo_tmp $grid_dir $orog_dir $topo_am + err=$? + if [ $err != 0 ]; then + exit $err + fi + fi tile=$(( $tile + 1 )) done fi @@ -235,6 +270,9 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then while [ $tile -le $ntiles ]; do cp $filter_dir/oro.C${res}.tile${tile}.nc $out_dir/C${res}_oro_data.tile${tile}.nc cp $grid_dir/C${res}_grid.tile${tile}.nc $out_dir/C${res}_grid.tile${tile}.nc + if [ $make_gsl_orog = true ]; then + cp $orog_dir/C${res}_oro_data_*.tile${tile}*.nc $out_dir/ # gsl drag suite oro_data files + fi tile=`expr $tile + 1 ` done @@ -450,6 +488,36 @@ elif [ $gtype = regional_gfdl ] || [ $gtype = regional_esg ]; then cp $grid_dir/C${res}_*mosaic.nc $out_dir + +#---------------------------------------------------------------------------------- +# Now that C${res}_grid.tile${tile}.halo0.nc has been created, we can use it +# to generate gsl drag suite oro_data files, which are generated only for halo0 +# Note: This is carried out only if $make_gsl_orog = true +#---------------------------------------------------------------------------------- + + if [ $make_gsl_orog = true ]; then + export halo_tmp="0" + ln -sf $out_dir/C${res}_grid.tile${tile}.halo0.nc $grid_dir/ + if [ $machine = WCOSS_C ]; then + echo $script_dir/fv3gfs_make_orog_gsl.sh $res $tile $halo_tmp $grid_dir $orog_dir $topo_am >>$TEMP_DIR/orog_gsl.file1 + aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TEMP_DIR/orog_gsl.file1 + err=$? + rm $TEMP_DIR/orog_gsl.file1 + else + set +x + echo + echo "............ Execute fv3gfs_make_orog_gsl.sh for tile $tile .................." + echo + set -x + $script_dir/fv3gfs_make_orog_gsl.sh $res $tile $halo_tmp $grid_dir $orog_dir $topo_am + err=$? + if [ $err != 0 ]; then + exit $err + fi + fi + cp $orog_dir/C${res}_oro_data_*.tile${tile}*.nc $out_dir/ # gsl drag suite oro_data files + fi + echo "Grid and orography files are now prepared for regional grid" fi diff --git a/ush/fv3gfs_make_orog_gsl.sh b/ush/fv3gfs_make_orog_gsl.sh new file mode 100755 index 000000000..40dfee059 --- /dev/null +++ b/ush/fv3gfs_make_orog_gsl.sh @@ -0,0 +1,70 @@ +#!/bin/bash +# +#----------------------------------------------------------------------- +# Script to run 'orog_gsl' executable, which generates 'oro_data' +# static topographic statistics files needed for GSL orographic +# drag suite. Source code is gsl_oro_data.f90. +#----------------------------------------------------------------------- +# + +set -eux + +nargv=$# + +if [ $nargv -ne 6 ]; then + echo "Number of arguments must be 6" + echo "Usage: $0 resolution tile griddir outdir script_dir hist_dir" + exit 1 +fi + +res=$1 +tile=$2 +halo=$3 +griddir=$4 +outdir=$5 +topo_am=$6 +workdir=$TEMP_DIR/C${res}/orog/tile$tile + + +executable=$exec_dir/orog_gsl +if [ ! -s $executable ]; then + echo "executable does not exist" + exit 1 +fi + +if [ ! -s $workdir ]; then mkdir -p $workdir ;fi +if [ ! -s $outdir ]; then mkdir -p $outdir ;fi + + +if [ $halo -eq "-999" ]; then + OUTGRID="C${res}_grid.tile${tile}.nc" +else + OUTGRID="C${res}_grid.tile${tile}.halo${halo}.nc" +fi + + +cd $workdir + +ln -sf ${griddir}/$OUTGRID . +ln -sf ${topo_am}/"HGT.Beljaars_filtered.lat-lon.30s_res.nc" . +ln -sf ${topo_am}/"geo_em.d01.lat-lon.2.5m.HGT_M.nc" . + +cp $executable . + +echo $tile > grid_info.dat +echo $res >> grid_info.dat +echo $halo >> grid_info.dat +time $executable < grid_info.dat + +if [ $? -ne 0 ]; then + echo "ERROR in running $executable " + exit 1 +else + mv ./C*oro_data_*.nc $outdir/ + echo "*oro_data_ls* and *oro_data_ss* files created" + echo "Successfully running $executable " + exit 0 +fi + +exit +