diff --git a/.gitignore b/.gitignore index f306e88a3c..079791b0f1 100644 --- a/.gitignore +++ b/.gitignore @@ -94,6 +94,7 @@ streamflow_obs_diag cam_dart_obs_preprocessor aether_to_dart dart_to_aether +create_geometry_file # Observation converter exectutables convert_aviso diff --git a/index.rst b/index.rst index 06b4686699..4b2ec13673 100644 --- a/index.rst +++ b/index.rst @@ -440,6 +440,7 @@ References models/9var/readme models/aether_lat-lon/readme + models/aether_cube_sphere/readme models/am2/readme models/bgrid_solo/readme models/cam-fv/readme diff --git a/models/README.rst b/models/README.rst index 932d60424e..9cd2ead3f9 100644 --- a/models/README.rst +++ b/models/README.rst @@ -6,6 +6,8 @@ Supported Models DART supported models: - :doc:`9var/readme` +- :doc:`aether_lat-lon/readme` +- :doc:`aether_cube_sphere/readme` - :doc:`am2/readme` - :doc:`bgrid_solo/readme` - :doc:`cam-fv/readme` diff --git a/models/aether_cube_sphere/aether_to_dart.f90 b/models/aether_cube_sphere/aether_to_dart.f90 new file mode 100644 index 0000000000..282ccaa2ed --- /dev/null +++ b/models/aether_cube_sphere/aether_to_dart.f90 @@ -0,0 +1,18 @@ +program aether_to_dart + +use utilities_mod, only : initialize_utilities, finalize_utilities +use transform_state_mod, only : initialize_transform_state_mod, model_to_dart, finalize_transform_state_mod + +implicit none + +call initialize_utilities(progname='aether_to_dart') + +call initialize_transform_state_mod() + +call model_to_dart() + +call finalize_transform_state_mod() + +call finalize_utilities('aether_to_dart') + +end program aether_to_dart diff --git a/models/aether_cube_sphere/create_geometry_file.f90 b/models/aether_cube_sphere/create_geometry_file.f90 new file mode 100644 index 0000000000..cfd4dbf392 --- /dev/null +++ b/models/aether_cube_sphere/create_geometry_file.f90 @@ -0,0 +1,625 @@ +program create_geometry_file + +use netcdf +use sort_mod, only : index_sort +use types_mod, only : i8, r8, PI, RAD2DEG, DEG2RAD +use location_mod, only : location_type, get_dist, get_location, set_location, VERTISHEIGHT +use utilities_mod, only : initialize_utilities, finalize_utilities, find_namelist_in_file, & + check_namelist_read +use netcdf_utilities_mod, only : nc_open_file_readonly, nc_get_variable_size, nc_get_variable, & + nc_create_file, nc_end_define_mode, nc_close_file +use transform_state_mod, only : file_type, integer_to_string, zero_fill +use quad_utils_mod, only : in_quad + +implicit none + +integer :: iunit, io + +! There will always be exactly 8 vertices on a cube sphere, regardless of how many blocks +integer, parameter :: nvertex_columns = 8 +integer, parameter :: nside_faces = 4 +integer, parameter :: nvertex_neighbors = 3 +integer, parameter :: nquad_neighbors = 4 +integer, dimension(nvertex_columns, 3) :: matrix_of_corner_indices +logical :: is_a_vertex +integer, dimension(3) :: corner_indices + +integer, dimension(3) :: variable_size + +integer :: nzs_per_centers_block, nys_per_centers_block, nxs_per_centers_block, & + nzs_per_corners_block, nys_per_corners_block, nxs_per_corners_block + +integer :: truncated_nys_per_centers_block, truncated_nxs_per_centers_block, & + truncated_nys_per_side_corners_block, truncated_nxs_per_side_corners_block, & + truncated_nys_per_top_bottom_corners_block, truncated_nxs_per_top_bottom_corners_block + +integer :: center_ncols, corner_ncols + +real(r8), allocatable, dimension(:, :, :) :: center_altitude_block + +real(r8), allocatable, dimension(:, :, :) :: center_lats_block, center_lons_block, & + corner_lats_block, corner_lons_block + +real(r8), allocatable, dimension(:, :, :) :: center_lats_truncated, center_lons_truncated, & + corner_lats_side_truncated, corner_lons_side_truncated, & + corner_lats_top_bottom_truncated, corner_lons_top_bottom_truncated + +real(r8), allocatable, dimension(:) :: center_latitude, center_longitude, & + quad_latitude, quad_longitude + +integer, allocatable, dimension(:, :) :: quad_neighbor_indices + +integer, dimension(nvertex_columns, nvertex_neighbors) :: vertex_neighbor_indices + +real(r8), dimension(nvertex_columns) :: vertex_latitude, vertex_longitude + +real(r8), dimension(3) :: location + +! Declare variables for the geometry_file that will be output by this program +real(r8) :: min_alt +real(r8) :: max_lon +real(r8) :: min_lon +real(r8) :: max_lat +real(r8) :: min_lat + +type :: vertex + type(location_type) :: loc + integer, dimension(3) :: column_index_of_neighbors +end type vertex + +type :: quad + type(location_type) :: loc + logical :: spans_prime_meridian + logical :: encloses_pole + integer, dimension(4) :: column_index_of_neighbors +end type quad + +! This is an intentionally unecessary nesting of a derived type just so that the lat and lon of each +! center_location is accessed using the same %loc attibute syntax +type :: center_location + type(location_type) :: loc +end type center_location + +type(vertex), dimension(nvertex_columns) :: vertexs +type(quad), allocatable, dimension(:) :: quads +type(center_location), allocatable, dimension(:) :: center_locations + +integer :: ivertex +integer :: iquad + +integer :: icol, ix, iy, nx, ny, inorth_pole_quad_column, isouth_pole_quad_column +real(r8), allocatable, dimension(:) :: distances_from_corners_to_centers + +character(len=256) :: restart_directory, grid_directory, filter_directory +namelist /directory_nml/ restart_directory, grid_directory, filter_directory + +character(len=256) :: grid_centers_file_prefix, grid_centers_file_suffix, grid_corners_file_prefix, grid_corners_file_suffix +namelist /grid_nml/ grid_centers_file_prefix, grid_centers_file_suffix, grid_corners_file_prefix, grid_corners_file_suffix + +integer :: nblocks, nhalos +character(len=256) :: restart_file_prefix, restart_file_middle, restart_file_suffix, filter_input_prefix, filter_input_suffix +namelist /transform_state_nml/ nblocks, nhalos,restart_file_prefix, restart_file_middle, restart_file_suffix, filter_input_prefix, filter_input_suffix + +! Assign initialized variables +max_lon = 360.0 +min_lon = 0.0 +max_lat = 90.0 +min_lat = -90.0 +corner_indices = [0, 0, 0] +location(:) = 0.0 +vertex_latitude(:) = 0.0 +vertex_longitude(:) = 0.0 +vertex_neighbor_indices(:, :) = 0 + +call initialize_utilities(progname='create_geometry_file') + +call initialize_program() + +call assign_triangles_and_quads() + +call sort_neighbors() + +call check_if_point_in_quad() + +call output_triangles_and_quads_to_netcdf() + +call finalize_utilities('create_geometry_file') + +contains + +subroutine initialize_program() + + call find_namelist_in_file('input.nml', 'directory_nml', iunit) + read(iunit, nml = directory_nml, iostat = io) + call check_namelist_read(iunit, io, 'directory_nml') + + call find_namelist_in_file('input.nml', 'grid_nml', iunit) + read(iunit, nml = grid_nml, iostat = io) + call check_namelist_read(iunit, io, 'grid_nml') + + call find_namelist_in_file('input.nml', 'transform_state_nml', iunit) + read(iunit, nml = transform_state_nml, iostat = io) + call check_namelist_read(iunit, io, 'transform_state_nml') + +end subroutine initialize_program + +subroutine assign_triangles_and_quads() + + type(file_type), allocatable, dimension(:) :: center_files + type(file_type), allocatable, dimension(:) :: corner_files + + integer :: iblock + + center_files = assign_grid_files_array(nblocks, grid_directory, grid_centers_file_prefix, grid_centers_file_suffix) + corner_files = assign_grid_files_array(nblocks, grid_directory, grid_corners_file_prefix, grid_corners_file_suffix) + + do iblock = 1, nblocks + + ! Center files + + center_files(iblock)%ncid = nc_open_file_readonly(center_files(iblock)%file_path) + + if (iblock == 1) then + ! Only call the nc_get_variable_size subroutine once for the center_files because the + ! Longitude and Latitude fields must be the same shape + call nc_get_variable_size(center_files(iblock)%ncid, 'Longitude', variable_size) + + nzs_per_centers_block = variable_size(1) + nys_per_centers_block = variable_size(2) + nxs_per_centers_block = variable_size(3) + + truncated_nys_per_centers_block = nys_per_centers_block-2*nhalos + truncated_nxs_per_centers_block = nxs_per_centers_block-2*nhalos + + allocate(center_altitude_block(nzs_per_centers_block, nys_per_centers_block, nxs_per_centers_block)) + center_altitude_block(:, :, :) = 0 + + allocate(center_lons_block(nzs_per_centers_block, nys_per_centers_block, nxs_per_centers_block)) + center_lons_block(:, :, :) = 0 + + allocate(center_lats_block(nzs_per_centers_block, nys_per_centers_block, nxs_per_centers_block)) + center_lats_block(:, :, :) = 0 + + allocate(center_lons_truncated(nblocks, truncated_nys_per_centers_block, truncated_nxs_per_centers_block)) + center_lons_truncated(:, :, :) = 0 + + allocate(center_lats_truncated(nblocks, truncated_nys_per_centers_block, truncated_nxs_per_centers_block)) + center_lats_truncated(:, :, :) = 0 + end if + + call nc_get_variable(center_files(iblock)%ncid, 'Altitude', center_altitude_block) + + call nc_get_variable(center_files(iblock)%ncid, 'Longitude', center_lons_block) + center_lons_truncated(iblock, :, :) = center_lons_block(1, nhalos+1:nys_per_centers_block-nhalos, nhalos+1:nxs_per_centers_block-nhalos) + + call nc_get_variable(center_files(iblock)%ncid, 'Latitude', center_lats_block) + center_lats_truncated(iblock, :, :) = center_lats_block(1, nhalos+1:nys_per_centers_block-nhalos, nhalos+1:nxs_per_centers_block-nhalos) + + ! Corner files + + corner_files(iblock)%ncid = nc_open_file_readonly(corner_files(iblock)%file_path) + + if (iblock == 1) then + ! Only call the nc_get_variable_size subroutine once for the corner_files because the + ! Longitude Corners and Latitude Corners fields must be the same shape + call nc_get_variable_size(corner_files(iblock)%ncid, 'Longitude Corners', variable_size) + + nzs_per_corners_block = variable_size(1) + nys_per_corners_block = variable_size(2) + nxs_per_corners_block = variable_size(3) + + truncated_nys_per_side_corners_block = nys_per_corners_block-(2*nhalos+2) + truncated_nxs_per_side_corners_block = nxs_per_corners_block-(2*nhalos+1) + + truncated_nys_per_top_bottom_corners_block = nys_per_corners_block-2*nhalos + truncated_nxs_per_top_bottom_corners_block = nxs_per_corners_block-2*nhalos + + allocate(corner_lons_block(nzs_per_corners_block, nys_per_corners_block, nxs_per_corners_block)) + corner_lons_block(:, :, :) = 0 + + allocate(corner_lats_block(nzs_per_corners_block, nys_per_corners_block, nxs_per_corners_block)) + corner_lats_block(:, :, :) = 0 + + allocate(corner_lons_side_truncated(nblocks, truncated_nys_per_side_corners_block, truncated_nxs_per_side_corners_block)) + corner_lons_side_truncated(:, :, :) = 0 + + allocate(corner_lats_side_truncated(nblocks, truncated_nys_per_side_corners_block, truncated_nxs_per_side_corners_block)) + corner_lats_side_truncated(:, :, :) = 0 + + allocate(corner_lons_top_bottom_truncated(nblocks, truncated_nys_per_top_bottom_corners_block, truncated_nxs_per_top_bottom_corners_block)) + corner_lons_top_bottom_truncated(:, :, :) = 0 + + allocate(corner_lats_top_bottom_truncated(nblocks, truncated_nys_per_top_bottom_corners_block, truncated_nxs_per_top_bottom_corners_block)) + corner_lats_top_bottom_truncated(:, :, :) = 0 + end if + + call nc_get_variable(corner_files(iblock)%ncid, 'Longitude Corners', corner_lons_block) + call nc_get_variable(corner_files(iblock)%ncid, 'Latitude Corners', corner_lats_block) + + if (iblock <= nside_faces) then + corner_lons_side_truncated(iblock, :, :) = corner_lons_block(1, nhalos+2:nys_per_corners_block-(nhalos+1), nhalos+2:nxs_per_corners_block-nhalos) + corner_lats_side_truncated(iblock, :, :) = corner_lats_block(1, nhalos+2:nys_per_corners_block-(nhalos+1), nhalos+2:nxs_per_corners_block-nhalos) + else + corner_lons_top_bottom_truncated(iblock, :, :) = corner_lons_block(1, nhalos+1:nys_per_corners_block-nhalos, nhalos+1:nxs_per_corners_block-nhalos) + corner_lats_top_bottom_truncated(iblock, :, :) = corner_lats_block(1, nhalos+1:nys_per_corners_block-nhalos, nhalos+1:nxs_per_corners_block-nhalos) + end if + + end do + + center_ncols = nblocks*truncated_nys_per_centers_block*truncated_nxs_per_centers_block + corner_ncols = nside_faces*truncated_nys_per_side_corners_block*truncated_nxs_per_side_corners_block + & + (nblocks-nside_faces)*truncated_nys_per_top_bottom_corners_block*truncated_nxs_per_top_bottom_corners_block + + allocate(center_locations(center_ncols)) + allocate(center_latitude(center_ncols)) + allocate(center_longitude(center_ncols)) + + min_alt = center_altitude_block(1, 1, 1) + + icol = 0 + do iblock = 1, nblocks + do iy = 1, truncated_nys_per_centers_block + do ix = 1, truncated_nxs_per_centers_block + icol = icol + 1 + ! lon and lat are private attributes of the location type stored internally as radians + ! and set_location requires the arguments to be in degrees. So, when invoking + ! set_location, the aether coordinates in radians must be converted to degrees for the + ! arguments before they are converted back to radians in the function + center_locations(icol)%loc = set_location(max(min_lon, min(max_lon, center_lons_truncated(iblock, iy, ix)*RAD2DEG)), max(min_lat, min(max_lat, center_lats_truncated(iblock, iy, ix)*RAD2DEG)), min_alt, VERTISHEIGHT) + + location = get_location(center_locations(icol)%loc) + center_longitude(icol) = location(1) + center_latitude(icol) = location(2) + + end do + end do + end do + + matrix_of_corner_indices = create_matrix_of_corner_indices(nblocks, nvertex_columns, & + truncated_nys_per_top_bottom_corners_block, & + truncated_nxs_per_top_bottom_corners_block) + + allocate(distances_from_corners_to_centers(center_ncols)) + allocate(quads(corner_ncols-nvertex_columns)) + allocate(quad_latitude(corner_ncols-nvertex_columns)) + allocate(quad_longitude(corner_ncols-nvertex_columns)) + allocate(quad_neighbor_indices(corner_ncols-nvertex_columns, nquad_neighbors)) + + quad_latitude(:) = 0.0 + quad_longitude(:) = 0.0 + quad_neighbor_indices(:, :) = 0 + + ivertex = 0 + iquad = 0 + do iblock = 1, nblocks + if (iblock <= nside_faces) then + ny = truncated_nys_per_side_corners_block + nx = truncated_nxs_per_side_corners_block + else + ny = truncated_nys_per_top_bottom_corners_block + nx = truncated_nxs_per_top_bottom_corners_block + end if + do iy = 1, ny + do ix = 1, nx + corner_indices(1) = iblock + corner_indices(2) = iy + corner_indices(3) = ix + is_a_vertex = is_corner_a_vertex(nvertex_columns, matrix_of_corner_indices, corner_indices) + if (is_a_vertex .eqv. .true.) then + ivertex = ivertex + 1 + if (iblock <= nside_faces) then + vertexs(ivertex)%loc = set_location(max(min_lon, min(max_lon, corner_lons_side_truncated(iblock, iy, ix)*RAD2DEG)), max(min_lat, min(max_lat, corner_lats_side_truncated(iblock, iy, ix)*RAD2DEG)), min_alt, VERTISHEIGHT) + else + vertexs(ivertex)%loc = set_location(max(min_lon, min(max_lon, corner_lons_top_bottom_truncated(iblock, iy, ix)*RAD2DEG)), max(min_lat, min(max_lat, corner_lats_top_bottom_truncated(iblock, iy, ix)*RAD2DEG)), min_alt, VERTISHEIGHT) + end if + + location = get_location(vertexs(ivertex)%loc) + vertex_longitude(ivertex) = location(1) + vertex_latitude(ivertex) = location(2) + + else + iquad = iquad + 1 + if (iblock <= nside_faces) then + quads(iquad)%loc = set_location(max(min_lon, min(max_lon, corner_lons_side_truncated(iblock, iy, ix)*RAD2DEG)), max(min_lat, min(max_lat, corner_lats_side_truncated(iblock, iy, ix)*RAD2DEG)), min_alt, VERTISHEIGHT) + else + quads(iquad)%loc = set_location(max(min_lon, min(max_lon, corner_lons_top_bottom_truncated(iblock, iy, ix)*RAD2DEG)), max(min_lat, min(max_lat, corner_lats_top_bottom_truncated(iblock, iy, ix)*RAD2DEG)), min_alt, VERTISHEIGHT) + end if + + location = get_location(quads(iquad)%loc) + quad_longitude(iquad) = location(1) + quad_latitude(iquad) = location(2) + + end if + end do + end do + end do + + do ivertex = 1, nvertex_columns + distances_from_corners_to_centers(:) = 0.0 + + do icol = 1, center_ncols + distances_from_corners_to_centers(icol) = get_dist(vertexs(ivertex)%loc, center_locations(icol)%loc) + end do + + call find_indices_of_n_smallest_elements_of_vector(distances_from_corners_to_centers, vertex_neighbor_indices(ivertex, :)) + + end do + + do iquad = 1, corner_ncols-nvertex_columns + distances_from_corners_to_centers(:) = 0.0 + + do icol = 1, center_ncols + distances_from_corners_to_centers(icol) = get_dist(quads(iquad)%loc, center_locations(icol)%loc) + end do + + call find_indices_of_n_smallest_elements_of_vector(distances_from_corners_to_centers, quad_neighbor_indices(iquad, :)) + + end do + +end subroutine assign_triangles_and_quads + +subroutine output_triangles_and_quads_to_netcdf() + + type(file_type) :: geometry_file + integer :: dim_id_center_altitudes + integer :: dim_id_vertex_columns, dim_id_quad_columns, dim_id_center_columns + integer :: dim_id_vertex_neighbors, dim_id_quad_neighbors + integer :: var_id_center_altitude + integer :: var_id_center_longitude, var_id_center_latitude + integer :: var_id_vertex_longitude, var_id_vertex_latitude, var_id_vertex + integer :: var_id_quad_longitude, var_id_quad_latitude, var_id_quad + integer :: xtype + + xtype = 5 + + geometry_file%file_path = trim(filter_directory) // 'geometry_file.nc' + + geometry_file%ncid = nc_create_file(geometry_file%file_path) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'center_altitudes', nzs_per_centers_block, dim_id_center_altitudes) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'vertex_columns', nvertex_columns, dim_id_vertex_columns) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'vertex_neighbors', nvertex_neighbors, dim_id_vertex_neighbors) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'quad_columns', (corner_ncols-nvertex_columns), dim_id_quad_columns) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'quad_neighbors', nquad_neighbors, dim_id_quad_neighbors) + + geometry_file%ncstatus = nf90_def_dim(geometry_file%ncid, 'center_columns', center_ncols, dim_id_center_columns) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'center_altitude', xtype, dim_id_center_altitudes, var_id_center_altitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'vertex_longitude', xtype, dim_id_vertex_columns, var_id_vertex_longitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'vertex_latitude', xtype, dim_id_vertex_columns, var_id_vertex_latitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'quad_longitude', xtype, dim_id_quad_columns, var_id_quad_longitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'quad_latitude', xtype, dim_id_quad_columns, var_id_quad_latitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'center_longitude', xtype, dim_id_center_columns, var_id_center_longitude) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'center_latitude', xtype, dim_id_center_columns, var_id_center_latitude) + + xtype = 4 + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'vertex_neighbor_indices', xtype, [dim_id_vertex_columns, dim_id_vertex_neighbors], var_id_vertex) + + geometry_file%ncstatus = nf90_def_var(geometry_file%ncid, 'quad_neighbor_indices', xtype, [dim_id_quad_columns, dim_id_quad_neighbors], var_id_quad) + + ! Add the indices of the north and south pole quad + geometry_file%ncstatus = nf90_put_att(geometry_file%ncid, NF90_GLOBAL, 'index_of_north_pole_quad_column', inorth_pole_quad_column) + geometry_file%ncstatus = nf90_put_att(geometry_file%ncid, NF90_GLOBAL, 'index_of_south_pole_quad_column', isouth_pole_quad_column) + + call nc_end_define_mode(geometry_file%ncid) + + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_center_altitude, center_altitude_block(:, 1, 1)) + + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_center_longitude, center_longitude) + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_center_latitude, center_latitude) + + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_vertex_longitude, vertex_longitude) + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_vertex_latitude, vertex_latitude) + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_vertex, vertex_neighbor_indices) + + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_quad_longitude, quad_longitude) + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_quad_latitude, quad_latitude) + geometry_file%ncstatus = nf90_put_var(geometry_file%ncid, var_id_quad, quad_neighbor_indices) + + call nc_close_file(geometry_file%ncid) + +end subroutine output_triangles_and_quads_to_netcdf + +function assign_grid_files_array(nblocks, grid_directory, grid_file_prefix, grid_file_suffix) result(grid_files) + + integer, intent(in) :: nblocks + character(len=*), intent(in) :: grid_directory + character(len=*), intent(in) :: grid_file_prefix + character(len=*), intent(in) :: grid_file_suffix + type(file_type), allocatable, dimension(:) :: grid_files + character(len=4) :: block_name + integer :: iblock + + print *, 'nblocks' + print *, nblocks + + allocate(grid_files(nblocks)) + + do iblock = 1, nblocks + block_name = zero_fill(integer_to_string(iblock-1), 4) + grid_files(iblock)%file_path = trim(grid_directory) // trim(grid_file_prefix) // & + block_name // trim(grid_file_suffix) + print *, grid_files(iblock)%file_path + end do + +end function assign_grid_files_array + +subroutine find_indices_of_n_smallest_elements_of_vector(vector, smallest_indices) + real(r8), dimension(:), intent(in) :: vector + integer, dimension(:), intent(inout) :: smallest_indices + integer :: iindex, nindices + logical, allocatable, dimension(:) :: vector_mask + + nindices = size(smallest_indices) + allocate(vector_mask(size(vector))) + vector_mask(:) = .true. + + do iindex=1, nindices + ! This peculiar (iindex:iindex) array indexing is required because minloc returns an array of dim 1 + smallest_indices(iindex:iindex) = minloc(vector, mask=vector_mask) + vector_mask(smallest_indices(iindex)) = .false. + end do + +end subroutine find_indices_of_n_smallest_elements_of_vector + +function create_matrix_of_corner_indices(nblocks, nvertex_columns, truncated_nys_per_corners_block, & + truncated_nxs_per_corners_block) result(matrix_of_corner_indices) + integer, intent(in) :: nblocks + integer, intent(in) :: nvertex_columns + integer, intent(in) :: truncated_nys_per_corners_block + integer, intent(in) :: truncated_nxs_per_corners_block + integer, dimension(nvertex_columns, 3) :: matrix_of_corner_indices + + matrix_of_corner_indices(:, :) = 0 + + matrix_of_corner_indices(1, :) = [nblocks-1, 1, 1] + matrix_of_corner_indices(2, :) = [nblocks-1, 1, truncated_nxs_per_corners_block] + matrix_of_corner_indices(3, :) = [nblocks-1, truncated_nys_per_corners_block, 1] + matrix_of_corner_indices(4, :) = [nblocks-1, truncated_nys_per_corners_block, truncated_nxs_per_corners_block] + matrix_of_corner_indices(5, :) = [nblocks, 1, 1] + matrix_of_corner_indices(6, :) = [nblocks, 1, truncated_nxs_per_corners_block] + matrix_of_corner_indices(7, :) = [nblocks, truncated_nys_per_corners_block, 1] + matrix_of_corner_indices(8, :) = [nblocks, truncated_nys_per_corners_block, truncated_nxs_per_corners_block] + +end function create_matrix_of_corner_indices + +function is_corner_a_vertex(nvertex_columns, matrix_of_corner_indices, corner_indices) result(is_a_vertex) + integer, intent(in) :: nvertex_columns + integer, dimension(nvertex_columns, 3), intent(in) :: matrix_of_corner_indices + integer, dimension(3), intent(in) :: corner_indices + logical :: is_a_vertex + integer :: ivertex + + is_a_vertex = .false. + do ivertex = 1, nvertex_columns + if (all(matrix_of_corner_indices(ivertex, :) == corner_indices)) then + is_a_vertex = .true. + exit + end if + end do + +end function is_corner_a_vertex + +subroutine check_if_point_in_quad() + + real(r8), dimension(4) :: x_neighbors, y_neighbors + logical :: inside + integer :: num_outside, num_inside, num_poles + integer :: icolumn, ineighbor + real(r8) :: tolerance + + tolerance = 0.00000314159_r8 + + num_outside = 0 + num_inside = 0 + num_poles = 0 + + do icolumn=1, (corner_ncols-nvertex_columns) + x_neighbors(:) = 0.0 + y_neighbors(:) = 0.0 + do ineighbor=1, nquad_neighbors + x_neighbors(ineighbor) = center_longitude(quad_neighbor_indices(icolumn, ineighbor)) + y_neighbors(ineighbor) = center_latitude(quad_neighbor_indices(icolumn, ineighbor)) + end do + + inside = in_quad(quad_longitude(icolumn), quad_latitude(icolumn), x_neighbors, y_neighbors) + + if (inside) then + num_inside = num_inside + 1 + else if (abs(-90.0-quad_latitude(icolumn)) < tolerance) then + isouth_pole_quad_column = icolumn + num_poles = num_poles + 1 + else if (abs(90.0-quad_latitude(icolumn)) < tolerance) then + inorth_pole_quad_column = icolumn + num_poles = num_poles + 1 + else + num_outside = num_outside + 1 + end if + end do + +end subroutine check_if_point_in_quad + +subroutine sort_neighbors() + + real(r8), dimension(nquad_neighbors) :: atan2_results + real(r8), dimension(nquad_neighbors) :: offset_lons, offset_lats + integer, dimension(nquad_neighbors) :: sorted_array + real(r8) :: central_lon, central_lat + integer(i8), dimension(4) :: indices + integer(i8) :: nneighbors + integer :: ineighbor + logical :: in_first_quadrant, in_fourth_quadrant + + ! Need to declare an additional integer equal to 4, since the interface defined for index sort, + ! when passed r8 values, requires an i8 length + nneighbors = nquad_neighbors + atan2_results(:) = 0 + + do iquad = 1, (corner_ncols-nvertex_columns) + + central_lon = quad_longitude(iquad) + central_lat = quad_latitude(iquad) + + ! Assign the offset_lons and offset_lats elements + do ineighbor = 1, nquad_neighbors + offset_lats(ineighbor) = center_latitude(quad_neighbor_indices(iquad, ineighbor)) + offset_lons(ineighbor) = center_longitude(quad_neighbor_indices(iquad, ineighbor)) + end do + + ! Check if the quad spans the prime meridian. It's definitely possible to perform this check + ! with only a single logical variable but using two makes the code more readable. + in_first_quadrant = .false. + in_fourth_quadrant = .false. + + do ineighbor = 1, nquad_neighbors + if (offset_lons(ineighbor) >= 0.0 .and. offset_lons(ineighbor) <= 90.0) then + in_first_quadrant = .true. + else if (offset_lons(ineighbor) >= 270.0 .and. offset_lons(ineighbor) <= 360.0) then + in_fourth_quadrant = .true. + end if + end do + + if (in_first_quadrant .and. in_fourth_quadrant) then + central_lon = central_lon + 180.0 + if (central_lon >= 360.0) then + central_lon = central_lon - 360.0 + end if + ! The quad spans the prime meridian + do ineighbor = 1, nquad_neighbors + offset_lons(ineighbor) = offset_lons(ineighbor) + 180.0 + if (offset_lons(ineighbor) >= 360.0) then + offset_lons(ineighbor) = offset_lons(ineighbor) - 360.0 + end if + end do + end if + + do ineighbor = 1, nquad_neighbors + offset_lons(ineighbor) = (central_lon-offset_lons(ineighbor))*DEG2RAD + offset_lats(ineighbor) = (central_lat-offset_lats(ineighbor))*DEG2RAD + atan2_results(ineighbor) = atan2(offset_lats(ineighbor), offset_lons(ineighbor)) + end do + + call index_sort(atan2_results, indices, nneighbors) + + do ineighbor = 1, nneighbors + sorted_array(ineighbor) = quad_neighbor_indices(iquad, indices(ineighbor)) + end do + + quad_neighbor_indices(iquad, :) = sorted_array(:) + + end do + +end subroutine sort_neighbors + +end program create_geometry_file diff --git a/models/aether_cube_sphere/dart_to_aether.f90 b/models/aether_cube_sphere/dart_to_aether.f90 new file mode 100644 index 0000000000..0fb06c94b4 --- /dev/null +++ b/models/aether_cube_sphere/dart_to_aether.f90 @@ -0,0 +1,18 @@ +program dart_to_aether + +use utilities_mod, only : initialize_utilities, finalize_utilities +use transform_state_mod, only : initialize_transform_state_mod, dart_to_model, finalize_transform_state_mod + +implicit none + +call initialize_utilities(progname='dart_to_aether') + +call initialize_transform_state_mod() + +call dart_to_model() + +call finalize_transform_state_mod() + +call finalize_utilities('dart_to_aether') + +end program dart_to_aether diff --git a/models/aether_cube_sphere/model_mod.f90 b/models/aether_cube_sphere/model_mod.f90 new file mode 100644 index 0000000000..f50d0c390d --- /dev/null +++ b/models/aether_cube_sphere/model_mod.f90 @@ -0,0 +1,835 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +module model_mod + +! This is a template showing the interfaces required for a model to be compliant +! with the DART data assimilation infrastructure. Do not change the arguments +! for the public routines. + +use netcdf + +use types_mod, only : r8, i8, MISSING_R8, vtablenamelength, DEG2RAD, radius => earth_radius + +use time_manager_mod, only : time_type, set_time + +use location_mod, only : location_type, get_close_type, & + loc_get_close_obs => get_close_obs, & + loc_get_close_state => get_close_state, & + set_location, set_location_missing, VERTISHEIGHT, query_location, & + get_location + +use utilities_mod, only : register_module, error_handler, & + E_ERR, E_MSG, & + nmlfileunit, do_output, do_nml_file, do_nml_term, & + find_namelist_in_file, check_namelist_read, to_upper, & + find_enclosing_indices + +use obs_kind_mod, only : get_index_for_quantity + +use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & + nc_add_global_creation_time, & + nc_begin_define_mode, nc_end_define_mode, nc_open_file_readonly, & + nc_close_file +use distributed_state_mod, only : get_state +use state_structure_mod, only : add_domain, get_dart_vector_index, get_domain_size, & + get_model_variable_indices, get_varid_from_kind + +use transform_state_mod, only : file_type + +use ensemble_manager_mod, only : ensemble_type + +! These routines are passed through from default_model_mod. +! To write model specific versions of these routines +! remove the routine from this use statement and add your code to +! this the file. +use default_model_mod, only : pert_model_copies, read_model_time, write_model_time, & + init_time => fail_init_time, & + init_conditions => fail_init_conditions, & + convert_vertical_obs, convert_vertical_state, adv_1step + +use quad_utils_mod, only : in_quad, quad_bilinear_interp + +implicit none +private + +! routines required by DART code - will be called from filter and other +! DART executables. +public :: get_model_size, & + get_state_meta_data, & + model_interpolate, & + end_model, & + static_init_model, & + nc_write_model_atts, & + get_close_obs, & + get_close_state, & + pert_model_copies, & + convert_vertical_obs, & + convert_vertical_state, & + read_model_time, & + adv_1step, & + init_time, & + init_conditions, & + shortest_time_between_assimilations, & + write_model_time + +integer :: iunit, io + +character(len=256), parameter :: source = "model_mod.f90" +logical :: module_initialized = .false. +integer :: dom_id ! used to access the state structure +type(time_type) :: assimilation_time_step + +real(r8), parameter :: roundoff = 1.0e-12_r8 + +! Geometry variables that are used throughout the module +integer :: ncenter_columns, ncenter_altitudes ! The number of center_columns and altitudes are read from the geometry file +integer, parameter :: nvertex_columns = 8 +integer, parameter :: nvertex_neighbors = 3 +integer :: nquad_columns ! The number of quad_columns is read from the geometry file +integer, parameter :: nquad_neighbors = 4 + +! Just like in cam-se, the aether cube_sphere filter input files are created to have a horizonal +! column dimension rather than being functions of latitude and longitude. +integer :: no_third_dimension = -99 + +integer :: inorth_pole_quad_column, isouth_pole_quad_column +real(r8), allocatable, dimension(:) :: center_latitude, center_longitude, center_altitude +integer, allocatable, dimension(:, :) :: quad_neighbor_indices, vertex_neighbor_indices + +! Error codes +integer, parameter :: INVALID_VERT_COORD_ERROR_CODE = 15 +integer, parameter :: INVALID_ALTITUDE_VAL_ERROR_CODE = 17 +integer, parameter :: UNKNOWN_OBS_QTY_ERROR_CODE = 20 + +! Example Namelist +! Use the namelist for options to be set at runtime. +character(len=256) :: template_file = 'model_restart.nc' +integer :: time_step_days = 0 +integer :: time_step_seconds = 3600 + +integer, parameter :: MAX_STATE_VARIABLES = 100 +integer, parameter :: NUM_STATE_TABLE_COLUMNS = 5 +character(len=vtablenamelength) :: variables(NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '' + +type :: var_type + integer :: count + character(len=64), allocatable :: names(:) + integer, allocatable :: qtys(:) + real(r8), allocatable :: clamp_values(:, :) + logical, allocatable :: updates(:) +end type var_type + +namelist /model_nml/ template_file, time_step_days, time_step_seconds, variables + +contains + +!------------------------------------------------------------------ +! +! Called to do one time initialization of the model. As examples, +! might define information about the model size or model timestep. +! In models that require pre-computed static data, for instance +! spherical harmonic weights, these would also be computed here. + +subroutine static_init_model() + +type(var_type) :: var + +module_initialized = .true. + +! Print module information to log file and stdout. +call register_module(source) + +call find_namelist_in_file('input.nml', 'model_nml', iunit) +read(iunit, nml = model_nml, iostat = io) +call check_namelist_read(iunit, io, 'model_nml') + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=model_nml) +if (do_nml_term()) write( * , nml=model_nml) + +! This time is both the minimum time you can ask the model to advance +! (for models that can be advanced by filter) and it sets the assimilation +! window. All observations within +/- 1/2 this interval from the current +! model time will be assimilated. If this is not settable at runtime +! feel free to hardcode it and remove from the namelist. +assimilation_time_step = set_time(time_step_seconds, & + time_step_days) + +var = assign_var(variables, MAX_STATE_VARIABLES) + +call read_geometry_file() + +! Define which variables are in the model state +dom_id = add_domain(template_file, var%count, var%names, var%qtys, & + var%clamp_values, var%updates) + +end subroutine static_init_model + +!------------------------------------------------------------------ +! Returns the number of items in the state vector as an integer. + +function get_model_size() + +integer(i8) :: get_model_size + +if ( .not. module_initialized ) call static_init_model + +get_model_size = get_domain_size(dom_id) + +end function get_model_size + + +!------------------------------------------------------------------ +! Given a state handle, a location, and a state quantity, +! interpolates the state variable fields to that location and returns +! the values in expected_obs. The istatus variables should be returned as +! 0 unless there is some problem in computing the interpolation in +! which case a positive istatus should be returned. +! +! For applications in which only perfect model experiments +! with identity observations (i.e. only the value of a particular +! state variable is observed), this can be a NULL INTERFACE. + +subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs, istatus) + + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: qty +real(r8), intent(out) :: expected_obs(ens_size) !< array of interpolated values +integer, intent(out) :: istatus(ens_size) + +character(len=512) :: error_string_1 + +! Location values stored in a vector +real(r8), dimension(3) :: lon_lat_alt + +! Logical needed for quad_bilinear_interp +logical :: cyclic + +! Vertical interpolation variables +integer(i8) :: state_index +integer :: below_index, above_index, enclosing_status, which_vertical +real(r8) :: fraction + +real(r8), dimension(nvertex_neighbors, ens_size) :: vertex_temp_values +real(r8), dimension(nquad_neighbors, ens_size) :: quad_temp_values +real(r8), dimension(ens_size) :: above_values, below_values + +! Actual variables needed for this routine + +integer :: icolumn, ineighbor, iens + +logical :: inside + +real(r8), dimension(nvertex_neighbors, 3) :: vertex_xyz +real(r8), dimension(3) :: r +real(r8), dimension(3) :: weights +real(r8), dimension(nquad_neighbors) :: x_neighbors_quad, y_neighbors_quad + +! Begin local variables for model_interpolate +integer :: varid + +! Begin variables for horiziontal interpolation + +! End variables for horizontal interpolation +cyclic = .true. + +lon_lat_alt = get_location(location) +which_vertical = nint(query_location(location)) + +! See if the state contains the obs quantity +varid = get_varid_from_kind(dom_id, qty) + +if (varid > 0) then + istatus = 0 +else + istatus = UNKNOWN_OBS_QTY_ERROR_CODE +endif + +if ( .not. module_initialized ) call static_init_model + +expected_obs(:) = MISSING_R8 +istatus(:) = 1 + +! Find the vertical levels +if ( which_vertical == VERTISHEIGHT ) then + call find_enclosing_indices(ncenter_altitudes, center_altitude, lon_lat_alt(3), below_index, & + above_index, fraction, enclosing_status) + if (enclosing_status /= 0) then + istatus(:) = INVALID_ALTITUDE_VAL_ERROR_CODE + end if +else + istatus(:) = INVALID_VERT_COORD_ERROR_CODE + write(error_string_1, *) 'unsupported vertical type: ', which_vertical + call error_handler(E_ERR, 'model_interpolate', error_string_1, source) +end if + +! If the vertical location is acceptable, then do the horizontal interpolation +if (istatus(1) == 1) then + + ! Find the enclosing triangle or quad + inside = .false. + do icolumn = 1, nvertex_columns + do ineighbor = 1, nvertex_neighbors + call latlon_to_xyz(center_latitude(vertex_neighbor_indices(icolumn, ineighbor)), & + center_longitude(vertex_neighbor_indices(icolumn, ineighbor)), & + vertex_xyz(ineighbor, 1), & + vertex_xyz(ineighbor, 2), & + vertex_xyz(ineighbor, 3)) + end do + + call latlon_to_xyz(lon_lat_alt(2), lon_lat_alt(1), r(1), r(2), r(3)) + + call inside_triangle(vertex_xyz(1, :), vertex_xyz(2, :), vertex_xyz(3, :), r, lon_lat_alt(2), lon_lat_alt(1), inside, weights) + + if (inside) then + + ! Do above level + do ineighbor = 1, nvertex_neighbors + state_index = get_dart_vector_index(vertex_neighbor_indices(icolumn, ineighbor), above_index, no_third_dimension, dom_id, varid) + vertex_temp_values(ineighbor, :) = get_state(state_index, state_handle) + end do + + above_values = barycentric_average(ens_size, weights, vertex_temp_values) + + ! Do below level + do ineighbor = 1, nvertex_neighbors + state_index = get_dart_vector_index(vertex_neighbor_indices(icolumn, ineighbor), below_index, no_third_dimension, dom_id, varid) + vertex_temp_values(ineighbor, :) = get_state(state_index, state_handle) + end do + + below_values = barycentric_average(ens_size, weights, vertex_temp_values) + + call vert_interp(ens_size, below_values, above_values, fraction, expected_obs) + + exit + end if + end do + + ! If the location is not inside any of the vertex triangles, do the quad loop + + if (.not. inside) then + + do icolumn = 1, nquad_columns + if (icolumn == inorth_pole_quad_column) then + inside = .true. + else if (icolumn == isouth_pole_quad_column) then + inside = .true. + else + do ineighbor = 1, nquad_neighbors + x_neighbors_quad(ineighbor) = center_longitude(quad_neighbor_indices(icolumn, ineighbor)) + y_neighbors_quad(ineighbor) = center_latitude(quad_neighbor_indices(icolumn, ineighbor)) + end do + inside = in_quad(lon_lat_alt(1), lon_lat_alt(2), x_neighbors_quad, y_neighbors_quad) + end if + + if (inside) then + + ! Get quad temp_values for the above level + do ineighbor = 1, nquad_neighbors + + state_index = get_dart_vector_index(quad_neighbor_indices(icolumn, ineighbor), above_index, no_third_dimension, dom_id, varid) + + quad_temp_values(ineighbor, :) = get_state(state_index, state_handle) + end do + + do iens = 1, ens_size + call quad_bilinear_interp(lon_lat_alt(1), lon_lat_alt(2), x_neighbors_quad, & + y_neighbors_quad, cyclic, quad_temp_values(:,iens), & + above_values(iens)) + enddo + + ! Get quad temp_values for the below level + do ineighbor =1, nquad_neighbors + state_index = get_dart_vector_index(quad_neighbor_indices(icolumn, ineighbor), below_index, no_third_dimension, dom_id, varid) + quad_temp_values(ineighbor, :) = get_state(state_index, state_handle) + end do + + do iens = 1, ens_size + call quad_bilinear_interp(lon_lat_alt(1), lon_lat_alt(2), x_neighbors_quad, & + y_neighbors_quad, cyclic, quad_temp_values(:,iens), & + below_values(iens)) + enddo + + call vert_interp(ens_size, below_values, above_values, fraction, expected_obs) + + exit + end if + + end do + + end if + + ! All good. + istatus(:) = 0 + +end if + +end subroutine model_interpolate + + +!------------------------------------------------------------------ +! Returns the smallest increment in time that the model is capable +! of advancing the state in a given implementation, or the shortest +! time you want the model to advance between assimilations. + +function shortest_time_between_assimilations() + +type(time_type) :: shortest_time_between_assimilations + +if ( .not. module_initialized ) call static_init_model + +shortest_time_between_assimilations = assimilation_time_step + +end function shortest_time_between_assimilations + + +!------------------------------------------------------------------ +! Given an integer index into the state vector, returns the +! associated location and optionally the physical quantity. + +subroutine get_state_meta_data(index_in, location, qty) + +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, intent(out), optional :: qty + +! Local variables + +integer :: lev_index, col_index +integer :: my_var_id, my_qty + +if ( .not. module_initialized ) call static_init_model + +call get_model_variable_indices(index_in, col_index, lev_index, no_third_dimension, & + var_id=my_var_id, kind_index=my_qty) + +! should be set to the actual location using set_location() +location = set_location(center_longitude(col_index), center_latitude(col_index), center_altitude(lev_index), VERTISHEIGHT) + +! should be set to the physical quantity, e.g. QTY_TEMPERATURE +if (present(qty)) qty = my_qty + +end subroutine get_state_meta_data + + +!------------------------------------------------------------------ +! Any model specific distance calcualtion can be done here +subroutine get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & + num_close, close_ind, dist, ens_handle) + +type(get_close_type), intent(in) :: gc ! handle to a get_close structure +integer, intent(in) :: base_type ! observation TYPE +type(location_type), intent(inout) :: base_loc ! location of interest +type(location_type), intent(inout) :: locs(:) ! obs locations +integer, intent(in) :: loc_qtys(:) ! QTYS for obs +integer, intent(in) :: loc_types(:) ! TYPES for obs +integer, intent(out) :: num_close ! how many are close +integer, intent(out) :: close_ind(:) ! incidies into the locs array +real(r8), optional, intent(out) :: dist(:) ! distances in radians +type(ensemble_type), optional, intent(in) :: ens_handle + +character(len=*), parameter :: routine = 'get_close_obs' + +call loc_get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & + num_close, close_ind, dist, ens_handle) + +end subroutine get_close_obs + + +!------------------------------------------------------------------ +! Any model specific distance calcualtion can be done here +subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & + num_close, close_ind, dist, ens_handle) + +type(get_close_type), intent(in) :: gc ! handle to a get_close structure +type(location_type), intent(inout) :: base_loc ! location of interest +integer, intent(in) :: base_type ! observation TYPE +type(location_type), intent(inout) :: locs(:) ! state locations +integer, intent(in) :: loc_qtys(:) ! QTYs for state +integer(i8), intent(in) :: loc_indx(:) ! indices into DART state vector +integer, intent(out) :: num_close ! how many are close +integer, intent(out) :: close_ind(:) ! indices into the locs array +real(r8), optional, intent(out) :: dist(:) ! distances in radians +type(ensemble_type), optional, intent(in) :: ens_handle + +character(len=*), parameter :: routine = 'get_close_state' + + +call loc_get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & + num_close, close_ind, dist, ens_handle) + + +end subroutine get_close_state + + +!------------------------------------------------------------------ +! Does any shutdown and clean-up needed for model. Can be a NULL +! INTERFACE if the model has no need to clean up storage, etc. + +subroutine end_model() + + +end subroutine end_model + + +!------------------------------------------------------------------ +! write any additional attributes to the output and diagnostic files + +subroutine nc_write_model_atts(ncid, domain_id) + +integer, intent(in) :: ncid ! netCDF file identifier +integer, intent(in) :: domain_id + +if ( .not. module_initialized ) call static_init_model + +! put file into define mode. + +call nc_begin_define_mode(ncid) + +call nc_add_global_creation_time(ncid) + +call nc_add_global_attribute(ncid, "model_source", source ) +call nc_add_global_attribute(ncid, "model", "template") + +call nc_end_define_mode(ncid) + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine nc_write_model_atts + +!----------------------------------------------------------------------- +! Parse the table of variables characteristics into arrays for easier access. + +function assign_var(variables, MAX_STATE_VARIABLES) result(var) + +character(len=vtablenamelength), intent(in) :: variables(:, :) +integer, intent(in) :: MAX_STATE_VARIABLES + +type(var_type) :: var +integer :: ivar +character(len=vtablenamelength) :: table_entry + +!----------------------------------------------------------------------- +! Codes for interpreting the NUM_STATE_TABLE_COLUMNS of the variables table +integer, parameter :: NAME_INDEX = 1 ! ... variable name +integer, parameter :: QTY_INDEX = 2 ! ... DART qty +integer, parameter :: MIN_VAL_INDEX = 3 ! ... minimum value if any +integer, parameter :: MAX_VAL_INDEX = 4 ! ... maximum value if any +integer, parameter :: UPDATE_INDEX = 5 ! ... update (state) or not + +! Loop through the variables array to get the actual count of the number of variables +do ivar = 1, MAX_STATE_VARIABLES + ! If the element is an empty string, the loop has exceeded the extent of the variables + if (variables(1, ivar) == '') then + var%count = ivar-1 + exit + endif +enddo + +! Allocate the arrays in the var derived type +allocate(var%names(var%count), var%qtys(var%count), var%clamp_values(var%count, 2), var%updates(var%count)) + +do ivar = 1, var%count + + var%names(ivar) = trim(variables(NAME_INDEX, ivar)) + + table_entry = variables(QTY_INDEX, ivar) + call to_upper(table_entry) + + var%qtys(ivar) = get_index_for_quantity(table_entry) + + if (variables(MIN_VAL_INDEX, ivar) /= 'NA') then + read(variables(MIN_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,1) + else + var%clamp_values(ivar,1) = MISSING_R8 + endif + + if (variables(MAX_VAL_INDEX, ivar) /= 'NA') then + read(variables(MAX_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,2) + else + var%clamp_values(ivar,2) = MISSING_R8 + endif + + table_entry = variables(UPDATE_INDEX, ivar) + call to_upper(table_entry) + + if (table_entry == 'UPDATE') then + var%updates(ivar) = .true. + else + var%updates(ivar) = .false. + endif + +enddo + +end function assign_var + +subroutine read_geometry_file() + + integer :: dimid, varid + character(len=256) :: name + + type(file_type) :: geometry_file + character(len=256) :: restart_directory, grid_directory, filter_directory + namelist /directory_nml/ restart_directory, grid_directory, filter_directory + + call find_namelist_in_file('input.nml', 'directory_nml', iunit) + read(iunit, nml = directory_nml, iostat = io) + call check_namelist_read(iunit, io, 'directory_nml') + + geometry_file%file_path = trim(filter_directory) // 'geometry_file.nc' + + geometry_file%ncid = nc_open_file_readonly(geometry_file%file_path) + + ! attributes + geometry_file%ncstatus = nf90_get_att(geometry_file%ncid, NF90_GLOBAL, 'index_of_north_pole_quad_column', inorth_pole_quad_column) + geometry_file%ncstatus = nf90_get_att(geometry_file%ncid, NF90_GLOBAL, 'index_of_south_pole_quad_column', isouth_pole_quad_column) + + ! dimensions + geometry_file%ncstatus = nf90_inq_dimid(geometry_file%ncid, 'center_altitudes', dimid) + geometry_file%ncstatus = nf90_inquire_dimension(geometry_file%ncid, dimid, name, ncenter_altitudes) + + geometry_file%ncstatus = nf90_inq_dimid(geometry_file%ncid, 'quad_columns', dimid) + geometry_file%ncstatus = nf90_inquire_dimension(geometry_file%ncid, dimid, name, nquad_columns) + + geometry_file%ncstatus = nf90_inq_dimid(geometry_file%ncid, 'center_columns', dimid) + geometry_file%ncstatus = nf90_inquire_dimension(geometry_file%ncid, dimid, name, ncenter_columns) + + ! allocate arrays + allocate(center_altitude(ncenter_altitudes)) + allocate(center_latitude(ncenter_columns)) + allocate(center_longitude(ncenter_columns)) + + allocate(vertex_neighbor_indices(nvertex_columns, nvertex_neighbors)) + allocate(quad_neighbor_indices(nquad_columns, nquad_neighbors)) + + ! variables + geometry_file%ncstatus = nf90_inq_varid(geometry_file%ncid, 'center_altitude', varid) + geometry_file%ncstatus = nf90_get_var(geometry_file%ncid, varid, center_altitude) + + geometry_file%ncstatus = nf90_inq_varid(geometry_file%ncid, 'center_longitude', varid) + geometry_file%ncstatus = nf90_get_var(geometry_file%ncid, varid, center_longitude) + + geometry_file%ncstatus = nf90_inq_varid(geometry_file%ncid, 'center_latitude', varid) + geometry_file%ncstatus = nf90_get_var(geometry_file%ncid, varid, center_latitude) + + geometry_file%ncstatus = nf90_inq_varid(geometry_file%ncid, 'vertex_neighbor_indices', varid) + geometry_file%ncstatus = nf90_get_var(geometry_file%ncid, varid, vertex_neighbor_indices) + + geometry_file%ncstatus = nf90_inq_varid(geometry_file%ncid, 'quad_neighbor_indices', varid) + geometry_file%ncstatus = nf90_get_var(geometry_file%ncid, varid, quad_neighbor_indices) + + call nc_close_file(geometry_file%ncid) + +end subroutine read_geometry_file + +! Barycentric procedures + +subroutine inside_triangle(t1, t2, t3, r, lat, lon, inside, weights) + + ! given 3 corners of a triangle and an xyz point, compute whether + ! the point is inside the triangle. this assumes r is coplanar + ! with the triangle - the caller must have done the lat/lon to + ! xyz conversion with a constant radius and then this will be + ! true (enough). sets inside to true/false, and returns the + ! weights if true. weights are set to 0 if false. + + real(r8), intent(in) :: t1(3), t2(3), t3(3) + real(r8), intent(in) :: r(3), lat, lon + logical, intent(out) :: inside + real(r8), intent(out) :: weights(3) + + ! check for degenerate cases first - is the test point located + ! directly on one of the vertices? (this case may be common + ! if we're computing on grid point locations.) + if (all(abs(r - t1) < roundoff)) then + inside = .true. + weights = (/ 1.0_r8, 0.0_r8, 0.0_r8 /) + return + else if (all(abs(r - t2) < roundoff)) then + inside = .true. + weights = (/ 0.0_r8, 1.0_r8, 0.0_r8 /) + return + else if (all(abs(r - t3) < roundoff)) then + inside = .true. + weights = (/ 0.0_r8, 0.0_r8, 1.0_r8 /) + return + endif + + ! not a vertex. compute the weights. if any are + ! negative, the point is outside. since these are + ! real valued computations define a lower bound for + ! numerical roundoff error and be sure that the + ! weights are not just *slightly* negative. + call get_3d_weights(r, t1, t2, t3, lat, lon, weights) + + if (any(weights < -roundoff)) then + inside = .false. + weights = 0.0_r8 + return + endif + + ! truncate barely negative values to 0 + inside = .true. + where (weights < 0.0_r8) weights = 0.0_r8 + return + +end subroutine inside_triangle + +subroutine get_3d_weights(p, v1, v2, v3, lat, lon, weights) + + ! Given a point p (x,y,z) inside a triangle, and the (x,y,z) + ! coordinates of the triangle corner points (v1, v2, v3), + ! find the weights for a barycentric interpolation. this + ! computation only needs two of the three coordinates, so figure + ! out which quadrant of the sphere the triangle is in and pick + ! the 2 axes which are the least planar: + ! (x,y) near the poles, + ! (y,z) near 0 and 180 longitudes near the equator, + ! (x,z) near 90 and 270 longitude near the equator. + ! (lat/lon are the coords of p. we could compute them here + ! but since in all cases we already have them, pass them + ! down for efficiency) + + real(r8), intent(in) :: p(3) + real(r8), intent(in) :: v1(3), v2(3), v3(3) + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: weights(3) + + real(r8) :: cxs(3), cys(3) + + ! above or below 45 in latitude, where -90 < lat < 90: + if (lat >= 45.0_r8 .or. lat <= -45.0_r8) then + cxs(1) = v1(1) + cxs(2) = v2(1) + cxs(3) = v3(1) + cys(1) = v1(2) + cys(2) = v2(2) + cys(3) = v3(2) + call get_barycentric_weights(p(1), p(2), cxs, cys, weights) + return + endif + + ! nearest 0 or 180 in longitude, where 0 < lon < 360: + if ( lon <= 45.0_r8 .or. lon >= 315.0_r8 .or. & + (lon >= 135.0_r8 .and. lon <= 225.0_r8)) then + cxs(1) = v1(2) + cxs(2) = v2(2) + cxs(3) = v3(2) + cys(1) = v1(3) + cys(2) = v2(3) + cys(3) = v3(3) + call get_barycentric_weights(p(2), p(3), cxs, cys, weights) + return + endif + + ! last option, nearest 90 or 270 in lon: + cxs(1) = v1(1) + cxs(2) = v2(1) + cxs(3) = v3(1) + cys(1) = v1(3) + cys(2) = v2(3) + cys(3) = v3(3) + call get_barycentric_weights(p(1), p(3), cxs, cys, weights) + +end subroutine get_3d_weights + +subroutine get_barycentric_weights(x, y, cxs, cys, weights) + + ! Computes the barycentric weights for a 2d interpolation point + ! (x,y) in a 2d triangle with the given (cxs,cys) corners. + + real(r8), intent(in) :: x, y, cxs(3), cys(3) + real(r8), intent(out) :: weights(3) + + real(r8) :: denom + + ! Get denominator + denom = (cys(2) - cys(3)) * (cxs(1) - cxs(3)) + & + (cxs(3) - cxs(2)) * (cys(1) - cys(3)) + + weights(1) = ((cys(2) - cys(3)) * (x - cxs(3)) + & + (cxs(3) - cxs(2)) * (y - cys(3))) / denom + + weights(2) = ((cys(3) - cys(1)) * (x - cxs(3)) + & + (cxs(1) - cxs(3)) * (y - cys(3))) / denom + + weights(3) = 1.0_r8 - weights(1) - weights(2) + + if (any(abs(weights) < roundoff)) then + where (abs(weights) < roundoff) weights = 0.0_r8 + where (abs(1.0_r8 - abs(weights)) < roundoff) weights = 1.0_r8 + endif + +end subroutine get_barycentric_weights + +subroutine latlon_to_xyz(lat, lon, x, y, z) + + ! Given a lat, lon in degrees, return the cartesian x,y,z coordinate + ! on the surface of a specified radius relative to the origin + ! at the center of the earth. + + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: x, y, z + + real(r8) :: rlat, rlon + + rlat = lat * deg2rad + rlon = lon * deg2rad + + x = radius * cos(rlon) * cos(rlat) + y = radius * sin(rlon) * cos(rlat) + z = radius * sin(rlat) + +end subroutine latlon_to_xyz + +!----------------------------------------------------------------------- +! interpolate in the vertical between 2 arrays of items. + +! vert_fracts: 0 is 100% of the first level and +! 1 is 100% of the second level + +subroutine vert_interp(nitems, levs1, levs2, vert_fract, out_vals) + + integer, intent(in) :: nitems + real(r8), intent(in) :: levs1(nitems) + real(r8), intent(in) :: levs2(nitems) + real(r8), intent(in) :: vert_fract + real(r8), intent(out) :: out_vals(nitems) + + out_vals(:) = (levs1(:) * (1.0_r8 - vert_fract)) + & + (levs2(:) * vert_fract ) + +end subroutine vert_interp + + +function barycentric_average(nitems, weights, vertex_temp_values) result (averaged_values) + + integer, intent(in) :: nitems + real(r8), dimension(nvertex_neighbors), intent(in) :: weights + real(r8), dimension(nvertex_neighbors, nitems), intent(in) :: vertex_temp_values + + real(r8), dimension(nitems) :: averaged_values + + integer :: iweight, iitem + + averaged_values(:) = 0 + + do iitem = 1, nitems + do iweight = 1, nvertex_neighbors + averaged_values(iitem) = averaged_values(iitem) + weights(iweight)*vertex_temp_values(iweight, iitem) + end do + end do + +end function barycentric_average + +!=================================================================== +! End of model_mod +!=================================================================== +end module model_mod diff --git a/models/aether_cube_sphere/readme.rst b/models/aether_cube_sphere/readme.rst new file mode 100644 index 0000000000..9d0cf12940 --- /dev/null +++ b/models/aether_cube_sphere/readme.rst @@ -0,0 +1,174 @@ +Aether Cube Sphere +================== + +This document describes the DART interface to the Aether ionosphere-thermosphere model in its cube +sphere implementation. + +In addition to the standard DART programs associated with a given model, this interface creates +three additional programs for the Aether cube sphere: ``aether_to_dart``, ``dart_to_aether``, and +``create_geometry_file``. + +Use `this link `__ to +download a set of Aether cube sphere restart files. + +aether_to_dart +-------------- + +This program takes six Aether restart files and combines them into a single filter input file for +DART. This program reads entries from two namelists in ``input.nml``: + +In ``&directory_nml``: + +- ``restart_directory`` specifies the path where the restart files are saved +- ``filter_directory`` specifies the path where the ``filter_input`` file will be created (this + entry can be the same as ``restart_directory`` but an additional entry is included for + flexibility) + +In ``&transform_state_nml``: + +- ``restart_file_prefix``, ``restart_file_middle``, and ``restart_file_suffix`` specify substrings + of the restart filename that are combined to create the full + filename. The ensemble member string is inserted in between ``restart_file_prefix`` and + ``restart_file_middle``. The grid face is inserted in between ``restart_file_middle`` and + ``restart_file suffix``. +- ``filter_input_prefix`` and ``filter_input_suffix`` specify substrings of the filter input + filename. The ensemble member string is inserted in between ``filter_input_prefix`` and + ``filter_input_suffix``. + +Using aether_to_dart +~~~~~~~~~~~~~~~~~~~~ + +When executing ``aether_to_dart`` the ensemble member must be specified as a command line argument: + +.. code-block:: + + ./aether_to_dart 0000 + +The program must be run once for each ensemble member, ``./aether_to_dart 0000``, +``./aether_to_dart 0001`` and so on. + +dart_to_aether +-------------- + +This program takes a single filter output file and inserts them back into the six initial Aether +restart files. This program uses all of the same namelist entries as ``aether_to_dart`` except for +instead of using ``filter_input_prefix`` and ``filter_input_suffix`` in ``&transform_state_nml`` it +uses: + +- ``filter_output_prefix`` and ``filter_output_suffix`` specify substrings of the filter output + filename. The ensemble member string is inserted in between ``filter_output_prefix`` and + ``filter_output_suffix``. + +Using dart_to_aether +~~~~~~~~~~~~~~~~~~~~ + +When executing ``aether_to_dart`` the ensemble member must be specified as a command line argument: + +.. code-block:: + + ./dart_to_aether 0000 + +The program must be run once for each ensemble member, ``./dart_to_aether 0000``, +``./dart_to_aether 0001`` and so on. + +create_geometry_file +-------------------- + +This program takes two sets of six Aether grid files: + +- ``grid_corners_g000?.nc``, where ``?`` corresponds to the six grid corner files, numbered 0-5, and +- ``grid_g000?.nc``, where ``?`` corresponds to the six grid center files, numbered 0-5. + +Use `this link `__ to +download a set of Aether cube sphere grid files. + +Aether's model fields are defined at the grid centers. The grid corners are at locations +approximately equidistant from either three (if the grid corner is at at a cube vertex) or four +(if the grid corner is on a cube face) grid centers. + +The program iterates through the grid corners and finds the nearest three or four grid centers in +order to create a single ``geometry_file.nc`` that saves the cube sphere grid topology so that the +center fields can be used for interpolated by DART. + +The two key compontents of the ``geometry_file.nc`` are: + +1. The locations of the eight grid corners that coincide with the cube vertices and the column + indices of the nearest three grid centers to each of the cube vertices. This structure defines a + triangle surrounding each cube vertex. +2. The locations of the remaining ``N`` grid corners that lie on the cube faces and the column + indices of the nearest four grid centers to each of the cube face grid corners. This structure + defines a quadrilateral surrouding each of the grid corners that lie on a cube face. + +.. note:: + + The ``in_quad`` function in DART's ``quad_utilities_mod`` does not implement a method to + determine whether a point is in a quad if the given quad encloses either the North or South Pole. + Thus ``create_geometry_file`` also determines which two grid corners coincide with the North and + South Poles and stores the index of these corners as global attributes in ``geometry_file.nc`` + that can be read and used by the model mod to account for the special case where a quad encloses + a pole. + +The header of the ``geometry_file.nc`` created by ``create_geometry_file`` should resemble the +following: + +.. code-block:: + + netcdf geometry_file { + dimensions: + center_altitudes = 44 ; + vertex_columns = 8 ; + vertex_neighbors = 3 ; + quad_columns = 1938 ; + quad_neighbors = 4 ; + center_columns = 1944 ; + variables: + float center_altitude(center_altitudes) ; + float vertex_longitude(vertex_columns) ; + float vertex_latitude(vertex_columns) ; + float quad_longitude(quad_columns) ; + float quad_latitude(quad_columns) ; + float center_longitude(center_columns) ; + float center_latitude(center_columns) ; + int vertex_neighbor_indices(vertex_neighbors, vertex_columns) ; + int quad_neighbor_indices(quad_neighbors, quad_columns) ; + + // global attributes: + :index_of_north_pole_quad_column = 1760 ; + :index_of_south_pole_quad_column = 1403 ; + } + +geometry_file.nc dimensions +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``vertex_columns`` and ``quad_columns`` correspond to all of the grid corners on the sphere. +Each of the ``vertex_columns`` are enclosed by three grid centers, which are referred to as +``vertex_neighbors``. Each of the ``quad_columns`` are enclosed by four grid centers, which are +refered to as ``quad_neighbors``. + +The dimensions of the grid centers are ``center_altitudes`` in the vertical and ``center_columns`` +in the horizontal. + +geometry_file.nc variables +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The longitudes and latitudes for each of grid corners corresponding to the eight cube vertices are +stored in the ``vertex_longitude`` and ``vertex_latitude`` fields, respectively. + +The longitudes and latitudes for each of the grid corners on the cube faces enclosed by +quadrilaterals are stored in the ``quad_longitude`` and ``quad_latitude`` fields, respectively. + +The altitudes, longitudes and latitudes for each of the grid centers are stored in the +``center_altitude``, ``center_longitude`` and ``center_latitude`` fields, respectively. + +.. important:: + + The key feature of the ``geometry_file.nc`` is the relationship between the grid corners and the + grid centers. This relationship is stored in the ``vertex_neighbor_indices`` and + ``quad_neighbor_indices`` integer fields. + +``vertex_neighbor_indices`` is a 3x8 integer array where each of the 8 columns corresponds to the +a grid corner coinciding with a cube vertex and each of the three rows corresponds to the indices of +the center columns that define a triangle enclosing the cube vertex. +``quad_neighbor_indices`` is a 4xN integer array where each of the N columns corresponds to +a grid corner on a cube face and each of the four rows corresponds to the indices of the center +columns that define a quad that encloses the grid corner. diff --git a/models/aether_cube_sphere/transform_state_mod.f90 b/models/aether_cube_sphere/transform_state_mod.f90 new file mode 100644 index 0000000000..1c97b591d6 --- /dev/null +++ b/models/aether_cube_sphere/transform_state_mod.f90 @@ -0,0 +1,444 @@ +module transform_state_mod + +use netcdf +use types_mod, only : r4, r8, varnamelength +use netcdf_utilities_mod, only : nc_open_file_readonly, nc_open_file_readwrite, nc_close_file, & + nc_create_file, nc_end_define_mode +use utilities_mod, only : open_file, close_file, find_namelist_in_file, & + check_namelist_read, error_handler, E_ERR, string_to_integer + +implicit none +private + +public :: initialize_transform_state_mod, & + finalize_transform_state_mod, & + model_to_dart, & + dart_to_model, & + integer_to_string, & + file_type, & + zero_fill + +integer :: iunit, io + +character(len=4) :: restart_ensemble_member, dart_ensemble_member + +type :: file_type +character(len=256) :: file_path +integer :: ncid, ncstatus, unlimitedDimId, nDimensions, nVariables, nAttributes, formatNum +end type file_type + +type(file_type), allocatable, dimension(:) :: block_files +type(file_type) :: filter_input_file, filter_output_file + +integer :: nblocks, nhalos +character(len=256) :: restart_file_prefix, restart_file_middle, restart_file_suffix, & + filter_input_prefix, filter_input_suffix, filter_output_prefix, & + filter_output_suffix +namelist /transform_state_nml/ nblocks, nhalos, restart_file_prefix, restart_file_middle, & + restart_file_suffix, filter_input_prefix, filter_input_suffix, & + filter_output_prefix, filter_output_suffix + +character(len=256) :: restart_directory, grid_directory, filter_directory +namelist /directory_nml/ restart_directory, grid_directory, filter_directory + +contains + +subroutine initialize_transform_state_mod() + + restart_ensemble_member = get_ensemble_member_from_command_line() + dart_ensemble_member = zero_fill(integer_to_string(string_to_integer(restart_ensemble_member)+1), 4) + + call find_namelist_in_file('input.nml', 'transform_state_nml', iunit) + read(iunit, nml = transform_state_nml, iostat = io) + call check_namelist_read(iunit, io, 'transform_state_nml') + + call find_namelist_in_file('input.nml', 'directory_nml', iunit) + read(iunit, nml = directory_nml, iostat = io) + call check_namelist_read(iunit, io, 'directory_nml') + + block_files = assign_block_files_array(nblocks, restart_ensemble_member, restart_directory, & + restart_file_prefix, restart_file_middle, & + restart_file_suffix) + + + +end subroutine initialize_transform_state_mod + +subroutine finalize_transform_state_mod() + + integer :: iblock + + ! Close all of the files + + do iblock = 1, nblocks + call nc_close_file(block_files(iblock)%ncid) + end do + +end subroutine finalize_transform_state_mod + +subroutine model_to_dart() + + integer :: iblock + integer :: dimid, dart_dimid + integer :: ix, iy, iz, icol + integer :: varid + character(len=NF90_MAX_NAME) :: name + character(len=NF90_MAX_NAME) :: attribute + integer :: length + integer :: xtype, nDimensions, nAtts + integer, dimension(NF90_MAX_VAR_DIMS) :: dimids + integer :: ntimes + integer :: nxs_per_block, nys_per_block, truncated_nxs_per_block, truncated_nys_per_block, total_truncated_ncols + integer :: nzs + + integer, dimension(3) :: time_lev_col_dims + integer, allocatable, dimension (:) :: dart_varids + + ! The time variable in the block files is a double + real(r8), allocatable, dimension(:) :: time_array + ! The other variables are floats + real(r4), allocatable, dimension(:, :, :) :: block_array + real(r4), allocatable, dimension(:) :: spatial_array + real(r4), allocatable, dimension(:, :, :) :: variable_array + + filter_input_file = assign_filter_file(dart_ensemble_member, filter_directory, filter_input_prefix, filter_input_suffix) + + ! The block files are read only + do iblock = 1, nblocks + block_files(iblock)%ncid = nc_open_file_readonly(block_files(iblock)%file_path) + end do + + ! The dart file is create + filter_input_file%ncid = nc_create_file(filter_input_file%file_path) + + ! The first set of nested loops iterates through all of the block files and all of the dimensions + ! of each block file and counts the lengths of each dimension. + + do iblock = 1, nblocks + ! There doesn't seem to be a helper procedure corresponding to nf90_inquire in + ! netcdf_utilities_mod so this uses the external function directly from the netcdf library + block_files(iblock)%ncstatus = nf90_inquire(block_files(iblock)%ncid, & + block_files(iblock)%nDimensions, & + block_files(iblock)%nVariables, & + block_files(iblock)%nAttributes, & + block_files(iblock)%unlimitedDimId, & + block_files(iblock)%formatNum) + + if (iblock == 1) then + do dimid = 1, block_files(iblock)%nDimensions + ! There doesn't seem to be a helper procedure corresponding to nf90_inquire_dimension that + ! assigns name and length in netcdf_utilities_mod so this uses the external function + ! directly from the netcdf library + block_files(iblock)%ncstatus = nf90_inquire_dimension(block_files(iblock)%ncid, dimid, name, length) + + if (trim(name) == 'time') then + ntimes = length + else if (trim(name) == 'x') then + truncated_nxs_per_block = length-2*nhalos + nxs_per_block = length + else if (trim(name) == 'y') then + truncated_nys_per_block = length-2*nhalos + nys_per_block = length + else if (trim(name) == 'z') then + nzs = length + end if + end do + end if + end do + + total_truncated_ncols = truncated_nxs_per_block*truncated_nys_per_block*nblocks + + ! All of the lengths have been counted properly, create each dimension in the filter_input_file and save + ! the dimensions to the time_x_y_z and x_y_z arrays used during variable definition + filter_input_file%ncstatus = nf90_def_dim(filter_input_file%ncid, 'time', NF90_UNLIMITED, dart_dimid) + time_lev_col_dims(3) = dart_dimid + + filter_input_file%ncstatus = nf90_def_dim(filter_input_file%ncid, 'z', nzs, dart_dimid) + time_lev_col_dims(2) = dart_dimid + + filter_input_file%ncstatus = nf90_def_dim(filter_input_file%ncid, 'col', total_truncated_ncols, dart_dimid) + time_lev_col_dims(1) = dart_dimid + + ! Allocate all of the storage arrays + allocate(time_array(ntimes)) + allocate(block_array(nzs, nys_per_block, nxs_per_block)) + allocate(spatial_array(total_truncated_ncols)) + allocate(variable_array(total_truncated_ncols, nzs, ntimes)) + allocate(dart_varids(block_files(1)%nVariables)) + + block_array(:, :, :) = 0 + spatial_array(:) = 0 + variable_array(:, :, :) = 0 + + ! The filter_input_file is still in define mode. Create all of the variables before entering data mode. + do varid = 1, block_files(1)%nVariables + block_files(1)%ncstatus = nf90_inquire_variable(block_files(1)%ncid, varid, name, xtype, nDimensions, dimids, nAtts) + if (trim(name) == 'time') then + filter_input_file%ncstatus = nf90_def_var(filter_input_file%ncid, name, xtype, time_lev_col_dims(3), dart_varids(varid)) + else if (trim(name) == 'z') then + ! Rename the 'z' variable as 'alt' so there isn't a dimension and a variable with the same name + filter_input_file%ncstatus = nf90_def_var(filter_input_file%ncid, 'alt', xtype, time_lev_col_dims(2), dart_varids(varid)) + else if ((trim(name) == 'lon') .or. (trim(name) == 'lat')) then + filter_input_file%ncstatus = nf90_def_var(filter_input_file%ncid, name, xtype, time_lev_col_dims(1), dart_varids(varid)) + else + filter_input_file%ncstatus = nf90_def_var(filter_input_file%ncid, name, xtype, time_lev_col_dims, dart_varids(varid)) + end if + + ! In the block files, time does not have units + if (trim(name) /= 'time') then + block_files(iblock)%ncstatus = nf90_get_att(block_files(1)%ncid, varid, 'units', attribute) + filter_input_file%ncstatus = nf90_put_att(filter_input_file%ncid, dart_varids(varid), 'units', attribute) + end if + + ! In the block files, only lon, lat and z have long_name + if ((trim(name) == 'lon') .or. (trim(name) == 'lat') .or. (trim(name) == 'z')) then + block_files(iblock)%ncstatus = nf90_get_att(block_files(1)%ncid, varid, 'long_name', attribute) + filter_input_file%ncstatus = nf90_put_att(filter_input_file%ncid, dart_varids(varid), 'long_name', attribute) + end if + + ! print *, 'name: ' // name + ! print *, 'dart_varids(varid): ' // integer_to_string(dart_varids(varid)) + + end do + + call nc_end_define_mode(filter_input_file%ncid) + + ! The second set of nested loops has a different loop order. The outer loop is all of the + ! variables while the inner loop is all of the blocks. The order is switched because all of the + ! ncid pointers to each of the block files have already been assigned and it is more + ! straightforward to assign all of the elements in the variable arrays if the blocks are the + ! inner loop. + + do varid = 1, block_files(1)%nVariables + icol = 0 + do iblock = 1, nblocks + block_files(iblock)%ncstatus = nf90_inquire_variable(block_files(iblock)%ncid, varid, name, xtype, nDimensions, dimids, nAtts) + + if (trim(name) == 'time') then + ! This is a 1-D time array + if (iblock == 1) then + block_files(iblock)%ncstatus = nf90_get_var(block_files(iblock)%ncid, varid, time_array) + filter_input_file%ncstatus = nf90_put_var(filter_input_file%ncid, dart_varids(varid), time_array) + end if + else if (trim(name) == 'z') then + if (iblock == 1) then + block_files(iblock)%ncstatus = nf90_get_var(block_files(iblock)%ncid, varid, block_array) + filter_input_file%ncstatus = nf90_put_var(filter_input_file%ncid, dart_varids(varid), block_array(:,1,1)) + end if + else + ! All of the variables besides time can be read into the block array + block_files(iblock)%ncstatus = nf90_get_var(block_files(iblock)%ncid, varid, block_array) + + if ((trim(name) == 'lon') .or. (trim(name) == 'lat')) then + do iy = 1, truncated_nys_per_block + do ix = 1, truncated_nxs_per_block + icol = icol + 1 + spatial_array(icol) = block_array(1, nhalos+iy, nhalos+ix) + end do + end do + + if (iblock == nblocks) then + filter_input_file%ncstatus = nf90_put_var(filter_input_file%ncid, dart_varids(varid), spatial_array) + end if + else + ! This is one of the other non-spatial variables + + do iy = 1, truncated_nys_per_block + do ix = 1, truncated_nxs_per_block + icol = icol + 1 + do iz = 1, nzs + variable_array(icol, iz, 1) = block_array(iz, nhalos+iy, nhalos+ix) + end do + end do + end do + + if (iblock == nblocks) then + filter_input_file%ncstatus = nf90_put_var(filter_input_file%ncid, dart_varids(varid), variable_array) + end if + end if + end if + end do + end do + + call nc_close_file(filter_input_file%ncid) + +end subroutine model_to_dart + +subroutine dart_to_model() + + integer :: iblock, icol, ix, iy, iz + integer :: ntimes, nxs, nys, nzs, ncols + integer :: filter_varid, block_varid, dimid + character(len=NF90_MAX_NAME) :: filter_name, block_name + integer :: filter_xtype, filter_nDimensions, filter_nAtts + integer, dimension(NF90_MAX_VAR_DIMS) :: filter_dimids + real(r4), allocatable, dimension(:, :, :) :: filter_array, block_array + + filter_output_file = assign_filter_file(dart_ensemble_member, filter_directory, filter_output_prefix, filter_output_suffix) + + ! The block files are read/write + do iblock = 1, nblocks + block_files(iblock)%ncid = nc_open_file_readwrite(block_files(iblock)%file_path) + end do + ! The dart file is read only + filter_output_file%ncid = nc_open_file_readonly(filter_output_file%file_path) + + ! Get the variables list from the filter_output_file + + filter_output_file%ncstatus = nf90_inquire(filter_output_file%ncid, & + filter_output_file%nDimensions, & + filter_output_file%nVariables, & + filter_output_file%nAttributes, & + filter_output_file%unlimitedDimId, & + filter_output_file%formatNum) + + filter_output_file%ncstatus = nf90_inq_dimid(filter_output_file%ncid, 'time', dimid) + filter_output_file%ncstatus = nf90_inquire_dimension(filter_output_file%ncid, dimid, filter_name, ntimes) + + filter_output_file%ncstatus = nf90_inq_dimid(filter_output_file%ncid, 'z', dimid) + filter_output_file%ncstatus = nf90_inquire_dimension(filter_output_file%ncid, dimid, filter_name, nzs) + + filter_output_file%ncstatus = nf90_inq_dimid(filter_output_file%ncid, 'col', dimid) + filter_output_file%ncstatus = nf90_inquire_dimension(filter_output_file%ncid, dimid, filter_name, ncols) + + allocate(filter_array(ncols, nzs, ntimes)) + filter_array(:, :, :) = 0 + + ! We need full blocks from the block files + do filter_varid = 1, filter_output_file%nVariables + icol = 0 + filter_output_file%ncstatus = nf90_inquire_variable(filter_output_file%ncid, filter_varid, filter_name, filter_xtype, filter_nDimensions, filter_dimids, filter_nAtts) + + if (filter_name /= 'time') then + + filter_output_file%ncstatus = nf90_get_var(filter_output_file%ncid, filter_varid, filter_array) + + do iblock = 1, nblocks + + if (filter_varid == 1 .and. iblock == 1) then + block_files(iblock)%ncstatus = nf90_inq_dimid(block_files(iblock)%ncid, 'x', dimid) + block_files(iblock)%ncstatus = nf90_inquire_dimension(block_files(iblock)%ncid, dimid, block_name, nxs) + + block_files(iblock)%ncstatus = nf90_inq_dimid(block_files(iblock)%ncid, 'y', dimid) + block_files(iblock)%ncstatus = nf90_inquire_dimension(block_files(iblock)%ncid, dimid, block_name, nys) + + allocate(block_array(nzs, nys, nxs)) + block_array(:, :, :) = 0 + end if + + block_files(iblock)%ncstatus = nf90_inq_varid(block_files(iblock)%ncid, filter_name, block_varid) + block_files(iblock)%ncstatus = nf90_get_var(block_files(iblock)%ncid, block_varid, block_array) + + do iy = 1, nys-2*nhalos + do ix = 1, nxs-2*nhalos + icol = icol + 1 + do iz = 1, nzs + block_array(iz, nhalos+iy, nhalos+ix) = filter_array(icol, iz, 1) + end do + end do + end do + + block_files(iblock)%ncstatus = nf90_put_var(block_files(iblock)%ncid, block_varid, block_array) + print *, block_files(iblock)%ncstatus + + end do + end if + end do + + call nc_close_file(filter_output_file%ncid) + +end subroutine dart_to_model + +function get_ensemble_member_from_command_line() result(ensemble_member) + ! Calls Fortran intrinsic subroutine get_command_argument and returns + ! a string with four characters + + character(len=4) :: ensemble_member + integer :: nargs + + nargs = command_argument_count() + + if (nargs /= 1) then + call error_handler(E_ERR, 'get_ensemble_member_from_command_line', & + 'ensemble member must be passed as a command line argument') + end if + + call get_command_argument(1, ensemble_member) + +end function get_ensemble_member_from_command_line + +function assign_block_files_array(nblocks, ensemble_member, restart_directory, & + restart_file_prefix, restart_file_middle, restart_file_suffix) & + result(block_files) + + integer, intent(in) :: nblocks + character(len=4), intent(in) :: ensemble_member + character(len=*), intent(in) :: restart_directory + character(len=*), intent(in) :: restart_file_prefix + character(len=*), intent(in) :: restart_file_middle + character(len=*), intent(in) :: restart_file_suffix + type(file_type), allocatable, dimension(:) :: block_files + character(len=4) :: block_name + integer :: iblock + + allocate(block_files(nblocks)) + + do iblock = 1, nblocks + block_name = zero_fill(integer_to_string(iblock-1), 4) + block_files(iblock)%file_path = trim(restart_directory) // trim(restart_file_prefix) // & + ensemble_member // trim(restart_file_middle) // & + block_name // trim(restart_file_suffix) + end do + +end function assign_block_files_array + +function assign_filter_file(ensemble_member, filter_directory, filter_input_prefix, filter_input_suffix) & + result(filter_file) + + character(len=4), intent(in) :: ensemble_member + character(len=*), intent(in) :: filter_directory + character(len=*), intent(in) :: filter_input_prefix + character(len=*), intent(in) :: filter_input_suffix + type(file_type) :: filter_file + + filter_file%file_path = trim(filter_directory) // trim(filter_input_prefix) // ensemble_member // trim(filter_input_suffix) + +end function assign_filter_file + +function integer_to_string(int) result(string) + + integer, intent(in) :: int + character(len=varnamelength) :: string + + write(string,'(I0)') int + string = trim(string) + +end function integer_to_string + +function zero_fill(string, desired_length) result(filled_string) + + character(len=*), intent(in) :: string + integer, intent(in) :: desired_length + + character(len=varnamelength) :: filled_string + integer :: length_of_string + integer :: string_index, difference_of_string_lengths + + filled_string = '' + length_of_string = len_trim(string) + difference_of_string_lengths = desired_length - length_of_string + + if (difference_of_string_lengths < 0) then + print *, 'Error: input string is longer than the desired output string.' + stop + else if (difference_of_string_lengths > 0) then + do string_index = 1, difference_of_string_lengths + filled_string(string_index:string_index) = '0' + end do + end if + + filled_string(difference_of_string_lengths+1:desired_length) = trim(string) + +end function zero_fill + +end module transform_state_mod diff --git a/models/aether_cube_sphere/work/input.nml b/models/aether_cube_sphere/work/input.nml new file mode 100644 index 0000000000..3b1f6dfecb --- /dev/null +++ b/models/aether_cube_sphere/work/input.nml @@ -0,0 +1,242 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + +&perfect_model_obs_nml + read_input_state_from_file = .true., + single_file_in = .true. + input_state_files = "perfect_input.nc" + + write_output_state_to_file = .true., + single_file_out = .true. + output_state_files = "perfect_output.nc" + output_interval = 1, + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_seq_in_file_name = "obs_seq.in", + obs_seq_out_file_name = "obs_seq.out", + init_time_days = 0, + init_time_seconds = 0, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + trace_execution = .false., + output_timestamps = .false., + print_every_nth_obs = -1, + output_forward_op_errors = .false., + silence = .false., + / + +&filter_nml + single_file_in = .false., + input_state_files = '' + input_state_file_list = 'filter_input_files.txt' + + stages_to_write = 'analysis' + + single_file_out = .false., + output_state_files = '' + output_state_file_list = 'filter_output_files.txt' + output_interval = 1, + output_members = .true. + num_output_state_members = 3, + output_mean = .true. + output_sd = .true. + write_all_stages_at_end = .false. + + ens_size = 3, + num_groups = 1, + perturb_from_single_instance = .false., + perturbation_amplitude = 0.2, + distributed_state = .true. + + async = 0, + adv_ens_command = "./advance_model.csh", + + obs_sequence_in_name = "obs_seq.out.sample", + obs_sequence_out_name = "obs_seq.final", + num_output_obs_members = 3, + init_time_days = 153130, + init_time_seconds = 84600, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + + inf_flavor = 0, 0, + inf_initial_from_restart = .false., .false., + inf_sd_initial_from_restart = .false., .false., + inf_deterministic = .true., .true., + inf_initial = 1.0, 1.0, + inf_lower_bound = 1.0, 1.0, + inf_upper_bound = 100.0, 1000000.0, + inf_damping = 1.0, 1.0, + inf_sd_initial = 0.0, 0.0, + inf_sd_lower_bound = 0.0, 0.0, + inf_sd_max_change = 1.05, 1.05, + + trace_execution = .false., + output_timestamps = .false., + output_forward_op_errors = .false., + silence = .false., + / + + +&ensemble_manager_nml + / + +&assim_tools_nml + cutoff = 1000000.0 + sort_obs_inc = .false., + spread_restoration = .false., + sampling_error_correction = .false., + adaptive_localization_threshold = -1, + distribute_mean = .false. + output_localization_diagnostics = .false., + localization_diagnostics_file = 'localization_diagnostics', + print_every_nth_obs = 0 + / + +&cov_cutoff_nml + select_localization = 1 + / + +®_factor_nml + select_regression = 1, + input_reg_file = "time_mean_reg", + save_reg_diagnostics = .false., + reg_diagnostics_file = "reg_diagnostics" + / + +&obs_sequence_nml + write_binary_obs_sequence = .false. + / + +&obs_kind_nml + assimilate_these_obs_types = 'AIRS_TEMPERATURE' + evaluate_these_obs_types = '' + / + +&model_nml + template_file = 'aether_restart_m0001.nc' + time_step_days = 0, + time_step_seconds = 3600 + variables = 'Temperature', 'QTY_TEMPERATURE', '0.0', 'NA', 'UPDATE', + 'O+', 'QTY_DENSITY_ION_OP', '0.0', 'NA', 'UPDATE' + / + +&utilities_nml + TERMLEVEL = 1, + module_details = .false., + logfilename = 'dart_log.out', + nmlfilename = 'dart_log.nml', + write_nml = 'none' + / + +&preprocess_nml + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_upper_atm_mod.f90', + '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + '../../../observations/forward_operators/obs_def_gps_mod.f90', + '../../../observations/forward_operators/obs_def_vortex_mod.f90', + '../../../observations/forward_operators/obs_def_gts_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90', + '../../../assimilation_code/modules/observations/space_quantities_mod.f90', + '../../../assimilation_code/modules/observations/chemistry_quantities_mod.f90' + / + +&obs_sequence_tool_nml + filename_seq = 'obs_seq.one', 'obs_seq.two', + filename_out = 'obs_seq.processed', + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + print_only = .false., + gregorian_cal = .false. + / + +&obs_diag_nml + obs_sequence_name = 'obs_seq.final', + bin_width_days = -1, + bin_width_seconds = -1, + init_skip_days = 0, + init_skip_seconds = 0, + Nregions = 3, + trusted_obs = 'null', + lonlim1 = 0.00, 0.00, 0.50 + lonlim2 = 1.01, 0.50, 1.01 + reg_names = 'whole', 'yin', 'yang' + create_rank_histogram = .true., + outliers_in_histogram = .true., + use_zero_error_obs = .false., + verbose = .false. + / + +&state_vector_io_nml + / + +&model_mod_check_nml + input_state_files = 'aether_restart_m0001.nc' + output_state_files = 'mmc_output.nc' + test1thru = 0 + run_tests = 1,2,3,4,5,6,7 + x_ind = 12 + loc_of_interest = 98.5, 85.5, 96344 + quantity_of_interest = 'QTY_TEMPERATURE' + interp_test_dlon = 4.0 + interp_test_lonrange = 0.0, 360.0 + interp_test_dlat = 4.0 + interp_test_latrange = -90.0, 90.0 + interp_test_dvert = 2000.0 + interp_test_vertrange = 100000.0, 106000.0 + interp_test_vertcoord = 'VERTISHEIGHT' + verbose = .false. + / + +&quality_control_nml + input_qc_threshold = 3.0, + outlier_threshold = -1.0, +/ + +&location_nml + / + +&directory_nml + restart_directory = './' + grid_directory = './' + filter_directory = './' + / + +&grid_nml + grid_centers_file_prefix = 'grid_g' + grid_centers_file_suffix = '.nc' + grid_corners_file_prefix = 'grid_corners_g' + grid_corners_file_suffix = '.nc' + / + +&transform_state_nml + nblocks = 6 + nhalos = 2 + restart_file_prefix = 'aether_restart_m' + restart_file_middle = '_g' + restart_file_suffix = '.nc' + filter_input_prefix = 'aether_restart_m' + filter_input_suffix = '.nc' + filter_output_prefix = 'analysis_member_' + filter_output_suffix = '.nc' + / diff --git a/models/aether_cube_sphere/work/obs_seq.out.sample b/models/aether_cube_sphere/work/obs_seq.out.sample new file mode 100644 index 0000000000..1cd14ccb15 --- /dev/null +++ b/models/aether_cube_sphere/work/obs_seq.out.sample @@ -0,0 +1,31 @@ + obs_sequence +obs_kind_definitions + 1 + 33 AIRS_TEMPERATURE + num_copies: 1 num_qc: 1 + num_obs: 2 max_num_obs: 2 +observation +Data QC + first: 1 last: 2 + OBS 1 + 271.330627441406 + 0.000000000000000E+000 + -1 2 -1 +obdef +loc3d + 3.406717740263719 0.5806184282903090 100000.0000000000 3 +kind + 33 +84601 153130 + 1.07229244766182 + OBS 2 + 27450.2966235645 + 0.000000000000000E+000 + 1 -1 -1 +obdef +loc3d + 3.484538383406885 0.5925166389933947 120000.0000000000 3 +kind + 33 +84601 153130 + 1.03153675838621 \ No newline at end of file diff --git a/models/aether_cube_sphere/work/quickbuild.sh b/models/aether_cube_sphere/work/quickbuild.sh new file mode 100755 index 0000000000..6f4e628fa0 --- /dev/null +++ b/models/aether_cube_sphere/work/quickbuild.sh @@ -0,0 +1,62 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL=aether_cube_sphere +LOCATION=threed_sphere + + +programs=( +closest_member_tool +filter +model_mod_check +perfect_model_obs +) + +serial_programs=( +create_fixed_network_seq +create_obs_sequence +fill_inflation_restart +integrate_model +obs_common_subset +obs_diag +obs_sequence_tool +) + +model_programs=( +) + +model_serial_programs=( +aether_to_dart +dart_to_aether +create_geometry_file +) + +# quickbuild arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build any NetCDF files from .cdl files +cdl_to_netcdf + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@" diff --git a/models/utilities/quad_utils_mod.f90 b/models/utilities/quad_utils_mod.f90 index b9d11e6334..d7173aeb33 100644 --- a/models/utilities/quad_utils_mod.f90 +++ b/models/utilities/quad_utils_mod.f90 @@ -78,7 +78,10 @@ module quad_utils_mod QUAD_LOCATED_CELL_CORNERS, & get_quad_grid_size, & get_quad_global, & - print_quad_handle ! debug + print_quad_handle, & ! debug + in_quad, & + quad_bilinear_interp, & + line_intercept ! version controlled file description for error handling, do not edit