From 46c85c3e5801075360910831cd5b6900ae62a149 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Mon, 4 Oct 2021 20:36:53 +0000 Subject: [PATCH 01/15] Updated the initialize module in location_mod.f90 to handle vertical localization coordinates --- .../threed_cartesian/location_mod.f90 | 186 ++++++++++++++++-- 1 file changed, 175 insertions(+), 11 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 75320ea98a..5956d6bbd2 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -45,6 +45,11 @@ module location_mod character(len = 129), parameter :: LocationStorageOrder = "X Y Z" character(len = 129), parameter :: LocationUnits = "none none none" +! The possible numeric values for the location_type%which_vert component. +! The numeric values are PRIVATE to this module. The parameter names are PUBLIC. +integer, parameter :: VERTISHEIGHT = 1 ! by height (in meters) +!--- Can eventually add other heights if necesary (JDL) + type location_type private real(r8) :: x, y, z @@ -96,22 +101,28 @@ module location_mod + +! Support more than a single cutoff distance. nt is the count of +! distinct cutoffs, which are selected per specific observation type. +! The map associates the incoming location type with the +! corresponding gtt index. There are only as many close_types +! as there are distinct cutoff distances; if more than one specific +! type has the same cutoff distances they share the type. type get_close_type private 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 get_close_type -! private -! integer :: num -! real(r8) :: maxdist -! type(box_type) :: box -!end type get_close_type - +! Some calls include information about the type or kind of the location. +! If the location refers to an observation it is possible to have both +! a specific type and a generic kind. If the location refers to a +! state vector item, it only has a generic kind. Some routines have +! a specific type and a generic kind as arguments; look carefully at +! the argument names before using. - Message stollen from 3d_Sphere + type(random_seq_type) :: ran_seq logical :: ran_seq_init = .false. logical, save :: module_initialized = .false. @@ -119,6 +130,28 @@ module location_mod character(len = 512) :: errstring character(len = 512) :: msgstring, msgstring1, msgstring2 +! JDL NEW +! Horizontal localization/cutoff values are passed in by the caller. +! The Vertical normalization values are globals; are specified by namelist +! here, and apply to all specific types unless a 'special list' is also specified +! that overrides the default values. + +logical :: has_special_vertical_norms = .false. +integer :: num_special_vert_norms = 0 +integer :: location_vertical_localization_coord = 0 + +! JDL NEW +! Global storage for vertical distance normalization factors +! The 4 below is for the 4 vertical units (pressure, level, height, +! scale height). undefined and surface don't need vert distances. +! NOTE: Code that uses VERT_TYPE_COUNT depends on pressure, level, +! height, and scale height having actual values between 1 and 4, or +! this code will break. +integer, parameter :: VERT_TYPE_COUNT = 1 ! JDL JUST WORKING ON HEIGHT +real(r8) :: vert_normalization(VERT_TYPE_COUNT) +real(r8), allocatable :: per_type_vert_norm(:,:) ! if doing per-type + + ! for sanity when i'm using arrays of length(3): integer, parameter :: IX = 1 integer, parameter :: IY = 2 @@ -173,12 +206,34 @@ module location_mod integer :: ny = 10 integer :: nz = 10 + +!--- JDL: Additional Items for Vertical Localization +!------------------------------------------------------------------------ +! Namelist with default values +! horiz_dist_only == .true. -> Only the great circle horizontal distance is +! computed in get_dist. +! horiz_dist_only == .false. -> Square root of sum of squared horizontal and +! normalized vertical dist computed in get_dist +! vert_normalization_height -> Number meters that give a distance equivalent +! to one radian in horizontal +! special_vert_normalization_obs_types -> Which obs types to modify the default vert +! normalization values +! special_vert_normalization_pressure -> must give all 4 values for each type listed + +logical :: horiz_dist_only = .true. ! JDL While DEFAULT - You'll want false +real(r8) :: vert_normalization_height = 10000.0_r8 +integer, parameter :: MAX_ITEMS = 500 +character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS) +real(r8) :: special_vert_normalization_heights(MAX_ITEMS) + namelist /location_nml/ & - compare_to_correct, output_box_info, print_box_level, & + compare_to_correct, output_box_info, print_box_level, & x_is_periodic, min_x_for_periodic, max_x_for_periodic, & y_is_periodic, min_y_for_periodic, max_y_for_periodic, & z_is_periodic, min_z_for_periodic, max_z_for_periodic, & - nx, ny, nz, debug + nx, ny, nz, horiz_dist_only, vert_normalization_height, & + special_vert_normalization_obs_types, & + special_vert_normalization_heights, debug !> @todo: @@ -230,12 +285,16 @@ subroutine initialize_module ! things which need doing exactly once. -integer :: iunit, io +integer :: iunit, io, i, k, typecount, type_index if (module_initialized) return module_initialized = .true. +! give these initial values before reading them from the namelist. +special_vert_normalization_obs_types(:) = 'null' +special_vert_normalization_heights(:) = missing_r8 + ! Read the namelist entry call find_namelist_in_file("input.nml", "location_nml", iunit) read(iunit, nml = location_nml, iostat = io) @@ -246,6 +305,111 @@ subroutine initialize_module if(do_nml_file()) write(nmlfileunit, nml=location_nml) if(do_nml_term()) write( * , nml=location_nml) +! Copy the normalization factors in the vertical into an array +vert_normalization(VERTISHEIGHT) = vert_normalization_height + +! If the user has set a namelist to make different vertical normalization factors +! when computing the localization distances, allocate and set that up here. +! It overrides the defaults above. + +if (special_vert_normalization_obs_types(1) /= 'null' .or. & + special_vert_normalization_heights(1) /= missing_r8 .or. & + + ! FIXME: add code to check for mismatched length lists. are we going to force + ! users to specify all 4 values for any obs type that is not using the defaults? + + typecount = get_num_types_of_obs() ! ignore function name, this is specific type count + allocate(per_type_vert_norm(VERT_TYPE_COUNT, typecount)) + + ! Set the defaults for all specific types not listed in the special list + per_type_vert_norm(VERTISHEIGHT, :) = vert_normalization_height + + ! Go through special-treatment observation kinds, if any. + num_special_vert_norms = 0 + k = 0 + do i = 1, MAX_ITEMS + if(special_vert_normalization_obs_types(i) == 'null') exit + k = k + 1 + enddo + num_special_vert_norms = k + + if (num_special_vert_norms > 0) has_special_vertical_norms = .true. + + do i = 1, num_special_vert_norms + type_index = get_index_for_type_of_obs(special_vert_normalization_obs_types(i)) + if (type_index < 0) then + write(msgstring, *) 'unrecognized TYPE_ in the special vertical normalization namelist:' + call error_handler(E_ERR,'location_mod:', msgstring, source, & + text2=trim(special_vert_normalization_obs_types(i))) + endif + + per_type_vert_norm(VERTISHEIGHT, type_index) = special_vert_normalization_heights(i) + + enddo + + if (any(per_type_vert_norm == missing_r8)) then + write(msgstring, *) 'one or more special vertical normalization values is uninitialized.' + call error_handler(E_ERR,'location_mod:', & + 'special vert normalization value namelist requires all 4 values per type', & + source, text2=trim(msgstring)) + endif + +endif +if (horiz_dist_only) then + call error_handler(E_MSG,'location_mod:', & + 'Ignoring vertical separation when computing distances; horizontal distances only', & + source) +else + call error_handler(E_MSG,'location_mod:', & + 'Including vertical separation when computing distances:', source) + write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height + call error_handler(E_MSG,'location_mod:',msgstring,source) + + if (allocated(per_type_vert_norm)) then + typecount = get_num_types_of_obs() ! ignore function name, this is specific type count + do i = 1, typecount + !if ((per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) .or. & + ! (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) .or. & + ! (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) .or. & + ! (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height)) then + + ! write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) + ! call error_handler(E_MSG,'location_mod:',msgstring,source) + ! if (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) then + ! write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', & + ! per_type_vert_norm(VERTISPRESSURE, i) + ! call error_handler(E_MSG,'location_mod:',msgstring,source) + ! endif + ! if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then + ! write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', & + ! per_type_vert_norm(VERTISHEIGHT, i) + ! call error_handler(E_MSG,'location_mod:',msgstring,source) + ! endif + ! if (per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) then + ! write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', & + ! per_type_vert_norm(VERTISLEVEL, i) + ! call error_handler(E_MSG,'location_mod:',msgstring,source) + ! endif + ! if (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height) then + ! write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', & + ! per_type_vert_norm(VERTISSCALEHEIGHT, i) + ! call error_handler(E_MSG,'location_mod:',msgstring,source) + ! endif + ! JDL - Strip Down to Just VERTISHEIGHT Right Now (CM1 is in meter coordinates only right now) + if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then + write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) + call error_handler(E_MSG,'location_mod:',msgstring,source) + if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then + write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz meter: ', & + per_type_vert_norm(VERTISHEIGHT, i) + call error_handler(E_MSG,'location_mod:',msgstring,source) + endif + endif + enddo + endif +endif + + end subroutine initialize_module !---------------------------------------------------------------------------- From e5d9b391d03a7f992afb4f83e19d7d0350bd0da1 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Tue, 5 Oct 2021 15:45:52 +0000 Subject: [PATCH 02/15] Updated Assimilation Code to handle vertical localization. This is the first version that was able to compile but likely contains bugs --- .../threed_cartesian/location_mod.f90 | 96 ++++++++++++++----- 1 file changed, 72 insertions(+), 24 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 5956d6bbd2..8f5d7a71e9 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -8,7 +8,7 @@ module location_mod ! Has interfaces to convert spherical lat/lon coords in degrees, ! plus a radius, into cartesian coords. -use types_mod, only : r8, i8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD +use types_mod, only : r8, i8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD, OBSTYPELENGTH use utilities_mod, only : error_handler, E_ERR, ascii_file_format, & E_MSG, open_file, close_file, set_output, & logfileunit, nmlfileunit, find_namelist_in_file, & @@ -313,7 +313,7 @@ subroutine initialize_module ! It overrides the defaults above. if (special_vert_normalization_obs_types(1) /= 'null' .or. & - special_vert_normalization_heights(1) /= missing_r8 .or. & + special_vert_normalization_heights(1) /= missing_r8) then ! FIXME: add code to check for mismatched length lists. are we going to force ! users to specify all 4 values for any obs type that is not using the defaults? @@ -416,6 +416,10 @@ end subroutine initialize_module function get_dist(loc1, loc2, type1, kind2) +! JDL I don't think I need to include the no_vert variable, but keep this +! in mind in case problems arise later on. Probbably also don't need kind2 since +! we are just looking at heights + ! returns the distance between 2 locations ! the names are correct here - the first location gets the corresponding @@ -432,6 +436,8 @@ function get_dist(loc1, loc2, type1, kind2) real(r8) :: diff(3), next_diff(3) logical :: below_L1(3), below_L2(3) real(r8) :: square_dist, this_dist +!logical :: comp_h_only ! Only comput the horizontal component + if ( .not. module_initialized ) call initialize_module @@ -442,13 +448,16 @@ function get_dist(loc1, loc2, type1, kind2) diff(IZ) = loc1%z - loc2%z ! hope for the simple case first +! JDL - In this instance you will want to pass +! on what the vertical normmalization component is... if (.not. any_periodic) then if (debug > 0) write(0,*) 'non-periodic distance: ' - get_dist = dist_3d(diff) + get_dist = dist_3d(diff,type1=type1) return + else if (debug > 0) write(0,*) 'periodic distance: ' - square_dist = dist_3d_sq(diff) + square_dist = dist_3d_sq(diff,type1=type1) if (debug > 0) write(0,*) 'non-periodic 0 sq_dist, dist: ', square_dist, sqrt(square_dist) endif @@ -495,7 +504,7 @@ function get_dist(loc1, loc2, type1, kind2) next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) endif endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic XYZ dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -528,7 +537,7 @@ function get_dist(loc1, loc2, type1, kind2) endif endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic XY dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -561,7 +570,7 @@ function get_dist(loc1, loc2, type1, kind2) next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) endif endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic XZ dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -593,7 +602,7 @@ function get_dist(loc1, loc2, type1, kind2) next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) endif endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic YZ dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -615,7 +624,8 @@ function get_dist(loc1, loc2, type1, kind2) else next_diff(IX) = loopy%offset(IX) + next_diff(IX) endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) + if (debug > 0) write(0,*) 'periodic X dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -637,7 +647,7 @@ function get_dist(loc1, loc2, type1, kind2) next_diff(IY) = loopy%offset(IY) + next_diff(IY) endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic Y dist: ', square_dist, this_dist if(this_dist < square_dist) then @@ -660,7 +670,7 @@ function get_dist(loc1, loc2, type1, kind2) else next_diff(IZ) = loopy%offset(IZ) + next_diff(IZ) endif - this_dist = dist_3d_sq(next_diff) + this_dist = dist_3d_sq(next_diff,type1=type1) if (debug > 0) write(0,*) 'periodic Y dist: ', square_dist, this_dist if(this_dist < square_dist) then square_dist = this_dist @@ -687,14 +697,31 @@ end function get_dist ! return the 3d distance given the separation along each axis !pure function dist_3d(separation) -function dist_3d(separation) result(val) +!function dist_3d(separation) result(val) +function dist_3d(separation,type1) result(val) real(r8), intent(in) :: separation(3) -real(r8) :: val +real(r8) :: val, vert_normal +integer, optional, intent(in) :: type1 ! JDL Addition +!--- JDL Update Function to account for Normalization in Vertical +!--- JDL Addition to Account for Vertical Normalization +if (allocated(per_type_vert_norm)) then + if (.not. present(type1)) then + write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' + call error_handler(E_MSG, 'get_dist', msgstring, source) + endif + print*,'JDL type1 = ',type1 + vert_normalization = separation(IZ)/per_type_vert_norm(VERTISHEIGHT, type1) +else + vert_normalization = separation(IZ)/vert_normalization(VERTISHEIGHT) +endif +!--- JDL End Addition for Vertical Normalization + val = sqrt(separation(IX)*separation(IX) + & separation(IY)*separation(IY) + & - separation(IZ)*separation(IZ) ) + vert_normal*vert_normal ) + !separation(IZ)*separation(IZ) ) if (debug > 0) write(0,*) 'dist_3d called, distance computed: ', val if (debug > 0) write(0,*) 'XYZ separations: ', separation @@ -706,14 +733,31 @@ end function dist_3d ! (saves doing a square root) !pure function dist_3d_sq(separation) -function dist_3d_sq(separation) result(val) +!function dist_3d_sq(separation) result(val) +function dist_3d_sq(separation,type1) result(val) real(r8), intent(in) :: separation(3) -real(r8) :: val +real(r8) :: val, vert_normal +integer, optional, intent(in) :: type1 ! JDL Addition + +!--- JDL Update Function to account for Normalization in Vertical +!--- JDL Addition to Account for Vertical Normalization +if (allocated(per_type_vert_norm)) then + if (.not. present(type1)) then + write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' + call error_handler(E_MSG, 'get_dist', msgstring, source) + endif + print*,'JDL type1 = ',type1 + vert_normal = separation(IZ) / per_type_vert_norm(VERTISHEIGHT, type1) +else + vert_normal = separation(IZ) / vert_normalization(VERTISHEIGHT) +endif +!--- JDL End Addition for Vertical Normalization val = separation(IX)*separation(IX) + & separation(IY)*separation(IY) + & - separation(IZ)*separation(IZ) + vert_normal * vert_normal + !separation(IZ)*separation(IZ) if (debug > 0) write(0,*) 'dist_3d_sq called, distance computed: ', val if (debug > 0) write(0,*) 'XYZ separations: ', separation @@ -1403,11 +1447,11 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & !> but should give the right answer. if(.true.) then if (present(dist)) then - call exhaustive_collect(gc%box(bt), base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, base_type, locs, & num_close, close_ind, dist) else allocate(cdist(size(locs))) - call exhaustive_collect(gc%box(bt), base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, base_type, locs, & num_close, close_ind, cdist) deallocate(cdist) endif @@ -1418,7 +1462,7 @@ 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%box(bt), base_loc, locs, & + call exhaustive_collect(gc%box(bt), base_loc, base_type, locs, & cnum_close, cclose_ind, cdist) endif @@ -1477,7 +1521,7 @@ subroutine get_close(gc, base_loc, base_type, locs, loc_qtys, & ! Only compute distance if dist is present if(present(dist)) then - this_dist = get_dist(base_loc, locs(t_ind)) + this_dist = get_dist(base_loc, locs(t_ind), type1=base_type) !write(0,*) 'this_dist = ', this_dist ! If this loc's distance is less than cutoff, add it in list if(this_dist <= this_maxdist) then @@ -1656,6 +1700,9 @@ subroutine find_nearest(box, base_loc, loc_list, nearest, rc) !write(0,*) 'x: ', start_x, end_x !write(0,*) 'y: ', start_y, end_y !write(0,*) 'z: ', start_z, end_z +! JDL WARNING STATEMENT IN CASE YOU HAPPEN TO CALL GET_DIST() FROM THIS SUBROUTINE +! YOU SHOULDN'T BE USING THIS SUBROUTINE +print*,'JDL IF YOUR HERE MAKE SURE TO WORK ON FIND_NEAREST SUBROUTINE' ! Next, loop through each box that is close to this box do i = start_x, end_x @@ -1673,7 +1720,7 @@ subroutine find_nearest(box, base_loc, loc_list, nearest, rc) t_ind = box%loc_box(st - 1 + l) !write(0,*) 'l, t_ind = ', l, t_ind - this_dist = get_dist(base_loc, loc_list(t_ind)) + this_dist = get_dist(base_loc, loc_list(t_ind)) !JDL Probably do not need to include base type here because this is never called during DA (OR IT SHOULDNT BE) !write(0,*) 'this_dist = ', this_dist ! If this loc's distance is less than current nearest, it's new nearest if(this_dist <= dist) then @@ -2215,13 +2262,14 @@ end function is_location_in_region !--------------------------------------------------------------------------- -subroutine exhaustive_collect(box, base_loc, loc_list, num_close, close_ind, close_dist) +subroutine exhaustive_collect(box, base_loc, base_type, loc_list, num_close, close_ind, close_dist) ! For validation, it is useful to be able to compare against exact ! exhaustive search type(box_type), intent(in) :: box type(location_type), intent(in) :: base_loc, loc_list(:) +integer, intent(in) :: base_type ! JDL Addition integer, intent(out) :: num_close integer, intent(out) :: close_ind(:) real(r8), intent(out) :: close_dist(:) @@ -2231,7 +2279,7 @@ subroutine exhaustive_collect(box, base_loc, loc_list, num_close, close_ind, clo num_close = 0 do i = 1, box%num - this_dist = get_dist(base_loc, loc_list(i)) + this_dist = get_dist(base_loc, loc_list(i), type1=base_type) if(this_dist <= box%maxdist) then ! Add this loc to correct list num_close = num_close + 1 From 9d72fe37a57435ae1004f652b7c3a467a1889468 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Thu, 21 Oct 2021 18:42:41 +0000 Subject: [PATCH 03/15] Added print statments to find out why interpolation crashes near boundary edges. Do not interpolate outside of the second unstaggered gridpoint --- models/cm1/model_mod.f90 | 45 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/models/cm1/model_mod.f90 b/models/cm1/model_mod.f90 index 5fde4101f0..3aa1be4ad7 100644 --- a/models/cm1/model_mod.f90 +++ b/models/cm1/model_mod.f90 @@ -885,7 +885,6 @@ subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, Q12(level) = zhalf(x_ind(1), y_ind(2), level) Q21(level) = zhalf(x_ind(2), y_ind(1), level) Q22(level) = zhalf(x_ind(2), y_ind(2), level) - enddo else ! on z half levels @@ -903,6 +902,20 @@ subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, Z(:) = bilinear_interpolation(nlevs, obs_loc_array(1), obs_loc_array(2), & x_val(1), x_val(2), y_val(1), y_val(2), Q11, Q12, Q21, Q22) +do level = 1, nlevs + if (isnan(Z(level))) then + print*,'JDL x_val(1)',x_val(1) + print*,'JDL x_val(2)',x_val(2) + print*,'JDL y_val(1)',y_val(1) + print*,'JDL y_val(2)',y_val(2) + print*,'JDL OBS ARRAY(1)',obs_loc_array(1) + print*,'JDL OBS ARRAY(2)',obs_loc_array(2) + print*,'JDL Q11',Q11(level) + print*,'JDL Q12',Q12(level) + print*,'JDL Q21',Q21(level) + print*,'JDL Q22',Q22(level) + endif +enddo ! Find out which level the point is in: call get_enclosing_coord(obs_loc_array(3), Z, z_ind, z_val) @@ -1083,6 +1096,24 @@ subroutine height_interpolate(obs_loc_array, varid, nlevs, nlevs_shrink, x_ind, x_val(1), x_val(2), y_val(1), y_val(2), Q11, Q12, Q21, Q22) +!do level = 1, nlevs +xloop: do level = 1, nlevs + if (isnan(Z(level))) then + print*,'JDL x_val(1)',x_val(1) + print*,'JDL x_val(2)',x_val(2) + print*,'JDL y_val(1)',y_val(1) + print*,'JDL y_val(2)',y_val(2) + print*,'JDL OBS ARRAY(1)',obs_loc_array(1) + print*,'JDL OBS ARRAY(2)',obs_loc_array(2) + print*,'JDL Q11',Q11(level) + print*,'JDL Q12',Q12(level) + print*,'JDL Q21',Q21(level) + print*,'JDL Q22',Q22(level) + print*,'**********' + exit + endif +enddo xloop + ! Find out which level the point is in: call get_enclosing_coord(obs_loc_array(3), Z, z_ind, z_val) @@ -1550,6 +1581,18 @@ subroutine get_enclosing_coord(x, xcoords, ind, val) enddo xloop if (ind(1) == -999) then + + print*,'JDL X = ',x + print*,'JDL COORD SIZE =',size(xcoords) + print*,'JDL LOC(4) = ',xcoords(4) + !do i = 2, size(xcoords) + ! if(isnan(xcoords(i))) THEN + ! print*,'JDL THIS IS A NAN' + ! endif + + ! print*,'JDL coord = ',xcoords(i) + !enddo + call error_handler(E_ERR, 'get_enclosing_coord', 'off the grid, unexpected 3', & source, revision, revdate) endif From d7bd4353c1dce98688da9b02289f002dc67b4876 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Thu, 21 Oct 2021 18:44:20 +0000 Subject: [PATCH 04/15] Updated input namelist to have vertical localization --- models/cm1/work/input.nml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index b3d14558ca..11049cba2a 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -103,6 +103,9 @@ nx = 10 ny = 10 nz = 10 + vert_normalization_height = 1 + special_vert_normalization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_SPECIFIC_HUMIDITY', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_vert_normalization_heights= 0.0075, 0.0075, 0.0075, 0.0075, 0.0800, 0.0800, 0.0833, 0.1250, 0.5000, 0.5000, / @@ -163,8 +166,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', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_SPECIFIC_HUMIDITY', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' - special_localization_cutoffs = 30000.,30000.,30000.,30000.,6000.,6000. + special_localization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_SPECIFIC_HUMIDITY', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_localization_cutoffs = 400000.00, 400000.00, 400000.00, 400000.00, 25000.0000, 25000.0000, 15000.0000, 10000.0000, 6000.0000, 6000.0000, / From 1110d824eec5d86d4c5ea3a11e98fa7840d161b8 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Thu, 21 Oct 2021 18:45:09 +0000 Subject: [PATCH 05/15] Updated the threed_cartesian module to include vertical localization --- .../threed_cartesian/location_mod.f90 | 61 ++++++++++--------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 8f5d7a71e9..7ba15c8011 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -210,6 +210,7 @@ module location_mod !--- JDL: Additional Items for Vertical Localization !------------------------------------------------------------------------ ! Namelist with default values +!--- horiz_dist_only is currently ignored (JDL) ! horiz_dist_only == .true. -> Only the great circle horizontal distance is ! computed in get_dist. ! horiz_dist_only == .false. -> Square root of sum of squared horizontal and @@ -220,7 +221,7 @@ module location_mod ! normalization values ! special_vert_normalization_pressure -> must give all 4 values for each type listed -logical :: horiz_dist_only = .true. ! JDL While DEFAULT - You'll want false +!logical :: horiz_dist_only = .true. ! JDL While DEFAULT - You'll want false real(r8) :: vert_normalization_height = 10000.0_r8 integer, parameter :: MAX_ITEMS = 500 character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS) @@ -231,7 +232,8 @@ module location_mod x_is_periodic, min_x_for_periodic, max_x_for_periodic, & y_is_periodic, min_y_for_periodic, max_y_for_periodic, & z_is_periodic, min_z_for_periodic, max_z_for_periodic, & - nx, ny, nz, horiz_dist_only, vert_normalization_height, & + !nx, ny, nz, horiz_dist_only, vert_normalization_height, & + nx, ny, nz, vert_normalization_height, & special_vert_normalization_obs_types, & special_vert_normalization_heights, debug @@ -355,19 +357,19 @@ subroutine initialize_module endif endif -if (horiz_dist_only) then - call error_handler(E_MSG,'location_mod:', & - 'Ignoring vertical separation when computing distances; horizontal distances only', & - source) -else - call error_handler(E_MSG,'location_mod:', & - 'Including vertical separation when computing distances:', source) - write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height - call error_handler(E_MSG,'location_mod:',msgstring,source) - - if (allocated(per_type_vert_norm)) then - typecount = get_num_types_of_obs() ! ignore function name, this is specific type count - do i = 1, typecount +!if (horiz_dist_only) then +! call error_handler(E_MSG,'location_mod:', & +! 'Ignoring vertical separation when computing distances; horizontal distances only', & +! source) +!else +call error_handler(E_MSG,'location_mod:', & + 'Including vertical separation when computing distances:', source) +write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height +call error_handler(E_MSG,'location_mod:',msgstring,source) + +if (allocated(per_type_vert_norm)) then + typecount = get_num_types_of_obs() ! ignore function name, this is specific type count + do i = 1, typecount !if ((per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) .or. & ! (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) .or. & ! (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) .or. & @@ -396,18 +398,18 @@ subroutine initialize_module ! call error_handler(E_MSG,'location_mod:',msgstring,source) ! endif ! JDL - Strip Down to Just VERTISHEIGHT Right Now (CM1 is in meter coordinates only right now) - if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then - write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) + if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then + write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) + call error_handler(E_MSG,'location_mod:',msgstring,source) + if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then + write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz meter: ', & + per_type_vert_norm(VERTISHEIGHT, i) call error_handler(E_MSG,'location_mod:',msgstring,source) - if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then - write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz meter: ', & - per_type_vert_norm(VERTISHEIGHT, i) - call error_handler(E_MSG,'location_mod:',msgstring,source) - endif endif - enddo - endif + endif + enddo endif +!endif end subroutine initialize_module @@ -710,10 +712,12 @@ function dist_3d(separation,type1) result(val) write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' call error_handler(E_MSG, 'get_dist', msgstring, source) endif - print*,'JDL type1 = ',type1 - vert_normalization = separation(IZ)/per_type_vert_norm(VERTISHEIGHT, type1) + vert_normal = separation(IZ)/per_type_vert_norm(VERTISHEIGHT, type1) + !print*,'JDL per_type non-periodic = ',per_type_vert_norm(VERTISHEIGHT, type1) else - vert_normalization = separation(IZ)/vert_normalization(VERTISHEIGHT) + vert_normal = separation(IZ)/vert_normalization(VERTISHEIGHT) + !print*,'JDL standard norm non-periodic = ',vert_normalization(VERTISHEIGHT) + endif !--- JDL End Addition for Vertical Normalization @@ -747,10 +751,11 @@ function dist_3d_sq(separation,type1) result(val) write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' call error_handler(E_MSG, 'get_dist', msgstring, source) endif - print*,'JDL type1 = ',type1 vert_normal = separation(IZ) / per_type_vert_norm(VERTISHEIGHT, type1) + !print*,'JDL per_type = ',per_type_vert_norm(VERTISHEIGHT, type1) else vert_normal = separation(IZ) / vert_normalization(VERTISHEIGHT) + !print*,'JDL standard norm = ',vert_normalization(VERTISHEIGHT) endif !--- JDL End Addition for Vertical Normalization From 6dcdb6d9f53050822a898b4571cadb2d4c19b8d0 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Mon, 31 Jan 2022 20:59:49 +0000 Subject: [PATCH 06/15] Added the ability to assimilate dewpoint temperature observations in addition to other common observations. Assimilating dewpoint temperature rather than water vapor mixing ratio appears to work better because the range of observed values is much smaller. --- models/cm1/work/input.nml | 62 +++++++++++-------- .../obs_def_dew_point_mod.f90 | 18 +++++- 2 files changed, 51 insertions(+), 29 deletions(-) diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index 11049cba2a..18b97307b1 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -15,7 +15,7 @@ first_obs_seconds = -1 last_obs_days = -1 last_obs_seconds = -1 - distributed_state = .false. + distributed_state = .true. ens_size = 40 num_output_obs_members = 40 num_output_state_members = 40 @@ -23,9 +23,7 @@ output_members = .true. output_mean = .true. output_sd = .true. - !stages_to_write = 'preassim', 'analysis', 'output' - !stages_to_write = 'forecast', 'preassim', 'postassim', 'analysis', 'output' - stages_to_write = 'output' + stages_to_write = 'preassim','analysis','output' write_all_stages_at_end = .true. distributed_state = .true. compute_posterior = .true. @@ -34,12 +32,13 @@ num_groups = 1 output_forward_op_errors = .false. - inf_flavor = 0, 4 - inf_initial_from_restart = .false., .false. - inf_sd_initial_from_restart = .false., .false. - inf_initial = 1.0, 0.98 + 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 @@ -59,30 +58,31 @@ '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' , 0.0000, 'NULL', 'UPDATE', 'air_temp' , 'QTY_TEMPERATURE' , 'NULL', 'NULL', 'UPDATE', - 'spec_hum' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', + 'dewpoint' , 'QTY_DEWPOINT' , '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', + 'qhl' , 'QTY_HAIL_MIXING_RATIO' , 0.0000, 'NULL', 'UPDATE', + 'ccw' , 'QTY_DROPLET_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'crw' , 'QTY_RAIN_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'cci' , 'QTY_ICE_NUMBER_CONCENTRATION', 0.0000, 'NULL', 'UPDATE', + 'csw' , 'QTY_SNOW_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'chw' , 'QTY_GRAUPEL_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'chl' , 'QTY_HAIL_NUMBER_CONCENTR' , 0.0000, 'NULL', 'UPDATE', + 'vhw' , 'QTY_GRAUPEL_VOLUME' , 0.0000, 'NULL', 'UPDATE', + 'vhl' , 'QTY_HAIL_VOLUME' , 0.0000, 'NULL', 'UPDATE', 'dbz' , 'QTY_RADAR_REFLECTIVITY' , 0.0000, 'NULL', 'UPDATE', 'fall_spd' , 'QTY_POWER_WEIGHTED_FALL_SPEED', 0.0000, 'NULL', 'UPDATE', / &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', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_DEWPOINT', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' evaluate_these_obs_types = '' / @@ -104,8 +104,8 @@ ny = 10 nz = 10 vert_normalization_height = 1 - special_vert_normalization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_SPECIFIC_HUMIDITY', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' - special_vert_normalization_heights= 0.0075, 0.0075, 0.0075, 0.0075, 0.0800, 0.0800, 0.0833, 0.1250, 0.5000, 0.5000, + special_vert_normalization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_DEWPOINT', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_vert_normalization_heights = 0.0075, 0.0075, 0.0075, 0.0075, 0.04, 0.04, 0.04, 0.04, 0.33, 0.33 / @@ -138,8 +138,10 @@ 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_radar_mod.f90' - quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' + '../../../observations/forward_operators/obs_def_radar_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' / @@ -166,8 +168,11 @@ distribute_mean = .true. output_localization_diagnostics = .false. localization_diagnostics_file = 'localization_diagnostics' - special_localization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_SPECIFIC_HUMIDITY', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' - special_localization_cutoffs = 400000.00, 400000.00, 400000.00, 400000.00, 25000.0000, 25000.0000, 15000.0000, 10000.0000, 6000.0000, 6000.0000, + adjust_obs_impact = .false. + obs_impact_filename = "NONE" + special_localization_obs_types = 'DROPSONDE_U_WIND_COMPONENT','DROPSONDE_V_WIND_COMPONENT', 'DROPSONDE_TEMPERATURE', 'DROPSONDE_SPECIFIC_HUMIDITY','RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE', 'RADIOSONDE_DEWPOINT', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_localization_cutoffs = 400000., 400000., 400000., 400000., 50000.,50000.,50000.,50000.,6000.,6000. + / @@ -301,6 +306,11 @@ verbose = .true. / +&obs_impact_tool_nml + input_filename = 'cross_correlations.txt' + output_filename = 'control_impact_runtime.txt' + debug = .false. + / &obs_def_radar_mod_nml apply_ref_limit_to_obs = .false., diff --git a/observations/forward_operators/obs_def_dew_point_mod.f90 b/observations/forward_operators/obs_def_dew_point_mod.f90 index 8d50a7e32f..8fc890628f 100644 --- a/observations/forward_operators/obs_def_dew_point_mod.f90 +++ b/observations/forward_operators/obs_def_dew_point_mod.f90 @@ -172,8 +172,8 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist call interpolate(state_handle, ens_size,location, QTY_VAPOR_MIXING_RATIO, qv, qv_istatus) call track_status(ens_size, qv_istatus, td, istatus, return_now) if (return_now) return - - +print*,'JDL P = ',p_Pa +print*,'JDL qv = ',qv where (qv < 0.0_r8 .or. qv >= 1.0_r8) istatus = 1 td = missing_r8 @@ -204,14 +204,26 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist ! The following expression can fail numerically for dewpoints very close to 0 C !td = t_kelvin + (243.5_r8 / ((17.67_r8 / log(e_mb/6.112_r8)) - 1.0_r8) ) + + + ! A numerically robust formula that avoids the failure near dewpoints of 0 C log_term = log(e_mb / 6.112_r8) td = t_kelvin + 243.5_r8 * log_term / (17.67_r8 - log_term) + !print *,'********************************' + !print *,'e_mb = ',e_mb(1) + !print *,'p_mb = ',p_mb(1) + !print *,'qv =',qv(1) + !print *,'td = ',td(1) + !print *,'t_kelvin = ',t_kelvin + !print *,'***********' + + elsewhere td = missing_r8 end where - +print*,'JDL TD = ',td end subroutine get_expected_dew_point !---------------------------------------------------------------------------- From 465c8e3f13e9342acf77bf0956287f877cbb98d8 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Tue, 22 Feb 2022 21:16:33 +0000 Subject: [PATCH 07/15] Cleaning up the code to remove some of the early error messages I originally created --- .../threed_cartesian/location_mod.f90 | 10 ++--- models/cm1/model_mod.f90 | 45 ------------------- .../obs_def_dew_point_mod.f90 | 9 ---- 3 files changed, 4 insertions(+), 60 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 245c1c278c..dde0346e45 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -193,7 +193,9 @@ module location_mod !--- JDL: Additional Items for Vertical Localization !------------------------------------------------------------------------ ! Namelist with default values -!--- horiz_dist_only is currently ignored (JDL) +! +! horiz_dist_only is currently ignored (JDL) +! ! horiz_dist_only == .true. -> Only the great circle horizontal distance is ! computed in get_dist. ! horiz_dist_only == .false. -> Square root of sum of squared horizontal and @@ -694,10 +696,8 @@ function dist_3d(separation,type1) result(val) call error_handler(E_MSG, 'get_dist', msgstring, source) endif vert_normal = separation(IZ)/per_type_vert_norm(VERTISHEIGHT, type1) - !print*,'JDL per_type non-periodic = ',per_type_vert_norm(VERTISHEIGHT, type1) else vert_normal = separation(IZ)/vert_normalization(VERTISHEIGHT) - !print*,'JDL standard norm non-periodic = ',vert_normalization(VERTISHEIGHT) endif !--- JDL End Addition for Vertical Normalization @@ -733,10 +733,8 @@ function dist_3d_sq(separation,type1) result(val) call error_handler(E_MSG, 'get_dist', msgstring, source) endif vert_normal = separation(IZ) / per_type_vert_norm(VERTISHEIGHT, type1) - !print*,'JDL per_type = ',per_type_vert_norm(VERTISHEIGHT, type1) else vert_normal = separation(IZ) / vert_normalization(VERTISHEIGHT) - !print*,'JDL standard norm = ',vert_normalization(VERTISHEIGHT) endif !--- JDL End Addition for Vertical Normalization @@ -1691,7 +1689,7 @@ subroutine find_nearest(gc,base_loc,loc_list,nearest, rc) !write(0,*) 'z: ', start_z, end_z ! JDL WARNING STATEMENT IN CASE YOU HAPPEN TO CALL GET_DIST() FROM THIS SUBROUTINE ! YOU SHOULDN'T BE USING THIS SUBROUTINE -print*,'JDL IF YOUR HERE MAKE SURE TO WORK ON FIND_NEAREST SUBROUTINE' +print*,'WARNING - IF YOUR HERE MAKE SURE TO WORK ON FIND_NEAREST SUBROUTINE' ! Next, loop through each box that is close to this box do i = start_x, end_x diff --git a/models/cm1/model_mod.f90 b/models/cm1/model_mod.f90 index d2d1416416..b7d89ffd47 100644 --- a/models/cm1/model_mod.f90 +++ b/models/cm1/model_mod.f90 @@ -894,20 +894,6 @@ subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, Z(:) = bilinear_interpolation(nlevs, obs_loc_array(1), obs_loc_array(2), & x_val(1), x_val(2), y_val(1), y_val(2), Q11, Q12, Q21, Q22) -do level = 1, nlevs - if (isnan(Z(level))) then - print*,'JDL x_val(1)',x_val(1) - print*,'JDL x_val(2)',x_val(2) - print*,'JDL y_val(1)',y_val(1) - print*,'JDL y_val(2)',y_val(2) - print*,'JDL OBS ARRAY(1)',obs_loc_array(1) - print*,'JDL OBS ARRAY(2)',obs_loc_array(2) - print*,'JDL Q11',Q11(level) - print*,'JDL Q12',Q12(level) - print*,'JDL Q21',Q21(level) - print*,'JDL Q22',Q22(level) - endif -enddo ! Find out which level the point is in: call get_enclosing_coord(obs_loc_array(3), Z, z_ind, z_val) @@ -1087,25 +1073,6 @@ subroutine height_interpolate(obs_loc_array, varid, nlevs, nlevs_shrink, x_ind, Z(:) = bilinear_interpolation(nlevs, obs_loc_array(1), obs_loc_array(2), & x_val(1), x_val(2), y_val(1), y_val(2), Q11, Q12, Q21, Q22) - -!do level = 1, nlevs -xloop: do level = 1, nlevs - if (isnan(Z(level))) then - print*,'JDL x_val(1)',x_val(1) - print*,'JDL x_val(2)',x_val(2) - print*,'JDL y_val(1)',y_val(1) - print*,'JDL y_val(2)',y_val(2) - print*,'JDL OBS ARRAY(1)',obs_loc_array(1) - print*,'JDL OBS ARRAY(2)',obs_loc_array(2) - print*,'JDL Q11',Q11(level) - print*,'JDL Q12',Q12(level) - print*,'JDL Q21',Q21(level) - print*,'JDL Q22',Q22(level) - print*,'**********' - exit - endif -enddo xloop - ! Find out which level the point is in: call get_enclosing_coord(obs_loc_array(3), Z, z_ind, z_val) @@ -1573,18 +1540,6 @@ subroutine get_enclosing_coord(x, xcoords, ind, val) enddo xloop if (ind(1) == -999) then - - print*,'JDL X = ',x - print*,'JDL COORD SIZE =',size(xcoords) - print*,'JDL LOC(4) = ',xcoords(4) - !do i = 2, size(xcoords) - ! if(isnan(xcoords(i))) THEN - ! print*,'JDL THIS IS A NAN' - ! endif - - ! print*,'JDL coord = ',xcoords(i) - !enddo - call error_handler(E_ERR, 'get_enclosing_coord', 'off the grid, unexpected 3', & source, revision, revdate) endif diff --git a/observations/forward_operators/obs_def_dew_point_mod.f90 b/observations/forward_operators/obs_def_dew_point_mod.f90 index 668aea7968..4635035e14 100644 --- a/observations/forward_operators/obs_def_dew_point_mod.f90 +++ b/observations/forward_operators/obs_def_dew_point_mod.f90 @@ -209,15 +209,6 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist log_term = log(e_mb / 6.112_r8) td = t_kelvin + 243.5_r8 * log_term / (17.67_r8 - log_term) - !print *,'********************************' - !print *,'e_mb = ',e_mb(1) - !print *,'p_mb = ',p_mb(1) - !print *,'qv =',qv(1) - !print *,'td = ',td(1) - !print *,'t_kelvin = ',t_kelvin - !print *,'***********' - - elsewhere td = missing_r8 end where From 440e635b2cd94e0a22c603a391192153efe149cd Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Wed, 23 Feb 2022 14:05:56 +0000 Subject: [PATCH 08/15] corrected a bug to perform vertical localization in mixed boundary condition cases --- .../location/threed_cartesian/location_mod.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index dde0346e45..be38b7822e 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -443,7 +443,7 @@ function get_dist(loc1, loc2, type1, kind2) ! everything below deals with the periodic options if (debug > 0) write(0,*) 'periodic distance: ' -square_dist = dist_3d_sq(diff) +square_dist = dist_3d_sq(diff,type1=type1) if (debug > 0) write(0,*) 'non-periodic 0 sq_dist, dist: ', square_dist, sqrt(square_dist) ! ok, not simple. @@ -692,7 +692,7 @@ function dist_3d(separation,type1) result(val) !--- JDL Addition to Account for Vertical Normalization if (allocated(per_type_vert_norm)) then if (.not. present(type1)) then - write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' + write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization 3d' call error_handler(E_MSG, 'get_dist', msgstring, source) endif vert_normal = separation(IZ)/per_type_vert_norm(VERTISHEIGHT, type1) @@ -729,7 +729,7 @@ function dist_3d_sq(separation,type1) result(val) !--- JDL Addition to Account for Vertical Normalization if (allocated(per_type_vert_norm)) then if (.not. present(type1)) then - write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization' + write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization 3d sq' call error_handler(E_MSG, 'get_dist', msgstring, source) endif vert_normal = separation(IZ) / per_type_vert_norm(VERTISHEIGHT, type1) From 8b5da4d55c58b9c74edef96d8f21103b089361f1 Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Wed, 23 Feb 2022 16:01:00 +0000 Subject: [PATCH 09/15] added additional documentation on the updates made to the cm1 code --- models/cm1/readme.rst | 105 ++++++++++++++++++++++++++++++++++++++ models/cm1/work/input.nml | 8 ++- 2 files changed, 112 insertions(+), 1 deletion(-) diff --git a/models/cm1/readme.rst b/models/cm1/readme.rst index c2e5017296..32b59c6a63 100644 --- a/models/cm1/readme.rst +++ b/models/cm1/readme.rst @@ -69,6 +69,7 @@ specified in the ``¶m11`` namelist. Observation locations are specified in meters relative to the domain origin as defined the ``iorigin`` setting of ``¶m2``. + About Testing CM1 and DART -------------------------- @@ -439,3 +440,107 @@ Description of each namelist entry 'air_temp' , 'QTY_TEMPERATURE' , 'NULL', 'NULL', 'UPDATE', 'spec_hum' , 'QTY_SPECIFIC_HUMIDITY' , 0.0000, 'NULL', 'UPDATE', + +Some other areas of interest in input.nml will also be discussed + +.. &assim_tools_nml:: + +The CM1 DART code was updated so that horizontal and vertical localization radii +can be defined for each assimilated observation type. The horizontal and vertical localization +radius is first defined by the ``cutoff`` variable in the &assim_tools_nml +code block of input.nml. However, you can use ``special_localization_obs_types`` and +``special_localization_cutoffs`` variables to define the observation type and corresponding +half localization radius, respectively. + +horizontal localization radius = 2 * special_localization_cutoffs[ob_type] + +An example of the code block is provided below. + + &assim_tools_nml + adaptive_localization_threshold = -1 + cutoff = 15000.0 + filter_kind = 1 + print_every_nth_obs = 100 + rectangular_quadrature = .true. + sampling_error_correction = .false. + sort_obs_inc = .false. + spread_restoration = .false. + gaussian_likelihood_tails = .false. + 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. ! Horizontal Localization radii (100km, 100km, 12km, 12km) + / + + + +.. &location_nml:: + +You can also define the vertical localization radius for each observation +type in this code block. ``special_vert_normalization_obs_types`` is +where you define the observation type and ``special_vert_normalization_heights`` is +where you define the normalization factor in the vertical. + +vertical localization radius = 2 * special_localization_cutoffs[ob_type]* special_vert_normalization_heights[ob_type] + +You can also define periodic boundary conditions in the &location_nml code block. Do not change ``nx``, ``ny``, and ``nz`` +they are defined search radii used by DART (not the domain size itself). + +An example of the &location_nml is provided below. + + &location_nml + x_is_periodic = .false. + y_is_periodic = .false. + z_is_periodic = .false. + min_x_for_periodic = 0.0 + max_x_for_periodic = 200000.0 + min_y_for_periodic = 0.0 + max_y_for_periodic = 200000.0 + min_z_for_periodic = 0.0 + max_z_for_periodic = 1.0 + compare_to_correct = .false. + output_box_info = .false. + print_box_level = 0 + nx = 10 + ny = 10 + nz = 10 + vert_normalization_height = 1 + special_vert_normalization_obs_types = 'RADIOSONDE_U_WIND_COMPONENT','RADIOSONDE_V_WIND_COMPONENT', 'RADAR_REFLECTIVITY', 'DOPPLER_RADIAL_VELOCITY' + special_vert_normalization_heights = 0.04, 0.04, 0.33, 0.33 + / + +.. &obs_def_radar_mod_nml:: +The block defines how radar observations are assimilated +e.g., what forward operator is used for reflectivity, radial velocity, +microphysical information, observation cutoff values. + +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. + +A sample of the &obs_def_radar_mod_nml code block is provided below: + + &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 + / diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index 39a3f55157..008b43dcd9 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -169,7 +169,7 @@ &location_nml x_is_periodic = .false. - y_is_periodic = .true. + y_is_periodic = .false. z_is_periodic = .false. min_x_for_periodic = 0.0 max_x_for_periodic = 200000.0 @@ -235,6 +235,12 @@ ! If performing vertical conversion, try ! 'distribute_mean = .false.' +! 'cutoff' initially defines the localization radius (horizontal and vertical +! for all observations. If you want to specify the localization +! radius for each observation type you can list the observation types in 'special_localization_obs_types' +! 'special_localization_cutoffs' defines the localization half-radius for each +! observation type + &assim_tools_nml adaptive_localization_threshold = -1 cutoff = 15000.0 From 56f01d87abb4c917a557953581c0e1efc962872d Mon Sep 17 00:00:00 2001 From: jonlabriola Date: Wed, 23 Feb 2022 16:14:07 +0000 Subject: [PATCH 10/15] Further improving documentations --- models/cm1/readme.rst | 6 +++++- models/cm1/work/input.nml | 8 +++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/models/cm1/readme.rst b/models/cm1/readme.rst index 32b59c6a63..365987bd8e 100644 --- a/models/cm1/readme.rst +++ b/models/cm1/readme.rst @@ -480,9 +480,13 @@ An example of the code block is provided below. You can also define the vertical localization radius for each observation type in this code block. ``special_vert_normalization_obs_types`` is where you define the observation type and ``special_vert_normalization_heights`` is -where you define the normalization factor in the vertical. +where you define the normalization factor in the vertical. If you just one set one normalization for all +observation types use ``vert_normalization_height`` (not typically recommended since horizontal localization +radii often switch between obs types). vertical localization radius = 2 * special_localization_cutoffs[ob_type]* special_vert_normalization_heights[ob_type] +or if no obs types are defined... +vertical localization radius = 2 * cutoff * vert_normalization_height You can also define periodic boundary conditions in the &location_nml code block. Do not change ``nx``, ``ny``, and ``nz`` they are defined search radii used by DART (not the domain size itself). diff --git a/models/cm1/work/input.nml b/models/cm1/work/input.nml index 008b43dcd9..1c734fdc14 100644 --- a/models/cm1/work/input.nml +++ b/models/cm1/work/input.nml @@ -160,12 +160,18 @@ ! Define the observation type in 'special_vert_normalization_obs_types' and define the vertical ! radius in 'special_vert_normalization_heights'. ! -! The vertical localization radius is defined as '2*special_vert_normalization_heights x special_localization_cutoffs' +! The vertical localization radius is defined as '2 x special_vert_normalization_heights x special_localization_cutoffs' ! For example if... special_vert_normalization_heights = 50000 m (for RADIOSONDE) ! special_localization_cutoffs = 0.04 (for RADIOSONDE) ! ! the horizontal localization radius will be 100000 m ! the vertical locaalization radius will be 4000 m +! +! If the observation type you are assimilating is not defined then CM1 will use vert_normalization_height to +! define the vertical localization radius. i.e., vertical localization = 2 x cutoff x vert_normalization_height +! +! You should also define boundary condition information in this code block. Do not change nx, ny, nor nz. +! These are specific search radii used by DART, not the forecast domain size. &location_nml x_is_periodic = .false. From 452e41931ebc7a592287b13aab28e2ed84580513 Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Wed, 23 Feb 2022 15:28:56 -0700 Subject: [PATCH 11/15] revert whitespace changes to main --- models/cm1/model_mod.f90 | 2 ++ observations/forward_operators/obs_def_dew_point_mod.f90 | 6 +++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/models/cm1/model_mod.f90 b/models/cm1/model_mod.f90 index b7d89ffd47..8d407350bc 100644 --- a/models/cm1/model_mod.f90 +++ b/models/cm1/model_mod.f90 @@ -877,6 +877,7 @@ subroutine height_interpolate_s_grid(obs_loc_array, varid, nlevs, nlevs_shrink, Q12(level) = zhalf(x_ind(1), y_ind(2), level) Q21(level) = zhalf(x_ind(2), y_ind(1), level) Q22(level) = zhalf(x_ind(2), y_ind(2), level) + enddo else ! on z half levels @@ -1073,6 +1074,7 @@ subroutine height_interpolate(obs_loc_array, varid, nlevs, nlevs_shrink, x_ind, Z(:) = bilinear_interpolation(nlevs, obs_loc_array(1), obs_loc_array(2), & x_val(1), x_val(2), y_val(1), y_val(2), Q11, Q12, Q21, Q22) + ! Find out which level the point is in: call get_enclosing_coord(obs_loc_array(3), Z, z_ind, z_val) diff --git a/observations/forward_operators/obs_def_dew_point_mod.f90 b/observations/forward_operators/obs_def_dew_point_mod.f90 index 4635035e14..8d50a7e32f 100644 --- a/observations/forward_operators/obs_def_dew_point_mod.f90 +++ b/observations/forward_operators/obs_def_dew_point_mod.f90 @@ -172,6 +172,8 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist call interpolate(state_handle, ens_size,location, QTY_VAPOR_MIXING_RATIO, qv, qv_istatus) call track_status(ens_size, qv_istatus, td, istatus, return_now) if (return_now) return + + where (qv < 0.0_r8 .or. qv >= 1.0_r8) istatus = 1 td = missing_r8 @@ -202,9 +204,6 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist ! The following expression can fail numerically for dewpoints very close to 0 C !td = t_kelvin + (243.5_r8 / ((17.67_r8 / log(e_mb/6.112_r8)) - 1.0_r8) ) - - - ! A numerically robust formula that avoids the failure near dewpoints of 0 C log_term = log(e_mb / 6.112_r8) td = t_kelvin + 243.5_r8 * log_term / (17.67_r8 - log_term) @@ -212,6 +211,7 @@ subroutine get_expected_dew_point(state_handle, ens_size, location, key, td, ist elsewhere td = missing_r8 end where + end subroutine get_expected_dew_point !---------------------------------------------------------------------------- From f900c4d50bc39cb0b519ea7f45360dfa487bc8ad Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Thu, 31 Mar 2022 14:53:10 -0600 Subject: [PATCH 12/15] comment tidy. removed unused variable location_vertical_localization_coord --- .../threed_cartesian/location_mod.f90 | 67 ++----------------- 1 file changed, 7 insertions(+), 60 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index be38b7822e..7b0797f568 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -45,10 +45,9 @@ module location_mod character(len = 129), parameter :: LocationStorageOrder = "X Y Z" character(len = 129), parameter :: LocationUnits = "none none none" -! The possible numeric values for the location_type%which_vert component. -! The numeric values are PRIVATE to this module. The parameter names are PUBLIC. +! The numeric value for the vertical component +! Used with CM1 when using normalized vertical localization integer, parameter :: VERTISHEIGHT = 1 ! by height (in meters) -!--- Can eventually add other heights if necesary (JDL) type location_type private @@ -118,7 +117,6 @@ module location_mod character(len = 512) :: errstring character(len = 512) :: msgstring, msgstring1, msgstring2 -! JDL NEW ! Horizontal localization/cutoff values are passed in by the caller. ! The Vertical normalization values are globals; are specified by namelist ! here, and apply to all specific types unless a 'special list' is also specified @@ -126,16 +124,9 @@ module location_mod logical :: has_special_vertical_norms = .false. integer :: num_special_vert_norms = 0 -integer :: location_vertical_localization_coord = 0 -! JDL NEW ! Global storage for vertical distance normalization factors -! The 4 below is for the 4 vertical units (pressure, level, height, -! scale height). undefined and surface don't need vert distances. -! NOTE: Code that uses VERT_TYPE_COUNT depends on pressure, level, -! height, and scale height having actual values between 1 and 4, or -! this code will break. -integer, parameter :: VERT_TYPE_COUNT = 1 ! JDL JUST WORKING ON HEIGHT +integer, parameter :: VERT_TYPE_COUNT = 1 ! only height supported real(r8) :: vert_normalization(VERT_TYPE_COUNT) real(r8), allocatable :: per_type_vert_norm(:,:) ! if doing per-type @@ -190,23 +181,15 @@ module location_mod integer :: nz = 10 -!--- JDL: Additional Items for Vertical Localization +! Additional Items for Vertical Localization !------------------------------------------------------------------------ ! Namelist with default values ! -! horiz_dist_only is currently ignored (JDL) -! -! horiz_dist_only == .true. -> Only the great circle horizontal distance is -! computed in get_dist. -! horiz_dist_only == .false. -> Square root of sum of squared horizontal and -! normalized vertical dist computed in get_dist ! vert_normalization_height -> Number meters that give a distance equivalent ! to one radian in horizontal ! special_vert_normalization_obs_types -> Which obs types to modify the default vert ! normalization values -! special_vert_normalization_pressure -> must give all 4 values for each type listed - -!logical :: horiz_dist_only = .true. ! JDL While DEFAULT - You'll want false +! special_vert_normalization_heights -> value for each obs type real(r8) :: vert_normalization_height = 10000.0_r8 integer, parameter :: MAX_ITEMS = 500 character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS) @@ -217,7 +200,6 @@ module location_mod x_is_periodic, min_x_for_periodic, max_x_for_periodic, & y_is_periodic, min_y_for_periodic, max_y_for_periodic, & z_is_periodic, min_z_for_periodic, max_z_for_periodic, & - !nx, ny, nz, horiz_dist_only, vert_normalization_height, & nx, ny, nz, vert_normalization_height, & special_vert_normalization_obs_types, & special_vert_normalization_heights, debug @@ -302,8 +284,7 @@ subroutine initialize_module if (special_vert_normalization_obs_types(1) /= 'null' .or. & special_vert_normalization_heights(1) /= missing_r8) then - ! FIXME: add code to check for mismatched length lists. are we going to force - ! users to specify all 4 values for any obs type that is not using the defaults? + ! FIXME: add code to check for mismatched length lists. typecount = get_num_types_of_obs() ! ignore function name, this is specific type count allocate(per_type_vert_norm(VERT_TYPE_COUNT, typecount)) @@ -355,34 +336,6 @@ subroutine initialize_module if (allocated(per_type_vert_norm)) then typecount = get_num_types_of_obs() ! ignore function name, this is specific type count do i = 1, typecount - !if ((per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) .or. & - ! (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) .or. & - ! (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) .or. & - ! (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height)) then - - ! write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) - ! call error_handler(E_MSG,'location_mod:',msgstring,source) - ! if (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) then - ! write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', & - ! per_type_vert_norm(VERTISPRESSURE, i) - ! call error_handler(E_MSG,'location_mod:',msgstring,source) - ! endif - ! if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then - ! write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', & - ! per_type_vert_norm(VERTISHEIGHT, i) - ! call error_handler(E_MSG,'location_mod:',msgstring,source) - ! endif - ! if (per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) then - ! write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', & - ! per_type_vert_norm(VERTISLEVEL, i) - ! call error_handler(E_MSG,'location_mod:',msgstring,source) - ! endif - ! if (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height) then - ! write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', & - ! per_type_vert_norm(VERTISSCALEHEIGHT, i) - ! call error_handler(E_MSG,'location_mod:',msgstring,source) - ! endif - ! JDL - Strip Down to Just VERTISHEIGHT Right Now (CM1 is in meter coordinates only right now) if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i)) call error_handler(E_MSG,'location_mod:',msgstring,source) @@ -394,7 +347,6 @@ subroutine initialize_module endif enddo endif -!endif end subroutine initialize_module @@ -403,7 +355,7 @@ end subroutine initialize_module function get_dist(loc1, loc2, type1, kind2) -! JDL I don't think I need to include the no_vert variable, but keep this +! @todo JDL I don't think I need to include the no_vert variable, but keep this ! in mind in case problems arise later on. Probbably also don't need kind2 since ! we are just looking at heights @@ -423,7 +375,6 @@ function get_dist(loc1, loc2, type1, kind2) real(r8) :: diff(3), next_diff(3) logical :: below_L1(3), below_L2(3) real(r8) :: square_dist, this_dist -!logical :: comp_h_only ! Only comput the horizontal component if ( .not. module_initialized ) call initialize_module @@ -681,8 +632,6 @@ end function get_dist ! return the 3d distance given the separation along each axis -!pure function dist_3d(separation) -!function dist_3d(separation) result(val) function dist_3d(separation,type1) result(val) real(r8), intent(in) :: separation(3) @@ -717,8 +666,6 @@ end function dist_3d ! return the square of the 3d distance given the separation along each axis ! (saves doing a square root) -!pure function dist_3d_sq(separation) -!function dist_3d_sq(separation) result(val) function dist_3d_sq(separation,type1) result(val) real(r8), intent(in) :: separation(3) From 22cdd0f81608826b8ac089dad4110eb45789747e Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Thu, 31 Mar 2022 15:09:40 -0600 Subject: [PATCH 13/15] doc: fix rst syntax --- models/cm1/readme.rst | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/models/cm1/readme.rst b/models/cm1/readme.rst index 2301a0f85e..58026d048f 100644 --- a/models/cm1/readme.rst +++ b/models/cm1/readme.rst @@ -447,7 +447,8 @@ Description of each namelist entry Some other areas of interest in input.nml will also be discussed -.. &assim_tools_nml:: +&assim_tools_nml +~~~~~~~~~~~~~~~~~~ The CM1 DART code was updated so that horizontal and vertical localization radii can be defined for each assimilated observation type. The horizontal and vertical localization @@ -456,9 +457,12 @@ code block of input.nml. However, you can use ``special_localization_obs_types` ``special_localization_cutoffs`` variables to define the observation type and corresponding half localization radius, respectively. -horizontal localization radius = 2 * special_localization_cutoffs[ob_type] -An example of the code block is provided below. +``horizontal localization radius = 2 * special_localization_cutoffs[ob_type]`` + +An example of the &assim_tools_nml using per-type radii is provided below. + +.. code-block:: text &assim_tools_nml adaptive_localization_threshold = -1 @@ -479,24 +483,30 @@ An example of the code block is provided below. -.. &location_nml:: +&location_nml +~~~~~~~~~~~~~~~ You can also define the vertical localization radius for each observation -type in this code block. ``special_vert_normalization_obs_types`` is +type. ``special_vert_normalization_obs_types`` is where you define the observation type and ``special_vert_normalization_heights`` is where you define the normalization factor in the vertical. If you just one set one normalization for all observation types use ``vert_normalization_height`` (not typically recommended since horizontal localization radii often switch between obs types). -vertical localization radius = 2 * special_localization_cutoffs[ob_type]* special_vert_normalization_heights[ob_type] + +``vertical localization radius = 2 * special_localization_cutoffs[ob_type]* special_vert_normalization_heights[ob_type]`` + or if no obs types are defined... -vertical localization radius = 2 * cutoff * vert_normalization_height + +``vertical localization radius = 2 * cutoff * vert_normalization_height`` You can also define periodic boundary conditions in the &location_nml code block. Do not change ``nx``, ``ny``, and ``nz`` they are defined search radii used by DART (not the domain size itself). An example of the &location_nml is provided below. +.. code-block:: text + &location_nml x_is_periodic = .false. y_is_periodic = .false. @@ -518,8 +528,11 @@ An example of the &location_nml is provided below. special_vert_normalization_heights = 0.04, 0.04, 0.33, 0.33 / -.. &obs_def_radar_mod_nml:: -The block defines how radar observations are assimilated + +&obs_def_radar_mod_nml +~~~~~~~~~~~~~~~~~~~~~~~~ + +The block defines how radar observations are assimilated e.g., what forward operator is used for reflectivity, radial velocity, microphysical information, observation cutoff values. @@ -531,7 +544,9 @@ 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. -A sample of the &obs_def_radar_mod_nml code block is provided below: +A sample &obs_def_radar_mod_nml is provided below: + +.. code-block:: text &obs_def_radar_mod_nml apply_ref_limit_to_obs = .false., From c5769a6296d13ef09aa955129e24c0b53c7f3ca0 Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Fri, 1 Apr 2022 09:02:32 -0600 Subject: [PATCH 14/15] default vert_normalization_height set to 1.0 This keeps the default behaviour for dist_3d remain the same (no change to vertical) --- .../location/threed_cartesian/location_mod.f90 | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/assimilation_code/location/threed_cartesian/location_mod.f90 b/assimilation_code/location/threed_cartesian/location_mod.f90 index 7b0797f568..32c712d165 100644 --- a/assimilation_code/location/threed_cartesian/location_mod.f90 +++ b/assimilation_code/location/threed_cartesian/location_mod.f90 @@ -190,7 +190,7 @@ module location_mod ! special_vert_normalization_obs_types -> Which obs types to modify the default vert ! normalization values ! special_vert_normalization_heights -> value for each obs type -real(r8) :: vert_normalization_height = 10000.0_r8 +real(r8) :: vert_normalization_height = 1.0_r8 integer, parameter :: MAX_ITEMS = 500 character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS) real(r8) :: special_vert_normalization_heights(MAX_ITEMS) @@ -637,8 +637,7 @@ function dist_3d(separation,type1) result(val) real(r8), intent(in) :: separation(3) real(r8) :: val, vert_normal integer, optional, intent(in) :: type1 ! JDL Addition -!--- JDL Update Function to account for Normalization in Vertical -!--- JDL Addition to Account for Vertical Normalization + if (allocated(per_type_vert_norm)) then if (.not. present(type1)) then write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization 3d' @@ -649,13 +648,10 @@ function dist_3d(separation,type1) result(val) vert_normal = separation(IZ)/vert_normalization(VERTISHEIGHT) endif -!--- JDL End Addition for Vertical Normalization - val = sqrt(separation(IX)*separation(IX) + & separation(IY)*separation(IY) + & vert_normal*vert_normal ) - !separation(IZ)*separation(IZ) ) if (debug > 0) write(0,*) 'dist_3d called, distance computed: ', val if (debug > 0) write(0,*) 'XYZ separations: ', separation @@ -672,8 +668,6 @@ function dist_3d_sq(separation,type1) result(val) real(r8) :: val, vert_normal integer, optional, intent(in) :: type1 ! JDL Addition -!--- JDL Update Function to account for Normalization in Vertical -!--- JDL Addition to Account for Vertical Normalization if (allocated(per_type_vert_norm)) then if (.not. present(type1)) then write(msgstring, *) 'obs type required in get_dist`() if doing per-type vertical normalization 3d sq' @@ -683,12 +677,10 @@ function dist_3d_sq(separation,type1) result(val) else vert_normal = separation(IZ) / vert_normalization(VERTISHEIGHT) endif -!--- JDL End Addition for Vertical Normalization val = separation(IX)*separation(IX) + & separation(IY)*separation(IY) + & vert_normal * vert_normal - !separation(IZ)*separation(IZ) if (debug > 0) write(0,*) 'dist_3d_sq called, distance computed: ', val if (debug > 0) write(0,*) 'XYZ separations: ', separation From 032baa625fb228c82357684e437a964199076f84 Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Fri, 1 Apr 2022 09:25:51 -0600 Subject: [PATCH 15/15] bump version and changelog for release --- CHANGELOG.rst | 6 ++++++ conf.py | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index e205edf7eb..10c0af119c 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,12 @@ individual files. The changes are now listed with the most recent at the top. +**April 1 2022 :: Per-obs-type localization for 3D Cartesian location_mod. Tag: v9.16.1** + +- Optional per-obs-type localization for 3D Cartesian location + +*Contributed by Jon Labriola for use with CM1* + **March 31 2022 :: MiTgcm-ocean NBLING. Tag: v9.16.0** - MITgcm-ocean interface updated to Manhattan. diff --git a/conf.py b/conf.py index 0569f38062..db800e96cc 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.16.0' +release = '9.16.1' master_doc = 'README' # -- General configuration ---------------------------------------------------