From d4f0db20cabfd22a143557fe8e400687f9010d92 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Fri, 17 Jan 2025 08:38:06 -0500 Subject: [PATCH] Update orography code for low-res grids (#1012) Fix logic that caused crash for very low-resolution grids. Update I/O routines to use F90 NetCDF calls. Add unit tests for 20 subroutines/functions. For complete details see #1000 and #1012. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 382 ++++++++++-------- .../orog.fd/mtnlm7_oclsm.F90 | 50 ++- .../orog_mask_tools.fd/orog.fd/orog_utils.F90 | 95 ++--- tests/orog/CMakeLists.txt | 86 +++- tests/orog/data/C12_grid.tile1.nc | Bin 0 -> 51102 bytes tests/orog/data/C48.mx500.tile1.nc | Bin 0 -> 55384 bytes tests/orog/data/landcover.umd.lowres.nc | Bin 0 -> 13774 bytes .../data/topography.antarctica.ramp.lowres.nc | Bin 0 -> 13633 bytes .../orog/data/topography.gmted2010.lowres.nc | Bin 0 -> 13565 bytes tests/orog/ftst_find_nearest_pole_pts.F90 | 156 +++++++ tests/orog/ftst_find_poles.F90 | 72 ++++ tests/orog/ftst_get_index.F90 | 123 ++++++ tests/orog/ftst_get_xnsum.F90 | 125 ++++++ tests/orog/ftst_get_xnsum2.F90 | 74 ++++ tests/orog/ftst_get_xnsum3.F90 | 76 ++++ tests/orog/ftst_inside_polygon.F90 | 223 ++++++++++ tests/orog/ftst_qc_orog_by_ramp.F90 | 98 +++++ tests/orog/ftst_read_global_mask.F90 | 38 ++ tests/orog/ftst_read_global_orog.F90 | 38 ++ tests/orog/ftst_read_mask.F90 | 82 ++++ tests/orog/ftst_read_mdl_dims.F90 | 32 ++ tests/orog/ftst_read_mdl_grid_file.F90 | 77 ++++ tests/orog/ftst_rm_isolated_pts.F90 | 217 ++++++++++ tests/orog/ftst_transpose.F90 | 90 +++++ util/gdas_init/set_fixed_files.sh | 6 + 25 files changed, 1886 insertions(+), 254 deletions(-) create mode 100644 tests/orog/data/C12_grid.tile1.nc create mode 100644 tests/orog/data/C48.mx500.tile1.nc create mode 100644 tests/orog/data/landcover.umd.lowres.nc create mode 100644 tests/orog/data/topography.antarctica.ramp.lowres.nc create mode 100644 tests/orog/data/topography.gmted2010.lowres.nc create mode 100644 tests/orog/ftst_find_nearest_pole_pts.F90 create mode 100644 tests/orog/ftst_find_poles.F90 create mode 100644 tests/orog/ftst_get_index.F90 create mode 100644 tests/orog/ftst_get_xnsum.F90 create mode 100644 tests/orog/ftst_get_xnsum2.F90 create mode 100644 tests/orog/ftst_get_xnsum3.F90 create mode 100644 tests/orog/ftst_inside_polygon.F90 create mode 100644 tests/orog/ftst_qc_orog_by_ramp.F90 create mode 100644 tests/orog/ftst_read_global_mask.F90 create mode 100644 tests/orog/ftst_read_global_orog.F90 create mode 100644 tests/orog/ftst_read_mask.F90 create mode 100644 tests/orog/ftst_read_mdl_dims.F90 create mode 100644 tests/orog/ftst_read_mdl_grid_file.F90 create mode 100644 tests/orog/ftst_rm_isolated_pts.F90 create mode 100644 tests/orog/ftst_transpose.F90 diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 51a646779..7ccb86765 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -39,6 +39,7 @@ module io_utils !! @param[in] lat Latitude of the first column of the model grid tile. !! @author Jordan Alpert NOAA/EMC GFDL Programmer subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolon, geolat, lon, lat) + use netcdf implicit none integer, intent(in):: im, jm, ntiles, tile real, intent(in) :: lon(im), lat(jm) @@ -56,7 +57,6 @@ subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolo integer :: id_oa1,id_oa2,id_oa3,id_oa4 integer :: id_ol1,id_ol2,id_ol3,id_ol4 integer :: id_theta,id_gamma,id_sigma,id_elvmax - include "netcdf.inc" if(ntiles > 1) then write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc' @@ -68,164 +68,165 @@ subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolo dim2=size(lat,1) !--- open the file - error = NF__CREATE(outfile, IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid) + error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), ncid, & + initialsize=inital, chunksize=fsize) call netcdf_err(error, 'Creating file '//trim(outfile) ) !--- define dimension - error = nf_def_dim(ncid, 'lon', dim1, dim_lon) + error = nf90_def_dim(ncid, 'lon', dim1, dim_lon) call netcdf_err(error, 'define dimension lon for file='//trim(outfile) ) - error = nf_def_dim(ncid, 'lat', dim2, dim_lat) + error = nf90_def_dim(ncid, 'lat', dim2, dim_lat) call netcdf_err(error, 'define dimension lat for file='//trim(outfile) ) !--- define field !---geolon - error = nf_def_var(ncid, 'geolon', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolon) + error = nf90_def_var(ncid, 'geolon', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolon) call netcdf_err(error, 'define var geolon for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "long_name", 9, "Longitude") + error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude") call netcdf_err(error, 'define geolon name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "units", 12, "degrees_east") + error = nf90_put_att(ncid, id_geolon, "units", "degrees_east") call netcdf_err(error, 'define geolon units for file='//trim(outfile) ) !---geolat - error = nf_def_var(ncid, 'geolat', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolat) + error = nf90_def_var(ncid, 'geolat', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolat) call netcdf_err(error, 'define var geolat for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "long_name", 8, "Latitude") + error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude") call netcdf_err(error, 'define geolat name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "units", 13, "degrees_north") + error = nf90_put_att(ncid, id_geolat, "units", "degrees_north") call netcdf_err(error, 'define geolat units for file='//trim(outfile) ) !---slmsk - error = nf_def_var(ncid, 'slmsk', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_slmsk) + error = nf90_def_var(ncid, 'slmsk', NF90_FLOAT,(/dim_lon,dim_lat/), id_slmsk) call netcdf_err(error, 'define var slmsk for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_slmsk, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat") call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) ) !--- land_frac - error = nf_def_var(ncid, 'land_frac', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_land_frac) + error = nf90_def_var(ncid, 'land_frac', NF90_FLOAT, (/dim_lon,dim_lat/), id_land_frac) call netcdf_err(error, 'define var land_frac for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_land_frac, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat") call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) ) !---orography - raw - error = nf_def_var(ncid, 'orog_raw', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_orog_raw) + error = nf90_def_var(ncid, 'orog_raw', NF90_FLOAT, (/dim_lon,dim_lat/), id_orog_raw) call netcdf_err(error, 'define var orog_raw for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_orog_raw, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_orog_raw, "coordinates", "geolon geolat") call netcdf_err(error, 'define orog_raw coordinates for file='//trim(outfile) ) !---orography - filtered - error = nf_def_var(ncid, 'orog_filt', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_orog_filt) + error = nf90_def_var(ncid, 'orog_filt', NF90_FLOAT, (/dim_lon,dim_lat/), id_orog_filt) call netcdf_err(error, 'define var orog_filt for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_orog_filt, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_orog_filt, "coordinates", "geolon geolat") call netcdf_err(error, 'define orog_filt coordinates for file='//trim(outfile) ) !---stddev - error = nf_def_var(ncid, 'stddev', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_stddev) + error = nf90_def_var(ncid, 'stddev', NF90_FLOAT, (/dim_lon,dim_lat/), id_stddev) call netcdf_err(error, 'define var stddev for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_stddev, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_stddev, "coordinates", "geolon geolat") call netcdf_err(error, 'define stddev coordinates for file='//trim(outfile) ) !---convexity - error = nf_def_var(ncid, 'convexity', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_convex) + error = nf90_def_var(ncid, 'convexity', NF90_FLOAT, (/dim_lon,dim_lat/), id_convex) call netcdf_err(error, 'define var convexity for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_convex, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_convex, "coordinates", "geolon geolat") call netcdf_err(error, 'define convexity coordinates for file='//trim(outfile) ) !---oa1 -> oa4 - error = nf_def_var(ncid, 'oa1', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa1) + error = nf90_def_var(ncid, 'oa1', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa1) call netcdf_err(error, 'define var oa1 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa1, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa1, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa1 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa2', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa2) + error = nf90_def_var(ncid, 'oa2', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa2) call netcdf_err(error, 'define var oa2 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa2, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa2, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa2 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa3', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa3) + error = nf90_def_var(ncid, 'oa3', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa3) call netcdf_err(error, 'define var oa3 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa3, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa3, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa3 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa4', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa4) + error = nf90_def_var(ncid, 'oa4', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa4) call netcdf_err(error, 'define var oa4 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa4, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa4, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa4 coordinates for file='//trim(outfile) ) !---ol1 -> ol4 - error = nf_def_var(ncid, 'ol1', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol1) + error = nf90_def_var(ncid, 'ol1', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol1) call netcdf_err(error, 'define var ol1 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol1, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol1, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol1 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol2', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol2) + error = nf90_def_var(ncid, 'ol2', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol2) call netcdf_err(error, 'define var ol2 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol2, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol2, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol2 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol3', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol3) + error = nf90_def_var(ncid, 'ol3', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol3) call netcdf_err(error, 'define var ol3 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol3, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol3, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol3 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol4', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol4) + error = nf90_def_var(ncid, 'ol4', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol4) call netcdf_err(error, 'define var ol4 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol4, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol4, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol4 coordinates for file='//trim(outfile) ) !---theta gamma sigma elvmax - error = nf_def_var(ncid, 'theta', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_theta) + error = nf90_def_var(ncid, 'theta', NF90_FLOAT, (/dim_lon,dim_lat/), id_theta) call netcdf_err(error, 'define var theta for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_theta, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_theta, "coordinates", "geolon geolat") call netcdf_err(error, 'define theta coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'gamma', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_gamma) + error = nf90_def_var(ncid, 'gamma', NF90_FLOAT, (/dim_lon,dim_lat/), id_gamma) call netcdf_err(error, 'define var gamma for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_gamma, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_gamma, "coordinates", "geolon geolat") call netcdf_err(error, 'define gamma coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'sigma', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_sigma) + error = nf90_def_var(ncid, 'sigma', NF90_FLOAT, (/dim_lon,dim_lat/), id_sigma) call netcdf_err(error, 'define var sigma for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_sigma, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_sigma, "coordinates", "geolon geolat") call netcdf_err(error, 'define sigma coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'elvmax', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_elvmax) + error = nf90_def_var(ncid, 'elvmax', NF90_FLOAT, (/dim_lon,dim_lat/), id_elvmax) call netcdf_err(error, 'define var elvmax for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_elvmax, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_elvmax, "coordinates", "geolon geolat") call netcdf_err(error, 'define elvmax coordinates for file='//trim(outfile) ) - error = nf__enddef(ncid, header_buffer_val,4,0,4) + error = nf90_enddef(ncid, header_buffer_val,4,0,4) call netcdf_err(error, 'end meta define for file='//trim(outfile) ) !--- write out data - error = nf_put_var_double( ncid, id_geolon, geolon(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2)) call netcdf_err(error, 'write var geolon for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_geolat, geolat(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2)) call netcdf_err(error, 'write var geolat for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_slmsk, slm(:dim1,:dim2)) + error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2)) call netcdf_err(error, 'write var slmsk for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_land_frac, land_frac(:dim1,:dim2)) + error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2)) call netcdf_err(error, 'write var land_frac for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_orog_raw, oro(:dim1,:dim2)) + error = nf90_put_var( ncid, id_orog_raw, oro(:dim1,:dim2)) call netcdf_err(error, 'write var orog_raw for file='//trim(outfile) ) ! We no longer filter the orog, so the raw and filtered records are the same. - error = nf_put_var_double( ncid, id_orog_filt, oro(:dim1,:dim2)) + error = nf90_put_var( ncid, id_orog_filt, oro(:dim1,:dim2)) call netcdf_err(error, 'write var orog_filt for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_stddev, hprime(:dim1,:dim2,1)) + error = nf90_put_var( ncid, id_stddev, hprime(:dim1,:dim2,1)) call netcdf_err(error, 'write var stddev for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_convex, hprime(:dim1,:dim2,2)) + error = nf90_put_var( ncid, id_convex, hprime(:dim1,:dim2,2)) call netcdf_err(error, 'write var convex for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa1, hprime(:dim1,:dim2,3)) + error = nf90_put_var( ncid, id_oa1, hprime(:dim1,:dim2,3)) call netcdf_err(error, 'write var oa1 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa2, hprime(:dim1,:dim2,4)) + error = nf90_put_var( ncid, id_oa2, hprime(:dim1,:dim2,4)) call netcdf_err(error, 'write var oa2 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa3, hprime(:dim1,:dim2,5)) + error = nf90_put_var( ncid, id_oa3, hprime(:dim1,:dim2,5)) call netcdf_err(error, 'write var oa3 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa4, hprime(:dim1,:dim2,6)) + error = nf90_put_var( ncid, id_oa4, hprime(:dim1,:dim2,6)) call netcdf_err(error, 'write var oa4 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol1, hprime(:dim1,:dim2,7)) + error = nf90_put_var( ncid, id_ol1, hprime(:dim1,:dim2,7)) call netcdf_err(error, 'write var ol1 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol2, hprime(:dim1,:dim2,8)) + error = nf90_put_var( ncid, id_ol2, hprime(:dim1,:dim2,8)) call netcdf_err(error, 'write var ol2 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol3, hprime(:dim1,:dim2,9)) + error = nf90_put_var( ncid, id_ol3, hprime(:dim1,:dim2,9)) call netcdf_err(error, 'write var ol3 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol4, hprime(:dim1,:dim2,10)) + error = nf90_put_var( ncid, id_ol4, hprime(:dim1,:dim2,10)) call netcdf_err(error, 'write var ol4 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_theta, hprime(:dim1,:dim2,11)) + error = nf90_put_var( ncid, id_theta, hprime(:dim1,:dim2,11)) call netcdf_err(error, 'write var theta for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_gamma, hprime(:dim1,:dim2,12)) + error = nf90_put_var( ncid, id_gamma, hprime(:dim1,:dim2,12)) call netcdf_err(error, 'write var gamma for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_sigma, hprime(:dim1,:dim2,13)) + error = nf90_put_var( ncid, id_sigma, hprime(:dim1,:dim2,13)) call netcdf_err(error, 'write var sigma for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_elvmax, hprime(:dim1,:dim2,14)) + error = nf90_put_var( ncid, id_elvmax, hprime(:dim1,:dim2,14)) call netcdf_err(error, 'write var elvmax for file='//trim(outfile) ) - error = nf_close(ncid) + error = nf90_close(ncid) call netcdf_err(error, 'close file='//trim(outfile) ) end subroutine write_netcdf @@ -236,18 +237,19 @@ end subroutine write_netcdf !! @param[in] string The NetCDF error message !! @author Jordan Alpert NOAA/EMC subroutine netcdf_err( err, string ) + use netcdf + implicit none integer, intent(in) :: err character(len=*), intent(in) :: string character(len=256) :: errmsg - include "netcdf.inc" - if( err.EQ.NF_NOERR )return - errmsg = NF_STRERROR(err) + if( err.EQ.NF90_NOERR )return + errmsg = NF90_STRERROR(err) print*, 'FATAL ERROR: ', trim(string), ': ', trim(errmsg) call abort return - end subroutine netcdf_err + end subroutine netcdf_err !> Write the land mask file !! @@ -262,6 +264,7 @@ end subroutine netcdf_err !! @author George Gayno NOAA/EMC subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat) + use netcdf implicit none integer, intent(in):: im, jm, ntiles, tile real, intent(in), dimension(im,jm) :: slm, geolon, geolat, land_frac @@ -273,7 +276,6 @@ subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geola integer :: dim_lon, dim_lat integer :: id_geolon,id_geolat integer :: id_slmsk,id_land_frac - include "netcdf.inc" if(ntiles > 1) then write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc' @@ -285,57 +287,58 @@ subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geola dim2=jm !--- open the file - error = NF__CREATE(outfile, IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid) + error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), ncid, & + initialsize=inital, chunksize=fsize) call netcdf_err(error, 'Creating file '//trim(outfile) ) !--- define dimension - error = nf_def_dim(ncid, 'lon', dim1, dim_lon) + error = nf90_def_dim(ncid, 'lon', dim1, dim_lon) call netcdf_err(error, 'define dimension lon for file='//trim(outfile) ) - error = nf_def_dim(ncid, 'lat', dim2, dim_lat) + error = nf90_def_dim(ncid, 'lat', dim2, dim_lat) call netcdf_err(error, 'define dimension lat for file='//trim(outfile) ) !--- define field !---geolon - error = nf_def_var(ncid, 'geolon', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolon) + error = nf90_def_var(ncid, 'geolon', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolon) call netcdf_err(error, 'define var geolon for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "long_name", 9, "Longitude") + error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude") call netcdf_err(error, 'define geolon name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "units", 12, "degrees_east") + error = nf90_put_att(ncid, id_geolon, "units", "degrees_east") call netcdf_err(error, 'define geolon units for file='//trim(outfile) ) !---geolat - error = nf_def_var(ncid, 'geolat', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolat) + error = nf90_def_var(ncid, 'geolat', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolat) call netcdf_err(error, 'define var geolat for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "long_name", 8, "Latitude") + error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude") call netcdf_err(error, 'define geolat name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "units", 13, "degrees_north") + error = nf90_put_att(ncid, id_geolat, "units", "degrees_north") call netcdf_err(error, 'define geolat units for file='//trim(outfile) ) !---slmsk - error = nf_def_var(ncid, 'slmsk', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_slmsk) + error = nf90_def_var(ncid, 'slmsk', NF90_FLOAT, (/dim_lon,dim_lat/), id_slmsk) call netcdf_err(error, 'define var slmsk for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_slmsk, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat") call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) ) !--- land_frac - error = nf_def_var(ncid, 'land_frac', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_land_frac) + error = nf90_def_var(ncid, 'land_frac', NF90_FLOAT, (/dim_lon,dim_lat/), id_land_frac) call netcdf_err(error, 'define var land_frac for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_land_frac, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat") call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) ) - error = nf__enddef(ncid, header_buffer_val,4,0,4) + error = nf90_enddef(ncid, header_buffer_val,4,0,4) call netcdf_err(error, 'end meta define for file='//trim(outfile) ) !--- write out data - error = nf_put_var_double( ncid, id_geolon, geolon(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2)) call netcdf_err(error, 'write var geolon for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_geolat, geolat(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2)) call netcdf_err(error, 'write var geolat for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_slmsk, slm(:dim1,:dim2)) + error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2)) call netcdf_err(error, 'write var slmsk for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_land_frac, land_frac(:dim1,:dim2)) + error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2)) call netcdf_err(error, 'write var land_frac for file='//trim(outfile) ) - error = nf_close(ncid) + error = nf90_close(ncid) call netcdf_err(error, 'close file='//trim(outfile) ) end subroutine write_mask_netcdf @@ -343,17 +346,18 @@ end subroutine write_mask_netcdf !> Read the land mask file !! !! @param[in] merge_file path -!! @param[in] slm Land-sea mask. -!! @param[in] land_frac Land fraction. -!! @param[in] lake_frac Lake fraction +!! @param[out] slm Land-sea mask. +!! @param[out] land_frac Land fraction. +!! @param[out] lake_frac Lake fraction !! @param[in] im 'i' dimension of a model grid tile. !! @param[in] jm 'j' dimension of a model grid tile. !! @author George Gayno NOAA/EMC subroutine read_mask(merge_file,slm,land_frac,lake_frac,im,jm) + use netcdf + implicit none - include "netcdf.inc" character(len=*), intent(in) :: merge_file @@ -363,30 +367,28 @@ subroutine read_mask(merge_file,slm,land_frac,lake_frac,im,jm) real, intent(out) :: lake_frac(im,jm) real, intent(out) :: slm(im,jm) - integer ncid, error, fsize, id_var - - fsize = 66536 + integer ncid, error, id_var print*,'- READ IN EXTERNAL LANDMASK FILE: ',trim(merge_file) - error=NF__OPEN(merge_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(merge_file,nf90_nowrite,ncid) call netcdf_err(error, 'Open file '//trim(merge_file) ) - error=nf_inq_varid(ncid, 'land_frac', id_var) + error=nf90_inq_varid(ncid, 'land_frac', id_var) call netcdf_err(error, 'inquire varid of land_frac') - error=nf_get_var_double(ncid, id_var, land_frac) + error=nf90_get_var(ncid, id_var, land_frac) call netcdf_err(error, 'inquire data of land_frac') - error=nf_inq_varid(ncid, 'slmsk', id_var) + error=nf90_inq_varid(ncid, 'slmsk', id_var) call netcdf_err(error, 'inquire varid of slmsk') - error=nf_get_var_double(ncid, id_var, slm) + error=nf90_get_var(ncid, id_var, slm) call netcdf_err(error, 'inquire data of slmsk') - error=nf_inq_varid(ncid, 'lake_frac', id_var) + error=nf90_inq_varid(ncid, 'lake_frac', id_var) call netcdf_err(error, 'inquire varid of lake_frac') - error=nf_get_var_double(ncid, id_var, lake_frac) + error=nf90_get_var(ncid, id_var, lake_frac) call netcdf_err(error, 'inquire data of lake_frac') - error = nf_close(ncid) + error = nf90_close(ncid) end subroutine read_mask @@ -398,33 +400,32 @@ end subroutine read_mask !! @author George Gayno NOAA/EMC subroutine read_mdl_dims(mdl_grid_file, im, jm) + use netcdf + implicit none - include "netcdf.inc" character(len=*), intent(in) :: mdl_grid_file integer, intent(out) :: im, jm - integer ncid, error, fsize, id_dim, nx, ny - - fsize = 66536 + integer ncid, error, id_dim, nx, ny print*, "- READ MDL GRID DIMENSIONS FROM= ", trim(mdl_grid_file) - error=NF__OPEN(mdl_grid_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(mdl_grid_file, nf90_nowrite, ncid) call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) ) - error=nf_inq_dimid(ncid, 'nx', id_dim) + error=nf90_inq_dimid(ncid, 'nx', id_dim) call netcdf_err(error, 'inquire dimension nx from file '// trim(mdl_grid_file) ) - error=nf_inq_dimlen(ncid,id_dim,nx) + error=nf90_inquire_dimension(ncid, id_dim, len=nx) call netcdf_err(error, 'inquire nx from file '//trim(mdl_grid_file) ) - error=nf_inq_dimid(ncid, 'ny', id_dim) + error=nf90_inq_dimid(ncid, 'ny', id_dim) call netcdf_err(error, 'inquire dimension ny from file '// trim(mdl_grid_file) ) - error=nf_inq_dimlen(ncid,id_dim,ny) + error=nf90_inquire_dimension(ncid, id_dim, len=ny) call netcdf_err(error, 'inquire ny from file '//trim(mdl_grid_file) ) - error=nf_close(ncid) + error=nf90_close(ncid) IM = nx/2 JM = ny/2 @@ -451,10 +452,11 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & geolon, geolon_c, geolat, geolat_c, dx, dy, & is_north_pole, is_south_pole) + use netcdf + use orog_utils, only : find_poles, find_nearest_pole_points implicit none - include "netcdf.inc" character(len=*), intent(in) :: mdl_grid_file @@ -470,12 +472,11 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & real, intent(out) :: dx(im,jm), dy(im,jm) integer :: i, j - integer :: ncid, error, fsize, id_var, nx, ny + integer :: ncid, error, id_var, nx, ny integer :: i_south_pole,j_south_pole integer :: i_north_pole,j_north_pole real, allocatable :: tmpvar(:,:) - fsize = 66536 nx = 2*im ny = 2*jm @@ -484,12 +485,12 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & print*, "- OPEN AND READ= ", trim(mdl_grid_file) - error=NF__OPEN(mdl_grid_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(mdl_grid_file, nf90_nowrite, ncid) call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) ) - error=nf_inq_varid(ncid, 'x', id_var) + error=nf90_inq_varid(ncid, 'x', id_var) call netcdf_err(error, 'inquire varid of x from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of x from file ' // trim(mdl_grid_file)) ! Adjust lontitude to be between 0 and 360. @@ -503,9 +504,9 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & geolon(1:IM,1:JM) = tmpvar(2:nx:2,2:ny:2) geolon_c(1:IM+1,1:JM+1) = tmpvar(1:nx+1:2,1:ny+1:2) - error=nf_inq_varid(ncid, 'y', id_var) + error=nf90_inq_varid(ncid, 'y', id_var) call netcdf_err(error, 'inquire varid of y from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of y from file ' // trim(mdl_grid_file)) geolat(1:IM,1:JM) = tmpvar(2:nx:2,2:ny:2) @@ -522,12 +523,12 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & allocate(tmpvar(nx,ny)) - error=nf_inq_varid(ncid, 'area', id_var) + error=nf90_inq_varid(ncid, 'area', id_var) call netcdf_err(error, 'inquire varid of area from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of area from file ' // trim(mdl_grid_file)) - error = nf_close(ncid) + error = nf90_close(ncid) do j = 1, jm do i = 1, im @@ -550,28 +551,48 @@ end subroutine read_mdl_grid_file subroutine read_global_orog(imn,jmn,glob) use orog_utils, only : transpose_orog + use netcdf implicit none - include 'netcdf.inc' - integer, intent(in) :: imn, jmn integer*2, intent(out) :: glob(imn,jmn) - integer :: ncid, error, id_var, fsize - - fsize=65536 + integer :: ncid, error, id_dim, id_var, idim, jdim print*,"- OPEN AND READ ./topography.gmted2010.30s.nc" - error=NF__OPEN("./topography.gmted2010.30s.nc", & - NF_NOWRITE,fsize,ncid) - call netcdf_err(error, 'Open file topography.gmted2010.30s.nc' ) - error=nf_inq_varid(ncid, 'topo', id_var) + error=nf90_open("./topography.gmted2010.30s.nc", & + nf90_nowrite, ncid) + call netcdf_err(error, 'Opening file topography.gmted2010.30s.nc' ) + + error=nf90_inq_dimid(ncid, 'idim', id_dim) + call netcdf_err(error, 'Inquire dimid of idim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=idim) + call netcdf_err(error, 'Reading idim' ) + + if (imn /= idim) then + print*,"FATAL ERROR: i-dimensions do not match." + endif + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=jdim) + call netcdf_err(error, 'Reading jdim' ) + + if (jmn /= jdim) then + print*,"FATAL ERROR: j-dimensions do not match." + endif + + error=nf90_inq_varid(ncid, 'topo', id_var) call netcdf_err(error, 'Inquire varid of topo') - error=nf_get_var_int2(ncid, id_var, glob) - call netcdf_err(error, 'Read topo') - error = nf_close(ncid) + + error=nf90_get_var(ncid, id_var, glob) + call netcdf_err(error, 'Reading topo') + + error = nf90_close(ncid) print*,"- MAX/MIN OF OROGRAPHY DATA ",maxval(glob),minval(glob) @@ -589,28 +610,48 @@ end subroutine read_global_orog subroutine read_global_mask(imn, jmn, mask) use orog_utils, only : transpose_mask + use netcdf implicit none - include 'netcdf.inc' - integer, intent(in) :: imn, jmn integer(1), intent(out) :: mask(imn,jmn) - integer :: ncid, fsize, id_var, error - - fsize = 65536 + integer :: ncid, id_var, id_dim, error, idim, jdim print*,"- OPEN AND READ ./landcover.umd.30s.nc" - error=NF__OPEN("./landcover.umd.30s.nc",NF_NOWRITE,fsize,ncid) - call netcdf_err(error, 'Open file landcover.umd.30s.nc' ) - error=nf_inq_varid(ncid, 'land_mask', id_var) + error=nf90_open("./landcover.umd.30s.nc",nf90_nowrite,ncid) + call netcdf_err(error, 'Opening file landcover.umd.30s.nc' ) + + error=nf90_inq_dimid(ncid, 'idim', id_dim) + call netcdf_err(error, 'Inquire dimid of idim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=idim) + call netcdf_err(error, 'Reading idim' ) + + if (imn /= idim) then + print*,"FATAL ERROR: i-dimensions do not match." + endif + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=jdim) + call netcdf_err(error, 'Reading jdim' ) + + if (jmn /= jdim) then + print*,"FATAL ERROR: j-dimensions do not match." + endif + + error=nf90_inq_varid(ncid, 'land_mask', id_var) call netcdf_err(error, 'Inquire varid of land_mask') - error=nf_get_var_int1(ncid, id_var, mask) + + error=nf90_get_var(ncid, id_var, mask) call netcdf_err(error, 'Inquire data of land_mask') - error = nf_close(ncid) + + error = nf90_close(ncid) call transpose_mask(imn,jmn,mask) @@ -626,48 +667,53 @@ end subroutine read_global_mask !! @author G. Gayno subroutine qc_orog_by_ramp(imn, jmn, zavg, zslm) - implicit none + use netcdf - include 'netcdf.inc' + implicit none integer, intent(in) :: imn, jmn integer, intent(inout) :: zavg(imn,jmn) integer, intent(inout) :: zslm(imn,jmn) - integer :: i, j, error, ncid, id_var, fsize + integer :: i, j, error, ncid, id_var, id_dim, jramp real(4), allocatable :: gice(:,:) - fsize = 65536 - - allocate (GICE(IMN+1,3601)) - ! Read 30-sec Antarctica RAMP data. Points scan from South ! to North, and from Greenwich to Greenwich. print*,"- OPEN/READ RAMP DATA ./topography.antarctica.ramp.30s.nc" - error=NF__OPEN("./topography.antarctica.ramp.30s.nc", & - NF_NOWRITE,fsize,ncid) + error=nf90_open("./topography.antarctica.ramp.30s.nc", & + nf90_nowrite, ncid) call netcdf_err(error, 'Opening RAMP topo file' ) - error=nf_inq_varid(ncid, 'topo', id_var) + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid, id_dim, len=jramp) + call netcdf_err(error, 'Reading jdim' ) + + allocate (GICE(IMN+1,jramp)) + + error=nf90_inq_varid(ncid, 'topo', id_var) call netcdf_err(error, 'Inquire varid of RAMP topo') - error=nf_get_var_real(ncid, id_var, GICE) + + error=nf90_get_var(ncid, id_var, GICE) call netcdf_err(error, 'Inquire data of RAMP topo') - error = nf_close(ncid) + + error = nf90_close(ncid) print*,"- QC GLOBAL OROGRAPHY DATA WITH RAMP." ! If RAMP values are valid, replace the global value with the RAMP ! value. Invalid values are less than or equal to 0 (0, -1, or -99). - do j = 1, 3601 + do j = 1, jramp do i = 1, IMN - if( GICE(i,j) .ne. -99. .and. GICE(i,j) .ne. -1.0 ) then - if ( GICE(i,j) .gt. 0.) then - ZAVG(i,j) = int( GICE(i,j) + 0.5 ) - ZSLM(i,j) = 1 - endif + if ( GICE(i,j) .gt. 0.) then + ZAVG(i,j) = int( GICE(i,j) + 0.5 ) + ZSLM(i,j) = 1 endif enddo enddo diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index d8e55c96a..57cf5f2b7 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -419,15 +419,13 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & real, intent(out) :: slm(im,jm) real, intent(out) :: land_frac(im,jm) - integer, parameter :: MAXSUM=20000000 - real, parameter :: D2R = 3.14159265358979/180. integer jst, jen real GLAT(JMN), GLON(IMN) real LONO(4),LATO(4),LONI,LATI real LONO_RAD(4), LATO_RAD(4) - integer JM1,i,j,nsum,nsum_all,ii,jj,numx,i2 + integer JM1,i,j,ii,jj,numx,i2 integer ilist(IMN) real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1 real XNSUM_ALL,XLAND_ALL,XWATR_ALL @@ -451,20 +449,18 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & ! ! (*j*) for hard wired zero offset (lambda s =0) for terr05 !$omp parallel do & -!$omp private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,lono, & +!$omp private (j,i,xnsum,xland,xwatr,xl1,xs1,xw1,lono, & !$omp lato,lono_rad,lato_rad,jst,jen,ilist,numx,jj,i2,ii,loni,lati, & -!$omp xnsum_all,xland_all,xwatr_all,nsum_all) +!$omp xnsum_all,xland_all,xwatr_all) ! DO J=1,JM DO I=1,IM XNSUM = 0.0 XLAND = 0.0 XWATR = 0.0 - nsum = 0 XNSUM_ALL = 0.0 XLAND_ALL = 0.0 XWATR_ALL = 0.0 - nsum_all = 0 LONO(1) = lon_c(i,j) LONO(2) = lon_c(i+1,j) @@ -485,12 +481,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XLAND_ALL = XLAND_ALL + FLOAT(ZSLM(ii,jj)) XWATR_ALL = XWATR_ALL + FLOAT(1-ZSLM(ii,jj)) XNSUM_ALL = XNSUM_ALL + 1. - nsum_all = nsum_all+1 - if(nsum_all > MAXSUM) then - print*, "FATAL ERROR: nsum_all is greater than MAXSUM," - print*, "increase MAXSUM." - call ABORT() - endif if(inside_a_polygon(LONI*D2R,LATI*D2R,4, & LONO_RAD,LATO_RAD))then @@ -498,12 +488,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XLAND = XLAND + FLOAT(ZSLM(ii,jj)) XWATR = XWATR + FLOAT(1-ZSLM(ii,jj)) XNSUM = XNSUM + 1. - nsum = nsum+1 - if(nsum > MAXSUM) then - print*, "FATAL ERROR: nsum is greater than MAXSUM," - print*, "increase MAXSUM." - call ABORT() - endif endif enddo ; enddo @@ -559,13 +543,12 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real, intent(out) :: oro(im,jm) real, intent(out) :: var(im,jm),var4(im,jm) - integer, parameter :: MAXSUM=20000000 real, parameter :: D2R = 3.14159265358979/180. real, dimension(:), allocatable :: hgt_1d, hgt_1d_all real GLAT(JMN), GLON(IMN) - integer JST, JEN + integer JST, JEN, maxsum real LONO(4),LATO(4),LONI,LATI real LONO_RAD(4), LATO_RAD(4) real HEIGHT @@ -576,8 +559,6 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real XL1_ALL,XS1_ALL,XW1_ALL,XW2_ALL,XW4_ALL print*,'- CREATE OROGRAPHY AND CONVEXITY.' - allocate(hgt_1d(MAXSUM)) - allocate(hgt_1d_all(MAXSUM)) !---- GLOBAL XLAT AND XLON ( DEGREE ) ! JM1 = JM - 1 @@ -590,7 +571,24 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & GLON(I) = 0. + (I-1) * DELXN + DELXN * 0.5 ENDDO -! land_frac(:,:) = 0.0 + MAXSUM=0 + DO J=1,JM + DO I=1,IM + LONO(1) = lon_c(i,j) + LONO(2) = lon_c(i+1,j) + LONO(3) = lon_c(i+1,j+1) + LONO(4) = lon_c(i,j+1) + LATO(1) = lat_c(i,j) + LATO(2) = lat_c(i+1,j) + LATO(3) = lat_c(i+1,j+1) + LATO(4) = lat_c(i,j+1) + call get_index(IMN,JMN,4,LONO,LATO,DELXN,jst,jen,ilist,numx) + MAXSUM=MAX(MAXSUM,(JEN-JST+1)*IMN) + ENDDO + ENDDO + print*,"- MAXSUM IS ", maxsum + allocate(hgt_1d(MAXSUM)) + allocate(hgt_1d_all(MAXSUM)) ! !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX ! @@ -649,7 +647,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & nsum_all = nsum_all+1 if(nsum_all > MAXSUM) then print*, "FATAL ERROR: nsum_all is greater than MAXSUM," - print*, "increase MAXSUM." + print*, "increase MAXSUM.", jst,jen call ABORT() endif hgt_1d_all(nsum_all) = HEIGHT_ALL @@ -668,7 +666,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & nsum = nsum+1 if(nsum > MAXSUM) then print*, "FATAL ERROR: nsum is greater than MAXSUM," - print*, "increase MAXSUM." + print*, "increase MAXSUM.", jst,jen call ABORT() endif hgt_1d(nsum) = HEIGHT diff --git a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 index cae1f2bec..b0e14e9d4 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 @@ -574,7 +574,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) integer :: iw, ie, wgta, is, ise integer :: in, ine, inw, isw - real :: slma, oroa, vara, var4a, xn, xs + real :: slma, oroa, vara, var4a real, allocatable :: oaa(:), ola(:) ! REMOVE ISOLATED POINTS @@ -586,9 +586,10 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) iso_loop : DO J=2,JM-1 JN=J-1 JS=J+1 - i_loop : DO I=1,IM - IW=MOD(I+IM-2,IM)+1 - IE=MOD(I,IM)+1 + i_loop : DO I=2,IM-1 +! Check the points to the 'west' and 'east'. + IW=I-1 + IE=I+1 SLMA=SLM(IW,J)+SLM(IE,J) OROA=ORO(IW,J)+ORO(IE,J) VARA=VAR(IW,J)+VAR(IE,J) @@ -599,60 +600,33 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OLA(K)=OL(IW,J,K)+OL(IE,J,K) ENDDO WGTA=2 - XN=(I-1)+1 - IF(ABS(XN-NINT(XN)).LT.1.E-2) THEN - IN=MOD(NINT(XN)-1,IM)+1 - INW=MOD(IN+IM-2,IM)+1 - INE=MOD(IN,IM)+1 - SLMA=SLMA+SLM(INW,JN)+SLM(IN,JN)+SLM(INE,JN) - OROA=OROA+ORO(INW,JN)+ORO(IN,JN)+ORO(INE,JN) - VARA=VARA+VAR(INW,JN)+VAR(IN,JN)+VAR(INE,JN) - VAR4A=VAR4A+VAR4(INW,JN)+VAR4(IN,JN)+VAR4(INE,JN) - DO K=1,4 - OAA(K)=OAA(K)+OA(INW,JN,K)+OA(IN,JN,K)+OA(INE,JN,K) - OLA(K)=OLA(K)+OL(INW,JN,K)+OL(IN,JN,K)+OL(INE,JN,K) - ENDDO - WGTA=WGTA+3 - ELSE - INW=INT(XN) - INE=MOD(INW,IM)+1 - SLMA=SLMA+SLM(INW,JN)+SLM(INE,JN) - OROA=OROA+ORO(INW,JN)+ORO(INE,JN) - VARA=VARA+VAR(INW,JN)+VAR(INE,JN) - VAR4A=VAR4A+VAR4(INW,JN)+VAR4(INE,JN) - DO K=1,4 - OAA(K)=OAA(K)+OA(INW,JN,K)+OA(INE,JN,K) - OLA(K)=OLA(K)+OL(INW,JN,K)+OL(INE,JN,K) - ENDDO - WGTA=WGTA+2 - ENDIF - XS=(I-1)+1 - IF(ABS(XS-NINT(XS)).LT.1.E-2) THEN - IS=MOD(NINT(XS)-1,IM)+1 - ISW=MOD(IS+IM-2,IM)+1 - ISE=MOD(IS,IM)+1 - SLMA=SLMA+SLM(ISW,JS)+SLM(IS,JS)+SLM(ISE,JS) - OROA=OROA+ORO(ISW,JS)+ORO(IS,JS)+ORO(ISE,JS) - VARA=VARA+VAR(ISW,JS)+VAR(IS,JS)+VAR(ISE,JS) - VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(IS,JS)+VAR4(ISE,JS) - DO K=1,4 - OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(IS,JS,K)+OA(ISE,JS,K) - OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(IS,JS,K)+OL(ISE,JS,K) - ENDDO - WGTA=WGTA+3 - ELSE - ISW=INT(XS) - ISE=MOD(ISW,IM)+1 - SLMA=SLMA+SLM(ISW,JS)+SLM(ISE,JS) - OROA=OROA+ORO(ISW,JS)+ORO(ISE,JS) - VARA=VARA+VAR(ISW,JS)+VAR(ISE,JS) - VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(ISE,JS) - DO K=1,4 - OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(ISE,JS,K) - OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(ISE,JS,K) - ENDDO - WGTA=WGTA+2 - ENDIF +! Check the points to the 'northwest', 'north' and 'northeast' + IN=I + INW=I-1 + INE=I+1 + SLMA=SLMA+SLM(INW,JN)+SLM(IN,JN)+SLM(INE,JN) + OROA=OROA+ORO(INW,JN)+ORO(IN,JN)+ORO(INE,JN) + VARA=VARA+VAR(INW,JN)+VAR(IN,JN)+VAR(INE,JN) + VAR4A=VAR4A+VAR4(INW,JN)+VAR4(IN,JN)+VAR4(INE,JN) + DO K=1,4 + OAA(K)=OAA(K)+OA(INW,JN,K)+OA(IN,JN,K)+OA(INE,JN,K) + OLA(K)=OLA(K)+OL(INW,JN,K)+OL(IN,JN,K)+OL(INE,JN,K) + ENDDO + WGTA=WGTA+3 +! Check the points to the 'southwest', 'south' and 'southeast' + IS=I + ISW=I-1 + ISE=I+1 + SLMA=SLMA+SLM(ISW,JS)+SLM(IS,JS)+SLM(ISE,JS) + OROA=OROA+ORO(ISW,JS)+ORO(IS,JS)+ORO(ISE,JS) + VARA=VARA+VAR(ISW,JS)+VAR(IS,JS)+VAR(ISE,JS) + VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(IS,JS)+VAR4(ISE,JS) + DO K=1,4 + OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(IS,JS,K)+OA(ISE,JS,K) + OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(IS,JS,K)+OL(ISE,JS,K) + ENDDO + WGTA=WGTA+3 +! Take the average of the surrounding the points. OROA=OROA/WGTA VARA=VARA/WGTA VAR4A=VAR4A/WGTA @@ -660,6 +634,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OAA(K)=OAA(K)/WGTA OLA(K)=OLA(K)/WGTA ENDDO +! Isolated water point. IF(SLM(I,J).EQ.0..AND.SLMA.EQ.WGTA) THEN PRINT '(" - SEA ",2F8.0," MODIFIED TO LAND",2F8.0, & " AT ",2I8)',ORO(I,J),VAR(I,J),OROA,VARA,I,J @@ -671,6 +646,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OA(I,J,K)=OAA(K) OL(I,J,K)=OLA(K) ENDDO +! Isolated land point. ELSEIF(SLM(I,J).EQ.1..AND.SLMA.EQ.0.) THEN PRINT '(" - LAND",2F8.0," MODIFIED TO SEA ",2F8.0, & " AT ",2I8)',ORO(I,J),VAR(I,J),OROA,VARA,I,J @@ -699,7 +675,8 @@ end subroutine remove_isolated_pts !! @param[in] jmn 'j' dimension of the high-resolution orography !! data set. !! @param[in] npts Number of vertices to describe the cubed-sphere point. -!! @param[in] lonO The longitudes of the cubed-sphere vertices. +!! @param[in] lonO The longitudes of the cubed-sphere vertices. Must +!! range from 0 - 360. !! @param[in] latO The latitudes of the cubed-sphere vertices. !! @param[in] delxn Resolution of the high-resolution orography !! data set. diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 863a4d47d..7c2d894e1 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -8,7 +8,31 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8") endif() -include_directories(${PROJECT_SOURCE_DIR}) +# Copy necessary test files from the source data directory to the +# build test directory. + +# The "read_global_mask" test needs this file. Note that the copy command changes the +# file name to "landcover.umd.30s.nc", which is the name expected by the "read_global_mask" routine. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/landcover.umd.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/landcover.umd.30s.nc) + +# The "read_global_orog" test needs this file. Note that the copy command changes the +# file name to "topography.gmted2010.30s.nc", which is the name expected by the "read_global_orog" routine. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.gmted2010.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.gmted2010.30s.nc) + +# The "read_mdl_dims" and "read_mdl_grid_file" tests expects this file. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/C12_grid.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C12_grid.tile1.nc) + +# The "read_mask" test expects this file. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/C48.mx500.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C48.mx500.tile1.nc) + +# The "qc_orog_by_ramp" test needs this file. Note that the copy command changes the +# file name to "topography.antarctica.ramp.30s.nc", which is the name expected by the "qc_orog_by_ramp" routine. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.antarctica.ramp.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.antarctica.ramp.30s.nc) add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) @@ -21,3 +45,63 @@ target_link_libraries(ftst_minmax orog_lib) add_executable(ftst_get_ll_angle ftst_get_ll_angle.F90) add_test(NAME orog-ftst_get_ll_angle COMMAND ftst_get_ll_angle) target_link_libraries(ftst_get_ll_angle orog_lib) + +add_executable(ftst_get_index ftst_get_index.F90) +add_test(NAME orog-ftst_get_index COMMAND ftst_get_index) +target_link_libraries(ftst_get_index orog_lib) + +add_executable(ftst_inside_polygon ftst_inside_polygon.F90) +add_test(NAME orog-ftst_inside_polygon COMMAND ftst_inside_polygon) +target_link_libraries(ftst_inside_polygon orog_lib) + +add_executable(ftst_transpose ftst_transpose.F90) +add_test(NAME orog-ftst_transpose COMMAND ftst_transpose) +target_link_libraries(ftst_transpose orog_lib) + +add_executable(ftst_find_poles ftst_find_poles.F90) +add_test(NAME orog-ftst_find_poles COMMAND ftst_find_poles) +target_link_libraries(ftst_find_poles orog_lib) + +add_executable(ftst_find_nearest_pole_pts ftst_find_nearest_pole_pts.F90) +add_test(NAME orog-ftst_find_nearest_pole_pts COMMAND ftst_find_nearest_pole_pts) +target_link_libraries(ftst_find_nearest_pole_pts orog_lib) + +add_executable(ftst_get_xnsum ftst_get_xnsum.F90) +add_test(NAME orog-ftst_get_xnsum COMMAND ftst_get_xnsum) +target_link_libraries(ftst_get_xnsum orog_lib) + +add_executable(ftst_get_xnsum2 ftst_get_xnsum2.F90) +add_test(NAME orog-ftst_get_xnsum2 COMMAND ftst_get_xnsum2) +target_link_libraries(ftst_get_xnsum2 orog_lib) + +add_executable(ftst_get_xnsum3 ftst_get_xnsum3.F90) +add_test(NAME orog-ftst_get_xnsum3 COMMAND ftst_get_xnsum3) +target_link_libraries(ftst_get_xnsum3 orog_lib) + +add_executable(ftst_rm_isolated_pts ftst_rm_isolated_pts.F90) +add_test(NAME orog-ftst_rm_isolated_pts COMMAND ftst_rm_isolated_pts) +target_link_libraries(ftst_rm_isolated_pts orog_lib) + +add_executable(ftst_read_global_mask ftst_read_global_mask.F90) +add_test(NAME orog-ftst_read_global_mask COMMAND ftst_read_global_mask) +target_link_libraries(ftst_read_global_mask orog_lib) + +add_executable(ftst_read_global_orog ftst_read_global_orog.F90) +add_test(NAME orog-ftst_read_global_orog COMMAND ftst_read_global_orog) +target_link_libraries(ftst_read_global_orog orog_lib) + +add_executable(ftst_read_mdl_dims ftst_read_mdl_dims.F90) +add_test(NAME orog-ftst_read_mdl_dims COMMAND ftst_read_mdl_dims) +target_link_libraries(ftst_read_mdl_dims orog_lib) + +add_executable(ftst_read_mdl_grid_file ftst_read_mdl_grid_file.F90) +add_test(NAME orog-ftst_read_mdl_grid_file COMMAND ftst_read_mdl_grid_file) +target_link_libraries(ftst_read_mdl_grid_file orog_lib) + +add_executable(ftst_read_mask ftst_read_mask.F90) +add_test(NAME orog-ftst_read_mask COMMAND ftst_read_mask) +target_link_libraries(ftst_read_mask orog_lib) + +add_executable(ftst_qc_orog_by_ramp ftst_qc_orog_by_ramp.F90) +add_test(NAME orog-ftst_qc_orog_by_ramp COMMAND ftst_qc_orog_by_ramp) +target_link_libraries(ftst_qc_orog_by_ramp orog_lib) diff --git a/tests/orog/data/C12_grid.tile1.nc b/tests/orog/data/C12_grid.tile1.nc new file mode 100644 index 0000000000000000000000000000000000000000..5bcfb11514828e109fee9a770c4112b06ef9356c GIT binary patch literal 51102 zcmeHw2S8NUw)Rjg2qG#-vw&hRAfRGF_8FQ5ODteV1XL730To;93ijU77<)JN&bG#0 zV((pJ@4fx!%wFHj@g%+^zMI^4@9(^KU9->Gd#!Knvf3VFdX1ovoH<-_m}k$f!{-xy zX8tUmc+rah{h1Z@Y6LYf&pFo2ys+V+g)P0Sl3RQ^Nc1Vu9OgPR<&B-4^5so-iden{T6kEluZBlI6gy@vS z#Q5a0ZITibQp=Y2a`Uj6rtq|or5WV4k4sKTOzL3})rZD=c27*|=RM9<`?Zn2(?LmEalZd5a@p=V4|oD4HEIXWq>b4s#jxpJPdsj<95NwqoB`#OYTF( z;^@rG?3IV+W@cue=q;#)y)u!x-3uA#)-k)zv9H`Ki(#QGa)^bBDk}z#qwx;;4FfE) z=&(7o)TXr3)k213n5)1rsbNFIOeqVmEJ{$#I>B`t)~r`IvSFYU7EC*$u(D?Ir@bIL zg^jXqwzMPa#->EawDk>&>m;kI-4li}*nz}AA>9^R)vXs;qxK)87TPYWBnCGG4+5o| zdV5$6g^^7&vGKA9I!A6w2%e7q}E(A_OjeDH4<&xzoG{@KN=?JWlW?TdFt9MG>@JcC*I zoy7|rUaTuA{=Xo9_o~T9Qa7XExU>JdRlE34t=cfM8eetfl-({L&fk1mA2aX&n}>ht z!M|_CoCzoBpFL;UJjFomIm`a9Hj@72o4;zs47TidRxGcs#mc|5K)>z$_b*umT}j&` z;lDj(SSrvzyJTI<{XatfR$}!8HM|o&rreg}# z=@8J_6c_Q5mB89AHnCG|3dN|(7F1v+cWxgm$J9~rHYJqVvCvtI&IyT0DeWUWC&mk{ z45;MpUgf;pD$9S^$T6U$sdwk3#1661DTuXm67%S;ZDJ!4vf5OZSJE&Ah9EE|P7Zlv zQ{s9?8HYL0u=vDw3ZnQPkx4SssDyT1I(aUAt#!ViWDA+aR%}0I4#a0>sE7@p zQ6xNoo?{nUe%Z(eD_#H6E>u=Y`G6`0K(X5bf+d-PE=m{WL$*aUbP=W8)9LHi5rY|* zLtrJf`r@8fQVIp0yu;chMRjg3=O=Pz(k?EgYfP+7G3A{x_V1bymy&E)Gvx#;gJNRa zCB?=@#zrNlj5zY(&_A;DDioBK;6TZA~$Wd_g(Egeu=PBD?<@1j(oMLmuk?WAp zfIt^7i1x^JcC}v_(x66t+1@LYtJrn=AT3@&Z6|}>mK(|vC?8PK`n^z^-h57Ei(nin z+=iQCvHZ|We|%eaeWMydviX(!Dd#SxL>Yn0;rz0YeyudGThWpm2C48;N^s7J(+ zwBhCVAq{0-I%pvz3a12#3boYL*QM&*$B6tQy!v|J?iWIr-O8s@-SE3VO#T7xBfUQwwDcP-^tBypMUub7}noZuo5dMwv8-Hte0_% ziAd8hh|LA^HjIv=r*7}1uRE{+b zb7DxE&Tg^}Nj6d*a$&|9aw|AT*wm}mF@O>HNU0>va5 zuJZX4xtXN&`%781a8pR1qfbdQS7X%`w;FU$adk_3Q3 zrQr}gP*zv@;Q2iIf}jV_pHiHIa0k!PA8`+!qdy`ZJV$?oK6pNifFS6>b4KA7q8~g5 ze-)7+341Ak`Z!_?12ENU}w;A|01K(!g+YEf0f&cvssDHeH%fa`QV=H%! z7P_YL4X*MTCWuZKjxYauy4UaF^Un}5AtQ~_ZkrYgqQCMR7NyLod_a|Q27be$sHtu% zwT!|+`&WO4fF>nIey4oIoIImB*(4tRt#fjwXMvQm<1ZQtEU?#?LCIM&P(wAozRsg;H@-BTy=B&1QKY z5;A}Lry3b7U2EzlPW5>K{bp`y?b#jcC>$F&oC8@`e3X+pDaO z@=;URUZwtde5u*QkcJI|b@&yEnS2@kozi=r`H5Z1NYh6_X46~z;xCuZOZf$3gW!$`Yg^Q-HF}ol|M#3V#i%Kgv3iO-el+-XykRbu;&-M0 zs#U=+CqOq-u^+m3rP?^!ST(fU(yb-h2kG}dXiunKy4pXpicN5RAr?pHo1eZA*<#g{buft#7 zQ)ErrXI|Y*9AUYYRr%X^3GHxsVy65c#_%yuNX1Q#nI%rQc|s3d`|07xI-;*zh%|~GiqjuJ0A0xR7|UxB`&>(fzS5<{!rfTWi4m?tW4f`!Ebq~lugzIwyO%eXPIhMT$Ma>5 z{mFx+=%RKMZ(oLG%$xLF_o6&2xcjR0h^ZA=tCmqY6WuB^KR=s(cInla*)Ct!>3D$P z-%k=$)}!y*VT13g^jKiMJkoi&Nsm4E)1Dqw>0#I3X8tgh9($~;8x>UPF>}Y8<`Y$V zwDHe6w}VQLnfZIYJ^9$2H6lH7TuQ9IWVr`RCOu~ObUE$ytUSw4ddv(il%+#KlOEX? zHd+vIEPxrpN|heF+5A~nfu{8sxo}F=N-8}zCR(T8R_Rf8=G>oWsr2xFGVNVyl^(e_ zMCO~O(&J(14|if!dbEui>TQ14v>tI);;pkSHLXV{%bJNppP1I8T=%XaBl4T{SoLgr zmbOQJyB-$$*tW%0dZZrRp1-k5kE$Oo`&hV}_SHzgMm^W5^yr`1@?8~`9$rt9BW9`e zIBGj$VQZBhB?4-EH{hB{k5&&;*Lp4TU|nfF>Xp&$IQYo49=%^*^I7g<(&OyqB~#N6 z1u#SXrmtyVty&_D3smW`ukm`1E-F2C*YbIpU!}(~-mT{Wl^%_P3mmMW(xZ}JztVG6 zdhGkT?}_FrJ(`8Ld@}XCNss1ti~lrczG*$yH#{)I>!E2q;^z!o^T5fZ$C;p;=l2}= z?Rwm6awxKeN{`Uk&4-7n^k`q|R@(|HJxaKgG`phGqx_i7lfzYd#Ck2WUZm1v`jH0D z!&Q28zFn#CniD2H9&NgC?$aF8dh~3XF?P}Y@~i-zj}c31+BbBp%qr12Z_=fdMV0-( z)uTzd0>x8QdL$**tTsob$EQy#285~f=n=4>&NG!BF$Ja!h*s&5dgS`lWhy<=yz;kg zpwffad9u3tVUr%`?VM^)nr&JSIRiXd>7Hplu4Ip^deFhN9;q&6M{V8vTRrwgS^qpl zrN_YUyS-nh(&I|aTstCFdMwOSYrn;Glil^~JCi4<^mw}d-n5@odgQ8MA6r+YM@)OZ za@Jmxy?U{D>$)McOnSVFvc7C_*R&po9jEw>u{W(p^7Wdd$N&6WJu00Za(;?R5APMH zbM9B^ad~&c>Pad+u3VouGOtRHy}O*c_fY8(W&JpAtxAts&vQ+xrP6~>&qxm1X<85S zL4kEfWtjA^_;IJpgYC;mEwE^%yx~Txims-|BIC?%L?NDm~7|{@Ul1N{{G- zhoiq&>9Nt$Q!1#^giMnjAKMo?S80Yx zkNwWhW{xU7y8PI~eM%nFdW>B)Z`a1%ztyAbJ+@(mN{?~PdM~-A(&NbTwiAb`^q83E z)Z!8dn!FncABgkuhQe<=rY;7RCYFdw_h0|Xixn)|9B_Ce~*H!7^l0Pd; z`yD<0>GjC}2kQ~~mFw}}az6gA*{cJ8V6R4h!Cw8h_>X_LKY#iM{`}CN^ymLA@zp;a z?SF~B^jD}iM6I)lA)=WWcK!KY_Syuj@-edaaM^E+ul%YrNedB@-ZnH%Q>@(~fYuUk|M;*lTQ`MPfk;z#z^ z_P6&8<}=Ez8hdwqFn^qNtoN?l!F)rXS+VC!hVWH0S0|5X8p6YSPFT{STL{-FpMIt` z-lIRpPa}Td4|$LeJizBlc!7UVRIMs2OENzF!pP0X&M@A)PvH;c`$#-ui@8U*qn;l= zK562Lxq5yo=eizoMFV;L+{@>MP7dUGRvbwx{Wg$WnYI2{CoG8PYwVXNVPp`mdwKJl z23vyol+4@pj-3hOxos`ng07j`c<)60bI|zQ#1H%-5AuNr_<$Gq>-^N@QtoJf-hOt_ zt6M)+=PoyN)2iPO;G=WasdjA`<978!LYwE6_=4uo?PrXU_^}8(+g_g}-Y7G-Z-dr) zUgTZ5{ij#z`TlwRJ70gS=l2Tj>g!x6kgqI%`(v@Hf!s}ruY&^lPt?ZyebgW0!4Lc) z5AuNr_<$Gq$IZ%LV`owo-s9uB;wu(c<)dz4IcjKsU1YuSEl zUx_#E`>ybUG>NBB8}HE{<3ABU@E=3+ARl;u4|sw9%*ZQ0r3|RVS6fX9Yd^=2X9>7` zbjhm9{NtGDp@-L3;WsnAq?QY+^6vXi9G^C%8qfN$q2wFp&u8Bl`SI;rf8IB9$D7(C zs`J*H*Js?<1@JLT-Z|xN9Kb`|)+fZI1@Iz<`cnX}L2bNmMg1`z{J_5l$%B000Y2ab z{`|e?x)vMf%ZGTkoZn|+1%7b+InR4(75O5MDou~4RN`GX9lTMzjvt?r?Y?s!r^gy6AGG*XgQ!l_9S) zweem@{V{$Y@taTlArJDG5T1R6Z!qBn{x$}CRi0=1T&|ox)SHj~_2j5hkv@FVi8|>! z%KLK5we8cMKK126d+VCDT~vVwP2W-Xdaa5)eXV1Mrspg24uyu+Oleq&H&{RS;=DDL z_@%IeO?SSn#N#ivir84xj|cSgx|-M5&)CL$^vC#N#1H%fNgm__5AXpm@QZU&mR}fK z(6#xm<#>molLlQa<;9;mGn>+vy!drLTbBs~%JY{)?#`KB(3^kA$n|9ULT_GgL+ZWx z9zMMG)@gly9PYzQ4)W~~dD(}b**mKJE5#GFaR_?c#t&j; zA>jpnp^qbf+|y~vuzgOvZT+cUQ{FrC$o3vPel6?Di{(1#{4O{@?;Y)#mb+mAUf7M> zRh0_zBd?tY_I57BKU$SO7<8x*-~O!l75Tr!&4;;lE+oae@l$uU=Xf>Pjjx~kw8Ql= zrZ(QAKgJ&;e&7#zkRL^OfDd?q|EytOW#e`wvz;HgFFWaB#oN#5)AxB=><&@jd!uJowRi80A4e@Bkn10zd69M*EHQaAtHqNRRxC_ABX8kkR>I zw7(hccf#q$=zNeKrZ$}qV}FdN^I_x=h$>#mHVUvRB6b7*F=k$e-+$Q9kehAMgUd(5E;fd(6l_ zGqTr=>^CF-!N|TdviFSaKO=v@$UiXhAB_Bmv5oiWkMZCK{^UQ5@_`5VfEW15AF^}8 zpEL58jQl4ff6mCiGV-^K{5d0k%*a18^4Dyd@aK&DIV1mWe2@MZUq|?JBmZZ@pEL63 zMm)er{+t0n#V3s76-Mz3qj-iOM~vbnM)4D)_=-_{Wo+X;`eQuB zS4RF6Um4{C5AXpm@Ke0VDBfih?=p&a8O587;!j5LE~EIAQM}41-enZevMVCqWfbok z+jx)u7{6b{yNu#pHd(~iM)?%)GKzPNc!8hdcSi9%qxhatyw529XOw?1$`2Uj3yksy zM)?P${DM*b!6^S=lz$lCQ~qHb4}RcJ`G-+H@Bkn10zc(bELG%#jPfl;`4^*nj8Q(w zC|_fgzcI?^80B}2@;yfRAftScQ9fvVKThP2#_`|>{*(_IM)e6s^$JGy3r6(}M)eX#^%6$)4@UJ6M)eX#^%7$n@6n&?B}RVW4|$Le zJirIMz)$rbM)fO3^&v*}B1ZKiM)f2{^(992D@OGvM)fF0^(jX6D`p7Y`h^(Vc#r-V zPxUJ!f5?M;-~m3WUoqfcJoe@8yHD~;Gl)(kI-Y0;qHBrHB084n1fqwCCK1gb+JxwO zqSJ`3COVxc0rKdr5`l6fdrSqXpCLg_QFV%~=aQvlAdwKTieB+J1H2ZXZ+gCCzkY#t? zw%Ht{eve#xcT`b4??C+)(YR*B=P>a-O>(-BT(|ZDM-bs!L^z##3fw;@1vO6GVlUYx z3!2_!XX!(Q9Hfy%+g0g#>0oyU=_pa}gP*$8%yf{N5PfmAw(IUE4pP=@g3dZo;m~dC ze9}px=|t}lElB-7(YQe3Q-kPKl9NGnAmJEC^e*9jKs1kQ;TqE|9i%Bg_pY6q#X<7! z=Uw8#GzY0($lLx4r#ncc<_~XdRV1IZ)G5MxhD;|9m}$0qd_F1Xkl9gXCghW>D_rU~ z?M6Q7=QlfiPRq2jd+7n!Z5^fQcUpIPCDS@Yx6$(v)K8*ud5F&h;yZxkgp%BfgrhOx znnO50x(M9*y;IT`ggQuji3Zsg+%nN7pHz$J#?lE-mhgPiMWVgCjP%GG;V6wHdSUs= zzLR!3N^^;ZJS#UP?Y)efXc*B@qMNDT4x)pI&k&+JNX}8BA%vqQ(SC$;7*U?yBFDn6 z4${MQ9UdR*<{}F;?wzi|xtk8wvb~&7+ zT@xl=bF+4m28LcJmOjZz+E(#c^N%tueX3#bwkuAOe`xr3w^S6*Pf@?gG;R*@=}&yW zBRP<}hj6qfTyF{IYQkN6!AkRKpsT*i2${XxQF3h{=<%^O zEsrj8l3oz)Sn=ePrQbPAONst`%&Gb1#m>?d4aIZxA5HzQ5}(^dFOs~6M6C!%UZTl_ zvpZ3J%^|h#Zgh}N#uq;rw%I|_SzoQu_iaAu(wu~yH{a!x=0DHuc)qlg^j-fqx{_s` zq?b`H2V;GmrJR*Ux~Eofmb{mpDEj)fv-F|<%2CVSI7=ID9DY}^jf>RAH)~oq6~*(7 z)GwLFeX1z-J!Ky=2o z5BYKpaFU!e13$KN)=K>}3aOP?rZY=Y~5<7UL{x*J|P@NiNfi zRoLkybx-Ya{>UyTDQd@(T@Nf=r0Ku*XjmtUi?propzQuHT%<6|)6>elbdib`+!i&U zyQ@@rMZcw|WqN+Opo8dnR_d3fo*1{0_~?nRAITX=a`%lDI93y`vV;@3ZzW@b&)0#op5tPoywuE(kr4X$8CvU zbIeup8ZYQkYTHr!C+fGK=yKw-f~XtG$wPE1;Q%d6I6a9@*qmCU+XDwFsmsj%E{_}} zkBg~sHwHLL9n1oTEgj@2rFk|Q=$gw}@_U^-&ow(|Dg6BEZi5nBq>7iK`sD8HBE{c% zP;f=6tMsf%gNC~CD!){!pqp=7neHzjXtT@W`BLgPfyQ|gp9tb>O>%xB zx#@%>8{uk7IF}Lb@D-2uXMgMRb$SD^XFx@MvIU;iiP|37{I;~W>$*dNP(ll`eq_NVt3?2jAaTuFBGKi&Q_@?)pI z;D6_lA8hnD`QJC>e{+)m%}H{uliY>me|^8;fB)0tpT8P^Mo^r%_bcPilN5hWq4;w) z_z>R)6n{>jxK|s0zM=SYEybVz>G6Lq%D<-l$^2_7kIkUe|r9_tv@(Yo*zVc{~y&KzNb3E%)eNFSn`GX!+EMd#8dsjm2lMjLjB=C zz5X*vQwPF+LR|vq8+Hoz7I_=|AN)At5yVx93(>}V^v8JkA@GMh$j3PaKJ1$*g#SmP zsAHiXgt`jq6Ug(CFC*_n{)RjV`2^x}=z%ugBW}TX_)+kOJjjQg0Y2cJL-?^CsKcS2 zg}M;xE2vW--$&kz{1bT`@*(6Eh|kf+d-TV6@B@FygM8otKAcMW>NAn!;1j64+i9P%RM7ii->;zo=IKk$b<$Oj(agWW-07W)cz2=zqNHx^6k(VOBLmThWAL9{6fM*Ehpe}&? z8g0Bse~bq|@JCz;`G|9X4|s7tpa<>&P)|o)8ueY&Nl`yU-469H)R9n+L0tp&0krWR z{V^W=z#sA;A9#Qdal;1MS5t}N9s~CQsH>wsjXE#trKmfieup|3>Pe`}puT}N-lIRp zgCF=q9^?ZL@FDJ5Liphiao>Xb4BQK#zK(h{>bj_pqRxmq9_nGJE1@2PHr}H@#)BXD zLmuP<5AXpm?yYc6aF2xh7~{PL-3Oq~j(RofzNnv~4vBgm>SCxbp^f+GkMZCK{*Z_J z3E%-f;O$B4q16NTNw}B6{Ra3#4&(v{>ZqtkqOOPaKpXGTALGFf{2>qWfd}}Y$M7%g zE8K_SUJ3UzSZAy|_5t<<>cFU{qArQ`KpXGTALGFf{2>qWfd}}ouLh7FIOn*B#C;y_ z#c;obdm7ji+#8_&jyg8#!KkaEK8ZHoqd&%jANV6LhkWFTzz4iI=dioDN5y?4?)7j# zhI=O5*Wlg*_XnuMqn?erFzTzQQ{p}PV?6kQKk5pQk31LnkaxlELJ!;n<31JllDOZ) zJsIwsaBqYA58NZ59*?><>cgnB;ywCfJnB&3kGcrt0}t>4FYeb7C*mF(_rbVV#r-7i z`EXx`dnerA;2s3`38>4XzKys9@6jLQQAY#+fmGLleBc2-;ML|IxQEAmHtvOSzlwWG z-1p($4EImC$H9FF?iEm4|NCpjstN!eg}cy6X0GR z_uHt`;l363mbm}JJsR$zaIb@V5w!6h{V^VOSn!8D$VZ(D_)s^QM|Kx_;CCGOJp_JN zfctsebK|}k_pZ3V#62Kz;$8}RppEzFkMXD@gFoa!KI&}1hq@2)TI?(Q4h6sG!0#gP z`vTn4>avgzJirIMh)1v<_#F*?kAmNI z;P(;uodNFaaqq3YKgN9~?geqbhBn^g-UQ>p5ByPAhJ4@wKH!}Rf8M-H{uga4>Tmb> zH8OXVO8VNtOLJa!_0yj^zhS5UQ9ph4R-IlK8(mqSKd{lw>26i@2?MWuj9ynoZ!`Sq zf&llb`pOM^S@%h+s(;&YGe5Jhsy_5ZuGeWFtLk@bv$&9EU(HBu?u+N>m(oIv`=!6& z1HO<0xw&r&9E0l$T);V1gS+>qg)jwHXbV9eM)%24$9e6%}mA`&i;I8Rm=c?=d{91N(jt|i1e>ihi9yg|U@tPL5 zq7Ktr4KG@$$Y7?A&K+B5=>(?VZ6he2qaVhBPrWLF?=L%r9LNO@;PSj9<`?UM^?{zo z{*lA$R_a(yKgBtZ!{Kbz^%-Rg3~ao;y8ijlPrbc!GX1mUL(7{rWcugx3s~np%Jg;? zo@*|@WBRS#p@z zC(s+?r%xODQ?}ep?=iQ0@a%M^U$x$xt1NGJ)e>u9lxQXI9`hA=w#=Q*^e86|;EFtIPL?IVA%4`+;062jg=LP!? z>xX#&?wW^oEO_>j>En9)ZoH$@>qn#?9un0;uP=Ij+p#wJWt>js?6-#p>OC_a=tgc0 z)JIuNx%1<(K>ep|g5o*)wZ0+7Jsc$XWa}^ZLJsr*4(wmx1a6#roFD9OtS9CLxEGa5 zUv$7fP(N4Ktn|^9f%<;+9PT>I4$^mTz4%3z`$75#2PO^oNeI^OscY3XB+W#J9Tm^f zuRtR)?vj<@1HO<0xxfKj!0FRd*g@DY*nwn8oXat(Vt-?wV4kEbLSCbYYegHm2k8S= z%&5E3BUryKcb+SEItA-(ENX4L`8HVJ?Zm?GBJ4wqG~r`S!*lezxLS-`-$3vIU&w)6 zoDbjvPT+?9hCSQ#RM-=oFYHt76U;Yo=W4zu{Zjo9{l)`@Zk}8aqW_S-uSj%QsJ`zC z?Mj9hb$v!nq3nQP>CIg1v_SfInHZMc8lHH`o)LH|$%i zKjsm5>J@r5ZPbSl{Y>Xadn^`*>YvV8^W1uGsJ`(mLGgT*i|7||ROkyn_OW9BLC(|J z!X7NoChQt;!GC17`P^SRt`YtPc6w!RVFzJ1aQ?85u^%wsRSOE-QT_844T-6t@35wK z!ODYc=nt(BbXRZj9Q~$^66*mz;CpbkIRB6fy9-?KkMQSdYlWYOpM>98U@7ePBJ<-0 zI|zFc(@>mC>~HK7%r|h4X@7o(+58%MQ0Rg6z&P;1d50Xx1rGQ@_;ti5e-%IEKrV0q zKX8`q^!xbH597cG^Pzy+MZ4f_eZ6!GeF z`w2Vtwe}PH1ojhl2lf+w0CpL2V7Fn%VLxH#VL#yq;1}R0zSe)io?u^}JtF)Vo@3pB z8-DPw`mf2UpXUwm%ka~PGvUYK*Wu^C6hB}eh(rTFdkvCrc-_<5{1_ASmI>c;zq=8zy+L$-w=l)-&wa< z#Hol|5yyfb<^}Ny&IkPT?){&~&#>coj(*tRh@Tx+i?|$dC*pSC04{C(j64AOBJzX* zU&vn&4A|C`!tpUoeU=OTYZ9*n#g`6cwkdi_!Uj`$CL5PlN=9?xOFVFzLNI;;`r z68jtb890DTdMWA*z^$!AKrgHh&Iis5&OPcUu)`I)ejX3Pzk|ZhqaXYZ>^JyAjzwcp zx554fE-mf|+sUu4RrcfOvsc|z*sBuXcE9D>>-j72`lX9>t#9VT-`iwcJN}9nuYJ8k zpPaqQ@*%4%R^_cwl5Z=PKIXtv55C~iRga`~Zahz7lbMIjocY2TsWXPovEy@U<~Z*? zQ^)b2fS@)jf1qX3wcc%6xfPSVTJ>tn4v6uL__5A{KMND`*aIP-5gvA4;A4cB*$e!v zWPG*brGhGPZ!hb;CF)h=&ztuw?HuLHf9{w4W7Rls-qmk|2@@Ub!iFMHe3fS;Ay-hA2a z+7$8mYQ=J@a%4bNqPlpnG zvq!E&afLnjtCMTvv$(tQ9mP5qI#A7tckWrQO{aAVk@11>fHF<&_zoUl4;%CJ(rYR-0dFI+ChL!Ihn zODkQ3?^$!W{f=D)d9MaFg0GiwiZk{a0Xa3Tn%#mG1xZSc} z;ug(w;|WEI_I!HSg?pEHSI3jtbCImA)gfvHQ-@{kEIH{%u?WIGotQfxn9!u&N=+3t#EqE}YbV2^Z zy!@Z4G^!eHO-O_eX%=WxC>Ot**n{zD|t&;WVeDrmdDY71Y z#dt>i*l@w0ofPs|TOpqj9#%u(V}zGA6Zl!rB);YQv8CkyZe91Yy;qjxBVK2`k33bJ zkB@uOsqvkn{NeLu-`zf1gqK^|kF{Rl&U5$;_KaLwke`iRRbcr8XYRbSv2N=_d-*@h zUhtvi`dmE8=iI)L)w1xFGq>+P*z|=_kG!)C?SMKK&jwk3-vSh5>m(_qKyD{Tt!lI`XmS~%XBpM6rKY^5a*eE){$;q?dQ=KEgQ zwzQp;gKx=fUf%Vn8DF~W)QPSSo*4C@c0jH{c?X8edSrBclssJ44_uLtcX1ehM0sW)4IlA#H+fR77KP$j_*0LoAX1MaOvO`Yj ztett`&9+%)w9d!Nzu7x&Q}4X|-d>$^(Lr|n^p^ISVe517YzH#?bPToR;kt>tt>tsE zrq!KsMO)m<%vPaw>oq%_{2#V7v;(Sjznr{T*27`*g#B(Z-6qB};>XGe{>)v-W6453 z`&r;&TLeC~PvB)<0zZ3l@VhP}CgtM|>lM3M>$(FUwfxDACvENd#Wp$hR%i2YtKAm1 z^9$zYH}+&X>%730cenJ&dwWw3-fY6`w8=kZD$Uc9o@gJrA+6uX#{{rEgNgwFvNBTm4C$S!k)+4}|=8M)NfYW@^dIV&M^$5WHKkhZA z+N)i!w_tr}JvglgUqkCj>ycTH)|=KN^SxLPwN7_$xK8*W~&WCY7()nO? zJ{X-3<9?>|!RUMh(D?|!dB9ZTywLgJbUrwp594{G^TFwSa5^8x^GfH#c%JEe@cMM# z$zHJ?WDm$*v1Vj1$X>DPWKUpUU~kA?F|tKRU&&x z_KJ5Sdr9_+&m((E_KJs+y(Rx)w8yZ|u-D{282JxI{=;bR$$v2NAB_A5gTKIfz@L!+ z;N(9z`46K%BLBh3e{k|2Mt=tX27gEXgSQ}mNdBBnCx1!)oc%)nl>9lXPW~4Dm;5;+ zf6kcj=ZyS$z^Z#qt4D5k@+W_uN&cLZKj)Px9zc8$Nbv&sbDl}@1o?B`m*Nfb=bZdG zr})YkuTXqtjAtmmViaF7im#0E5XD!F;wwh+6+=9Q^+3Eu@fD}|ic@@LjMpf>;uK$T zim#0E9>rJ2c#z^NPVp|Ic-I(DQoPG3-enZ;vN(!IDc)r#C|;#_m)%)oh<6#{Tj+uK zm*QPc@h)FX@iN7`d=SOc6z}qi6mL_!%hM?yr+Am=p?ID04`Vz}`G+yyr~HFa{=q2! zFy;%Ce=y2F808-f`3B`5{>Vot|KOB=aLPZ7`3&VBobnG&`G+wdqWr^{FH!!%{V1P$ zMfuealy6Z!$S5CVln=7Tl&?`f$i`7VNBJP5e2^jkqkM2H@#F8b58lWG2f^BoKt?zDL*&X3n)J~))Oc{=Tt9YR4-xqsUAV~5*AMN3aXbd zs+TaTm#|q>@1S}KvlR6bhWZH32kIwOFX2=#;b~NFp?V4Tp?VC}OL!F3Yp7nra|wK$ z>Lr}&C7kM4#(EIduZ;B~s$VgxUoom*8S70{zhYFsVpP9ks8>M`)U&95#i@S9seWaw zhf)2CQ~iol{mNKRqxzMx-bVE+ZjjS2oDIA?D%84LI755s+>I+1s@nVdZTw+hchmbL zQClAl=%H#a$a(Z=FH<}8`P0fHzc;meTIbx8)Zf&u-sGvL=b&)<_YlLG=S+R>K10Hp zg#tjvQ@U?RWW}^_X80Gc<#wg!A@w8DP2)#7blLD+ZogLk0l&Vc=imK&DKtujd(d8& z)fOt;=j(6i<(_Klf9qYHkV7gxuJ8Z(rFn1D`=))Rw=-3GEbO-}OO#3v`wD(m9|oHG z7qeWo(OcF7^H_RC=%7nOP3`Y1-3)s#>xK7Q{M9K1-!oGnjJ(@){KafM#!|;!h%j*R_`sp4wPE_H3>vL$rc~!ey z@0stntMtfSeN5%bDn0UykA4-P(qro4hr?H?^r+c5cS;?V9#`XMWalb9?ru-qvq#ng zcyVxXVKqx%;Ycb~m5AcbKZZb**hiL)AVil+!|* zuhJtt=Yv`UReE&y9A9&vN{>!Q5=g?&)VK9J)S&$zpbfC55HPZ+ANj#8P;di zhxV1a%6j4a$b0((zE2f1l^$Dnj~ecyIxoFDI_4V}uCVO>Jm(9DK%6=_6xBgI~O0ORIhi}oT^k}&zyF+jDzV?IG1q##UA1E?fO=e%1U(U)66)RO#_1=u%M^l|9V9+5Oy2)p?v% zaBJ`wSr5pc^~?4w*H!D&YkR3_yJWrazU}EJMJucH8a>>7V1%lTe!qIIjbE+8nQQ*l z8v9l48dp1QQ`@7mnLRh0SLso!(YTlQReEGeOuiMY(j(v45zd=bdaPg4$ttf(kDuA# zP!E+JYYL8&6nplb_WzTn1=m+l>1A`Y*7z`08+tTfQlZI86|TfRxpr(;wGY;rSmdLs z9az5T-k~boea79|)Izn7(q7Km7_ZV}oK@?P@hUwAcXAG%quOWl4xg)(s?y76Rfh{^ zD!pb@e(WElYOiawe?dJ}8+r^`U0OP*!gcz6o^B^p?J=DiI-99*c6Tg&T#b9nv!8wL zs@^|}YMC-qrAMBbKMs#l=~30j@1TQ9kMX~rD5KV^-;M)fymRPB$y9x0birRS%+ou=nkwXaxhELBL=h93WJ zsvHA#oP>4wtGHfh^fSNrx&9V+KA)d+n)y4cnP2VvYu87+et9(OYvK6$`kOm^zCYYF z`v>~=)a)-lGwoNjhiUe6kY>MY&x7{7?56W1{{8Q8mY_Lrq~F2U0$O_Ek;3uvn8%q-*kx37UK) zT9dEXY4Vw$Dc^y8N~C-U_6zwE&JS?)*W_2(H2K$YO@8)BlfOOGl`7-ip^q;E9x4UTa@gz;Y-j?!t*sHvh z?_>T;YU%-jntH(kswZgcCCfGShu)g{L@`bM!c|k>D5$A_XzL@|`pJDwedU~{{&G@N zpBbR3-vn#wJDATtR1ZQskEUMaNA)DUFHiL*w7sbwh4w~Gy=sW2p0$JOU9eBbsUC*? mf0pWHnD5=1dRh|I+fXk_?i3Xt9~m8&6dfOH!0`XqAO8;keCB5W literal 0 HcmV?d00001 diff --git a/tests/orog/data/C48.mx500.tile1.nc b/tests/orog/data/C48.mx500.tile1.nc new file mode 100644 index 0000000000000000000000000000000000000000..51b5fc0a40eaee21145fa6b33a3dc215d82f0648 GIT binary patch literal 55384 zcmeI13vg7`8OQHtlMq6H@J7nxij0K;HX#HEA$u=jvmszeBtbxE37hPOEbJpD3mD`P zM@vPhP_T%V5zz`Z0y^a(uNIh5d{Y4jkSQ%Q7DfjZP#9}%rRUsxKlZ^s*xhq)vYqe3 z?sw1Y`~K&59(!}{s*SxNa%e%yV*1wo*`k->PZlB}ElU9gpae21=Rw zRHY2s1lsf@PDiU`X43C=diu?A#j@czXim}@U(ifFmqz|2b6lOh(c^At;<}N{56LHq zd}i3J{A7!Z=kxdxBlAa&%FoT?bFyYX ztap3YTz>auxEMa6u1y{*AQFf{UuYz>Wew~HTT$Tp1T z2bZ{=wtCAv8*k+6Jf8aetSncJ+t65RAPtA{Cb!k6WS-qqx3Jn!<95rR4lNjjYl1{$osI~}|(_YR-~LN1=`mB#g2LGDtdVb!P0!`l02QZct)He8Rs zQUiTEE^sCE4SXn6npwom7bf_ZK`4#nbX=q(hd@2`(i;-=D=s{_c1)rpSW?cT2m*_Tyy*nsTc@c7ysLNVKiFP7CJoGBd@0q5_ z9OVT$U@Onj@nMJ_+A{8SN;I`6m8wgYwsdl*ceq`9l}VY}TP{gBTWB2YrzQ{%Iebcng z>BhaICz{GjC5_Y>xuknym7(Q+R5HLP62Y z>DKLhs4sPbd}VirJG4ju2_OL^fCP{L5qBG-_8>ZdcSAjDl*_6+FBjxPgxqI2S8HX*Dmust3IRfNBDWpLuq(Ld92d5ohJgvfP))CsU zw}9maL4ShVnI-=D*-8>D(Dbi`DI|adkN^@u0!RP}Ab~_B5c6}R!@tMyya*Pb8~yoU z3JD+qB!C2v01`j~iB2Hq=f-vZJ%{H-u=w2Q&j(XT00|%gB!C2v01`-a0x>@~?(y$A zJTHR9=SF`%m_hk1tCRw1}lnjAn5~@lvsA^kA|5u?a?XVzBW# z;oQZsT}yHTSkpaJ{&YEetxf9x`I=3}f0v#Y#-5%k_P=BaIQR7Ge{6f((Z8eJo5mMs zye<57(-?8s*Gs&NuXy!e-BL=%!mEE(adN=euo?GgyllRIhj3|YsknR0^4pqyHU6jB zCOX}_g)!M?anF}4{ObI3zEhmrua{WyrBe`FCXS0HE~~nGyLR`}ohj{pS#dH`YUdAW5R##+pVMj+D$h4(|`G(SaPC+IC$w(ny%?nU$M~-Yvt4CYBHYc zULz`gDtO}PZv^L!OTvcZQ<=;k(0_5M`ueBZa{cykh|R(dY*Jr;%|q(x-?C7B{guV) z>7U)8zJ6n#diqDFyTtboY!wc?l_$89N2v0fHg>3b`hC~XZ}-!^js0aAppE|W-RkK7 zsJk}$5A9M%|Jn}P=)cgSj{ZOEwb6fQyE^)7|2vZT?^pFJeFe*oE`o98WMM|rGNDuM zMxl4n1XZ4ep8pqn{f}HJ2v`fJ@$5ZX z(m4J0Z^+h3IZU55tCy>6L9q1|ZQg;ePBPyQL-!+Eb)oAoKk=Gz zh5d%{qT|Pe)bFmRYtw(ds-F(cGujffJF)eD>r8C@iyeurzuT$@#D|`$6R&>a&~P0@ zJIf>2zj;0}_fP2AjI)WYzw>Jo1n2q7!Vm5X8vG6bHMOkTzXxjR2EU1|f6X6tLVjj% z;mVZ3D(y+YNnQJR|6X^c|JM3K@q^FIBDdW8&My2jh`PGLd)JS9-j)9Qzd0-YjeYTj z9OKP*HVWO1y~Uy<)-Y`$g1WZW`R+J=-_yq2vt5Of@y$Z&l_SE2(@El=RvN_Yl-Z(b zwfw&Ex%0AqUp3b=*Nvm^6zRJ~`i{{%9=^3&_x{Vk7CbH+yFb;hE0V$F$l!sD}E4*GYM{Qf!`kFnRCKcttq>EI&q*K;i5 z-};-xA3F15?$Yms`7IlSozJ`<^xRxYKle?B@tM=a<+K)m;* zgTk5Pe~|i9+al|K#ZBkDeCfBZHFnzYhEO=y6?9CX9&5XMczFHZHQ6jaZ&)T?t}=wC};uhm*N#EVxqvR=oUh1AVyh7@W-Mi~sffK*EE>bLiCsx- zEyHOERz{-)%ZI$!RK(s2;@+(dra$)Rp<4#i1GR&d#ghkqS=wNFz<;o^c=Et63mcD) zXxg&S8q7yD8se!NEH;b0cr-=RmW5XE5lviHbs;tj%Sux;epu)RAJN2RRTpBju&gvi z=11Euu)db*KngV#gQT^*IPi%Jp)$K3{X* zwRZR77*R?Al>`VV4T?~!f-I4?sah$ZN>ypZ$B!To2_oX-2c$wRKOzy7!XF@1YI$#F z-udct!r6)~B=*^NZ+G_1o0<1H_jZo=#S*I;H#F8YG%$h>B>bvFv;^gYjsAZ+Lw&K4 zx>YZzbz69&xTZEwl(?~uDcG`R4Srwu;WuM@1!CG(yil12fAA=6X4}(_ zYZtOETLbeE)H|saWp3@)C6(-V?q1MhI zej>VSIBrM4%<_fvY4X_@tts1ZpuMq*WpKwlEG55o&M$)eh;fY!a8BW1Fxc8242F2; zcH)BUYUEO85#%hw!~%KMz@i|96YdUmceTRz4zVC&gZ$=U(>+HZ$Od~$9^pB%Z4aVg z`2+zoZJIGZsMRObIv5?|Dt$~EfoZtWSVA9~H0%uM+tI<29IzCMecNgfW@P z8+wl}Ez!(c0cR`#9b5INm8*D?FSz5SLxHR{(FR2L0f)K+29of?)VWuL$1+7(hmCcL zqNws|8Lva{x@+^Mvcoksu%;)0%Fn%5Z4}-+(;_-PBm3Xe%sZ%lMiFB!);@-$Dolx6 zcr>{?mEyv{Xn^(|>)#(AO!W^9rc=?~$+!#?9e(V_*tkFp*`V95mSJMhaI?9wNUShX z$W;Xj-eu)-1vD#m>ip~%+1n@Y;d|BNmb4mM<%T+xS0u5McH?4u(YyO23lGy*z7J6rpi4KG)#j&w^{ ztxOERk{^$NRzTDkgO;%kU&fu+}4Cp*>p zm4<4q$*d>c3ml?GtqpAx>zb=I-%U+-kq-AgNVW)z(i&PjxMmgnY4gZ8su~y5u=I+H zu=r>9{k&Rn^>Vs|YgJsAj2%m^xJMhs3>T$19{=Z2+&Aw#e_}Cus_L`!itA#>bASI) zS#b~G6zU+yQEDs4NjbOjy@n|d_n0E8H!y!aUaH&h{wa|G?$t!EZ^_P-Ic|~f;z+8p zb&#ka-KJ%tidvSPE0`JBLVY-0K{Dpv1DyyX1Du0)zsN5#b(YM6>+l4Fu8PE2V3)`J zqD+?1)aS1cbD8fDedRV6wwn30nVAq<xqOnf^b>scu>*1^ z6Z-Ev8T}zh$_TwIvXVtU0Bi8C7U~u4`0P7-8`dILm4N)6BNxg@&$kG(yLr!r=y^VI z2m#+%h0GW;1!7uK2*kL{40V& zDPfr#7pTCoIW>l5Q6TzV%9g@Sk>!*>Y9OqBsm@V))A|#khOr@S-&$ z)kzgAs1(h@xlfP2m>B6B=8bhch{}ibvlV_-aG;+VJc?lpdVB>?(W|0q-b--ca!li|A-%%vnr013D8FptLe=P_9E6ono|oD$ zM^a%eW!Az9%I+$Rt-M(i_B8d(wJu`r;boBkrYhPdH_KGbcg--mb`WEU3K6N+eJWWB zVT?btPsQuXd!+f-nb)`OjSeMwoZ>GKC}xShFOaZzq`eVOvYF?lrI9&OJJFL}QhQ3; z9vLXL6CbivYS$<{<0JqclUqS?=Y4w;@Fd_#z>|O{0Z#&+1Uw0N67VG8Nx+kUCjn0a z3z0x&1%m3d^1UoC1*`1BR>AM9#P$<(ar^Q8AHC+xHzRs=h-GM%?UJ?e<(Sd$fBT!~ z_zQ*4Q8gxrZw*J{iTFr-Fd9$uw}zci9>D?%9$<6~clPyR`BV~Xz~`bWqVr>7R8b*( zL7p`rslGyLkPWA=KCx>)*~r2T(~9(7^LEXU*M(qo!udIe8P5Ikp-*zzBx&SYD^urLD@3vWFTA?EQm2`;kT|8(>;CcwyHe_>Wt9$I?Q*!iZgN9_((xM zJS|F9L&ZnlQ!fwdQk_xpQDw4To>HZXqT-_=R^Y*k^hCxvOunc1s6MH}Yk2ZN#Yff2 zMrcio^Bm+1ijS)5DxB~j^HF?M$ZUjG*=4R=am^lh67VG8Nx+kUCxQPz35c_G{2R~N z9X|%J?%`f82*5|D0PcRc_tHfGxQ`(KKK~{Fw6`Mwe(7ZZdqA{11HDLo(H-c$@)kg8 vp!cuW0gem^zDEao^(_EEAHn`W&-nuV=m7wr&y4{9Js*E+pckig4Ep^CnyjEv literal 0 HcmV?d00001 diff --git a/tests/orog/data/topography.antarctica.ramp.lowres.nc b/tests/orog/data/topography.antarctica.ramp.lowres.nc new file mode 100644 index 0000000000000000000000000000000000000000..1023adf6e562af9431c42ebcfedb8a023ead7a22 GIT binary patch literal 13633 zcmeHNdu)@}6+ez05)u;^vY1FAwxp=@+dHmxmmmFRO~bXukkt(sa@?5}B6I<=ct@y|pz)UvgMopbL!Nr;m* z$wD@1-%0$P``!0FzjN=o-?{f-n?JCmY-L$dX(=#WQ{-)wG4(3p)5i|)>~8aiik5ha zifToqd9g~1SINLfimzHG-^+@?M4`$`dLA&-6QETxux_A$8Laqc0qD)doV|HHmyGF9 ziT%MSBn^u2Nwrva5;sE$lu}wgJCuth;3iDjhCj*`^7_4@us5Q4yCU9@FVgPQI=#VQ zdsnL#40Yex;)`f?A#Z1JgBIxycDIJS!M3|Kzc=EwZl(&%P%f(;%w_LS#Pnn~BdTAj zs?R0*hf>j;mWt||Q`wC5Gvoq@ zYp8E%+^%I}59GDm;#xAEOxqkb%%EO+2l373aM1hxl-S^My6DlwF_eQNw`qEIFx#Jt z4h{_4qZvJ#!*0i-_FObQXm{4-g_=^Rc?feA z3%xRzNxeF@P8NJZmA`qFC}1uI*5%3Ex{0!yu}eY}z7B6VEQCRNN}A+v?`-J`w|95- zhJD_S7Nx*%D+L~)(nwSdsZE7a;I4!oi}$(w$#gPapu(a=g?L`~aOQ(6x_xb*=Uy+G zw%2mB>qOrH!d#1(VbFf^!@p5A)_deEDOMS_ir*# zcQFtSw?MNpR3@tlE>pJkzr6BDLH21PZa8X=dN@^GQeB8zx#|(~nWA=e8r0@}x_w+K zdpIZ5#40-lYE`#C(J&1Zvx~m*z&YUT`}_S8ILBw3IL(c#HIWBg<>0!85-6@9<7eKe z1c_2x>+Ka!$Z(%khCk;q(A;@>$nTBT*6(Ni>~60P+C)+N+n&1;A`BGY0@a z`u^F0LIBysF`pZN7ZiX&C4vP3kVGjFSg+Q_tZH7yE>kPP6afBiE~l z74obP9MAQ$VsJdy!*bZA%sIas+>n4F0Yd_Y1Plon5-=oSNWhSQApt`Ih6D@=7!sJd z1jtp*2c7ICRjxjk2~!xY93hVE*~H&z~(5&hw7>JmVzah3&ULyFArcSExZ) zQ`}-bovMaf0xhAIE?-NpI91((Rgi-u9@V4ygg$<}I(D50N-14>J+vuLU03pPifZx7 zraycz|6?V8OkYv@Pg9MxPPI%9bu)Zt0dYP0$M(^wdPM$3q8l&$b#0 zTZ=y$p1L);qrryv>+m+<4wrqKy}>pR?`z!Qusd9KXZ&-w9}cwnn}Ki5Qp>MgZi<@N zMF|$V`N@t*lZ87|LeY5rlUFL=tZ}v*Z;D+PtdDw(Z8y% zD}KZ+>p!91PpWo)tm03p_MTViFQ|53Qu$s{?f=yMvi6Lpzvq-^|HIFDp!tU$`0F=4 zaQu)5_6~VK>+rytT2E^5n)}(|tL|t2Ug@cPsn)Zqvf2~wEAqVE`JQ{%lRt1b1)Xls N>(zJ2G5mL){{mapfdK#j literal 0 HcmV?d00001 diff --git a/tests/orog/data/topography.gmted2010.lowres.nc b/tests/orog/data/topography.gmted2010.lowres.nc new file mode 100644 index 0000000000000000000000000000000000000000..80c71b1a5da21da3c9c19feebc3efb7a9f05d7f4 GIT binary patch literal 13565 zcmeHNeQaA-6+bU_KI*z@y0s0BHuhFkS~iU9*iQ2yRXFuW+{R7PIGNU|#dGZ!$HTFM z{n9oFt!NV3N-3xmRV-thn)0!ATD1wLfm&@7(k9r(A43xof(#)LAJSBYU}#XmIOpDb zoHkB2CnOOP?<9Wr`Mz_{J@5SPx##-adomh|Y;d$VY8o1V=>tXo^2-{JdhpE?ry2*M zp`n@$_u6WgG?mF*-fPl>Y;(2l?s!^lmXP!$MHj4Bk|t2-n)M|9P)bFJo@m! z;6U$C&tUYR9_ktHaqc888?2Y@CyIspjg*-!2K0~1CuiRhQ~qeiiW z<3?|Hx4C!etHGw-+E687ieRiK9v8+S52n!`+Pgp8AKyFBpNI#0VqxXle^RbJ zMy-)Z4Eb84a_xS@Or=Ntq3n1zUE$oKMWeXha2eTCL>-eA472r6OjJyE)$fT*~3+CLTqsd}6 zc_e4JbtWI&^y!w1mCz}}D4@`^RYEu2D&Jy0boND!uuSM?)aY9lbi%?5K{x#1tknO% zyskcQeh7bV-;XPSW2fY@THstFssYRgZpD4>9*mJh%L0xDJbGsZ8Sgg4BJp?_4yh2X zIW_nju|0F);2$d37pYi%)ONvH*Fz1JsMRZtxU3en7Z*Wo#ocX_Lu$Az)WrP10BZKh z%P%bg#nP%*AGi)Unf}X*61cXH4srUYaJD82fPX#sKR^rAHd64j-KB;!g!ldkt$dy6w}#!5|?T>#w(1*GOWl&sp?0C!B^MI z?|D`2jFD&?Y=kBfJ{LPWmIeey=P3tKmm%~fmG{hn|Yfh#O=XM#cRn(u%XU&p`jpRnu zGHqmvhEaOA{m{QFeQuEuw?ya;DCl(Tj%ZDwBa_O-!$TWHufWjB+CZzU#ifI!TE2mk zXHQR)bdjwMJ!6;ge*F4AOIzR6K}pyw;TWq#r!e={QDC>9INy_4sl$n zNppz+epUgzs*K%hnq2oBmK;F0_a$FxwzrAZ8xX(3`avyOMDc_)2)}(B6wZ7m-kQI z@qNKnM_e@JAlDWtV~%U&{ddM9ag+u87V2eG&)g@a9R}ezYC;Ai+NxGHG;a`bIey^W zTYH9M5vXmyGkP1)kE>?*#sjfdxQMq~;mcPOoj^mp1v(zyF$EF9@Y82!Ps1_6;Q!|@ z&%(=idn=q7FHHlzcMF`Zzx53u^_$_POzurMEf`*j&%S~;7X5EoK;X~Q0=hVA*E~0% zB%tLXE(ZWlE^27dR47P}iSlAb!d@0FhZKhtHEPtr1U2?v+ji|7OS8%~u|RKdASPxB zHCytW!pNLcX!980`U47Q9t>RXPzdv=;QE)<4UY$|XZbuZxSr*(5$;nfG5&6qmIW*e zSQfA>U|GPjfMo&80+t0V3s@GgEMQr{vOx6~pj0()bn>|!r)o}FU>p2Q?LvKth2G6K zD__oJXGcZyJhqa1oRqsbd{6Ybh2E8u#LgV|=%HQJp>QNT6z&g(6Jl3&7m6UxbkawLfMT1cZ?l*X%X zqi5%3c?w@S{QUJA{X*?2?O=@#F>2NBg>xQ@OUifgUO45ZuK!LssZ+_=71F@@d3l^~ zBdpZ;U=6sOa;TS5MNw^FajpL9>!r2PBK9i3pG0{Q?CxfazawGvq9~2ZjM1_3Bh8XCyT{ zcVWZxxr~Bl=R~#xYZEE4J;yGr**T4A!$T2lF3rwKOb6DgCOo1wJ7+fSv}wgArPs8` zvVdg)%L0}KEDQWUTj0-EkzN}NWJZy&e0vt@CxZbW(d&bOr!OL5{dOWOryoJW=TB07 zu`h7tGE%88@CK3D7x)#mpX>{~LvoMy1%z#n_XTunkNkd;i~J!X_%7bfz EPSILON) stop 6 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum2 diff --git a/tests/orog/ftst_get_xnsum3.F90 b/tests/orog/ftst_get_xnsum3.F90 new file mode 100644 index 000000000..75a0a8fe8 --- /dev/null +++ b/tests/orog/ftst_get_xnsum3.F90 @@ -0,0 +1,76 @@ + program check_get_xnsum3 + +! Unit test for routine get_xnsum3, which counts the +! number of high-resolution orography points within a +! model grid box, and counts the number of high-resolution +! points higher than a critical height. The critical +! height is passed into routine get_xnsum3, whereas in +! get_xnsum2 it is computed. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : get_xnsum3 + + implicit none + + integer, parameter :: imn=360 ! i-dimension of high-res grid + integer, parameter :: jmn=181 ! j-dimension of high-res grid + + integer :: j + integer :: zavg(imn,jmn) ! High-res grid terrain + + real :: delxn=360.0/float(imn) ! Resolution of high-res grid in degrees. + real :: glat(jmn) ! Latitude of each high-res grid row in degrees. + real :: lon1,lat1,lon2,lat2,hc + real :: xnsum1,xnsum2 + +! Variables holding the expected test results. + + integer :: expected_xnsum1, expected_xnsum2 + + data expected_xnsum1 /4/ ! Expected number of high-res pts in model + ! grid box that are above the critical height. + data expected_xnsum2 /16/ ! Expected total number of high-res pts in model grid box. + + print*,"Begin test of routine get_xnsum3." + +! The high-res grid is a global one-degree lat/lon grid. Point (1,1) +! is the south pole/greenwich. + + do j = 1, jmn + glat(j) = -90.0 + float(j-1) * delxn + enddo + +! Bounds of model grid box - straddles equator/greenwich. +! There are 16 high-res points located within the model +! box. The i/j range of these 16 points is i=359,360,1,2 +! and j=0,91,92,93. + + lon1 = -2.5 ! in degrees. + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + +! Initialize high-res orography. Half of the points +! within the model grid box are at sea level, one quarter +! are at 500 meters and one quarter are at 1000 meters. +! The critical height was chosen as 700 meters. So, +! four points should be above the critical value. + + zavg = -999 ! Flag value for sea level. + zavg(359,90:93) = 1000 + zavg(360,90:93) = 500 + + hc = 700.0 ! Critical height in meters. + + call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn, & + glat,zavg,delxn,xnsum1,xnsum2,hc) + + if (nint(xnsum1) /= expected_xnsum1) stop 2 + if (nint(xnsum2) /= expected_xnsum2) stop 4 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum3 diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 new file mode 100644 index 000000000..c9f12e99e --- /dev/null +++ b/tests/orog/ftst_inside_polygon.F90 @@ -0,0 +1,223 @@ + program inside_polygon + +! Unit test for function inside_a_polygon, which determines +! whether a test point is located within an area specified +! by a polygon. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : inside_a_polygon + + implicit none + + integer, parameter :: npts=4 + + real, parameter :: D2R = 3.14159265358979/180. + + logical :: inside + + real :: lon1, lat1 + real :: lon2(npts), lat2(npts) + + print*,"Starting test of inside_a_polygon." + +! The first part of the function checks if the test point is well outside +! the neighborhood of the polygon - i.e., a gross check. There +! are six separate checks. The first six tests are designed +! so that the polygon is far enough from the test point that +! each check is tripped. + +! Test to trip the first 'if' gross check. Is point well outside the +! neighborhood in the plus 'x' direction? + + print*, "Test point 1" + + lon1 = 90.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 94.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 94.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 95.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 95.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 2 ! Test point should be outside polygon. + +! Test to trip the second 'if' gross check. Is point well outside the +! neighborhood in the minus 'x' direction? + + print*, "Test point 2" + + lon1 = 90.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 84.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 84.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 85.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 85.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 4 ! Test point should be outside polygon. + +! Test to trip the third 'if' gross check. Is point well outside the +! neighborhood in the plus 'y' direction? + + print*, "Test point 3" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 354.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 354.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 355.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 355.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 6 ! Test point should be outside polygon. + +! Test to trip the fourth 'if' gross check. Is point well outside the +! neighborhood in the minus 'y' direction? + + print*, "Test point 4" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 4.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 4.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 5.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 5.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 8 ! Test point should be outside polygon. + +! Test to trip the fifth 'if' gross check. Is point well outside the +! neighborhood in the plus 'z' direction? + + print*, "Test point 5" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 359.0 * D2R ! Polygon. + lat2(1) = -6.0 * D2R + lon2(2) = 359.0 * D2R + lat2(2) = -5.0 * D2R + lon2(3) = 1.0 * D2R + lat2(3) = -5.0 * D2R + lon2(4) = 1.0 * D2R + lat2(4) = -6.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 10 ! Test point should be outside polygon. + +! Test to trip the sixth 'if' gross check. Is point well outside the +! neighborhood in the minus 'z' direction? + + print*, "Test point 6" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 359.0 * D2R ! Polygon. + lat2(1) = 5.0 * D2R + lon2(2) = 359.0 * D2R + lat2(2) = 6.0 * D2R + lon2(3) = 1.0 * D2R + lat2(3) = 6.0 * D2R + lon2(4) = 1.0 * D2R + lat2(4) = 5.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 12 ! Test point should be outside polygon. + +! Test the case when the test point coincides with a corner point +! of the polygon. Result should be 'inside'. + + print*, "Test point 7" + + lon1 = 90.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (.not.inside) stop 14 ! Test point should be inside polygon. + +! Test the case when the test point is inside the polygon. +! Here, the test point is exactly in the middle of a square +! polygon. That means the computed angle to each corner +! is 90 degrees. + + print*, "Test point 8" + + lon1 = 90.5 * D2R ! Test point. + lat1 = 0.5 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (.not.inside) stop 16 ! Test point should be inside polygon. + +! Test the case when the test point just outside the polygon. + + print*, "Test point 9" + + lon1 = 89.5 * D2R ! Test point. + lat1 = 0.5 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 18 ! Test point should be outside polygon. + + print*,"OK" + print*,"SUCCESS" + + end program inside_polygon diff --git a/tests/orog/ftst_qc_orog_by_ramp.F90 b/tests/orog/ftst_qc_orog_by_ramp.F90 new file mode 100644 index 000000000..06c4ad4cf --- /dev/null +++ b/tests/orog/ftst_qc_orog_by_ramp.F90 @@ -0,0 +1,98 @@ + program qc_orog_ramp + +! Test routine "qc_orog_by_ramp", which adjusts +! the global terrain in the vicinity of Antarctica +! using 'RAMP' data. +! +! In OPS, the global data is 30-sec with dimensions +! 43200 x 21600. The RAMP data is 30-sec with dimension +! 43201 x 3601 and only covers Antarctica. For this test, +! reduced versions of both grids are used: +! global - 9x7; RAMP 10x5. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : qc_orog_by_ramp + + implicit none + +! Dimensions of the global grid. + + integer, parameter :: imn = 9 + integer, parameter :: jmn = 7 + + integer :: i, j + +! The terrain height (zavg) and land mask (zslm) +! of the global grid. + + integer :: zavg(imn,jmn) + integer :: zslm(imn,jmn) + +! The expected values for a successful test. + + integer :: zavg_expected(imn,jmn) + integer :: zslm_expected(imn,jmn) + + print*,'- Starting test of qc_orog_by_ramp.' + +! Initialize the global grid to all ocean. + + zslm = 0 ! water mask + zavg = 0 ! sea level + +! These global grid points are outside the 'ramp' grid, +! so they should not change. + + zslm_expected(:,6:7) = 0 + zavg_expected(:,6:7) = 0 + +! These global grid points are located within the 'ramp' +! grid. For this test, the first two rows of the 'ramp' +! data have non-zero terrain. So, rows 1 and 2 of the global +! grid will become land, located above sea level. Rows +! 3,4 and 5 of the RAMP data are ocean. So, rows 3,4 and +! 5 of the global grid will remain ocean. + + zslm_expected(:,3:5) = 0 ! ocean mask + zavg_expected(:,3:5) = 0 ! terrain height + + zslm_expected(:,1:2) = 1 ! becomes land + + zavg_expected(1,1) = 5 ! acquires non-zero terrain. + zavg_expected(2,1) = 5 + zavg_expected(3,1) = 5 + zavg_expected(4,1) = 5 + zavg_expected(5,1) = 5 + zavg_expected(6,1) = 4 + zavg_expected(7,1) = 4 + zavg_expected(8,1) = 3 + zavg_expected(9,1) = 3 + + zavg_expected(1,2) = 2 + zavg_expected(2,2) = 2 + zavg_expected(3,2) = 3 + zavg_expected(4,2) = 2 + zavg_expected(5,2) = 2 + zavg_expected(6,2) = 2 + zavg_expected(7,2) = 1 + zavg_expected(8,2) = 1 + zavg_expected(9,2) = 0 ! Note: this 'ramp' point has non-zero terrain of + ! 0.14, which rounds down to zero. + +! Note: The path/name of the RAMP data is set in the routine. + + call qc_orog_by_ramp(imn, jmn, zavg, zslm) + + do i = 1, imn + do j = 1, jmn + if (zavg(i,j) /= zavg_expected(i,j)) stop 4 + if (zslm(i,j) /= zslm_expected(i,j)) stop 8 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program qc_orog_ramp diff --git a/tests/orog/ftst_read_global_mask.F90 b/tests/orog/ftst_read_global_mask.F90 new file mode 100644 index 000000000..9672a9391 --- /dev/null +++ b/tests/orog/ftst_read_global_mask.F90 @@ -0,0 +1,38 @@ + program read_gbl_mask + +! Test routine "read_global_mask" using a reduced-size +! version (6 x 3 vs 43200 x 21600) of the umd land mask file. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : read_global_mask + + implicit none + + integer, parameter :: im=6 + integer, parameter :: jm=3 + + integer :: i, j + + integer(kind=1) :: mask(im,jm) + integer(kind=1) :: mask_expected(im,jm) + +! Note the routine flips the i and j directions. + data mask_expected /1,1,1,0,0,1, & + 1,1,1,0,0,0, & + 1,1,1,0,0,0/ + + print*,"- Begin test of read_global_mask" + + call read_global_mask(im, jm, mask) + + do j = 1, jm + do i = 1, im + if (mask(i,j) /= mask_expected(i,j)) stop 3 + enddo + enddo + + print*,"OK" + print*,"SUCCESS" + + end program read_gbl_mask diff --git a/tests/orog/ftst_read_global_orog.F90 b/tests/orog/ftst_read_global_orog.F90 new file mode 100644 index 000000000..0f196243a --- /dev/null +++ b/tests/orog/ftst_read_global_orog.F90 @@ -0,0 +1,38 @@ + program read_gbl_orog + +! Test routine "read_global_orog" using a reduced-size +! version (6 x 3 vs 43200 x 21600) of the gmted 2010 orog file. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : read_global_orog + + implicit none + + integer, parameter :: im=6 + integer, parameter :: jm=3 + + integer :: i, j + + integer(kind=2) :: glob(im,jm) + integer(kind=2) :: glob_expected(im,jm) + +! Note the routine flips the i and j directions. + data glob_expected /161,162,163,167,162,162, & + 157,153,148,166,165,162, & + 169,163,155,169,170,171/ + + print*,"- Begin test of read_global_orog" + + call read_global_orog(im, jm, glob) + + do j = 1, jm + do i = 1, im + if (glob(i,j) /= glob_expected(i,j)) stop 3 + enddo + enddo + + print*,"OK" + print*,"SUCCESS" + + end program read_gbl_orog diff --git a/tests/orog/ftst_read_mask.F90 b/tests/orog/ftst_read_mask.F90 new file mode 100644 index 000000000..e7d9f552f --- /dev/null +++ b/tests/orog/ftst_read_mask.F90 @@ -0,0 +1,82 @@ + program read_mask_test + +! Test routine 'read_mask' by reading a sample C48 global +! uniform tile file. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : read_mask + + implicit none + + character(len=20) :: merge_file + + integer, parameter :: im=48 + integer, parameter :: jm=48 + integer, parameter :: chk_pts = 4 + integer :: i, j, ij + + real, parameter :: EPSILON = 0.001 + real :: land_frac(im,jm) + real :: lake_frac(im,jm) + real :: slm(im,jm) + real :: land_frac_chk(chk_pts) + real :: lake_frac_chk(chk_pts) + real :: slm_chk(chk_pts) + + type :: indices + integer :: i + integer :: j + end type indices + + type(indices) :: idx(chk_pts) + +! The i/j indicies of the check points. + + data idx%i /1,29,47,48/ + data idx%j /1,25,23,48/ + +! The expected values of the check points. + + data land_frac_chk /0.58463, 0.8377, 0.0, 0.415374/ + data lake_frac_chk /0.0, 0.0, 1.0, 0.0/ + data slm_chk /1.0, 1.0, 0.0, 0.0/ + + print*,"- Starting test of routine read_mask." + + merge_file = "./C48.mx500.tile1.nc" + +! Initialize output fields to flag value. + + land_frac = -999.9 + lake_frac = -999.9 + slm = -999.9 + + call read_mask(merge_file,slm,land_frac,lake_frac,im,jm) + +! All flag values should have been replaced. All fields +! should have values greater than zero. + + do j = 1, jm + do i = 1, im + if (land_frac(i,j) < 0.0) stop 3 + if (lake_frac(i,j) < 0.0) stop 5 + if (slm(i,j) < 0.0) stop 7 + enddo + enddo + +! Check some points against expected values. + + do ij = 1, chk_pts + i = idx(ij)%i + j = idx(ij)%j + if (abs(land_frac(i,j)-land_frac_chk(ij)) > EPSILON) stop 12 + if (abs(lake_frac(i,j)-lake_frac_chk(ij)) > EPSILON) stop 15 + if (abs(slm(i,j)-slm_chk(ij)) > EPSILON) stop 18 + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program read_mask_test diff --git a/tests/orog/ftst_read_mdl_dims.F90 b/tests/orog/ftst_read_mdl_dims.F90 new file mode 100644 index 000000000..a072e9652 --- /dev/null +++ b/tests/orog/ftst_read_mdl_dims.F90 @@ -0,0 +1,32 @@ + program read_model_dims + +! Test routine "read_mdl_dims", which gets the +! i/j dimensions of the model tile from the +! 'grid' file. This test uses a global uniform +! C12 'grid' file. The dimensions returned from +! the routine should be i=12, j=12. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : read_mdl_dims + + implicit none + + character(len=19) :: mdl_grid_file + + integer :: im, jm + + print*,"- Begin test of routine read_mdl_dims." + + mdl_grid_file="./C12_grid.tile1.nc" + + call read_mdl_dims(mdl_grid_file,im,jm) + + if (im /= 12) stop 3 + if (jm /= 12) stop 6 + + print*,"OK" + + print*,"SUCCESS" + + end program read_model_dims diff --git a/tests/orog/ftst_read_mdl_grid_file.F90 b/tests/orog/ftst_read_mdl_grid_file.F90 new file mode 100644 index 000000000..952212be5 --- /dev/null +++ b/tests/orog/ftst_read_mdl_grid_file.F90 @@ -0,0 +1,77 @@ + program read_model_grid_file + +! Test routine 'read_mdl_grid_file' by reading a +! C12 'grid' file from a global uniform grid. + +! Author George Gayno NCEP/EMC + + use io_utils, only : read_mdl_grid_file + + implicit none + + character(len=19) :: mdl_grid_file + + integer, parameter :: im = 12 + integer, parameter :: jm = 12 + integer :: i, j + + logical :: is_north_pole(im,jm) + logical :: is_south_pole(im,jm) + + real, parameter :: EPSILON=0.01 + real :: geolat(im,jm) + real :: geolat_c(im+1,jm+1) + real :: geolon(im,jm) + real :: geolon_c(im+1,jm+1) + real :: dx(im,jm), dy(im,jm) + + print*,"- Begin test of routine read_mdl_grid_file." + + mdl_grid_file="./C12_grid.tile1.nc" + +! Initialize all output variables to flag values. + + is_north_pole=.true. + is_south_pole=.true. + + geolat = -999.9 + geolat_c = -999.9 + geolon = -999.9 + geolon_c = -999.9 + dx = -999.9 + dy = -999.9 + + call read_mdl_grid_file(mdl_grid_file, im, jm, & + geolon, geolon_c, geolat, geolat_c, dx, dy, & + is_north_pole, is_south_pole) + +! Check values at selected points. + + if (abs(geolon_c(1,1)-305.0) > EPSILON) stop 2 + if (abs(geolon_c(13,13)-35.0) > EPSILON) stop 4 + if (abs(geolat_c(13,1)-(-35.2644)) > EPSILON) stop 6 + if (abs(geolat_c(1,13)-35.2644) > EPSILON) stop 8 + if (abs(geolon(5,5)-337.656) > EPSILON) stop 10 + if (abs(geolon(10,10)-17.9245) > EPSILON) stop 12 + if (abs(geolat(7,3)-(-27.84)) > EPSILON) stop 14 + if (abs(geolat(2,9)-16.8678) > EPSILON) stop 16 + if (abs(dx(1,1)-631596.076) > EPSILON) stop 18 + if (abs(dx(12,12)-631596.076) > EPSILON) stop 20 + +! There is no pole on this tile, so the pole variables +! should be .false. The routine sets dx equal to dy, +! so they should match. + + do j = 1, jm + do i = 1, im + if (abs(dx(i,j)-dy(i,j)) > EPSILON) stop 22 + if (is_north_pole(i,j)) stop 24 + if (is_south_pole(i,j)) stop 26 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program read_model_grid_file diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 new file mode 100644 index 000000000..04375fd6b --- /dev/null +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -0,0 +1,217 @@ + program rm_isolated_pts + +! Unit test for subroutine remove_isolated_pts. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : remove_isolated_pts + + implicit none + + integer :: im, jm, i, j, k + + real, parameter :: EPSILON=0.001 + + real, allocatable :: slm(:,:), oro(:,:), var(:,:), & + var4(:,:), oa(:,:,:), ol(:,:,:) + + real, allocatable :: slm_expected(:,:), oro_expected(:,:), var_expected(:,:), & + var4_expected(:,:), oa_expected(:,:,:), ol_expected(:,:,:) + + print*,"- Starting test of remove_isolated_pts." + +! A 5x4 grid. + + im = 5 + jm = 4 + + allocate (slm(im,jm)) + allocate (oro(im,jm)) + allocate (var(im,jm)) + allocate (var4(im,jm)) + allocate (oa(im,jm,4)) + allocate (ol(im,jm,4)) + +! Test Point 1. + +! This is an isolated island. The island should be +! removed (slm set to 0.0) and all other fields +! should be the average of the surrounding points (which +! for this test case is zero). All surrounding points +! should be unchanged. + + print*,'- Test point 1.' + +! Initialize grid to all ocean. + + slm = 0.0 ! land/sea mask + oro = 0.0 ! orography + var = 0.0 ! standard deviation of orography + var4 = 0.0 ! convexity + oa = 0.0 ! orographic asymetry + ol = 0.0 ! orographic length scale + +! Initialize an island in the middle of the grid. + + slm(2,2) = 1.0 + oro(2,2) = 50.0 + var(2,2) = 10.0 + var4(2,2) = 5.0 + + oa(2,2,1) = -1.0 + oa(2,2,2) = -0.5 + oa(2,2,3) = 0.5 + oa(2,2,4) = 1.0 + + ol(2,2,1) = 0.1 + ol(2,2,2) = 0.25 + ol(2,2,3) = 0.5 + ol(2,2,4) = 1.0 + + allocate (slm_expected(im,jm)) + allocate (oro_expected(im,jm)) + allocate (var_expected(im,jm)) + allocate (var4_expected(im,jm)) + allocate (oa_expected(im,jm,4)) + allocate (ol_expected(im,jm,4)) + +! All points should be ocean (all fields zero). + + slm_expected = 0.0 + oro_expected = 0.0 + var_expected = 0.0 + var4_expected = 0.0 + oa_expected = 0.0 + ol_expected = 0.0 + + call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) + + do j = 1, jm + do i = 1, im + if (abs(slm(i,j)-slm_expected(i,j)) > EPSILON) stop 2 + if (abs(oro(i,j)-oro_expected(i,j)) > EPSILON) stop 4 + if (abs(var(i,j)-var_expected(i,j)) > EPSILON) stop 6 + if (abs(var4(i,j)-var4_expected(i,j)) > EPSILON) stop 8 + do k = 1, 4 + if (abs(oa(i,j,k)-oa_expected(i,j,k)) > EPSILON) stop 10 + if (abs(ol(i,j,k)-ol_expected(i,j,k)) > EPSILON) stop 12 + enddo + enddo + enddo + + deallocate (slm, oro, var, var4, oa, ol) + deallocate (slm_expected, oro_expected, var_expected, var4_expected, oa_expected, ol_expected) + +! Test Point 2 + +! Now remove an isolated water point. The point should be changed +! to land (slm=1.0) and all other fields should be the average +! of the eight surrounding points. + + print*,'- Test point 2.' + +! A 3x3 grid. + + im = 3 + jm = 3 + + allocate (slm(im,jm)) + allocate (oro(im,jm)) + allocate (var(im,jm)) + allocate (var4(im,jm)) + allocate (oa(im,jm,4)) + allocate (ol(im,jm,4)) + + slm = 1.0 ! Initialize grid to all land. + slm(2,2) = 0.0 ! Set interior point to water. + + oro(1,1) = 5.0 + oro(2,1) = 10.0 + oro(3,1) = 17.0 + oro(1,2) = 22.0 + oro(2,2) = 0.0 ! water point. + oro(3,2) = 100.0 + oro(1,3) = 34.0 + oro(2,3) = 55.0 + oro(3,3) = 68.0 + + var = 0.5 * oro + var4 = 0.25 * oro + + ol(1,1,1) = 0.1 + ol(2,1,1) = 0.2 + ol(3,1,1) = 0.3 + ol(1,2,1) = 0.4 + ol(2,2,1) = 0.0 ! water point + ol(3,2,1) = 0.6 + ol(1,3,1) = 0.7 + ol(2,3,1) = 0.8 + ol(3,3,1) = 0.9 + + do j = 1, jm + do i = 1, im + ol(i,j,2) = ol(i,j,1) * .75 + ol(i,j,3) = ol(i,j,1) * .5 + ol(i,j,4) = ol(i,j,1) * .25 + enddo + enddo + + oa = -(ol) + +! All points should be land. The former isolated ocean +! point should have values that are the average of the +! surrounding land points. All other points should remain +! unchanged. + + allocate (slm_expected(im,jm)) + allocate (oro_expected(im,jm)) + allocate (var_expected(im,jm)) + allocate (var4_expected(im,jm)) + allocate (oa_expected(im,jm,4)) + allocate (ol_expected(im,jm,4)) + + slm_expected = slm + oro_expected = oro + var_expected = var + var4_expected = var4 + ol_expected = ol + oa_expected = oa + +! Former ocean point should be average of the surrounding +! ocean points. + + slm_expected(2,2) = 1.0 ! land point + oro_expected(2,2) = 38.875 + var_expected(2,2) = 0.5 * oro_expected(2,2) + var4_expected(2,2) = 0.25 * oro_expected(2,2) + ol_expected(2,2,1) = 0.5 + ol_expected(2,2,2) = 0.375 + ol_expected(2,2,3) = 0.25 + ol_expected(2,2,4) = 0.125 + oa_expected(2,2,1) = -0.5 + oa_expected(2,2,2) = -0.375 + oa_expected(2,2,3) = -0.25 + oa_expected(2,2,4) = -0.125 + + call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) + + do j = 1, jm + do i = 1, im + if (abs(slm(i,j)-slm_expected(i,j)) > EPSILON) stop 22 + if (abs(oro(i,j)-oro_expected(i,j)) > EPSILON) stop 24 + if (abs(var(i,j)-var_expected(i,j)) > EPSILON) stop 26 + if (abs(var4(i,j)-var4_expected(i,j)) > EPSILON) stop 28 + do k = 1, 4 + if (abs(oa(i,j,k)-oa_expected(i,j,k)) > EPSILON) stop 30 + if (abs(ol(i,j,k)-ol_expected(i,j,k)) > EPSILON) stop 32 + enddo + enddo + enddo + + deallocate (slm, oro, var, var4, oa, ol) + deallocate (slm_expected, oro_expected, var_expected, var4_expected, oa_expected, ol_expected) + + print*,"OK" + print*,"SUCCESS" + + end program rm_isolated_pts diff --git a/tests/orog/ftst_transpose.F90 b/tests/orog/ftst_transpose.F90 new file mode 100644 index 000000000..31f8a8d4d --- /dev/null +++ b/tests/orog/ftst_transpose.F90 @@ -0,0 +1,90 @@ + program transpose + +! Unit test for routines transpose_mask and transpose_orog. +! They are essentially the same routine - one takes +! one byte integer data, and one takes two byte integer +! data. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : transpose_mask, transpose_orog + + implicit none + + integer, parameter :: imn = 360 + integer, parameter :: jmn = 181 + + integer(1) :: mask(imn,jmn) + integer(2) :: mask2(imn,jmn) + integer :: i, j, jj + + print*,"Starting test of transpose routines." + +! Transpose is from S to N to the NCEP standard N to S, +! and, in the E/W direction, from the dateline to the +! NCEP standard Greenwich. + +! Test N/S transpose. Set up a one-degree global mask. +! Although mask is a yes/no flag, for this test set each +! row to the latitude to simplify checking the answer. + + jj=0 + do j = -90, 90 ! row 1 is South Pole. + jj = jj + 1 + mask(:,jj) = j + mask2(:,jj) = j + enddo + + print*,"Call transpose routines to test N/S flip." + + call transpose_mask(imn, jmn, mask) + call transpose_orog(imn, jmn, mask2) + + jj=0 + do j = 90, -90, -1 ! row 1 is North Pole. + jj = jj + 1 + do i = 1, imn + if (mask(i,jj) /= j) stop 2 + if (mask2(i,jj) /= j) stop 3 + enddo + enddo + +! Now test the transpose in the E/W direction. +! Here, the East half of the domain is a flag value +! of minus 1 and the West half is plus 1. + + do i = 1, 180 + mask(i,:) = -1 + mask2(i,:) = -1 + enddo + do i = 181, 360 + mask(i,:) = +1 + mask2(i,:) = +1 + enddo + + print*,"Call transpose routines to test E/W flip." + + call transpose_mask(imn, jmn, mask) + call transpose_orog(imn, jmn, mask2) + +! After the transpose, the East half should be plus 1 +! and the West half should be minus 1. + + do i = 1, 180 + do j = 1, jmn + if (mask(i,j) /= 1) stop 4 + if (mask2(i,j) /= 1) stop 5 + enddo + enddo + do i = 181, 360 + do j = 1, jmn + if (mask(i,j) /= -1) stop 6 + if (mask2(i,j) /= -1) stop 7 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program transpose diff --git a/util/gdas_init/set_fixed_files.sh b/util/gdas_init/set_fixed_files.sh index 0c3e1b4fe..451dd0bed 100755 --- a/util/gdas_init/set_fixed_files.sh +++ b/util/gdas_init/set_fixed_files.sh @@ -17,6 +17,12 @@ elif [ ${CTAR} == 'C768' ]; then OCNRES='025' elif [ ${CTAR} == 'C1152' ]; then OCNRES='025' +elif [ ${CTAR} == 'C12' ]; then + OCNRES='900' +elif [ ${CTAR} == 'C18' ]; then + OCNRES='900' +elif [ ${CTAR} == 'C24' ]; then + OCNRES='900' else OCNRES='025' fi