diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 267fb5bd63..5756765af8 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,17 @@ individual files. The changes are now listed with the most recent at the top. +**February 7 2022 :: CM1 and 3D Cartesian location_mod updates. Tag v9.13.1** + +*Contributed by Jon Labriola* + +- Updated CM1 model_mod to use mixed-case boundary conditions, for example + periodic in the x-direction but non-periodic in the y-direction. +- Added capability to CM1 model_mod to interpolate 3D fields such as reflectivity. +- Added capability to use multiple localization radii to threed_cartisian + location_mod. +- Bug-fix for threed_cartesian location_mod for periodic boundaries. + **February 3 2022 :: CLM with SWE repartitioning. Tag: v9.13.0** - Updated Community Land Model (CLM) model_mod, scripting, and diagnostics. diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 8011b9c4a5..2f532a2fe5 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -16,6 +16,7 @@ module location_mod do_nml_term, is_longitude_between use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform use mpi_utilities_mod, only : my_task_id, task_count +use obs_kind_mod, only : get_num_types_of_obs, get_name_for_type_of_obs, get_index_for_type_of_obs use ensemble_manager_mod, only : ensemble_type use default_location_mod, only : has_vertical_choice, vertical_localization_on, & get_vertical_localization_coord, & @@ -69,9 +70,11 @@ module location_mod type box_type private - integer, pointer :: loc_box(:) ! (nloc); List of loc indices in boxes - integer, pointer :: count(:, :, :) ! (nx, ny, nz); # of locs in each box - integer, pointer :: start(:, :, :) ! (nx, ny, nz); Start of list of locs in this box + integer :: num + real(r8) :: maxdist ! furthest seperation between "close" locations + integer, allocatable :: loc_box(:) ! (nloc); List of loc indices in boxes + integer, allocatable :: count(:, :, :) ! (nx, ny, nz); # of locs in each box + integer, allocatable :: start(:, :, :) ! (nx, ny, nz); Start of list of locs in this box real(r8) :: bot_x, top_x ! extents in x, y, z real(r8) :: bot_y, top_y real(r8) :: bot_z, top_z @@ -79,20 +82,23 @@ module location_mod real(r8) :: nboxes_x, nboxes_y, nboxes_z ! based on maxdist how far to search end type box_type -! Type to facilitate efficient computation of observations close to a given location + + type get_close_type private - integer :: num - real(r8) :: maxdist - type(box_type) :: box + integer :: nt ! The number of distinct cutoffs + !real(r8),allocatable :: type_to_cutoff_map(:) ! mapping of types to index + integer, allocatable :: type_to_cutoff_map(:) ! mapping of types to index + type(box_type),allocatable :: box(:) ! Array of box types end type get_close_type + type(random_seq_type) :: ran_seq logical :: ran_seq_init = .false. logical, save :: module_initialized = .false. character(len = 512) :: errstring - +character(len = 512) :: msgstring, msgstring1, msgstring2 ! for sanity when i'm using arrays of length(3): integer, parameter :: IX = 1 @@ -126,11 +132,6 @@ module location_mod type(periodic_info_type) :: loopy -!> @todo: this code doesn't do the by-type variable maxdist computations. -!> need to lift it from the threed_sphere version. makes 'box' inside -!> get_close_type an array. (good thing i didn't remove it yet.) -!> also makes 'maxdist' an array. - !----------------------------------------------------------------- ! Namelist with default values @@ -252,17 +253,18 @@ function get_dist(loc1, loc2, type1, kind2) diff(IY) = loc1%y - loc2%y diff(IZ) = loc1%z - loc2%z -! hope for the simple case first +! the simple non-periodic case if (.not. any_periodic) then if (debug > 0) write(0,*) 'non-periodic distance: ' get_dist = dist_3d(diff) return -else - if (debug > 0) write(0,*) 'periodic distance: ' - square_dist = dist_3d_sq(diff) - if (debug > 0) write(0,*) 'non-periodic 0 sq_dist, dist: ', square_dist, sqrt(square_dist) endif +! everything below deals with the periodic options +if (debug > 0) write(0,*) 'periodic distance: ' +square_dist = dist_3d_sq(diff) +if (debug > 0) write(0,*) 'non-periodic 0 sq_dist, dist: ', square_dist, sqrt(square_dist) + ! ok, not simple. !> @todo: can we separate this? @@ -272,7 +274,6 @@ function get_dist(loc1, loc2, type1, kind2) !> for all combinations of location 1 and location 2 where !> they aren't both above or both below the midpoint add the !> offset, compute the distance, and keep the minimium - below_L1 = find_my_quadrant(loc1) if (debug > 0) write(0,*) 'below_L1, find my quadrant: ', below_L1 below_L2 = find_my_quadrant(loc2) @@ -282,13 +283,31 @@ function get_dist(loc1, loc2, type1, kind2) if (debug > 0) write(0,*) 'in the xyz_periodic case' - if ((below_L1(IX) .neqv. below_L2(IX)) .and. & - (below_L1(IY) .neqv. below_L2(IY)) .and. & + if ((below_L1(IX) .neqv. below_L2(IX)) .or. & + (below_L1(IY) .neqv. below_L2(IY)) .or. & (below_L1(IZ) .neqv. below_L2(IZ))) then next_diff = diff - next_diff(IX) = next_diff(IX) + loopy%offset(IX) - next_diff(IY) = next_diff(IY) + loopy%offset(IY) - next_diff(IZ) = next_diff(IZ) + loopy%offset(IZ) + if (below_L1(IX) .neqv. below_L2(IX)) then ! Adjust x + if (loc1%x > loopy%midline(IX)) then + next_diff(IX) = loopy%offset(IX) - next_diff(IX) + else + next_diff(IX) = loopy%offset(IX) + next_diff(IX) + endif + endif + if (below_L1(IY) .neqv. below_L2(IY)) then ! Adjust y + if (loc1%y > loopy%midline(IY)) then + next_diff(IY) = loopy%offset(IY) - next_diff(IY) + else + next_diff(IY) = loopy%offset(IY) + next_diff(IY) + endif + endif + if (below_L1(IZ) .neqv. below_L2(IZ)) then ! Adjust z + if (loc1%z > loopy%midline(IZ)) then + next_diff(IZ) = loopy%offset(IZ) - next_diff(IZ) + else + next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) + endif + endif this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic XYZ dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -304,11 +323,24 @@ function get_dist(loc1, loc2, type1, kind2) if (debug > 0) write(0,*) 'in the xy_periodic case' - if ((below_L1(IX) .neqv. below_L2(IX)) .and. & + if ((below_L1(IX) .neqv. below_L2(IX)) .or. & (below_L1(IY) .neqv. below_L2(IY))) then next_diff = diff - next_diff(IX) = next_diff(IX) + loopy%offset(IX) - next_diff(IY) = next_diff(IY) + loopy%offset(IY) + if (below_L1(IX) .neqv. below_L2(IX)) then + if (loc1%x > loopy%midline(IX)) then + next_diff(IX) = loopy%offset(IX) - next_diff(IX) + else + next_diff(IX) = loopy%offset(IX) + next_diff(IX) + endif + endif + if (below_L1(IY) .neqv. below_L2(IY)) then + if (loc1%y > loopy%midline(IY)) then + next_diff(IY) = loopy%offset(IY) - next_diff(IY) + else + next_diff(IY) = loopy%offset(IY) + next_diff(IY) + endif + endif + this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic XY dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -320,15 +352,28 @@ function get_dist(loc1, loc2, type1, kind2) endif + if (xz_periodic) then if (debug > 0) write(0,*) 'in the xz_periodic case' - if ((below_L1(IX) .neqv. below_L2(IX)) .and. & + if ((below_L1(IX) .neqv. below_L2(IX)) .or. & (below_L1(IZ) .neqv. below_L2(IZ))) then next_diff = diff - next_diff(IX) = next_diff(IX) + loopy%offset(IX) - next_diff(IZ) = next_diff(IZ) + loopy%offset(IZ) + if (below_L1(IX) .neqv. below_L2(IX)) then ! Adjust x + if (loc1%x > loopy%midline(IX)) then + next_diff(IX) = loopy%offset(IX) - next_diff(IX) + else + next_diff(IX) = loopy%offset(IX) + next_diff(IX) + endif + endif + if (below_L1(IZ) .neqv. below_L2(IZ)) then ! Adjust z + if (loc1%z > loopy%midline(IZ)) then + next_diff(IZ) = loopy%offset(IZ) - next_diff(IZ) + else + next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) + endif + endif this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic XZ dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -344,11 +389,23 @@ function get_dist(loc1, loc2, type1, kind2) if (debug > 0) write(0,*) 'in the yz_periodic case' - if ((below_L1(IY) .neqv. below_L2(IY)) .and. & + if ((below_L1(IY) .neqv. below_L2(IY)) .or. & (below_L1(IZ) .neqv. below_L2(IZ))) then next_diff = diff - next_diff(IY) = next_diff(IY) + loopy%offset(IY) - next_diff(IZ) = next_diff(IZ) + loopy%offset(IZ) + if (below_L1(IY) .neqv. below_L2(IY)) then ! Adjust y + if (loc1%y > loopy%midline(IY)) then + next_diff(IY) = loopy%offset(IY) - next_diff(IY) + else + next_diff(IY) = loopy%offset(IY) + next_diff(IY) + endif + endif + if (below_L1(IZ) .neqv. below_L2(IZ)) then ! Adjust z + if (loc1%z > loopy%midline(IZ)) then + next_diff(IZ) = loopy%offset(IZ) - next_diff(IZ) + else + next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) + endif + endif this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic YZ dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -366,7 +423,11 @@ function get_dist(loc1, loc2, type1, kind2) if (below_L1(IX) .neqv. below_L2(IX)) then next_diff = diff - next_diff(IX) = next_diff(IX) + loopy%offset(IX) + if (loc1%x > loopy%midline(IX)) then + next_diff(IX) = loopy%offset(IX) - next_diff(IX) + else + next_diff(IX) = loopy%offset(IX) + next_diff(IX) + endif this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic X dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -379,13 +440,18 @@ function get_dist(loc1, loc2, type1, kind2) endif if (y_periodic) then - if (debug > 0) write(0,*) 'in the y_periodic case' if (below_L1(IY) .neqv. below_L2(IY)) then next_diff = diff - next_diff(IY) = next_diff(IY) + loopy%offset(IY) + if (loc1%y > loopy%midline(IY)) then + next_diff(IY) = loopy%offset(IY) - next_diff(IY) + else + next_diff(IY) = loopy%offset(IY) + next_diff(IY) + endif + this_dist = dist_3d_sq(next_diff) + if (debug > 0) write(0,*) 'periodic Y dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -402,7 +468,11 @@ function get_dist(loc1, loc2, type1, kind2) if (below_L1(IZ) .neqv. below_L2(IZ)) then next_diff = diff - next_diff(IZ) = next_diff(IZ) + loopy%offset(IZ) + if (loc1%z > loopy%midline(IZ)) then + next_diff(IZ) = loopy%offset(IZ) - next_diff(IZ) + else + next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) + endif this_dist = dist_3d_sq(next_diff) if (debug > 0) write(0,*) 'periodic Y dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -805,36 +875,97 @@ end subroutine interactive_location ! Initializes get_close accelerator subroutine get_close_init(gc, num, maxdist, locs, maxdist_list) - type(get_close_type), intent(inout) :: gc integer, intent(in) :: num real(r8), intent(in) :: maxdist type(location_type), intent(in) :: locs(:) real(r8), intent(in), optional :: maxdist_list(:) -integer :: i, j, k, cum_start, l +integer :: i, j, k, cum_start, l, n integer :: x_box(num), y_box(num), z_box(num) integer :: tstart(nx, ny, nz) + +integer :: typecount, distcount +real(r8), allocatable :: distlist(:) + + if ( .not. module_initialized ) call initialize_module +typecount = get_num_types_of_obs() +allocate(gc%type_to_cutoff_map(typecount)) + +if (present(maxdist_list)) then + if (size(maxdist_list) .ne. typecount) then + write(msgstring,'(A,I8,A,I8)')'maxdist_list len must equal number of specific types, ', & + size(maxdist_list), ' /= ', typecount + call error_handler(E_ERR, 'get_close_init', msgstring, source) + endif + + allocate(distlist(typecount)) + call distinct_values(maxdist_list, distcount, distlist, gc%type_to_cutoff_map) + gc%nt = distcount + if (gc%nt <= 0) then + write(msgstring,'(A)')'error getting count of distinct cutoff dists; should not happen' + call error_handler(E_ERR, 'get_close_init', msgstring, source) + endif +else + gc%nt = 1 + gc%type_to_cutoff_map(:) = 1 +endif + +allocate(gc%box(gc%nt)) + +if (present(maxdist_list)) then + do i=1, gc%nt + gc%box(i)%maxdist = distlist(i) + enddo +else + ! no per-type settings, everyone uses same distance + gc%box(1)%maxdist = maxdist +endif + +if (present(maxdist_list)) deallocate(distlist) + +! Allocate the storage for the grid dependent boxes +do i=1, gc%nt + allocate(gc%box(i)%count(nx, ny, nz), gc%box(i)%start(nx, ny, nz)) + gc%box(i)%count = -1 + gc%box(i)%start = -1 +enddo + +! store the location counts in all derived types +do i=1, gc%nt + gc%box(i)%num = num +enddo + + +! If there are no locs to operate on, no point in going any further. +if (num == 0) return + ! Set the maximum localization distance -gc%maxdist = maxdist +!gc%maxdist = maxdist ! Allocate storage for number dependent part -allocate(gc%box%loc_box(num)) -gc%box%loc_box(:) = -1 +!allocate(gc%box%loc_box(num)) +!gc%box%loc_box(:) = -1 ! Allocate the storage for the grid dependent boxes -allocate(gc%box%count(nx,ny,nz), gc%box%start(nx,ny,nz)) -gc%box%count = -1 -gc%box%start = -1 +!allocate(gc%box%count(nx,ny,nz), gc%box%start(nx,ny,nz)) +!gc%box%count = -1 +!gc%box%start = -1 + +do i=1, gc%nt + ! Allocate storage for locs number dependent part + allocate(gc%box(i)%loc_box(num)) + gc%box(i)%loc_box(:) = -1 +enddo ! Set the value of num_locs in the structure -gc%num = num +!gc%num = num ! If num == 0, no point in going any further. -if (num == 0) return +!if (num == 0) return !! FIXME: compute nx, ny, nz from nboxes? or put in namelist !nx = nint(real(nboxes, r8)**0.33333) ! roughly cube root @@ -860,59 +991,60 @@ subroutine get_close_init(gc, num, maxdist, locs, maxdist_list) !> and nboxes*2 is the end of the halo on the upper side. ! Determine where the boxes should be for this set of locs and maxdist -call find_box_ranges(gc, locs, num) - -! Begin by computing the number of locations in each box in x,y,z -gc%box%count = 0 -do i = 1, num - -!write(0,*) i, locs(i)%x, locs(i)%y, locs(i)%z - x_box(i) = floor((locs(i)%x - gc%box%bot_x) / gc%box%x_width) + 1 - if(x_box(i) > nx) x_box(i) = nx - if(x_box(i) < 1) x_box(i) = 1 - - y_box(i) = floor((locs(i)%y - gc%box%bot_y) / gc%box%y_width) + 1 - if(y_box(i) > ny) y_box(i) = ny - if(y_box(i) < 1) y_box(i) = 1 - - z_box(i) = floor((locs(i)%z - gc%box%bot_z) / gc%box%z_width) + 1 - if(z_box(i) > nz) z_box(i) = nz - if(z_box(i) < 1) z_box(i) = 1 - - gc%box%count(x_box(i), y_box(i), z_box(i)) = gc%box%count(x_box(i), y_box(i), z_box(i)) + 1 -!write(0,*) 'adding count to box ', x_box(i), y_box(i), z_box(i), & -! gc%box%count(x_box(i), y_box(i), z_box(i)) -end do - -! Figure out where storage for each boxes members should begin -cum_start = 1 -do i = 1, nx - do j = 1, ny - do k = 1, nz - gc%box%start(i, j, k) = cum_start - cum_start = cum_start + gc%box%count(i, j, k) +do n=1, gc%nt + ! Determine where the boxes should be for this set of locs and maxdist + call find_box_ranges(gc%box(n), locs, num) + + ! Begin by computing the number of locations in each box in x,y,z + gc%box(n)%count = 0 + do i = 1, num + + !write(0,*) i, locs(i)%x, locs(i)%y, locs(i)%z + x_box(i) = floor((locs(i)%x - gc%box(n)%bot_x) / gc%box(n)%x_width) + 1 + if(x_box(i) > nx) x_box(i) = nx + if(x_box(i) < 1) x_box(i) = 1 + + y_box(i) = floor((locs(i)%y - gc%box(n)%bot_y) / gc%box(n)%y_width) + 1 + if(y_box(i) > ny) y_box(i) = ny + if(y_box(i) < 1) y_box(i) = 1 + + z_box(i) = floor((locs(i)%z - gc%box(n)%bot_z) / gc%box(n)%z_width) + 1 + if(z_box(i) > nz) z_box(i) = nz + if(z_box(i) < 1) z_box(i) = 1 + + gc%box(n)%count(x_box(i), y_box(i), z_box(i)) = gc%box(n)%count(x_box(i), y_box(i), z_box(i)) + 1 + !write(0,*) 'adding count to box ', x_box(i), y_box(i), z_box(i), & + ! gc%box%count(x_box(i), y_box(i), z_box(i)) + end do + + ! Figure out where storage for each boxes members should begin + cum_start = 1 + do i = 1, nx + do j = 1, ny + do k = 1, nz + gc%box(n)%start(i, j, k) = cum_start + cum_start = cum_start + gc%box(n)%count(i, j, k) + end do end do end do -end do - -! Now we know how many are in each box, get a list of which are in each box -tstart = gc%box%start -do i = 1, num - gc%box%loc_box(tstart(x_box(i), y_box(i), z_box(i))) = i - tstart(x_box(i), y_box(i), z_box(i)) = tstart(x_box(i), y_box(i), z_box(i)) + 1 -end do -do i = 1, nx - do j = 1, ny - do k = 1, nz -!if (gc%box%count(i,j,k) > 0) write(0,*) i,j,k, gc%box%count(i,j,k), gc%box%start(i,j,k) - do l=1, gc%box%count(i,j,k) -!write(0,*) l, gc%box%loc_box(l) - enddo + ! Now we know how many are in each box, get a list of which are in each box + tstart = gc%box(n)%start + do i = 1, num + gc%box(n)%loc_box(tstart(x_box(i), y_box(i), z_box(i))) = i + tstart(x_box(i), y_box(i), z_box(i)) = tstart(x_box(i), y_box(i), z_box(i)) + 1 + end do + do i = 1, nx + do j = 1, ny + do k = 1, nz + !if (gc%box%count(i,j,k) > 0) write(0,*) i,j,k, gc%box%count(i,j,k), gc%box%start(i,j,k) + do l=1, gc%box(n)%count(i,j,k) + !write(0,*) l, gc%box%loc_box(l) + enddo + end do end do end do end do - ! info on how well the boxes are working. by default print nothing. ! set print_box_level to higher values to get more and more detail. ! user info should be level 1; 2 and 3 should be for debug only. @@ -922,12 +1054,12 @@ subroutine get_close_init(gc, num, maxdist, locs, maxdist_list) ! if print level > 2, set all tasks to print and call print. ! then reset the status to off again. if (do_output()) then - call print_get_close_type(gc, print_box_level) + call print_get_close_type(gc, 1, print_box_level) else if (print_box_level >= 2 .or. print_box_level < 0) then ! print status was false, but turn on temporarily ! to output box info from all tasks. call set_output(.true.) - call print_get_close_type(gc, print_box_level) + call print_get_close_type(gc, 1, print_box_level) call set_output(.false.) endif endif @@ -937,10 +1069,16 @@ end subroutine get_close_init !---------------------------------------------------------------------------- subroutine get_close_destroy(gc) - type(get_close_type), intent(inout) :: gc +integer :: i -deallocate(gc%box%loc_box, gc%box%count, gc%box%start) +do i = 1, gc%nt + if (allocated(gc%box(i)%loc_box)) deallocate(gc%box(i)%loc_box) + deallocate(gc%box(i)%count, gc%box(i)%start) +enddo +deallocate(gc%type_to_cutoff_map) +deallocate(gc%box) +!deallocate(gc%box%loc_box, gc%box%count, gc%box%start) ! ORIG end subroutine get_close_destroy @@ -948,7 +1086,6 @@ end subroutine get_close_destroy subroutine get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & num_close, close_ind, dist, ens_handle) - ! The specific type of the base observation, plus the generic kinds list ! for either the state or obs lists are available if a more sophisticated ! distance computation is needed. @@ -969,7 +1106,6 @@ end subroutine get_close_obs subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & num_close, close_ind, dist, ens_handle) - ! The specific type of the base observation, plus the generic kinds list ! for either the state or obs lists are available if a more sophisticated ! distance computation is needed. @@ -991,7 +1127,6 @@ end subroutine get_close_state subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & num_close, close_ind, dist, ens_handle) - type(get_close_type), intent(in) :: gc type(location_type), intent(in) :: base_loc, locs(:) integer, intent(in) :: base_type, loc_qtys(:) @@ -1002,7 +1137,7 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & ! If dist is NOT present, just find everybody in a box, put them in the list, ! but don't compute any distances -integer :: x_box, y_box, z_box, i, j, k, l +integer :: x_box, y_box, z_box, i, j, k, l, bt integer :: start_x, end_x, start_y, end_y, start_z, end_z integer :: n_in_box, st, t_ind real(r8) :: this_dist, this_maxdist @@ -1020,30 +1155,72 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & if(present(dist)) dist = -1e38_r8 ! big but negative this_dist = 1e38_r8 ! something big and positive. + +! handle identity obs correctly -- only support them if there are no +! per-obs-type cutoffs. identity obs don't have a specific type, they +! only have a generic kind based on what index into the state vector is +! specified. if this is absolutely needed by someone, i can rejigger the +! code to try to use the default cutoff for identity obs, but it's tricky +! to identify which obs type is using the default. in theory it's possible +! to specify a default cutoff and then specify a per-type cutoff for every +! other type in the system. in that case i don't have a map entry for the +! default, and it's a pain to construct one. so i'm punting for now. +if (base_type < 0) then + if (gc%nt > 1) then + write(msgstring, '(A)') 'no support for identity obs if per-obs-type cutoffs are specified' + write(msgstring1, '(A)') 'contact dart support if you have a need for this combination' + call error_handler(E_ERR, 'get_close', msgstring, source, text2=msgstring1) + endif + bt = 1 +else + ! map from type index to gtt index + if (base_type < 1 .or. base_type > size(gc%type_to_cutoff_map)) then + write(msgstring,'(A,I8)')'base_type out of range, is ', base_type + write(msgstring1,'(A,2I8)')'must be between ', 1, size(gc%type_to_cutoff_map) + call write_location (0, base_loc, charstring=msgstring2) + call error_handler(E_ERR, 'get_close', msgstring, source, & + text2=msgstring1, text3=msgstring2) + endif + bt = gc%type_to_cutoff_map(base_type) + if (bt < 1 .or. bt > gc%nt) then + write(msgstring,'(A,I8)')'mapped type index out of range, is ', bt + write(msgstring1,'(A,2I8)')'must be between ', 1, gc%nt + write(msgstring2, '(A)')'internal error, should not happen. Contact DART Support' + call error_handler(E_ERR, 'get_close', msgstring, source, & + text2=msgstring1, text3=msgstring2) + endif +endif + + + ! the list of locations in the loc() argument must be the same ! as the list of locations passed into get_close_init(), so ! gc%num and size(loc) better be the same. if the list changes, ! you have to destroy the old gc and init a new one. -if (size(locs) /= gc%num) then - write(errstring,*)'locs() array must match one passed to get_close_init()' - call error_handler(E_ERR, 'get_close_boxes', errstring, source) +if (size(locs) /= gc%box(bt)%num) then + !write(errstring,*)'locs() array must match one passed to get_close_init()' + !call error_handler(E_ERR, 'get_close_boxes', errstring, source) + write(msgstring,'(A)')'locs() array must match one passed to get_close_init()' + write(msgstring1,'(A,2I8)')'init size, current list size: ', gc%box(bt)%num, size(locs) + write(msgstring2,'(A,I8)')'bt = ', bt + call error_handler(E_ERR, 'get_close', msgstring, source, & + text2=msgstring1, text3=msgstring2) endif ! If num == 0, no point in going any further. -if (gc%num == 0) return +if (gc%box(bt)%num == 0) return -this_maxdist = gc%maxdist +this_maxdist = gc%box(bt)%maxdist !> @todo this is doing an exhaustive search each time. expensive !> but should give the right answer. - if(.true.) then if (present(dist)) then - call exhaustive_collect(gc, base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, locs, & num_close, close_ind, dist) else allocate(cdist(size(locs))) - call exhaustive_collect(gc, base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, locs, & num_close, close_ind, cdist) deallocate(cdist) endif @@ -1054,15 +1231,15 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & ! exhaustive search if(compare_to_correct) then allocate(cclose_ind(size(locs)), cdist(size(locs))) - call exhaustive_collect(gc, base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, locs, & cnum_close, cclose_ind, cdist) endif ! Begin by figuring out which box the base loc is in -x_box = floor((base_loc%x - gc%box%bot_x) / gc%box%x_width) + 1 -y_box = floor((base_loc%y - gc%box%bot_y) / gc%box%y_width) + 1 -z_box = floor((base_loc%z - gc%box%bot_z) / gc%box%z_width) + 1 +x_box = floor((base_loc%x - gc%box(bt)%bot_x) / gc%box(bt)%x_width) + 1 +y_box = floor((base_loc%y - gc%box(bt)%bot_y) / gc%box(bt)%y_width) + 1 +z_box = floor((base_loc%z - gc%box(bt)%bot_z) / gc%box(bt)%z_width) + 1 ! If it is not in any box, then it is more than the maxdist away from everybody if(x_box > nx .or. x_box < 1 .or. x_box < 0) return @@ -1075,19 +1252,19 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & !write(0,*) 'nboxes x, y, z = ', gc%box%nboxes_x, gc%box%nboxes_y, gc%box%nboxes_z !write(0,*) 'base_loc in box ', x_box, y_box, z_box -start_x = x_box - gc%box%nboxes_x +start_x = x_box - gc%box(bt)%nboxes_x if (start_x < 1) start_x = 1 -end_x = x_box + gc%box%nboxes_x +end_x = x_box + gc%box(bt)%nboxes_x if (end_x > nx) end_x = nx -start_y = y_box - gc%box%nboxes_y +start_y = y_box - gc%box(bt)%nboxes_y if (start_y < 1) start_y = 1 -end_y = y_box + gc%box%nboxes_y +end_y = y_box + gc%box(bt)%nboxes_y if (end_y > ny) end_y = ny -start_z = z_box - gc%box%nboxes_z +start_z = z_box - gc%box(bt)%nboxes_z if (start_z < 1) start_z = 1 -end_z = z_box + gc%box%nboxes_z +end_z = z_box + gc%box(bt)%nboxes_z if (end_z > nz) end_z = nz !write(0,*) 'looping from ' @@ -1101,14 +1278,14 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & do k = start_z, end_z ! Box to search is i,j,k - n_in_box = gc%box%count(i, j, k) - st = gc%box%start(i,j,k) + n_in_box = gc%box(bt)%count(i, j, k) + st = gc%box(bt)%start(i,j,k) ! Loop to check how close all loc in the box are; add those that are close do l = 1, n_in_box - t_ind = gc%box%loc_box(st - 1 + l) + t_ind = gc%box(bt)%loc_box(st - 1 + l) !write(0,*) 'l, t_ind = ', l, t_ind ! Only compute distance if dist is present @@ -1144,7 +1321,7 @@ end subroutine get_close !-------------------------------------------------------------------------- -subroutine find_box_ranges(gc, locs, num) +subroutine find_box_ranges(box, locs, num) ! Finds boundaries for x,y,z boxes. ! FIXME: ways boxes could be divided: @@ -1153,58 +1330,58 @@ subroutine find_box_ranges(gc, locs, num) ! on each side of the dividing plane. ! - about 100 other schemes -type(get_close_type), intent(inout) :: gc +type(box_type), intent(inout) :: box integer, intent(in) :: num type(location_type), intent(in) :: locs(num) !logical :: old_out if (x_is_periodic) then - gc%box%bot_x = min_x_for_periodic - gc%box%top_x = max_x_for_periodic + box%bot_x = min_x_for_periodic + box%top_x = max_x_for_periodic else - gc%box%bot_x = minval(locs(:)%x - gc%maxdist) - gc%box%top_x = maxval(locs(:)%x + gc%maxdist) + box%bot_x = minval(locs(:)%x - box%maxdist) + box%top_x = maxval(locs(:)%x + box%maxdist) endif if (y_is_periodic) then - gc%box%bot_y = min_x_for_periodic - gc%box%top_y = max_x_for_periodic + box%bot_y = min_x_for_periodic + box%top_y = max_x_for_periodic else - gc%box%bot_y = minval(locs(:)%y - gc%maxdist) - gc%box%top_y = maxval(locs(:)%y + gc%maxdist) + box%bot_y = minval(locs(:)%y - box%maxdist) + box%top_y = maxval(locs(:)%y + box%maxdist) endif if (z_is_periodic) then - gc%box%bot_z = min_x_for_periodic - gc%box%top_z = max_x_for_periodic + box%bot_z = min_x_for_periodic + box%top_z = max_x_for_periodic else - gc%box%bot_z = minval(locs(:)%z - gc%maxdist) - gc%box%top_z = maxval(locs(:)%z + gc%maxdist) + box%bot_z = minval(locs(:)%z - box%maxdist) + box%top_z = maxval(locs(:)%z + box%maxdist) endif if (debug > 0) print *, 'nx/ny/nz: ', nx, ny, nz -if (debug > 0) print *, 'bots: ', gc%box%bot_x, gc%box%bot_y, gc%box%bot_z -if (debug > 0) print *, 'tops: ', gc%box%top_x, gc%box%top_y, gc%box%top_z +if (debug > 0) print *, 'bots: ', box%bot_x, box%bot_y, box%bot_z +if (debug > 0) print *, 'tops: ', box%top_x, box%top_y, box%top_z -gc%box%x_width = max(1.0_r8, (gc%box%top_x - gc%box%bot_x) / nx) -gc%box%y_width = max(1.0_r8, (gc%box%top_y - gc%box%bot_y) / ny) -gc%box%z_width = max(1.0_r8, (gc%box%top_z - gc%box%bot_z) / nz) +box%x_width = max(1.0_r8, (box%top_x - box%bot_x) / nx) +box%y_width = max(1.0_r8, (box%top_y - box%bot_y) / ny) +box%z_width = max(1.0_r8, (box%top_z - box%bot_z) / nz) -if (debug > 0) print *, 'widths = ', gc%box%x_width, gc%box%y_width, gc%box%z_width +if (debug > 0) print *, 'widths = ', box%x_width, box%y_width, box%z_width ! FIXME: compute a sphere of radius maxdist and see how ! many boxes in x, y, z that would include. -if (gc%box%x_width <= 0.0_r8) & +if (box%x_width <= 0.0_r8) & call error_handler(E_ERR, 'find_box_ranges', 'x_width <= 0', source) -if (gc%box%y_width <= 0.0_r8) & +if (box%y_width <= 0.0_r8) & call error_handler(E_ERR, 'find_boy_ranges', 'y_width <= 0', source) -if (gc%box%z_width <= 0.0_r8) & +if (box%z_width <= 0.0_r8) & call error_handler(E_ERR, 'find_boz_ranges', 'z_width <= 0', source) -gc%box%nboxes_x = aint((gc%maxdist + (gc%box%x_width-1)) / gc%box%x_width) -gc%box%nboxes_y = aint((gc%maxdist + (gc%box%y_width-1)) / gc%box%y_width) -gc%box%nboxes_z = aint((gc%maxdist + (gc%box%z_width-1)) / gc%box%z_width) +box%nboxes_x = aint((box%maxdist + (box%x_width-1)) / box%x_width) +box%nboxes_y = aint((box%maxdist + (box%y_width-1)) / box%y_width) +box%nboxes_z = aint((box%maxdist + (box%z_width-1)) / box%z_width) !if(compare_to_correct) then ! old_out = do_output() @@ -1222,8 +1399,10 @@ end subroutine find_box_ranges !---------------------------------------------------------------------------- -subroutine find_nearest(gc, base_loc, loc_list, nearest, rc) - type(get_close_type), intent(in), target :: gc +!subroutine find_nearest(box, base_loc, loc_list, nearest, rc) + !type(box_type), intent(in), target :: box +subroutine find_nearest(gc,base_loc,loc_list,nearest, rc) + type(get_close_type), intent(in), target :: gc type(location_type), intent(in) :: base_loc type(location_type), intent(in) :: loc_list(:) integer, intent(out) :: nearest @@ -1233,32 +1412,33 @@ subroutine find_nearest(gc, base_loc, loc_list, nearest, rc) integer :: x_box, y_box, z_box, i, j, k, l integer :: start_x, end_x, start_y, end_y, start_z, end_z -integer :: n_in_box, st, t_ind +integer :: n_in_box, st, t_ind +integer :: box_num real(r8) :: this_dist, dist ! First, set the intent out arguments to a missing value nearest = -99 rc = -1 dist = 1e38_r8 ! something big and positive. - +box_num = 1 ! JDL - HARD CODE TO JUST GET THE FIRST BOX OBJECT ! the list of locations in the loc() argument must be the same ! as the list of locations passed into get_close_init(), so ! gc%num and size(loc) better be the same. if the list changes, ! you have to destroy the old gc and init a new one. -if (size(loc_list) /= gc%num) then +if (size(loc_list) /= gc%box(box_num)%num) then write(errstring,*)'loc() array must match one passed to get_close_init()' call error_handler(E_ERR, 'find_nearest_boxes', errstring, source) endif ! If num == 0, no point in going any further. -if (gc%num == 0) return +if (gc%box(box_num)%num == 0) return !-------------------------------------------------------------- ! Begin by figuring out which box the base loc is in -x_box = floor((base_loc%x - gc%box%bot_x) / gc%box%x_width) + 1 -y_box = floor((base_loc%y - gc%box%bot_y) / gc%box%y_width) + 1 -z_box = floor((base_loc%z - gc%box%bot_z) / gc%box%z_width) + 1 +x_box = floor((base_loc%x - gc%box(box_num)%bot_x) / gc%box(box_num)%x_width) + 1 +y_box = floor((base_loc%y - gc%box(box_num)%bot_y) / gc%box(box_num)%y_width) + 1 +z_box = floor((base_loc%z - gc%box(box_num)%bot_z) / gc%box(box_num)%z_width) + 1 ! FIXME: this should figure out if it's > n or < 0 and ! set to n or 0 and always return something. @@ -1273,19 +1453,19 @@ subroutine find_nearest(gc, base_loc, loc_list, nearest, rc) !write(0,*) 'nboxes x, y, z = ', gc%box%nboxes_x, gc%box%nboxes_y, gc%box%nboxes_z !write(0,*) 'base_loc in box ', x_box, y_box, z_box -start_x = x_box - gc%box%nboxes_x +start_x = x_box - gc%box(box_num)%nboxes_x if (start_x < 1) start_x = 1 -end_x = x_box + gc%box%nboxes_x +end_x = x_box + gc%box(box_num)%nboxes_x if (end_x > nx) end_x = nx -start_y = y_box - gc%box%nboxes_y +start_y = y_box - gc%box(box_num)%nboxes_y if (start_y < 1) start_y = 1 -end_y = y_box + gc%box%nboxes_y +end_y = y_box + gc%box(box_num)%nboxes_y if (end_y > ny) end_y = ny -start_z = z_box - gc%box%nboxes_z +start_z = z_box - gc%box(box_num)%nboxes_z if (start_z < 1) start_z = 1 -end_z = z_box + gc%box%nboxes_z +end_z = z_box + gc%box(box_num)%nboxes_z if (end_z > nz) end_z = nz !write(0,*) 'looping from ' @@ -1299,14 +1479,14 @@ subroutine find_nearest(gc, base_loc, loc_list, nearest, rc) do k = start_z, end_z ! Box to search is i,j,k - n_in_box = gc%box%count(i, j, k) - st = gc%box%start(i,j,k) + n_in_box = gc%box(box_num)%count(i, j, k) + st = gc%box(box_num)%start(i,j,k) ! Loop to check how close all loc in the box are; add those that are close do l = 1, n_in_box - t_ind = gc%box%loc_box(st - 1 + l) + t_ind = gc%box(box_num)%loc_box(st - 1 + l) !write(0,*) 'l, t_ind = ', l, t_ind this_dist = get_dist(base_loc, loc_list(t_ind)) @@ -1420,7 +1600,7 @@ subroutine recompute_periodic() loopy%midline(IY) = (min_y_for_periodic + max_y_for_periodic ) / 2.0_r8 loopy%midline(IZ) = (min_z_for_periodic + max_z_for_periodic ) / 2.0_r8 -if (debug > 0) write(0,*) 'setting midlines: ', loopy%midline + if (debug > 0) write(0,*) 'setting midlines: ', loopy%midline ! update these in the get_close init routines loopy%midbox(IX) = -1 loopy%midbox(IY) = -1 @@ -1429,7 +1609,7 @@ subroutine recompute_periodic() loopy%offset(IX) = max_x_for_periodic - min_x_for_periodic loopy%offset(IY) = max_y_for_periodic - min_y_for_periodic loopy%offset(IZ) = max_z_for_periodic - min_z_for_periodic -if (debug > 0) write(0,*) 'setting offsets: ', loopy%offset + if (debug > 0) write(0,*) 'setting offsets: ', loopy%offset else loopy%midline(:) = MISSING_R8 loopy%midbox(:) = MISSING_I @@ -1441,6 +1621,61 @@ subroutine recompute_periodic() end subroutine recompute_periodic + +!--------------------------------------------------------------------------- +subroutine distinct_values(in_list, mycount, values, map) +!> parse an input list of values and return: +!> 1) the count of distinct values +!> 2) the list of unique values +!> 3) the mapping of the input list to the value list +!> the values and map list should already be allocated, and be the same +!> length as the incoming list length. + +real(r8), intent(in) :: in_list(:) !< incoming list of all values +integer, intent(out) :: mycount !< count of distinct values +real(r8), intent(inout) :: values(:) !< list of distinct values +integer, intent(inout) :: map(:) !< mapping of in_list to values + +integer :: i, j, listsize, nextslot +logical :: foundnew +real(r8) :: newval + +! set return values now; if we error out then we can +! just return. +mycount = 0 +values(:) = -1.0_r8 +map(:) = -1 + +listsize = size(in_list) +if (listsize <= 0) return +if (size(values) /= size(in_list)) return +if (size(map) /= size(in_list)) return + +! set values() with only the unique distances. +! when done, the valid values are only 'count' long, +! not 'listsize'. +OUTER: do i=1, listsize + newval = in_list(i) + foundnew = .true. + INNER: do j=1, listsize + nextslot = j + if (values(j) < 0.0_r8) exit INNER + if (abs(values(j) - newval) <= epsilon(newval)) then + foundnew = .false. + map(i) = j + exit INNER + endif + enddo INNER + if (foundnew) then + values(nextslot) = newval + map(i) = nextslot + mycount = nextslot + endif +enddo OUTER + +end subroutine distinct_values + + !--------------------------------------------------------------------------- function get_maxdist(gc, obs_type) @@ -1448,23 +1683,28 @@ function get_maxdist(gc, obs_type) integer, optional, intent(in) :: obs_type real(r8) :: get_maxdist -get_maxdist = gc%maxdist +integer :: bt + +bt = gc%type_to_cutoff_map(obs_type) +get_maxdist = gc%box(bt)%maxdist + +!get_maxdist = gc%maxdist end function get_maxdist !---------------------------------------------------------------------------- -subroutine print_get_close_type(gc, amount) +subroutine print_get_close_type(gc, tt, amount) ! print out debugging statistics, or optionally print out a full ! dump from all mpi tasks in a format that can be plotted with matlab. +type(get_close_type), intent(in) :: gc +integer, intent(in), optional :: tt +integer, intent(in), optional :: amount -type(get_close_type), intent(in), target :: gc -integer, intent(in), optional :: amount - -integer :: i, j, k, l, first, index, mytask, alltasks +integer :: i, j, k, l, first, index, mytask, alltasks, whichtt integer :: sample, nfull, nempty, howmuch, total, maxcount, maxi, maxj, maxk -logical :: tickmark(gc%num), iam0 +logical :: tickmark(gc%box(1)%num), iam0 real(r8) :: x_cen, y_cen, z_cen logical, save :: write_now = .true. @@ -1475,6 +1715,14 @@ subroutine print_get_close_type(gc, amount) ! cumulative times through this routine been_called = been_called + 1 +! second arg is optional, defaults to 1, and selects which +! of the cutoff structs to print +if (present(tt)) then + whichtt = tt +else + whichtt = 1 +endif + ! second arg is now an int, not logical, and means: ! 0 = very terse, only box summary (default). ! 1 = structs and first part of arrays. @@ -1506,7 +1754,7 @@ subroutine print_get_close_type(gc, amount) ! locations from the state vector in one set of boxes, but just a few ! locations from the observations in another. this lets you turn off ! the debugging level for the large set and leave it on for the small. -!if (gc%num > 100) howmuch = 0 +!if (gc%box(whichtt)%num > 100) howmuch = 0 ! print the get_close_type derived type values @@ -1514,19 +1762,19 @@ subroutine print_get_close_type(gc, amount) write(errstring,*) 'get_close_type values:' call error_handler(E_MSG, 'loc', errstring) - write(errstring,*) ' num = ', gc%num + write(errstring,*) ' num = ', gc%box(whichtt)%num call error_handler(E_MSG, 'loc', errstring) write(errstring,*) ' nx, ny, nz = ', nx, ny, nz call error_handler(E_MSG, 'loc', errstring) - write(errstring,"(A,F12.6)") ' maxdist = ', gc%maxdist + write(errstring,"(A,F12.6)") ' maxdist = ', gc%box(whichtt)%maxdist call error_handler(E_MSG, 'loc', errstring) - write(errstring, "(A,3(F12.6))") ' x_box: bot, top, width = ', gc%box%bot_x, gc%box%top_x, gc%box%x_width + write(errstring, "(A,3(F12.6))") ' x_box: bot, top, width = ', gc%box(whichtt)%bot_x, gc%box(whichtt)%top_x, gc%box(whichtt)%x_width call error_handler(E_MSG, 'loc', errstring) - write(errstring, "(A,3(F12.6))") ' y_box: bot, top, width = ', gc%box%bot_y, gc%box%top_y, gc%box%y_width + write(errstring, "(A,3(F12.6))") ' y_box: bot, top, width = ', gc%box(whichtt)%bot_y, gc%box(whichtt)%top_y, gc%box(whichtt)%y_width call error_handler(E_MSG, 'loc', errstring) - write(errstring, "(A,3(F12.6))") ' z_box: bot, top, width = ', gc%box%bot_z, gc%box%top_z, gc%box%z_width + write(errstring, "(A,3(F12.6))") ' z_box: bot, top, width = ', gc%box(whichtt)%bot_z, gc%box(whichtt)%top_z, gc%box(whichtt)%z_width call error_handler(E_MSG, 'loc', errstring) endif @@ -1534,19 +1782,19 @@ subroutine print_get_close_type(gc, amount) ! this one can be very large. print only the first nth unless ! instructed otherwise. (print n+1 because 1 more value fits on ! the line because it prints ( i ) and not ( i, j ) like the others.) -if (associated(gc%box%loc_box)) then - i = size(gc%box%loc_box,1) - if (i/= gc%num) then - write(errstring,*) ' warning: size of loc_box incorrect, nlocs, i =', gc%num, i +if (allocated(gc%box(whichtt)%loc_box)) then + i = size(gc%box(whichtt)%loc_box,1) + if (i/= gc%box(whichtt)%num) then + write(errstring,*) ' warning: size of loc_box incorrect, nlocs, i =', gc%box(whichtt)%num, i call error_handler(E_MSG, 'locations_mod', errstring) endif if (howmuch > 1) then ! DEBUG - write(errstring,"(A,I8,A,36(I8,1X))") ' loc_box(',i,') =', gc%box%loc_box(1:min(i,36)) ! (nlocs) + write(errstring,"(A,I8,A,36(I8,1X))") ' loc_box(',i,') =', gc%box(whichtt)%loc_box(1:min(i,36)) ! (nlocs) !write(errstring,*) ' loc_box(',i,') =', gc%box%loc_box ! (nlocs) call error_handler(E_MSG, 'locations_mod', errstring) else if(howmuch > 0) then - write(errstring,*) ' loc_box(',i,') =', gc%box%loc_box(1:min(i,sample+1)) + write(errstring,*) ' loc_box(',i,') =', gc%box(whichtt)%loc_box(1:min(i,sample+1)) call error_handler(E_MSG, 'locations_mod', errstring) write(errstring,*) ' ' call error_handler(E_MSG, 'locations_mod', errstring) @@ -1560,10 +1808,10 @@ subroutine print_get_close_type(gc, amount) ! like loc_box, this one can be very large. print only the first nth unless ! instructed otherwise -if (associated(gc%box%start)) then - i = size(gc%box%start,1) - j = size(gc%box%start,2) - k = size(gc%box%start,3) +if (allocated(gc%box(whichtt)%start)) then + i = size(gc%box(whichtt)%start,1) + j = size(gc%box(whichtt)%start,2) + k = size(gc%box(whichtt)%start,3) if ((i /= nx) .or. (j /= ny) .or. (k /= nz)) then write(errstring,*) ' warning: size of start incorrect, nx, ny, nz, i, j, k =', nx, ny, nz, i, j, k call error_handler(E_MSG, 'locations_mod', errstring) @@ -1572,11 +1820,11 @@ subroutine print_get_close_type(gc, amount) write(errstring,*) ' start(',i,j,k,') =' ! (nx, ny, nz) call error_handler(E_MSG, 'locations_mod', errstring) do l=1, j - write(errstring,"(36(I8,1X))") gc%box%start(1:min(i,36), l, 1) + write(errstring,"(36(I8,1X))") gc%box(whichtt)%start(1:min(i,36), l, 1) call error_handler(E_MSG, 'locations_mod', errstring) enddo else if (howmuch > 0) then - write(errstring,*) ' start(',i,j,k,') =', gc%box%start(1:min(i,sample), 1, 1) + write(errstring,*) ' start(',i,j,k,') =', gc%box(whichtt)%start(1:min(i,sample), 1, 1) call error_handler(E_MSG, 'locations_mod', errstring) write(errstring,*) ' ' call error_handler(E_MSG, 'locations_mod', errstring) @@ -1589,10 +1837,10 @@ subroutine print_get_close_type(gc, amount) endif ! as above, print only first n unless second arg is .true. -if (associated(gc%box%count)) then - i = size(gc%box%count,1) - j = size(gc%box%count,2) - k = size(gc%box%count,3) +if (allocated(gc%box(whichtt)%count)) then + i = size(gc%box(whichtt)%count,1) + j = size(gc%box(whichtt)%count,2) + k = size(gc%box(whichtt)%count,3) if ((i /= nx) .or. (j /= ny) .or. (k /= nz)) then write(errstring,*) ' warning: size of count incorrect, nx, ny, nz, i, j, k =', nx, ny, nz, i, j, k call error_handler(E_MSG, 'locations_mod', errstring) @@ -1601,11 +1849,11 @@ subroutine print_get_close_type(gc, amount) write(errstring,*) ' count(',i,j,k,') =' ! (nx, ny, nz) call error_handler(E_MSG, 'locations_mod', errstring) do l=1, j - write(errstring,"(36(I8,1X))") gc%box%count(1:min(i,36), l, 1) + write(errstring,"(36(I8,1X))") gc%box(whichtt)%count(1:min(i,36), l, 1) call error_handler(E_MSG, 'locations_mod', errstring) enddo else if (howmuch > 0) then - write(errstring,*) ' count(',i,j,k,') =', gc%box%count(1:min(i,sample), 1, 1) + write(errstring,*) ' count(',i,j,k,') =', gc%box(whichtt)%count(1:min(i,sample), 1, 1) call error_handler(E_MSG, 'locations_mod', errstring) write(errstring,*) ' ' call error_handler(E_MSG, 'locations_mod', errstring) @@ -1629,10 +1877,10 @@ subroutine print_get_close_type(gc, amount) do i=1, nx do j=1, ny do k=1, nz - first = gc%box%start(i, j, k) - do l=1, gc%box%count(i, j, k) + first = gc%box(whichtt)%start(i, j, k) + do l=1, gc%box(whichtt)%count(i, j, k) index = first + l - 1 - if ((index < 1) .or. (index > gc%num)) then + if ((index < 1) .or. (index > gc%box(whichtt)%num)) then write(errstring, *) 'exiting at first bad value; could be more' call error_handler(E_MSG, 'locations_mod', errstring) write(errstring, *) 'bad locs list index, in box: ', index, i, j, k @@ -1651,7 +1899,7 @@ subroutine print_get_close_type(gc, amount) enddo enddo -do i=1, gc%num +do i=1, gc%box(whichtt)%num if (.not. tickmark(i)) then write(errstring, *) 'exiting at first bad value; could be more' call error_handler(E_MSG, 'locations_mod', errstring) @@ -1684,24 +1932,24 @@ subroutine print_get_close_type(gc, amount) do i=1, nx if (howmuch == -8) then - x_cen = gc%box%bot_x + ((i-1)*gc%box%x_width) + (gc%box%x_width/2.0) + x_cen = gc%box(whichtt)%bot_x + ((i-1)*gc%box(whichtt)%x_width) + (gc%box(whichtt)%x_width/2.0) write(funit, '(A,I2,A,I4,A,F12.9,A)') 'xlocs(', i, ',', mytask+1, ') = ', x_cen, ';' endif do j=1, ny if (howmuch == -8 .and. i==1) then - y_cen = gc%box%bot_y + ((j-1)*gc%box%y_width) + (gc%box%y_width/2.0) + y_cen = gc%box(whichtt)%bot_y + ((j-1)*gc%box(whichtt)%y_width) + (gc%box(whichtt)%y_width/2.0) write(funit, '(A,I2,A,I4,A,F12.9,A)') 'ylocs(', j, ',', mytask+1, ') = ', y_cen, ';' endif do k=1, nz if (howmuch == -8 .and. i==1) then - z_cen = gc%box%bot_z + ((j-1)*gc%box%z_width) + (gc%box%z_width/2.0) + z_cen = gc%box(whichtt)%bot_z + ((j-1)*gc%box(whichtt)%z_width) + (gc%box(whichtt)%z_width/2.0) write(funit, '(A,I2,A,I4,A,F12.9,A)') 'zlocs(', k, ',', mytask+1, ') = ', z_cen, ';' endif - if (gc%box%count(i, j, k) > 0) then + if (gc%box(whichtt)%count(i, j, k) > 0) then nfull = nfull + 1 - total = total + gc%box%count(i, j, k) - if (gc%box%count(i, j, k) > maxcount) then - maxcount = gc%box%count(i, j, k) + total = total + gc%box(whichtt)%count(i, j, k) + if (gc%box(whichtt)%count(i, j, k) > maxcount) then + maxcount = gc%box(whichtt)%count(i, j, k) maxi = i maxj = j maxk = k @@ -1712,7 +1960,7 @@ subroutine print_get_close_type(gc, amount) ! output for grid boxes; in matlab-friendly format if (howmuch == -8) then write(funit, '(3(A,I2),A,I4,A,I8,A)') 'boxes(', i, ', ', j, ', ', k, & - ',', mytask+1, ') = ', gc%box%count(i, j, k), ';' + ',', mytask+1, ') = ', gc%box(whichtt)%count(i, j, k), ';' endif enddo enddo @@ -1728,7 +1976,7 @@ subroutine print_get_close_type(gc, amount) call error_handler(E_MSG, 'locations_mod', errstring) write(errstring, '(a,i9)') " Total boxes (nx * ny * nz): ", nfull + nempty call error_handler(E_MSG, 'locations_mod', errstring) -write(errstring, '(a,i9)') " Total items to put in boxes: ", gc%num +write(errstring, '(a,i9)') " Total items to put in boxes: ", gc%box(whichtt)%num call error_handler(E_MSG, 'locations_mod', errstring) if (howmuch > 0) then write(errstring, '(a,i9)') " Total boxes with 1+ items: ", nfull @@ -1783,12 +2031,12 @@ end function is_location_in_region !--------------------------------------------------------------------------- -subroutine exhaustive_collect(gc, base_loc, loc_list, num_close, close_ind, close_dist) +subroutine exhaustive_collect(box, base_loc, loc_list, num_close, close_ind, close_dist) ! For validation, it is useful to be able to compare against exact ! exhaustive search -type(get_close_type), intent(in) :: gc +type(box_type), intent(in) :: box type(location_type), intent(in) :: base_loc, loc_list(:) integer, intent(out) :: num_close integer, intent(out) :: close_ind(:) @@ -1798,9 +2046,9 @@ subroutine exhaustive_collect(gc, base_loc, loc_list, num_close, close_ind, clos integer :: i num_close = 0 -do i = 1, gc%num +do i = 1, box%num this_dist = get_dist(base_loc, loc_list(i)) - if(this_dist <= gc%maxdist) then + if(this_dist <= box%maxdist) then ! Add this loc to correct list num_close = num_close + 1 close_ind(num_close) = i diff --git a/conf.py b/conf.py index 367d25c394..337c1669a5 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '9.13.0' +release = '9.13.1' master_doc = 'README' # -- General configuration --------------------------------------------------- diff --git a/developer_tests/location/run_tests.csh b/developer_tests/location/run_tests.csh index 06b46fca9f..df0225b7b2 100755 --- a/developer_tests/location/run_tests.csh +++ b/developer_tests/location/run_tests.csh @@ -45,9 +45,7 @@ foreach i ( $LOCLIST ) cd $i/test - # The threed_sphere location_mod actually needs an obs_kind_mod.f90 - # Consequently, we need to run preprocess to generate the file. - # It is the only location module that needs it, AFAIK + # Need to run preprocess for location mods that use obs_kind_mod switch ( $i ) case threed_sphere \rm -rf Makefile @@ -55,6 +53,12 @@ foreach i ( $LOCLIST ) make >> $LOGDIR/buildlog.$i.preprocess.out ./preprocess > $LOGDIR/runlog.$i.preprocess.out breaksw + case threed_cartesian + \rm -rf Makefile + ./mkmf_preprocess > $LOGDIR/buildlog.$i.preprocess.out + make >> $LOGDIR/buildlog.$i.preprocess.out + ./preprocess > $LOGDIR/runlog.$i.preprocess.out + breaksw default breaksw endsw diff --git a/developer_tests/location/threed_cartesian/test/input.nml b/developer_tests/location/threed_cartesian/test/input.nml index 7ea001b4aa..6f2ead6f9c 100644 --- a/developer_tests/location/threed_cartesian/test/input.nml +++ b/developer_tests/location/threed_cartesian/test/input.nml @@ -1,3 +1,15 @@ + +# pick a random set of inputs +&preprocess_nml + overwrite_output = .true. + 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' + 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' + obs_type_files = '../../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90' + quantity_files = '../../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' + / + &location_nml nx = 10 ny = 10 diff --git a/developer_tests/location/threed_cartesian/test/path_names_location_test b/developer_tests/location/threed_cartesian/test/path_names_location_test index 02299a810f..e88f721c3d 100644 --- a/developer_tests/location/threed_cartesian/test/path_names_location_test +++ b/developer_tests/location/threed_cartesian/test/path_names_location_test @@ -1,5 +1,6 @@ assimilation_code/location/threed_cartesian/location_mod.f90 assimilation_code/location/utilities/default_location_mod.f90 +assimilation_code/modules/observations/obs_kind_mod.f90 assimilation_code/modules/utilities/ensemble_manager_mod.f90 assimilation_code/modules/utilities/netcdf_utilities_mod.f90 assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 diff --git a/developer_tests/location/threed_cartesian/test/path_names_location_test3 b/developer_tests/location/threed_cartesian/test/path_names_location_test3 index 175689c045..88f93851a3 100644 --- a/developer_tests/location/threed_cartesian/test/path_names_location_test3 +++ b/developer_tests/location/threed_cartesian/test/path_names_location_test3 @@ -1,6 +1,7 @@ assimilation_code/location/threed_cartesian/location_mod.f90 assimilation_code/location/utilities/default_location_mod.f90 assimilation_code/location/utilities/location_io_mod.f90 +assimilation_code/modules/observations/obs_kind_mod.f90 assimilation_code/modules/utilities/ensemble_manager_mod.f90 assimilation_code/modules/utilities/netcdf_utilities_mod.f90 assimilation_code/modules/utilities/null_mpi_utilities_mod.f90 diff --git a/models/cm1/model_mod.f90 b/models/cm1/model_mod.f90 index 93f378cce5..8d407350bc 100644 --- a/models/cm1/model_mod.f90 +++ b/models/cm1/model_mod.f90 @@ -148,7 +148,6 @@ module model_mod integer :: kind_list(MAX_STATE_VARIABLES) = MISSING_I real(r8) :: clamp_vals(MAX_STATE_VARIABLES,2) = MISSING_R8 - ! indices associated with variable_table columns integer, parameter :: VT_VARNAME_INDEX = 1 integer, parameter :: VT_DARTQTY_INDEX = 2 @@ -385,7 +384,6 @@ end function get_model_size !> given an index into the state vector, return its location and !> if given, the var kind. despite the name, var_type is a generic !> kind, like those in obs_kind/obs_kind_mod.f90, starting with QTY_ - subroutine get_state_meta_data(index_in, location, var_type) integer(i8), intent(in) :: index_in @@ -401,15 +399,16 @@ subroutine get_state_meta_data(index_in, location, var_type) if ( .not. module_initialized ) call static_init_model() ! get the local indicies and type from dart index. + call get_model_variable_indices(index_in, x_index, y_index, z_index, var_id=var_id) xlocation = 0.0_r8 ylocation = 0.0_r8 zlocation = 0.0_r8 -if (on_vgrid(var_id)) then +if (on_vgrid(var_id)) then - xlocation = xh(x_index) + xlocation = xh(x_index) ylocation = yf(y_index) if (y_index == 1 .or. y_index > nj ) then ! off height grid @@ -434,7 +433,7 @@ subroutine get_state_meta_data(index_in, location, var_type) endif -elseif (on_ugrid(var_id)) then +elseif (on_ugrid(var_id)) then xlocation = xf(x_index) ylocation = yh(y_index) @@ -461,7 +460,7 @@ subroutine get_state_meta_data(index_in, location, var_type) endif -else ! on sgrid +else ! on sgrid xlocation = xh(x_index) ylocation = yh(y_index) @@ -692,7 +691,9 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp integer :: hstatus ! for non-sgrid height interpolation ! need to know -integer :: nlevs ! number of levels (zhalf and zfull have different number?) +integer :: nlevs,nlevs_shrink ! number of levels (zhalf and zfull have different number?) +! nlevs_shrink is necesary when working with three-d data because only lowest two gridpoints +! are saved otherwise integer :: varid integer :: ndims real(r8), allocatable :: Q11_ens(:, :), Q12_ens(:, :), Q21_ens(:, :), Q22_ens(:, :) @@ -728,7 +729,6 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp istatus(:) = 11 return ! this kind isn't found in the state vector endif - nlevs = get_z_axis_length(varid) ndims = get_num_dims(domid, varid) if (debug > 99) then @@ -737,24 +737,23 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp endif ! 2d vs. 3d variable test if (ndims == 3) then - nlevs = 2 + nlevs_shrink = 2 else if (ndims == 2) then ! 2D, surface obs - nlevs = 1 + nlevs_shrink = 1 else ! should this be an error? write(string1, *) 'ndims not 3 or 2, unexpected? is ', ndims call say(string1) - nlevs = 1 !? just a guess + nlevs_shrink = 1 !? just a guess call write_location(0,location,charstring=string1) write(string1, *) 'task ', my_task_id(), ' kind ', obs_kind, ' (', & trim(get_name_for_quantity(obs_kind)), ') loc: ', trim(string1) call say(string1) - write(string1, *) 'varid, nlevs, ndims = ', varid, nlevs, ndims + write(string1, *) 'varid, nlevs, nlevs_shrink, ndims = ', varid, nlevs, nlevs_shrink, ndims call say(string1) endif - if( .not. observation_on_grid(obs_loc_array, ndims) ) then ! no need to interpolate istatus(:) = 12 @@ -764,8 +763,9 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp if (debug > 99) then write(string1, *) 'nlevs', nlevs call say(string1) + !write(string1, *) 'nlevs_shrink', nlevs_shrink + !call say(string1) endif - ! Interpolate the height field (z) to get the height at each level at ! the observation location. This allows us to find which level an observation ! is in. @@ -774,41 +774,35 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp ! Need grid from kind call get_x_axis(varid, axis, axis_length) call get_enclosing_coord(obs_loc_array(1), axis(1:axis_length), x_ind, x_val) - call get_y_axis(varid, axis, axis_length) call get_enclosing_coord(obs_loc_array(2), axis(1:axis_length), y_ind, y_val) - ! wrap the indicies if the observation is near the boundary if (periodic_x) call wrap_x( obs_loc_array(1), x_ind, x_val ) if (periodic_y) call wrap_y( obs_loc_array(2), y_ind, y_val ) - -if (nlevs == 2) then ! you need to find which 2 levels you are between +if (nlevs_shrink == 2) then ! you need to find which 2 levels you are between ! If variable is on ni, nj grid: if (is_on_s_grid(varid)) then - call height_interpolate_s_grid(obs_loc_array, varid, nlevs, z_ind, z_val) + call height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, z_ind, z_val) else - call height_interpolate(obs_loc_array, varid, nlevs, x_ind, y_ind, z_ind, x_val, y_val, z_val, hstatus) + call height_interpolate(obs_loc_array, varid, nlevs, nlevs_shrink, x_ind, y_ind, z_ind, x_val, y_val, z_val, hstatus) if (hstatus /= 0) return endif - if (debug > 99) print*, 'x y enclosing', x_ind, y_ind, x_val, y_val - if (debug > 99) print*, 'nlevs: ', nlevs, ', num_dims:', get_num_dims(domid, varid) - if (debug > 99) print*, 'z_ind', z_ind, 'varid', varid endif ! top and bottom, or just one value if 2d variable -allocate(Q11_ens(ens_size, nlevs)) -allocate(Q12_ens(ens_size, nlevs)) -allocate(Q21_ens(ens_size, nlevs)) -allocate(Q22_ens(ens_size, nlevs)) +allocate(Q11_ens(ens_size, nlevs_shrink)) +allocate(Q12_ens(ens_size, nlevs_shrink)) +allocate(Q21_ens(ens_size, nlevs_shrink)) +allocate(Q22_ens(ens_size, nlevs_shrink)) ! interpolated value -allocate(P_ens(ens_size, nlevs)) +allocate(P_ens(ens_size, nlevs_shrink)) -do i = 1, nlevs +do i = 1, nlevs_shrink indx = get_dart_vector_index(x_ind(1), y_ind(1), z_ind(i), domid, varid) Q11_ens(:, i) = get_state(indx, state_ens_handle) indx = get_dart_vector_index(x_ind(1), y_ind(2), z_ind(i), domid, varid) @@ -823,7 +817,7 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp ! P_ens is the interpolated value at the level below and above the obs for each ensemble member. ! P_ens is (ens_size, nlevs) -P_ens(:, :) = bilinear_interpolation_ens(ens_size, nlevs, obs_loc_array(1), obs_loc_array(2), & +P_ens(:, :) = bilinear_interpolation_ens(ens_size, nlevs_shrink, obs_loc_array(1), obs_loc_array(2), & x_val(1), x_val(2),y_val(1) , y_val(2), Q11_ens, Q12_ens, Q21_ens, Q22_ens) deallocate(Q11_ens, Q12_ens, Q21_ens, Q22_ens) @@ -831,7 +825,7 @@ subroutine model_interpolate(state_ens_handle, ens_size, location, obs_type, exp ! Interpolate between P to get the expected value ! If 2d don't call this. -if (nlevs == 2) then +if (nlevs_shrink == 2) then expected_obs = linear_interpolation(ens_size, obs_loc_array(3), z_val(1), z_val(2), P_ens(:, 1), P_ens(:, 2) ) else expected_obs = P_ens(:,1) @@ -848,11 +842,11 @@ end subroutine model_interpolate !> Interpolate height from corners of the bounding box locations to the !> observation location. -subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, z_ind, z_val) +subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, z_ind, z_val) real(r8), intent(in) :: obs_loc_array(3) integer, intent(in) :: varid -integer, intent(in) :: nlevs +integer, intent(in) :: nlevs, nlevs_shrink ! JDL nlevs_shrink is an addition to show number of gridpoints to collapse upon integer, intent(out) :: z_ind(2) real(r8), intent(out) :: z_val(2) @@ -916,11 +910,11 @@ end subroutine height_interpolate_s_grid !> of the observation. Then interpolate the height at each corner of the !> bounding box to the observation location -subroutine height_interpolate(obs_loc_array, varid, nlevs, x_ind, y_ind, z_ind, x_val, y_val, z_val, istatus) +subroutine height_interpolate(obs_loc_array, varid, nlevs, nlevs_shrink, x_ind, y_ind, z_ind, x_val, y_val, z_val, istatus) real(r8), intent(in) :: obs_loc_array(3) integer, intent(in) :: varid -integer, intent(in) :: nlevs +integer, intent(in) :: nlevs, nlevs_shrink integer, intent(in) :: x_ind(2) ! x indices on on the variable grid integer, intent(in) :: y_ind(2) ! y indices on on the variable grid integer, intent(out) :: z_ind(2) @@ -1127,6 +1121,37 @@ function observation_on_grid(obs_location, ndims) return ! exit early endif +elseif ( periodic_x) then + ! require that the point is contained within the staggered grid for the + ! y - direction since you cannot wrap-around + if ( (obs_location(1) < xf(1)) .or. (obs_location(1) > xf(nip1)) .or. & + (obs_location(2) < yh(1)) .or. (obs_location(2) > yh(nj)) ) then + + if (debug > 0) then + print *, 'periodic boundary conditions - in x-direction' + print *, 'OBSERVATION_x,y at ', obs_location(1:2), ' is off x,y grid' + print *, 'x_min, x_max ', xf(1), xf(nip1) + print *, 'y_min, y_max ', yh(1), yh(nj) + endif + + return ! exit early + endif +elseif ( periodic_y) then + ! require that the point is contained within the staggered grid for the + ! x - direction since you cannot wrap-arround + if ( (obs_location(1) < xh(1)) .or. (obs_location(1) > xh(ni)) .or. & + (obs_location(2) < yf(1)) .or. (obs_location(2) > yf(njp1)) ) then + + if (debug > 0) then + print *, 'periodic boundary conditions - in y-direction' + print *, 'OBSERVATION_x,y at ', obs_location(1:2), ' is off x,y grid' + print *, 'x_min, x_max ', xh(1), xh(ni) + print *, 'y_min, y_max ', yf(1), yf(njp1) + endif + + return ! exit early + endif + elseif ( (.not. periodic_x) .and. (.not. periodic_y) ) then ! require that the point is contained within the staggered grid for now. ! you could extrapolate for values that are within xf-xh, yf-yh @@ -1147,8 +1172,8 @@ function observation_on_grid(obs_location, ndims) endif else - write(string1,*) 'only grids with periodic x and y grids, or non-periodic ' - write(string2,*) 'boundary conditions supported' + write(string1,*) ' Unsupported combination of periodic boundary conditions.' + write(string2,*) 'Contact DART support for more help.' call error_handler(E_ERR, 'observation_on_grid', string1, & source, revision, revdate, text2=string2) endif @@ -2182,7 +2207,6 @@ subroutine parse_variable_input( state_variables, ngood ) end subroutine parse_variable_input - !=================================================================== ! End of model_mod !=================================================================== diff --git a/models/cm1/readme.rst b/models/cm1/readme.rst index f2e9f09783..a1312b5d12 100644 --- a/models/cm1/readme.rst +++ b/models/cm1/readme.rst @@ -43,7 +43,7 @@ only contain output at the analysis time (which requires setting Here is an example configuration of the ``¶m9`` namelist in ``namelist.input``: -.. code-block:: fortran +.. code-block:: text ¶m9 restart_format = 2 restart needs to be netCDF @@ -61,7 +61,7 @@ Here is an example configuration of the ``¶m9`` namelist in Additional state variables that have been tested within DART include: -``ua, va, wa, ppi, u0, v0, u10, v10, t2, th2, tsk, q2, psfc, qv, qc, qr, qi qs, & qg``. +``ua, va, wa, prs, u0, v0, u10, v10, t2, th2, tsk, q2, psfc, qv, qc, qr, qi qs, qg, nc, nr, ns, ni, nh, ng``. At present, observation times are evaluated relative to the date and time specified in the ``¶m11`` namelist. @@ -284,7 +284,7 @@ namelists start with an ampersand ``&`` and terminate with a slash ``/``. Character strings that contain a ``/`` must be enclosed in quotes to prevent them from prematurely terminating the namelist. -.. code-block:: fortran +.. code-block:: text &model_nml assimilation_period_days = 0 @@ -373,9 +373,11 @@ Description of each namelist entry .. note:: The values above are the default values. A more realistic example is shown - below and closely matches the values in the default ``input.nml``. + below and closely matches the values in the default ``input.nml``. The example input block is + for a case run using the Morrison Microphysics scheme. Any changes in microphysics will require the + user to update the hydrometeor state variables. -.. code-block:: fortran +.. code-block:: text &model_nml assimilation_period_days = 0 @@ -390,7 +392,33 @@ Description of each namelist entry 'va' , 'QTY_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', 'wa' , 'QTY_VERTICAL_VELOCITY' , 'NULL', 'NULL', 'UPDATE', 'theta', 'QTY_POTENTIAL_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', - 'ppi' , 'QTY_PRESSURE' , 'NULL', 'NULL', 'UPDATE', + 'prs' , 'QTY_PRESSURE' , 'NULL', 'NULL', 'UPDATE', + 'qv' , 'QTY_VAPOR_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qc' , 'QTY_CLOUD_LIQUID_WATER' , 0.0000, 'NULL', 'UPDATE', + 'qr' , 'QTY_RAINWATER_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qi' , 'QTY_CLOUD_ICE' , 0.0000, 'NULL', 'UPDATE', + 'qs' , 'QTY_SNOW_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qg' , 'QTY_GRAUPEL_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'ncr' , 'QTY_RAIN_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'nci' , 'QTY_ICE_NUMBER_CONCENTRATION', 0.0000, 'NULL', 'UPDATE', + 'ncs' , 'QTY_SNOW_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'ncg' , 'QTY_GRAUPEL_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'rho' , 'QTY_DENSITY' , 0.0000, 'NULL', 'UPDATE', + 'dbz' , 'QTY_RADAR_REFLECTIVITY' , 0.0000, 'NULL', 'UPDATE', + + + / + +.. note:: **From Jon Labriola on additional model variables that could be considered:** + + There are other model variables that can be included if you use a very specific forecast configuration. + The following model variables output by CM1 when you run with a simulation with + a surface model and set "output_sfcdiags = 1". These variables can be used to + simulated surface observations including TEMPERATURE_2M, U_WIND_10, V_WIND_10, SPECIFIC_HUMIDITY_2M, + and SURFACE_PRESSURE. + + .. code-block:: text + 'u10' , 'QTY_10M_U_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', 'v10' , 'QTY_10M_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', 't2' , 'QTY_2M_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', @@ -398,10 +426,20 @@ Description of each namelist entry 'tsk' , 'QTY_SURFACE_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', 'q2' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', 'psfc' , 'QTY_SURFACE_PRESSURE' , 0.0000, 'NULL', 'UPDATE', - 'qv' , 'QTY_VAPOR_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', - 'qc' , 'QTY_CLOUD_LIQUID_WATER' , 0.0000, 'NULL', 'UPDATE', - 'qr' , 'QTY_RAINWATER_MIXING_RATIO', 0.0000, 'NULL', 'UPDATE', - 'qi' , 'QTY_CLOUD_ICE' , 0.0000, 'NULL', 'UPDATE', - 'qs' , 'QTY_SNOW_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', - 'qg' , 'QTY_GRAUPEL_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE' - / + + If you want to assimilate pseudo "near-surface" observations but are not using a surface model + in your forecast, I recommend defining a single radionsonde or dropsonde observation that is located + near the surface (but at or above the lowest model level. The DART forward operator will perform + 3D interpolation to obtain the near-surface observation. + + Also be warned - from what I can tell DART forward operators are unable to calculate air temperature and + specific humidity from CM1 model fields. To simulate these fields (e.g., DROPSONDE_TEMPERATURE, RADIOSONDE_SPECIFIC_HUMIDITY,...) + you will have to manually go into the CM1 code (./src/writeout_nc.F and ./src/restart.F) and update + the model to output 3D air temperature (air_temp) and specific humidity (spec_hum) fields. Then you + can add the following fields to model_variables. + + .. code-block:: text + + 'air_temp' , 'QTY_TEMPERATURE' , 'NULL', 'NULL', 'UPDATE', + 'spec_hum' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', + diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index 1de01db7dc..94e409c15f 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -15,12 +15,11 @@ first_obs_seconds = -1 last_obs_days = -1 last_obs_seconds = -1 - distributed_state = .false. - ens_size = 20 - num_output_obs_members = 20 - num_output_state_members = 20 + ens_size = 40 + num_output_obs_members = 40 + num_output_state_members = 40 output_interval = 1 - output_members = .true. + output_members = .true. output_mean = .true. output_sd = .true. stages_to_write = 'preassim', 'analysis', 'output' @@ -32,66 +31,135 @@ num_groups = 1 output_forward_op_errors = .false. - inf_flavor = 2, 0 - inf_initial_from_restart = .false., .false. - inf_sd_initial_from_restart = .false., .false. - inf_initial = 1.0, 1.0 + + inf_flavor = 2, 0, + inf_initial_from_restart = .true., .false., + inf_sd_initial_from_restart = .true., .false., + inf_deterministic = .true., .true., + inf_initial = 1.0, 0.90, inf_lower_bound = 1.0, 1.0 - inf_upper_bound = 1000000.0, 1000000.0 + inf_upper_bound = 100.0, 100.0 inf_damping = 0.9, 0.9 inf_sd_initial = 0.6, 0.6 inf_sd_lower_bound = 0.6, 0.6 inf_sd_max_change = 1.05, 1.05 + / +! The model_nml namelist is currently made for a simulation run with the Morrison Microphysics scheme. +! The hydrometeor fields will need to be updated if a different microphysics option is used +! +! Several surface observations can only be assimilated if you run the CM1 +! with a surace parameterization and set "output_sfcdiags = 1". +! The following output variables include: +! +! Model Output Variables -> Assimilated Observations +! 1.) t2 -> TEMPERATURE_2M +! 2.) u10 -> U_WIND_10 +! 3.) v10 -> V_WIND_10 +! 4.) q2 -> SPECIFIC_HUMIDITY_2M +! 5.) psfc -> SURFACE_PRESSURE +! +! RADIOSONDE_TEMPERATURE and RADIOSONDE_SPECIFIC_HUMIDITY observations require +! air temperature and specific humidity as inputs. To do this you will have to +! go into CM1 code and tell the model to output three-dimensional air temperature and +! specific humidity variables in "./src/writeout_nc.F" and "./src/restart.F" +! + &model_nml assimilation_period_days = 0 assimilation_period_seconds = 60 calendar = 'Gregorian' cm1_template_file = 'cm1out_rst_000001.nc' - periodic_x = .true. - periodic_y = .true. + periodic_x = .false. + periodic_y = .false. debug = 0 - model_variables = 'ua' , 'QTY_U_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', - 'va' , 'QTY_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', - 'wa' , 'QTY_VERTICAL_VELOCITY' , 'NULL', 'NULL', 'UPDATE', - 'theta', 'QTY_POTENTIAL_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', - 'ppi' , 'QTY_PRESSURE' , 'NULL', 'NULL', 'UPDATE', - 'u10' , 'QTY_10M_U_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', - 'v10' , 'QTY_10M_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', - 't2' , 'QTY_2M_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', - 'th2' , 'QTY_POTENTIAL_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', - 'tsk' , 'QTY_SKIN_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', - 'q2' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', - 'qv' , 'QTY_VAPOR_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', - 'qc' , 'QTY_CLOUD_LIQUID_WATER' , 0.0000, 'NULL', 'UPDATE', - 'qr' , 'QTY_RAINWATER_MIXING_RATIO', 0.0000, 'NULL', 'UPDATE', - 'qi' , 'QTY_CLOUD_ICE' , 0.0000, 'NULL', 'UPDATE', - 'qs' , 'QTY_SNOW_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', - 'qg' , 'QTY_GRAUPEL_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', - 'psfc' , 'QTY_SURFACE_PRESSURE' , 0.0000, 'NULL', 'UPDATE' + model_variables = 'ua' , 'QTY_U_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', + 'va' , 'QTY_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', + 'wa' , 'QTY_VERTICAL_VELOCITY' , 'NULL', 'NULL', 'UPDATE', + 'theta', 'QTY_POTENTIAL_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', + 'prs' , 'QTY_PRESSURE' , 'NULL', 'NULL', 'UPDATE', + 'qv' , 'QTY_VAPOR_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qc' , 'QTY_CLOUD_LIQUID_WATER' , 0.0000, 'NULL', 'UPDATE', + 'qr' , 'QTY_RAINWATER_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qi' , 'QTY_CLOUD_ICE' , 0.0000, 'NULL', 'UPDATE', + 'qs' , 'QTY_SNOW_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'qg' , 'QTY_GRAUPEL_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'ncr' , 'QTY_RAIN_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'nci' , 'QTY_ICE_NUMBER_CONCENTRATION', 0.0000, 'NULL', 'UPDATE', + 'ncs' , 'QTY_SNOW_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'ncg' , 'QTY_GRAUPEL_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'rho' , 'QTY_DENSITY' , 0.0000, 'NULL', 'UPDATE', + 'dbz' , 'QTY_RADAR_REFLECTIVITY' , 0.0000, 'NULL', 'UPDATE', + + ! Model variables that are output when you run with a simulation with + ! a surface model and set "output_sfcdiags = 1" + ! + !'u10' , 'QTY_10M_U_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', + !'v10' , 'QTY_10M_V_WIND_COMPONENT' , 'NULL', 'NULL', 'UPDATE', + !'t2' , 'QTY_2M_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', + !'q2' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', + !'th2' , 'QTY_POTENTIAL_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', + !'tsk' , 'QTY_SKIN_TEMPERATURE' , 0.0000, 'NULL', 'UPDATE', + !'psfc' , 'QTY_SURFACE_PRESSURE' , 0.0000, 'NULL', 'UPDATE', + + ! Model variables that must be output by CM1 to assimilate RADIOSONDE TEMPERATURE + ! and RADIOSONDE_SPECIFIC_HUMIDITY observations + ! + !'air_temp' , 'QTY_TEMPERATURE' , 'NULL', 'NULL', 'UPDATE', + !'spec_hum' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', / +! +! This is a non-exhaustive list of variables that could be assimilated into CM1 forecasts +! +! Remember: TEMPERATURE_2M, U_WIND_10M, V_WIND_10M, SURFACE_PRESSURE, and SPECIFIC_HUMIDITY_2M +! require that CM1 is run with a surface model and outputs surface diagnostic variables +! +! If you are running CM1 simulations without a surface model you could theoretically assimilate a single +! RADIOSONDE observation that is near the surface to act as a "surface observation". DART will then ingest a +! 3D model field rather than referring to the 2D surface model output. +! +! Also remember RADIOSONDE_TEMPERATURE and RADIOSONDE_SPECIFIC_HUMIDITY observations +! requires CM1 to output three-dimensional arrays of air tempeature (i.e., air_temp) +! and specific humidity (spec_hum) +! &obs_kind_nml - assimilate_these_obs_types = 'TEMPERATURE_2M', - 'U_WIND_10M', - 'V_WIND_10M', - 'SURFACE_PRESSURE', - 'SPECIFIC_HUMIDITY_2M' + + assimilate_these_obs_types = 'RADIOSONDE_U_WIND_COMPONENT', + 'RADIOSONDE_V_WIND_COMPONENT', + 'RADAR_REFLECTIVITY', + 'DOPPLER_RADIAL_VELOCITY', + + ! Observations that require CM1 to output 3D fields of air temperature and specific humidity + ! This will require you to make minor modifications to the code + ! 'RADIOSONDE_TEMPERATURE', + ! 'RADIOSONDE_SPECIFIC_HUMIDITY' + + ! Observations that require CM1 simulations to be run with a surface model + ! Then set "output_sfcdiags = 1" + ! 'TEMPERATURE_2M', + ! 'U_WIND_10M', + ! 'V_WIND_10M', + ! 'SURFACE_PRESSURE', + ! 'SPECIFIC_HUMIDITY_2M' + + + evaluate_these_obs_types = '' / &location_nml - x_is_periodic = .true. + x_is_periodic = .false. y_is_periodic = .true. z_is_periodic = .false. min_x_for_periodic = 0.0 - max_x_for_periodic = 128000.0 + max_x_for_periodic = 200000.0 min_y_for_periodic = 0.0 - max_y_for_periodic = 128000.0 + max_y_for_periodic = 200000.0 min_z_for_periodic = 0.0 max_z_for_periodic = 1.0 compare_to_correct = .false. @@ -131,7 +199,8 @@ 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' obs_type_files = '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', - '../../../observations/forward_operators/obs_def_surface_mod.f90' + '../../../observations/forward_operators/obs_def_surface_mod.f90', + '../../../observations/forward_operators/obs_def_radar_mod.f90' quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' / @@ -148,7 +217,7 @@ &assim_tools_nml adaptive_localization_threshold = -1 - cutoff = 4000.0 + cutoff = 15000.0 filter_kind = 1 print_every_nth_obs = 100 rectangular_quadrature = .true. @@ -159,6 +228,8 @@ distribute_mean = .true. output_localization_diagnostics = .false. localization_diagnostics_file = 'localization_diagnostics' + special_localization_obs_types = 'RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_localization_cutoffs = 50000., 50000., 6000., 6000. / @@ -291,3 +362,32 @@ interp_test_zrange = 100.0, 101.0 verbose = .true. / + +! This block of code is for assimilating radar observations +! Rather than diagnosing radar reflectivity directly from model output +! I feed the DART system the 3D radar reflectivity variable output by cm1. +! If you decide to follow this same technique set microphysics_type = 5. +! +! To calculate radar radial velocity you need to set allow_dbztowt_conv=.true. +! this will allow the DART to diagnose particle fall speed via reflecvtivity +! which is used for the radial velocity calculation + +&obs_def_radar_mod_nml + apply_ref_limit_to_obs = .false., + reflectivity_limit_obs = -10.0, + lowest_reflectivity_obs = -10.0, + apply_ref_limit_to_fwd_op = .false., + reflectivity_limit_fwd_op = -10.0, + lowest_reflectivity_fwd_op = -10.0, + max_radial_vel_obs = 1000000, + allow_wet_graupel = .false., + microphysics_type = 5, + allow_dbztowt_conv = .true., + dielectric_factor = 0.224, + n0_rain = 8.0e6, + n0_graupel = 4.0e6, + n0_snow = 3.0e6, + rho_rain = 1000.0, + rho_graupel = 400.0, + rho_snow = 100.0 + /