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 000000000..5bcfb1151 Binary files /dev/null and b/tests/orog/data/C12_grid.tile1.nc differ diff --git a/tests/orog/data/C48.mx500.tile1.nc b/tests/orog/data/C48.mx500.tile1.nc new file mode 100644 index 000000000..51b5fc0a4 Binary files /dev/null and b/tests/orog/data/C48.mx500.tile1.nc differ diff --git a/tests/orog/data/landcover.umd.lowres.nc b/tests/orog/data/landcover.umd.lowres.nc new file mode 100644 index 000000000..5652f438f Binary files /dev/null and b/tests/orog/data/landcover.umd.lowres.nc differ 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 000000000..1023adf6e Binary files /dev/null and b/tests/orog/data/topography.antarctica.ramp.lowres.nc differ diff --git a/tests/orog/data/topography.gmted2010.lowres.nc b/tests/orog/data/topography.gmted2010.lowres.nc new file mode 100644 index 000000000..80c71b1a5 Binary files /dev/null and b/tests/orog/data/topography.gmted2010.lowres.nc differ diff --git a/tests/orog/ftst_find_nearest_pole_pts.F90 b/tests/orog/ftst_find_nearest_pole_pts.F90 new file mode 100644 index 000000000..4c6844790 --- /dev/null +++ b/tests/orog/ftst_find_nearest_pole_pts.F90 @@ -0,0 +1,156 @@ + program find_nearest_pole_pts + +! Unit test for routine find_nearest_pole_points. +! +! The routine is passed the supergrid indices of +! the pole (from the 'grid' files) and returns the +! nearest pole points (using a logical array) on +! the model tile. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : find_nearest_pole_points + + implicit none + + integer, parameter :: im=48 + integer, parameter :: jm=48 + + integer :: i, j + integer :: i_north_pole, j_north_pole + integer :: i_south_pole, j_south_pole + + logical :: is_north_pole(im,jm) + logical :: is_south_pole(im,jm) + + print*,"Starting test of find_nearest_pole_points routine." + +! Test 1 - C48 uniform tile containing north pole. + + print*,"Test number 1." + + i_north_pole = 49 ! Supergrid index + j_north_pole = 49 + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 24 .and. j == 24) .or. & + (i == 24 .and. j == 25) .or. & + (i == 25 .and. j == 24) .or. & + (i == 25 .and. j == 25)) then + if (.not.is_north_pole(i,j)) stop 2 + else + if (is_north_pole(i,j)) stop 4 + endif + if (is_south_pole(i,j)) stop 8 + enddo + enddo + +! Test 2 - C48 uniform tile containing south pole. + + print*,"Test number 2." + + i_north_pole = 0 + j_north_pole = 0 + i_south_pole = 49 ! Supergrid index. + j_south_pole = 49 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 24 .and. j == 24) .or. & + (i == 24 .and. j == 25) .or. & + (i == 25 .and. j == 24) .or. & + (i == 25 .and. j == 25)) then + if (.not.is_south_pole(i,j)) stop 12 + else + if (is_south_pole(i,j)) stop 14 + endif + if (is_north_pole(i,j)) stop 18 + enddo + enddo + +! Test 3 - C48 uniform tile containing no pole. + + print*,"Test number 3." + + i_north_pole = 0 ! Zero indicates no pole. + j_north_pole = 0 + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if (is_south_pole(i,j)) stop 24 + if (is_north_pole(i,j)) stop 26 + enddo + enddo + +! Test 4 - C48 stretched grid tile containing south pole. + + print*,"Test number 4." + + i_north_pole = 0 + j_north_pole = 0 + i_south_pole = 10 ! Supergrid index. + j_south_pole = 49 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 5 .and. j == 24) .or. & + (i == 5 .and. j == 25)) then + if (.not.is_south_pole(i,j)) stop 32 + else + if (is_south_pole(i,j)) stop 34 + endif + if (is_north_pole(i,j)) stop 38 + enddo + enddo + +! Test 5 - C48 stretched grid tile containing north pole. + + print*,"Test number 5." + + i_north_pole = 62 ! Supergrid index. + j_north_pole = 49 + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 31 .and. j == 24) .or. & + (i == 31 .and. j == 25)) then + if (.not.is_north_pole(i,j)) stop 32 + else + if (is_north_pole(i,j)) stop 34 + endif + if (is_south_pole(i,j)) stop 38 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program find_nearest_pole_pts diff --git a/tests/orog/ftst_find_poles.F90 b/tests/orog/ftst_find_poles.F90 new file mode 100644 index 000000000..b0b9e953f --- /dev/null +++ b/tests/orog/ftst_find_poles.F90 @@ -0,0 +1,72 @@ + program ftst_find_poles + +! Unit test for routine find_poles. +! +! Checks if a model tile contains a 'pole' +! point based on the latitude of each point. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : find_poles + + implicit none + + integer, parameter :: nx=9 + integer, parameter :: ny=9 + + integer :: i_north_pole, j_north_pole + integer :: i_south_pole, j_south_pole + + real :: geolat(nx+1,ny+1) + + print*,"Starting test of find_poles." + +! Point 1 - North Pole at (4,5) + + print*,"North pole test." + + geolat = 89.0 + geolat(4,5) = 89.95 ! North pole + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 4) stop 2 + if (j_north_pole /= 5) stop 4 + if (i_south_pole /= 0) stop 6 + if (j_south_pole /= 0) stop 8 + +! Point 2 - South Pole at (2,8) + + print*,"South pole test." + + geolat = -89.0 + geolat(2,8) = -89.95 ! South pole + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 0) stop 12 + if (j_north_pole /= 0) stop 14 + if (i_south_pole /= 2) stop 16 + if (j_south_pole /= 8) stop 18 + +! Point 3 - No pole points. + + print*,"No pole test." + + geolat = -89.3 + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 0) stop 22 + if (j_north_pole /= 0) stop 24 + if (i_south_pole /= 0) stop 26 + if (j_south_pole /= 0) stop 28 + + print*,"OK" + + print*,"SUCCESS" + + end program ftst_find_poles diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 new file mode 100644 index 000000000..3ebf86d42 --- /dev/null +++ b/tests/orog/ftst_get_index.F90 @@ -0,0 +1,123 @@ + program test_get_index + +! Unit test for routine "get_index", which finds the +! i/j location of the model point on the high-resolution +! mask and orography grids. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : get_index + + implicit none + + integer, parameter :: imn=360*120 + integer, parameter :: jmn=180*120 + integer, parameter :: npts=4 + + integer :: i, ii, jst, jen, ilist(imn), numx + + real :: lono(npts), lato(npts) + real, parameter :: delxn=360.0/float(imn) + + print*,"Starting test of get_index." + +! The get_index routine contains a main 'if/elseif/else' statement. The test +! points are designed to check each branch. + +! Point 1 - At equator. Western edge of grid cell at +! Greenwich. This tests the 'else' branch. + + print*,"Test point 1." + + lato(1) = 0.0; lono(1) = 0.0 + lato(2) = 1.0; lono(2) = 0.5 + lato(3) = 0.0; lono(3) = 1.0 + lato(4) = -1.0; lono(4) = 0.5 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 10676) stop 2 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 4 ! Ending 'j' index on hi-res grid. + if (numx /= 121) stop 6 ! Number of 'i' indices on the hi-res grid. + do i = 1, numx + if (ilist(i) /= i) stop 8 ! List of 'i' indices on the hi-res grid. + enddo + +! Point 2 - At equator. Grid cell crosses Greenwich. +! This tests the 'elseif' branch. + + print*,"Test point 2." + + lato(1) = 0.0; lono(1) = 359.5 + lato(2) = 1.0; lono(2) = 0.0 + lato(3) = 0.0; lono(3) = 0.5 + lato(4) = -1.0; lono(4) = 0.0 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 10676) stop 12 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 14 ! Ending 'j' index on hi-res grid. + if (numx /= 121) stop 16 ! Number of 'i' indices on the hi-res grid. + ii = 1 + do i = 43141, 43200 ! East of greenwich + if (ilist(ii) /= i) stop 18 ! List of 'i' indices on the hi-res grid. + ii = ii + 1 + enddo + do i = 1, 61 ! West of greenwich + if (ilist(ii) /= i) stop 18 + ii = ii + 1 + enddo + +! Point 3 - At equator. Grid cell centered at +! the dateline. This tests the 'else' branch. + + print*,"Test point 3." + + lato(1) = -1.0; lono(1) = 179.0 + lato(2) = 1.0; lono(2) = 179.0 + lato(3) = 1.0; lono(3) = 181.0 + lato(4) = -1.0; lono(4) = 181.0 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 10676) stop 22 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 24 ! Ending 'j' index on hi-res grid. + if (numx /= 241) stop 26 ! Number of 'i' indices on the hi-res grid. + ii = 1 + do i = 21481, 21721 + if (ilist(ii) /= i) stop 28 ! List of 'i' indices on the hi-res grid. + ii = ii + 1 + enddo + +! Point 4 - At North Pole. Grid cell centered at +! Greenwich. This tests the 'if' branch. + + print*,"Test point 4." + + lato(1) = 89.0; lono(1) = 359.5 + lato(2) = 90.0; lono(2) = 359.5 + lato(3) = 90.0; lono(3) = 0.5 + lato(4) = 89.0; lono(4) = 0.5 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 21476) stop 32 ! Starting 'j' index on hi-res grid. + if (jen /= 21600) stop 34 ! Ending 'j' index on hi-res grid. + if (numx /= 43200) stop 36 ! Number of 'i' indices on the hi-res grid. + do i = 1, numx + if (ilist(i) /= i) stop 38 ! List of 'i' indices on the hi-res grid. + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program test_get_index diff --git a/tests/orog/ftst_get_xnsum.F90 b/tests/orog/ftst_get_xnsum.F90 new file mode 100644 index 000000000..8ad59a27e --- /dev/null +++ b/tests/orog/ftst_get_xnsum.F90 @@ -0,0 +1,125 @@ + program check_get_xnsum + +! Unit test for function get_xnsum, which computes the +! the number of high-resolution orography data points +! within a model grid box that are above the average +! terrain height. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : get_xnsum + + 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 + integer :: zslm(imn,jmn) ! High-res grid land mask + + 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 + real :: xnsum + + print*,"Begin test of function get_xnsum" + +! 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 approximately 16 high-res points located within +! the model box. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + +! First test point. The high-res grid is all ocean. Since all points in +! the model grid box will be sea level, there will be no points +! higher than the average of 0 meters. + + print*,"Test point 1." + + zslm = 0 ! high-res grid all water + zavg = -999 ! high-res grid all sea level + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 0) stop 2 + +! Second test point. All high-res points are land with an elevation +! of 50 meters. There will be no points higher than the average. + + print*,"Test point 2." + + zslm = 1 ! high-res grid all land + zavg = 50 ! constant elevation of 50 meters. + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 0) stop 4 + +! Third test point. All high-res points are land. One point is +! 100 meters, the rest are 50 meters. Therefore, there should +! be one point higher than the average. + + print*,"Test point 3." + + zslm = 1 ! all land + zavg = 50 ! initialize elevation to 50 meters. + zavg(359,91) = 100 ! one point within the model box is 100 meters. + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 1) stop 6 + +! Fourth test point. All high-res points are water. One point is +! 20 meters, which represents an inland lake. The other points +! are ocean. Therefore, there should be one point higher than the average. + + print*,"Test point 4." + + zslm = 0 ! all water + zavg = -999 ! initialize to sea level. + zavg(2,93) = 20 ! this represents an inland lake above sea level. + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 1) stop 8 + +! Fifth test point. Half of the high-res points within the grid +! box are land and half are water. Four high-res points are +! above the average. + + print*,"Test point 5." + + zslm = 0 ! Initialize to all water + zavg = -999 ! Initialize to sea level. + zslm(1:2,90:93) = 1 ! Half the points in the grid box are land. + zavg(1,90:93) = 25 ! Set the elevation at the + zavg(2,90) = 100 ! land points so that + zavg(2,91) = 110 ! four points are above the average. + zavg(2,92) = 107 + zavg(2,93) = 207 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 4) stop 10 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum diff --git a/tests/orog/ftst_get_xnsum2.F90 b/tests/orog/ftst_get_xnsum2.F90 new file mode 100644 index 000000000..3ad78d395 --- /dev/null +++ b/tests/orog/ftst_get_xnsum2.F90 @@ -0,0 +1,74 @@ + program check_get_xnsum2 + +! Unit test for routine get_xnsum2, 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 a function of the standard deviation of height. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : get_xnsum2 + + 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. + + real, parameter :: EPSILON=0.001 + real :: expected_hc + integer :: expected_xnsum1, expected_xnsum2 + + data expected_hc /677.2/ ! Expected critical height in meters. + data expected_xnsum1 /8/ ! 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_xnsum2." + +! 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 and the +! other half at 1000 meters. + + zavg = -999 ! Flag value for sea level. + zavg(359:360,90:93) = 1000 + + call get_xnsum2(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 + if (abs(hc-expected_hc) > 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