diff --git a/assimilation_code/modules/observations/atmosphere_quantities_mod.f90 b/assimilation_code/modules/observations/atmosphere_quantities_mod.f90 index f580cd4f0..f01bba771 100644 --- a/assimilation_code/modules/observations/atmosphere_quantities_mod.f90 +++ b/assimilation_code/modules/observations/atmosphere_quantities_mod.f90 @@ -104,6 +104,7 @@ ! ! QTY_2M_SPECIFIC_HUMIDITY ! QTY_2M_TEMPERATURE +! QTY_2M_POTENTIAL_TEMPERATURE ! QTY_10M_U_WIND_COMPONENT ! QTY_10M_V_WIND_COMPONENT ! diff --git a/assimilation_code/modules/observations/default_quantities_mod.f90 b/assimilation_code/modules/observations/default_quantities_mod.f90 index 931e298a9..019ce87dc 100644 --- a/assimilation_code/modules/observations/default_quantities_mod.f90 +++ b/assimilation_code/modules/observations/default_quantities_mod.f90 @@ -11,6 +11,7 @@ ! QTY_2D_PARAMETER ! QTY_2M_SPECIFIC_HUMIDITY ! QTY_2M_TEMPERATURE +! QTY_2M_POTENTIAL_TEMPERATURE ! QTY_2M_VAPOR_MIXING_RATIO ! QTY_3D_PARAMETER ! QTY_ABSOLUTE_HUMIDITY diff --git a/index.rst b/index.rst index bf4596f72..ead1ed223 100644 --- a/index.rst +++ b/index.rst @@ -467,6 +467,7 @@ References models/wrf/readme models/wrf/WRF_DART_utilities/replace_wrf_fields models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess + models/wrf_unified/readme models/utilities/default_model_mod diff --git a/models/README.rst b/models/README.rst index 5fd8a609d..67f89f9ca 100644 --- a/models/README.rst +++ b/models/README.rst @@ -50,5 +50,6 @@ DART supported models: - :doc:`tiegcm/readme` - :doc:`wrf_hydro/readme` - :doc:`wrf/readme` +- :doc:`WRF CHEM ` If you are interested in creating a DART interface for a new model, see :ref:`Using new models` and :ref:`Porting new models`. diff --git a/models/wrf_unified/model_mod.f90 b/models/wrf_unified/model_mod.f90 new file mode 100644 index 000000000..f4b1ea2e1 --- /dev/null +++ b/models/wrf_unified/model_mod.f90 @@ -0,0 +1,4079 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +module model_mod + +use types_mod, only : r8, i8, MISSING_R8, digits12, & + gas_constant, gas_constant_v, ps0, gravity + +use time_manager_mod, only : time_type, set_time, GREGORIAN, set_date, & + set_calendar_type + +use location_mod, only : location_type, get_close_type, & + set_location, set_location_missing, & + set_vertical_localization_coord, set_vertical, & + VERTISHEIGHT, VERTISLEVEL, VERTISPRESSURE, & + VERTISSURFACE, VERTISUNDEF, VERTISSCALEHEIGHT, & + loc_get_close => get_close, get_location, & + query_location, is_vertical, vertical_localization_on, & + get_dist + +use utilities_mod, only : register_module, error_handler, & + E_ERR, E_MSG, & + nmlfileunit, do_output, do_nml_file, do_nml_term, & + find_namelist_in_file, check_namelist_read, & + to_upper + +use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & + nc_add_global_creation_time, & + nc_begin_define_mode, nc_end_define_mode, & + NF90_MAX_NAME, nc_get_variable_size, & + nc_get_variable, nc_close_file, nc_check, & + nc_open_file_readonly, nc_get_variable_size, & + nc_get_global_attribute, nc_get_dimension_size + +use state_structure_mod, only : add_domain, get_domain_size, get_model_variable_indices, & + get_dim_name, get_num_dims, get_dart_vector_index, & + get_varid_from_kind, get_varid_from_varname, state_structure_info + +use distributed_state_mod, only : get_state_array, get_state + +use obs_kind_mod, only : get_index_for_quantity, & + get_name_for_quantity, & + QTY_U_WIND_COMPONENT, & + QTY_v_WIND_COMPONENT, & + QTY_10M_U_WIND_COMPONENT, & + QTY_10M_V_WIND_COMPONENT, & + QTY_DENSITY, & + QTY_GEOPOTENTIAL_HEIGHT, & + QTY_PRESSURE, & + QTY_SURFACE_TYPE, & + QTY_SURFACE_ELEVATION, & + QTY_LANDMASK, & + QTY_SURFACE_PRESSURE, & + QTY_VAPOR_MIXING_RATIO, & + QTY_TEMPERATURE, & + QTY_POTENTIAL_TEMPERATURE, & + QTY_DENSITY, & + QTY_VERTICAL_VELOCITY, & + QTY_SPECIFIC_HUMIDITY, & + QTY_SURFACE_PRESSURE, & + QTY_VORTEX_LAT, & + QTY_VORTEX_LON, & + QTY_VORTEX_PMIN,QTY_VORTEX_WMAX, & + QTY_SKIN_TEMPERATURE, & + QTY_SURFACE_TYPE, & + QTY_2M_TEMPERATURE, & + QTY_2M_POTENTIAL_TEMPERATURE, & + QTY_2M_SPECIFIC_HUMIDITY, & + QTY_RAINWATER_MIXING_RATIO, & + QTY_GRAUPEL_MIXING_RATIO, & + QTY_HAIL_MIXING_RATIO, & + QTY_SNOW_MIXING_RATIO, & + QTY_ICE_MIXING_RATIO, & + QTY_CLOUDWATER_MIXING_RATIO, & + QTY_DROPLET_NUMBER_CONCENTR, & + QTY_ICE_NUMBER_CONCENTRATION, & + QTY_SNOW_NUMBER_CONCENTR, & + QTY_RAIN_NUMBER_CONCENTR, & + QTY_GRAUPEL_NUMBER_CONCENTR, & + QTY_HAIL_NUMBER_CONCENTR, & + QTY_SOIL_TEMPERATURE, & + QTY_SOIL_MOISTURE, & + QTY_SOIL_LIQUID_WATER, & + QTY_CLOUDWATER_DE, & + QTY_CLOUD_ICE_DE + +use ensemble_manager_mod, only : ensemble_type + +use default_model_mod, only : write_model_time, & + init_time => fail_init_time, & + init_conditions => fail_init_conditions, & + adv_1step + +use map_utils, only : latlon_to_ij, & + proj_info, & + map_set, & + map_init, & + gridwind_to_truewind, & + PROJ_LATLON, & + PROJ_LC, & + PROJ_PS, & + PROJ_PS_WGS84, & + PROJ_MERC, & + PROJ_GAUSS, & + PROJ_CYL, & + PROJ_CASSINI, & + PROJ_ALBERS_NAD83, & + PROJ_ROTLL + +use netcdf ! no get_char in netcdf_utilities_mod + +implicit none + +private + +! routines required by DART code - will be called from filter and other +! DART executables. +public :: get_model_size, & + get_state_meta_data, & + model_interpolate, & + end_model, & + static_init_model, & + nc_write_model_atts, & + get_close_obs, & + get_close_state, & + pert_model_copies, & + convert_vertical_obs, & + convert_vertical_state, & + read_model_time, & + adv_1step, & + init_time, & + init_conditions, & + shortest_time_between_assimilations, & + write_model_time + +! for wrf_dart_obs_preprocess.f90 +public :: get_domain_info + + +! module variables +character(len=256), parameter :: source = "wrf/model_mod.f90" +logical :: module_initialized = .false. + +integer, parameter :: MAX_STATE_VARIABLES = 100 +integer, parameter :: NUM_STATE_TABLE_COLUMNS = 4 +integer, parameter :: NUM_BOUNDS_TABLE_COLUMNS = 3 + +integer, allocatable :: wrf_dom(:) ! This needs a better name, it is the id from add_domain + ! for each wrf_domain added to the state + +!-- Namelist with default values -- +character(len=NF90_MAX_NAME) :: wrf_state_variables(MAX_STATE_VARIABLES*NUM_STATE_TABLE_COLUMNS) = 'NULL' +character(len=NF90_MAX_NAME) :: wrf_state_bounds(NUM_BOUNDS_TABLE_COLUMNS,MAX_STATE_VARIABLES) = 'NULL' +integer :: num_domains = 1 +logical :: chemistry_separate_file = .false. +character(len=NF90_MAX_NAME) :: chem_state_variables(MAX_STATE_VARIABLES*NUM_STATE_TABLE_COLUMNS) = 'NULL' +character(len=NF90_MAX_NAME) :: chem_state_bounds(NUM_BOUNDS_TABLE_COLUMNS,MAX_STATE_VARIABLES) = 'NULL' +integer :: calendar_type = GREGORIAN +integer :: assimilation_period_seconds = 21600 +! Max height a surface obs can be away from the actual model surface +! and still be accepted (in meters) +real (kind=r8) :: sfc_elev_max_diff = -1.0_r8 ! could be something like 200.0_r8 + +integer :: vert_localization_coord = VERTISHEIGHT + +! Allow observations above the surface but below the lowest sigma level. +logical :: allow_obs_below_vol = .false. + +! Do the interpolation of pressure values only after taking the log (.true.) +! vs doing a linear interpolation directly in pressure units (.false.) +! HK @todo these are not in the namelist in main:wrf. Should they be? +logical :: log_vert_interp = .true. +logical :: log_horz_interpM = .false. +logical :: log_horz_interpQ = .false. + +! wrf options, apply to domain 1 only +logical :: polar = .false. +logical :: periodic_x = .false. +logical :: periodic_y = .false. + +logical :: allow_perturbed_ics = .false. +!------------------------------- + +logical, parameter :: restrict_polar = .false. !HK what is this for? Hardcoded in original code + +namelist /model_nml/ & +wrf_state_variables, & +wrf_state_bounds, & +chemistry_separate_file, & +chem_state_variables, & +chem_state_bounds, & +num_domains, & +calendar_type, & +assimilation_period_seconds, & +sfc_elev_max_diff, & +vert_localization_coord, & +allow_perturbed_ics, & +allow_obs_below_vol, & +log_horz_interpM, & +log_horz_interpQ + +type grid_ll + integer :: map_proj ! MAP_PROJ in wrf netcdf file + type(proj_info) :: proj ! wrf map projection structure + real(r8) :: dx + real(r8) :: truelat1, truelat2, stand_lon + real(r8), dimension(:,:), allocatable :: latitude, latitude_u, latitude_v + real(r8), dimension(:,:), allocatable :: longitude, longitude_u, longitude_v + integer :: we, sn ! west-east, south-north number of grid points + integer :: wes, sns ! west-east staggered, south-north staggered number of grid points + integer :: bt ! bottom-top number of grid points + integer :: bts ! staggered bottom-top number of grid points + real(r8) :: dt ! time step + + ! wrf options, apply to domain 1 only. + logical :: polar = .false. + logical :: periodic_x = .false. + logical :: periodic_y = .false. + +end type grid_ll + +type static_data + real(r8), allocatable :: phb(:,:,:) ! base-state geopotential + real(r8), allocatable :: mub(:,:) ! base state dry air mass in column + real(r8), allocatable :: hgt(:,:) ! Terrain Height + real(r8), allocatable :: dnw(:) ! d(eta) values between full (w) level + real(r8), allocatable :: land(:,:) ! land mask (1 for land, 2 for water) + real(r8), allocatable :: zs(:) ! depths of center of soil layers + real(r8) :: p_top ! Pressure top of the model +end type static_data + +! need grid for each domain +type(grid_ll), allocatable :: grid(:) +type(static_data), allocatable :: stat_dat(:) + +! Physical constants +real(r8), parameter :: rd_over_rv = gas_constant / gas_constant_v +real(r8), parameter :: cpovcv = 1.4_r8 ! cp / (cp - gas_constant) +real(r8), parameter :: ts0 = 300.0_r8 ! Base potential temperature for all levels. +real(r8), parameter :: kappa = 2.0_r8/7.0_r8 ! gas_constant / cp + +contains + +!------------------------------------------------------------------ +subroutine static_init_model() + +integer :: iunit, io + +character(len=NF90_MAX_NAME) :: varname(MAX_STATE_VARIABLES) +integer :: state_qty(MAX_STATE_VARIABLES) +logical :: update_var(MAX_STATE_VARIABLES) +real(r8) :: bounds(MAX_STATE_VARIABLES, 2) ! lower, upper +real(r8) :: lower(MAX_STATE_VARIABLES), upper(MAX_STATE_VARIABLES) +character(len=9) :: in_domain(MAX_STATE_VARIABLES) ! assumes <=9 or 999 + +integer :: nfields, n +logical, allocatable :: domain_mask(:) +integer :: i, field ! loop indices +character (len=1) :: idom ! assumes <=9 + +character(len=256) :: errstring + +module_initialized = .true. + +call find_namelist_in_file("input.nml", "model_nml", iunit) +read(iunit, nml = model_nml, iostat = io) +call check_namelist_read(iunit, io, "model_nml") + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=model_nml) +if (do_nml_term()) write( * , nml=model_nml) + +call set_calendar_type(calendar_type) + +allocate(wrf_dom(num_state_domains()), grid(num_domains), stat_dat(num_domains)) + +call verify_state_variables(nfields, varname, state_qty, update_var, in_domain) +allocate(domain_mask(nfields)) +bounds(:,:) = MISSING_R8 ! default to no clamping + +do i = 1, num_domains + + do field = 1, nfields + domain_mask(field) = variable_is_on_domain(in_domain(field), i) + call get_variable_bounds(varname(field),lower(field),upper(field)) + end do + + n = count(domain_mask) + bounds(1:n, 1) = pack(lower(1:nfields),domain_mask) + bounds(1:n, 2) = pack(upper(1:nfields),domain_mask) + + write( idom , '(I1)') i + wrf_dom(i) = add_domain('wrfinput_d0'//idom, & + num_vars = n, & + var_names = pack(varname(1:nfields), domain_mask), & + kind_list = pack(state_qty(1:nfields), domain_mask), & + clamp_vals = bounds(1:n,:), & + update_list = pack(update_var(1:nfields), domain_mask) ) + +enddo + +deallocate(domain_mask) + +if (chemistry_separate_file) then + call verify_chem_state_variables(nfields, varname, state_qty, update_var, in_domain) + allocate(domain_mask(nfields)) + bounds(:,:) = MISSING_R8 ! default to no clamping + + do i = 1, num_domains + + do field = 1, nfields + domain_mask(field) = variable_is_on_domain(in_domain(field), i) + call get_chem_variable_bounds(varname(field),lower(field),upper(field)) + end do + + n = count(domain_mask) + bounds(1:n, 1) = pack(lower(1:nfields),domain_mask) + bounds(1:n, 2) = pack(upper(1:nfields),domain_mask) + + write( idom , '(I1)') i + wrf_dom(num_domains + i) = add_domain('wrfchem_d0'//idom, & + num_vars = n, & + var_names = pack(varname(1:nfields), domain_mask), & + kind_list = pack(state_qty(1:nfields), domain_mask), & + clamp_vals = bounds(1:n,:), & + update_list = pack(update_var(1:nfields), domain_mask) ) + + enddo + + deallocate(domain_mask) + +endif + +call read_grid() +call read_static_data() + +call set_vertical_localization_coord(vert_localization_coord) + +write(errstring,*) ' wrf model size is ', get_model_size() +call error_handler(E_MSG, 'static_init_model', errstring) + + +end subroutine static_init_model + +!------------------------------------------------------------------ +! Returns the number of items in the state vector as an i8 integer. +function get_model_size() + +integer(i8) :: get_model_size +integer :: i + +if ( .not. module_initialized ) call static_init_model + +get_model_size = 0 +do i = 1, num_state_domains() + get_model_size = get_model_size + get_domain_size(wrf_dom(i)) +enddo + +end function get_model_size +!------------------------------------------------------------------ + +function num_state_domains() + +integer :: num_state_domains + +if ( .not. module_initialized ) call static_init_model + +if (chemistry_separate_file) then + num_state_domains = num_domains * 2 +else + num_state_domains = num_domains +endif + +end function num_state_domains + +!------------------------------------------------------------------ +subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_obs, istatus) + + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: qty_in +real(r8), intent(out) :: expected_obs(ens_size) ! array of interpolated values +integer, intent(out) :: istatus(ens_size) + +integer, parameter :: FAILED_BOUNDS_CHECK = 44 +integer, parameter :: CANNOT_INTERPOLATE_QTY = 55 +integer, parameter :: POLAR_RESTRICTED = 10 ! polar observation while restrict_polar = .true. +integer, parameter :: NOT_IN_ANY_DOMAIN = 11 +integer, parameter :: VERTICAL_LOCATION_FAIL = 66 +integer, parameter :: SURFACE_ELEVATION_DIFF_FAIL = 77 +real(r8) :: lon_lat_vert(3) +real(r8) :: xloc, yloc ! WRF i,j in the grid +integer :: i, j ! grid +real(r8) :: dx, dxm, dy, dym ! grid fractions +integer :: ll(2), ul(2), lr(2), ur(2) !(x,y) of four corners +integer :: rc ! return code getCorners +integer :: id +integer :: k(ens_size) ! level +integer :: which_vert ! vertical coordinate of the observation +real(r8) :: zloc(ens_size) ! vertical location of the obs for each ens member +real(r8) :: fld_k1(ens_size), fld_k2(ens_size) ! value at level k and k+1 +real(r8) :: fld1(ens_size), fld2(ens_size) ! value at level k and k+1 non-negative if required +logical :: fail +integer :: qty + +if ( .not. module_initialized ) call static_init_model + +expected_obs(:) = MISSING_R8 + +istatus(:) = 1 + +which_vert = nint(query_location(location)) +lon_lat_vert = get_location(location) +call get_domain_info(lon_lat_vert(1),lon_lat_vert(2),id,xloc,yloc) ! mass points + +if (id == 0) then + istatus(:) = NOT_IN_ANY_DOMAIN + return +endif + +qty = update_qty_if_location_is_surface(qty_in, location) + +if (.not. able_to_interpolate_qty(id, qty) ) then + istatus(:) = CANNOT_INTERPOLATE_QTY + return +endif + +! horizontal location mass point +call toGrid(xloc,i,dx,dxm) +call toGrid(yloc,j,dy,dym) + +if ( .not. within_bounds_horizontal(i, j, id, qty) ) then + istatus(:) = FAILED_BOUNDS_CHECK + return +endif + +call getCorners(i, j, id, qty, ll, ul, lr, ur, rc) + +! vertical location +if (surface_qty(qty) .or. which_vert==VERTISSURFACE) then + zloc(:) = 1.0_r8 + k = 1 + fail = surface_elevation_difference_too_big(id, ll, ul, lr, ur, dx, dy, dxm, dym, lon_lat_vert(3), sfc_elev_max_diff) + if (fail) then + istatus(:) = SURFACE_ELEVATION_DIFF_FAIL + return + endif +else + call get_level_below_obs(which_vert, id, lon_lat_vert, ens_size, state_handle, ll, ul, lr, ur, dx, dy, dxm, dym, k, zloc, fail) + if (fail) then + istatus(:) = VERTICAL_LOCATION_FAIL + return + endif +endif + + +select case (qty) + case (QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT ) + fld_k1 = wind_interpolate(ens_size, state_handle, qty, id, k, xloc, yloc, i, j, dxm, dx, dy, dym, lon_lat_vert(1)) + fld_k2 = wind_interpolate(ens_size, state_handle, qty, id, k+1, xloc, yloc, i, j, dxm, dx, dy, dym, lon_lat_vert(1)) + case (QTY_10M_U_WIND_COMPONENT, QTY_10M_V_WIND_COMPONENT) + fld_k1(:) = wind_interpolate(ens_size, state_handle, qty, id, k, xloc, yloc, i, j, dxm, dx, dy, dym, lon_lat_vert(1)) + case (QTY_TEMPERATURE) + fld_k1 = temperature_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2 = temperature_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_POTENTIAL_TEMPERATURE) + fld_k1 = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + ts0 + fld_k2 = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + ts0 + case (QTY_DENSITY) + fld_k1 = density_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2 = density_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_VERTICAL_VELOCITY) + zloc(:) = zloc(:) + 0.5_r8 ! Adjust zloc for staggered + k(:) = max(1,int(zloc(:))) ! Adjust corresponding level k + fld_k1(:) = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2(:) = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_SPECIFIC_HUMIDITY) + fld_k1 = specific_humidity_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2 = specific_humidity_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_PRESSURE) + fld_k1 = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2 = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_SURFACE_PRESSURE) + fld_k1 = surface_pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, dxm, dx, dy, dym) + case (QTY_VORTEX_LAT, QTY_VORTEX_LON, QTY_VORTEX_PMIN, QTY_VORTEX_WMAX) + call vortex() + case (QTY_GEOPOTENTIAL_HEIGHT) + zloc(:) = zloc(:) + 0.5_r8 ! Adjust zloc for staggered + k(:) = max(1,int(zloc(:))) ! Adjust corresponding level k + fld_k1 = geopotential_height_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2 = geopotential_height_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + case (QTY_SURFACE_ELEVATION) + fld_k1 = surface_elevation_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + case (QTY_LANDMASK) + fld_k1 = landmask_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + case (QTY_SURFACE_TYPE) + fld_k1(:) = surface_type_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + case (QTY_SKIN_TEMPERATURE, QTY_2M_TEMPERATURE, QTY_2M_SPECIFIC_HUMIDITY) + fld_k1(:) = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + case default ! simple interpolation + ! HK todo 2D variables + fld_k1(:) = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + fld_k2(:) = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) +end select + +if (surface_qty(qty) .or. which_vert == VERTISSURFACE) then + expected_obs(:) = fld_k1(:) +else + fld1 = force_non_negative_if_required(ens_size, qty, fld_k1) + fld2 = force_non_negative_if_required(ens_size, qty, fld_k2) + expected_obs(:) = vertical_interpolation(ens_size, zloc, fld1, fld2) +endif + +call update_units_if_required(qty, expected_obs) + +istatus(:) = 0 + +end subroutine model_interpolate + +!------------------------------------------------------------------ +! Returns the smallest increment in time that the model is capable +! of advancing the state in a given implementation, or the shortest +! time you want the model to advance between assimilations. +function shortest_time_between_assimilations() + +type(time_type) :: shortest_time_between_assimilations +integer :: model_dt + +if ( .not. module_initialized ) call static_init_model + +model_dt = nint(grid(1)%dt) ! model time step in seconds lowest res domain + +if (assimilation_period_seconds < model_dt ) then + shortest_time_between_assimilations = set_time(model_dt) +else + shortest_time_between_assimilations = set_time(assimilation_period_seconds) +endif + +end function shortest_time_between_assimilations + +!------------------------------------------------------------------ +! Given an integer index into the state vector, returns the +! associated location and optionally the physical quantity. +subroutine get_state_meta_data(index_in, location, qty_out) + +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, optional, intent(out) :: qty_out + +integer :: i, j, k, id, var_id, state_id, qty + +if ( .not. module_initialized ) call static_init_model + +! wrf domain may not equal state_id +call get_model_variable_indices(index_in, i, j, k, var_id=var_id, dom_id=state_id, kind_index=qty) + +location = convert_indices_to_lon_lat_lev(i, j, k, var_id, state_id) + +! return DART variable qty if requested +if(present(qty_out)) qty_out = qty + +end subroutine get_state_meta_data + +!------------------------------------------------------------------ +! observations have a type and qty +! observation type not taken in to account for wrf get close calculations +subroutine get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & + num_close, close_ind, dist, state_handle) + +type(get_close_type), intent(in) :: gc ! handle to a get_close structure +integer, intent(in) :: base_type ! observation TYPE +type(location_type), intent(inout) :: base_loc ! location of interest +type(location_type), intent(inout) :: locs(:) ! obs/state locations +integer, intent(in) :: loc_qtys(:) ! QTYS for obs +integer, intent(in) :: loc_types(:) ! types for obs +integer, intent(out) :: num_close ! how many are close +integer, intent(out) :: close_ind(:) ! indices into the locs array +real(r8), optional, intent(out) :: dist(:) ! distances in radians +type(ensemble_type), optional, intent(in) :: state_handle + +character(len=*), parameter :: routine = 'get_close_obs' +integer :: istatus(1), loc_qtys_ar(1), loc_types_ar(1), i, t_ind +type(location_type) :: loc_ar(1) + +if (vertical_localization_on()) then + loc_ar(1) = base_loc + ! dummy qty 1 not used in convert_vertical_obs + call convert_vertical_obs(state_handle, 1, loc_ar, (/1/), (/base_type/), vert_localization_coord, istatus) + if (istatus(1) /= 0) then + num_close = 0 + return + endif + base_loc = loc_ar(1) +endif + +call loc_get_close(gc, base_loc, base_type, locs, loc_qtys, & + num_close, close_ind) + +if (.not. present(dist)) return + +do i = 1, num_close + t_ind = close_ind(i) + loc_ar(1) = locs(t_ind) + loc_qtys_ar(1) = loc_qtys(t_ind) ! HK not used in convert_vertical_obs todo: use dummy instead? + loc_types_ar(1) = loc_types(t_ind) ! HK not used in convert_vertical_obs todo: use dummy instead? + + istatus(1) = 0 + if (vertical_localization_on()) then + call convert_vertical_obs(state_handle, 1, loc_ar, loc_qtys_ar, loc_types_ar, vert_localization_coord, istatus) + endif + if (istatus(1) == 0) then + dist(i) = get_dist(base_loc, loc_ar(1), base_type, loc_qtys(t_ind)) + else + dist(i) = 1.0e9 + endif +enddo + + +end subroutine get_close_obs + +!------------------------------------------------------------------ +! state only has qty +subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & + num_close, close_ind, dist, state_handle) + +type(get_close_type), intent(in) :: gc +type(location_type), intent(inout) :: base_loc, locs(:) +integer, intent(in) :: base_type, loc_qtys(:) +integer(i8), intent(in) :: loc_indx(:) +integer, intent(out) :: num_close, close_ind(:) +real(r8), optional, intent(out) :: dist(:) +type(ensemble_type), optional, intent(in) :: state_handle + +character(len=*), parameter :: routine = 'get_close_state' +integer :: istatus(1), loc_qtys_ar(1), i, t_ind +integer(i8) :: loc_indx_ar(1) +type(location_type) :: loc_ar(1) + +loc_ar(1) = base_loc + +if (vertical_localization_on()) then + ! dummy qty 1 not used in convert_vertical_obs + call convert_vertical_obs(state_handle, 1, loc_ar, (/1/), (/base_type/), vert_localization_coord, istatus) + if (istatus(1) /= 0) then + num_close = 0 + return + endif + base_loc = loc_ar(1) +endif + +call loc_get_close(gc, base_loc, base_type, locs, loc_qtys, num_close, close_ind) + +if (.not. present(dist)) return + +do i = 1, num_close + t_ind = close_ind(i) + loc_ar(1) = locs(t_ind) + loc_qtys_ar(1) = loc_qtys(t_ind) ! HK not used in convert_vertical_state todo: use dummy instead? + loc_indx_ar(1) = loc_indx(t_ind) + + if (vertical_localization_on()) then + if (nint(query_location(loc_ar(1))) /= vert_localization_coord) then + ! convert_vertical_state always returns istatus = 0 + call convert_vertical_state(state_handle, 1, loc_ar, loc_qtys_ar, loc_indx_ar, vert_localization_coord, istatus(1)) + endif + endif + + dist(i) = get_dist(base_loc, loc_ar(1), base_type, loc_qtys(t_ind)) + +enddo + +end subroutine get_close_state + + +!------------------------------------------------------------------ +! write any additional attributes to netcdf files +! HK todo nc_write_model_atts +subroutine nc_write_model_atts(ncid, domain_id) + +integer, intent(in) :: ncid ! netCDF file identifier +integer, intent(in) :: domain_id + +if ( .not. module_initialized ) call static_init_model + +! put file into define mode. + +call nc_begin_define_mode(ncid) + +call nc_add_global_creation_time(ncid) + +call nc_add_global_attribute(ncid, "model_source", source ) +call nc_add_global_attribute(ncid, "model", "template") + +call nc_end_define_mode(ncid) + +! Flush the buffer and leave netCDF file open +call nc_synchronize_file(ncid) + +end subroutine nc_write_model_atts + +!------------------------------------------------------------------ +subroutine pert_model_copies(ens_handle, ens_size, dummy_pert_amp, interf_provided) + +type(ensemble_type), intent(inout) :: ens_handle +integer, intent(in) :: ens_size +real(r8), intent(in) :: dummy_pert_amp ! not used +logical, intent(out) :: interf_provided + +if (.not. allow_perturbed_ics) then +call error_handler(E_ERR,'pert_model_copies', & + 'starting WRF model from a single vector requires additional steps', & + text2='see comments in wrf/model_mod.f90::pert_model_copies()') +endif + +print*, 'not done' + +end subroutine pert_model_copies + +!------------------------------------------------------------------ +function read_model_time(filename) + +character(len=*), intent(in) :: filename +type(time_type) :: read_model_time + +character(len=*), parameter :: routine = 'read_model_time' +integer :: ncid, ret, var_id +integer :: dim_size(2) ! wrf netcdf: char Times(Time, DateStrLen) +character(len=19) :: timestring ! e.g. 2007-04-26_00:00:00 +integer :: year, month, day, hour, minute, second + +ncid = nc_open_file_readonly(filename, routine) + +call nc_get_variable_size(ncid, 'Times', dim_size, routine) + +ret = nf90_inq_varid(ncid, "Times", var_id) +call nc_check(ret, routine, 'inq_varid Times') + +! last slice of Time dimension +ret = nf90_get_var(ncid, var_id, timestring, start = (/ 1, dim_size(2) /)) +call nc_check(ret, routine, 'get_var Times') + +call get_wrf_date(timestring, year, month, day, hour, minute, second) +read_model_time = set_date(year, month, day, hour, minute, second) + +call nc_close_file(ncid, routine) + +end function read_model_time + +!------------------------------------------------------------------ +subroutine end_model() + +deallocate(wrf_dom, grid, stat_dat) + +end subroutine end_model + +!------------------------------------------------------------------ +! end of public routines +!------------------------------------------------------------------ +function force_non_negative_if_required(n, qty, fld) + +integer, intent(in) :: n +integer, intent(in) :: qty +real(r8), intent(in) :: fld(n) + +real(r8) :: force_non_negative_if_required(n) + +select case (qty) + case (QTY_RAINWATER_MIXING_RATIO, & + QTY_GRAUPEL_MIXING_RATIO, & + QTY_HAIL_MIXING_RATIO, & + QTY_SNOW_MIXING_RATIO, & + QTY_ICE_MIXING_RATIO, & + QTY_CLOUDWATER_MIXING_RATIO, & + QTY_DROPLET_NUMBER_CONCENTR, & + QTY_ICE_NUMBER_CONCENTRATION, & + QTY_SNOW_NUMBER_CONCENTR, & + QTY_RAIN_NUMBER_CONCENTR, & + QTY_GRAUPEL_NUMBER_CONCENTR, & + QTY_HAIL_NUMBER_CONCENTR ) + force_non_negative_if_required = max(0.0_r8, fld) + case default + force_non_negative_if_required = fld +end select + +end function force_non_negative_if_required +!------------------------------------------------------------------ +subroutine update_units_if_required(qty, expected_obs) + +integer, intent(in) :: qty +real(r8), intent(inout) :: expected_obs(:) + +select case (qty) + case (QTY_CLOUDWATER_DE, & + QTY_CLOUD_ICE_DE) + ! convert to diameter (*2) in micrometer (*1e6) + expected_obs = expected_obs * 2.0e6_r8 + case default + ! no unit conversion required +end select + +end subroutine update_units_if_required + +!------------------------------------------------------------------ +! 2D surface variables, or soil variables (z = soil_layers_stag) +function surface_qty(qty) + +integer, intent(in) :: qty +logical :: surface_qty + +!var_id get_varid_from_kind(wrf_dom(id), qty) +! HK todo soil_layers_sta +!do i = 1, get_num_dims(var_id) +!if () then +!get_dim_name +!endif + +select case (qty) + + case (QTY_2M_TEMPERATURE, & + QTY_2M_SPECIFIC_HUMIDITY, & + QTY_10M_U_WIND_COMPONENT, & + QTY_10M_V_WIND_COMPONENT, & + QTY_SURFACE_TYPE, & + QTY_SURFACE_PRESSURE, & + QTY_SKIN_TEMPERATURE, & + QTY_SURFACE_ELEVATION, & + QTY_LANDMASK) + surface_qty = .true. + case default + surface_qty = .false. + +end select + +end function surface_qty + +!------------------------------------------------------------------ +! WRF has separate model qtys for surface variables +function update_qty_if_location_is_surface(qty_in, location) result(qty) + +integer, intent(in) :: qty_in +type(location_type), intent(in) :: location +integer :: qty + +if (.not. is_vertical(location,"SURFACE")) then + qty = qty_in + return +endif + +select case (qty_in) + + case (QTY_U_WIND_COMPONENT); qty = QTY_10M_U_WIND_COMPONENT + case (QTY_V_WIND_COMPONENT); qty = QTY_10M_V_WIND_COMPONENT + case (QTY_POTENTIAL_TEMPERATURE); qty = QTY_2M_POTENTIAL_TEMPERATURE + case (QTY_SPECIFIC_HUMIDITY); qty = QTY_2M_SPECIFIC_HUMIDITY + case (QTY_VAPOR_MIXING_RATIO); qty = QTY_2M_SPECIFIC_HUMIDITY ! Vapor Mixing Ratio (QV, Q2) + case (QTY_PRESSURE); qty = QTY_SURFACE_PRESSURE + case (QTY_TEMPERATURE); qty = QTY_2M_TEMPERATURE + case default; qty = qty_in + +end select + +end function update_qty_if_location_is_surface + +!------------------------------------------------------------------ +function simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: qty +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym + +real(r8) :: simple_interpolation(ens_size) + +integer :: e ! loop variable +! lower left, upper left, lower right, upper right +integer(i8), dimension(ens_size) :: ill, iul, ilr, iur ! dart index at four corners +real(r8), dimension(ens_size) :: x_ill, x_iul, x_ilr, x_iur ! state value at four corners +integer :: var_id + + +var_id = get_varid_from_kind(wrf_dom(id), qty) + +do e = 1, ens_size + ! x, y, z, domain, variable + ill(e) = get_dart_vector_index(ll(1), ll(2), k(e), wrf_dom(id), var_id) + iul(e) = get_dart_vector_index(ul(1), ul(2), k(e), wrf_dom(id), var_id) + ilr(e) = get_dart_vector_index(lr(1), lr(2), k(e), wrf_dom(id), var_id) + iur(e) = get_dart_vector_index(ur(1), ur(2), k(e), wrf_dom(id), var_id) + +enddo + +call get_state_array(x_ill, ill, state_handle) +call get_state_array(x_iul, iul, state_handle) +call get_state_array(x_ilr, ilr, state_handle) +call get_state_array(x_iur, iur, state_handle) + +simple_interpolation = dym*( dxm*x_ill(:) + dx*x_ilr(:) ) + dy*( dxm*x_iul(:) + dx*x_iur(:) ) + +end function simple_interpolation + +!------------------------------------------------------------------ +function vertical_interpolation(ens_size, zloc, fld1, fld2) + +integer, intent(in) :: ens_size +real(r8), intent(in) :: zloc(ens_size) +real(r8), intent(in) :: fld1(ens_size), fld2(ens_size) + +real(r8) :: vertical_interpolation(ens_size) + +real(r8) :: dz, dzm +integer :: z ! level +integer :: e + +do e = 1, ens_size + call toGrid(zloc(e), z, dz, dzm) + ! comment from original code: + ! If you get here and zloc < 1.0, then z will be 0, and + ! we should extrapolate. fld(1,:) and fld(2,:) where computed + ! at levels 1 and 2. + + if (z >= 1) then + ! Linearly interpolate between levels + vertical_interpolation(e) = dzm*fld1(e) + dz*fld2(e) + else + ! Extrapolate below first level. + vertical_interpolation(e) = fld1(e) - (fld2(e)-fld1(e))*dzm + endif +enddo + +end function vertical_interpolation + +!------------------------------------------------------------------ +subroutine get_wrf_date(tstring, year, month, day, hour, minute, second) + +character(len=19), intent(in) :: tstring ! YYYY-MM-DD_hh:mm:ss +integer, intent(out) :: year, month, day, hour, minute, second + +read(tstring( 1: 4),'(i4)') year +read(tstring( 6: 7),'(i2)') month +read(tstring( 9:10),'(i2)') day +read(tstring(12:13),'(i2)') hour +read(tstring(15:16),'(i2)') minute +read(tstring(18:19),'(i2)') second + +end subroutine get_wrf_date + +!------------------------------------------------------------------ +! Verify that the namelist was filled in correctly +! and check there are valid entries for the dart qty. + +subroutine verify_state_variables_generic(nvar, varname, qty, update, in_domain, variable_table) +integer, intent(out) :: nvar +character(len=NF90_MAX_NAME), intent(out) :: varname(MAX_STATE_VARIABLES) +integer, intent(out) :: qty(MAX_STATE_VARIABLES) +logical, intent(out) :: update(MAX_STATE_VARIABLES) +character(len=9), intent(out) :: in_domain(MAX_STATE_VARIABLES) ! assumes <=9 or 999 +character(len=NF90_MAX_NAME), intent(in) :: variable_table(MAX_STATE_VARIABLES*NUM_STATE_TABLE_COLUMNS) + +integer :: i +character(len=NF90_MAX_NAME) :: qty_str, update_str +character(len=256) :: string1, string2 + +if ( .not. module_initialized ) call static_init_model + +nvar = 0 +varloop: do i = 1, MAX_STATE_VARIABLES + if ( variable_table(NUM_STATE_TABLE_COLUMNS*i -3) == 'NULL ' ) exit varloop ! Found end of list. + + varname(i) = trim(variable_table(NUM_STATE_TABLE_COLUMNS*i -3)) + qty_str = trim(variable_table(NUM_STATE_TABLE_COLUMNS*i -2)) + update_str = trim(variable_table(NUM_STATE_TABLE_COLUMNS*i -1)) + in_domain(i) = trim(variable_table(NUM_STATE_TABLE_COLUMNS*i )) + + call to_upper(update_str) + + ! Make sure DART qty is valid + qty(i) = get_index_for_quantity(qty_str) + if( qty(i) < 0 ) then + write(string1,'(''there is no QTY <'',a,''> in obs_kind_mod.f90'')') trim(qty_str) + call error_handler(E_ERR,'verify_state_variables',string1) + endif + + ! Force QTY_TEMPERATURE to QTY_POTENTIAL_TEMPERATURE + if (qty(i) == QTY_TEMPERATURE) qty(i) = QTY_POTENTIAL_TEMPERATURE + + select case (varname(i)) !HK do we need to worry about case sensitivity? + case ('MU'); qty(i) = QTY_PRESSURE !HK this is a hack to avoid 2 QTY_SURFACE_PRESSUREs + case ('PSFC'); qty(i) = QTY_SURFACE_PRESSURE + case ('T2'); qty(i) = QTY_2M_TEMPERATURE + case ('TH2'); qty(i) = QTY_2M_POTENTIAL_TEMPERATURE + case ('Q2'); qty(i) = QTY_2M_SPECIFIC_HUMIDITY !Q2 is actually a mixing ratio, not a specific humidity + end select + + ! Make sure the update variable has a valid name + select case (update_str) + case ('UPDATE') + update(i) = .true. + case ('NO_COPY_BACK') + update(i) = .false. + case default + write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' + write(string2,'(6A)') 'you provided : ', trim(varname(i)), ', ', trim(qty_str), ', ', trim(update_str) + call error_handler(E_ERR,'verify_state_variables',string1, text2=string2) + end select + + nvar = nvar + 1 +enddo varloop +end subroutine verify_state_variables_generic + +! Wrapper for state variables +subroutine verify_state_variables(nvar, varname, qty, update, in_domain) + integer, intent(out) :: nvar + character(len=NF90_MAX_NAME), intent(out) :: varname(MAX_STATE_VARIABLES) + integer, intent(out) :: qty(MAX_STATE_VARIABLES) + logical, intent(out) :: update(MAX_STATE_VARIABLES) + character(len=9), intent(out) :: in_domain(MAX_STATE_VARIABLES) + call verify_state_variables_generic(nvar, varname, qty, update, in_domain, wrf_state_variables) +end subroutine verify_state_variables + +! Wrapper for chemistry variables +subroutine verify_chem_state_variables(nvar, varname, qty, update, in_domain) + integer, intent(out) :: nvar + character(len=NF90_MAX_NAME), intent(out) :: varname(MAX_STATE_VARIABLES) + integer, intent(out) :: qty(MAX_STATE_VARIABLES) + logical, intent(out) :: update(MAX_STATE_VARIABLES) + character(len=9), intent(out) :: in_domain(MAX_STATE_VARIABLES) + call verify_state_variables_generic(nvar, varname, qty, update, in_domain, chem_state_variables) +end subroutine verify_chem_state_variables + +!------------------------------------------------------------------ +! matches WRF variable name in bounds table to input name, and assigns +! the bounds and if the variable is in the state +subroutine get_variable_bounds_generic(var_name, lower, upper, bounds_table) + +character(len=*), intent(in) :: var_name +real(r8), intent(out) :: lower,upper ! bounds +character(len=NF90_MAX_NAME), intent(in) :: bounds_table(NUM_BOUNDS_TABLE_COLUMNS,MAX_STATE_VARIABLES) +character(len=50) :: bound_trim +integer :: ivar + +! default to no bounds +lower = MISSING_R8 +upper = MISSING_R8 + +ivar = 1 +do while ( trim(bounds_table(1,ivar)) /= 'NULL' ) + + if ( trim(bounds_table(1,ivar)) == trim(var_name) ) then + + bound_trim = trim(bounds_table(2,ivar)) + if ( bound_trim /= 'NULL' ) then + read(bound_trim,'(d16.8)') lower + else + lower = MISSING_R8 + endif + + bound_trim = trim(bounds_table(3,ivar)) + if ( bound_trim /= 'NULL' ) then + read(bound_trim,'(d16.8)') upper + else + upper = MISSING_R8 + endif + + endif + + ivar = ivar + 1 + +enddo + +end subroutine get_variable_bounds_generic + +! Wrapper for state variables +subroutine get_variable_bounds(var_name, lower, upper) +character(len=*), intent(in) :: var_name +real(r8), intent(out) :: lower, upper +call get_variable_bounds_generic(var_name, lower, upper, wrf_state_bounds) +end subroutine get_variable_bounds + +! Wrapper for chemistry variables +subroutine get_chem_variable_bounds(var_name, lower, upper) +character(len=*), intent(in) :: var_name +real(r8), intent(out) :: lower, upper +call get_variable_bounds_generic(var_name, lower, upper, chem_state_bounds) +end subroutine get_chem_variable_bounds + + +!------------------------------------------------------------------ +! This is assuming that the number of domains <=9 +! do while loop could be replaced with intrinsic scan +logical function variable_is_on_domain(domain_id_string, id) + +integer, intent(in) :: id +character(len=*), intent(in) :: domain_id_string + +integer :: domain_int, i + +variable_is_on_domain = .false. + +! if '999' then counts all domains +if ( trim(domain_id_string) == '999' ) then + variable_is_on_domain = .true. +else +i = 1 + do while ( domain_id_string(i:i) /= ' ' ) + read(domain_id_string(i:i),'(i1)') domain_int + if ( domain_int == id ) variable_is_on_domain = .true. + i = i+1 + enddo +endif + +end function variable_is_on_domain + +!------------------------------------------------------------------ +subroutine read_grid() + +integer :: ncid, i +character (len=1) :: idom ! assumes <=9 +character(len=*), parameter :: routine = 'read_grid' + +integer :: dim_size(3) ! (west_east, south_north, Time) + +! This is assuming the unlimited dimension is length 1 +! Should we be reading the latest time slice instead? + +do i = 1, num_domains + + write( idom , '(I1)') i + ncid = nc_open_file_readonly('wrfinput_d0'//idom, routine) + + call nc_get_variable_size(ncid, 'XLONG', dim_size) + allocate(grid(i)%longitude(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLONG', grid(i)%longitude, routine) + grid(i)%we = dim_size(1); grid(i)%sn = dim_size(2) + + call nc_get_variable_size(ncid, 'XLONG_U', dim_size) + allocate(grid(i)%longitude_u(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLONG_U', grid(i)%longitude_u, routine) + grid(i)%wes = dim_size(1) + + call nc_get_variable_size(ncid, 'XLONG_V', dim_size) + allocate(grid(i)%longitude_v(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLONG_V', grid(i)%longitude_v, routine) + grid(i)%sns = dim_size(2) + + call nc_get_variable_size(ncid, 'XLAT', dim_size) + allocate(grid(i)%latitude(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLAT', grid(i)%latitude, routine) + + call nc_get_variable_size(ncid, 'XLAT_U', dim_size) + allocate(grid(i)%latitude_u(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLAT_U', grid(i)%latitude_u, routine) + + call nc_get_variable_size(ncid, 'XLAT_V', dim_size) + allocate(grid(i)%latitude_v(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLAT_V', grid(i)%latitude_v, routine) + + call nc_get_global_attribute(ncid, 'MAP_PROJ', grid(i)%map_proj) + call nc_get_global_attribute(ncid, 'DX', grid(i)%dx) + call nc_get_global_attribute(ncid, 'TRUELAT1', grid(i)%truelat1) + call nc_get_global_attribute(ncid, 'TRUELAT2', grid(i)%truelat2) + call nc_get_global_attribute(ncid, 'STAND_LON', grid(i)%stand_lon) + call nc_get_global_attribute(ncid, 'DT', grid(i)%dt) + + grid(i)%bt = nc_get_dimension_size(ncid, 'bottom_top', routine) + grid(i)%bts = nc_get_dimension_size(ncid, 'bottom_top_stag', routine) + + call nc_close_file(ncid, routine) + + call setup_map_projection(i) + + if (i == 1) then + grid(i)%periodic_x = periodic_x + grid(i)%periodic_y = periodic_y + grid(i)%polar = polar + else + grid(i)%periodic_x = .false. + grid(i)%periodic_y = .false. + grid(i)%polar = .false. + endif + +enddo + +end subroutine read_grid + +!------------------------------------------------------------------ +subroutine read_static_data() + +integer :: ncid, i +character (len=1) :: idom ! assumes <=9 +character(len=*), parameter :: routine = 'read_static_data' + +integer :: dim_size(4) ! (west_east, south_north, bottom_top{_stag}, Time) + +do i = 1, num_domains + + write( idom , '(I1)') i + ncid = nc_open_file_readonly('wrfinput_d0'//idom, routine) + + call nc_get_variable_size(ncid, 'PHB', dim_size) + allocate(stat_dat(i)%phb(dim_size(1), dim_size(2), dim_size(3))) + call nc_get_variable(ncid, 'PHB', stat_dat(i)%phb, routine) + + call nc_get_variable_size(ncid, 'MUB', dim_size) + allocate(stat_dat(i)%mub(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'MUB', stat_dat(i)%mub, routine) + + call nc_get_variable_size(ncid, 'HGT', dim_size) + allocate(stat_dat(i)%hgt(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'HGT', stat_dat(i)%hgt, routine) + + call nc_get_variable_size(ncid, 'DNW', dim_size) + allocate(stat_dat(i)%dnw(dim_size(1))) + call nc_get_variable(ncid, 'DNW', stat_dat(i)%dnw, routine) + + call nc_get_variable_size(ncid, 'XLAND', dim_size) + allocate(stat_dat(i)%land(dim_size(1), dim_size(2))) + call nc_get_variable(ncid, 'XLAND', stat_dat(i)%land, routine) + + call nc_get_variable_size(ncid, 'ZS', dim_size) + allocate(stat_dat(i)%zs(dim_size(1))) ! soil_layers_stag + call nc_get_variable(ncid, 'ZS', stat_dat(i)%zs, routine) + + call nc_get_variable(ncid, 'P_TOP', stat_dat(i)%p_top, routine) + + + call nc_close_file(ncid, routine) + +end do + +end subroutine read_static_data + +!------------------------------------------------------------------ +pure function compute_geometric_height(geopot, lat) + +real(r8), intent(in) :: geopot +real(r8), intent(in) :: lat +real(r8) :: compute_geometric_height + + +real(digits12) :: pi2, latr +real(digits12) :: semi_major_axis, semi_minor_axis, grav_polar, grav_equator +real(digits12) :: earth_omega, grav_constant, flattening, somigliana +real(digits12) :: grav_ratio, sin2, termg, termr, grav, eccentricity + +! Parameters below from WGS-84 model software inside GPS receivers. +parameter(semi_major_axis = 6378.1370d3) ! (m) +parameter(semi_minor_axis = 6356.7523142d3) ! (m) +parameter(grav_polar = 9.8321849378) ! (m/s2) +parameter(grav_equator = 9.7803253359) ! (m/s2) +parameter(earth_omega = 7.292115d-5) ! (rad/s) +parameter(grav_constant = 3.986004418d14) ! (m3/s2) +parameter(grav = 9.80665d0) ! (m/s2) WMO std g at 45 deg lat +parameter(eccentricity = 0.081819d0) ! unitless +parameter(pi2 = 3.14159265358979d0/180.d0) + +! Derived geophysical constants +parameter(flattening = (semi_major_axis-semi_minor_axis) / semi_major_axis) + +parameter(somigliana = (semi_minor_axis/semi_major_axis)*(grav_polar/grav_equator)-1.d0) + +parameter(grav_ratio = (earth_omega*earth_omega * & + semi_major_axis*semi_major_axis * semi_minor_axis)/grav_constant) + +! To use geopotential height uncomment the following two lines: +!compute_geometric_height = geopot +!return + +latr = lat * (pi2) ! in radians +sin2 = sin(latr) * sin(latr) +termg = grav_equator * ( (1.d0+somigliana*sin2) / & + sqrt(1.d0-eccentricity*eccentricity*sin2) ) +termr = semi_major_axis / (1.d0 + flattening + grav_ratio - 2.d0*flattening*sin2) + +compute_geometric_height = (termr*geopot) / ( (termg/grav) * termr - geopot ) + + +end function compute_geometric_height + +!------------------------------------------------------------------ +! This is mass level +subroutine get_level_below_obs(which_vert, id, lon_lat_vert, ens_size, state_handle, & + ll, ul, lr, ur, dx, dy, dxm, dym, & + level_below, zloc, fail) + +integer, intent(in) :: which_vert +integer, intent(in) :: id +real(r8), intent(in) :: lon_lat_vert(3) +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, dimension(2), intent(in) :: ll, ul, lr, ur ! (x,y) of each corner +real(r8), intent(in) :: dx, dxm, dy, dym ! grid fractions to obs +integer, intent(out) :: level_below(ens_size) +real(r8), intent(out) :: zloc(ens_size) ! vertical location of the obs for each ens member +logical, intent(out) :: fail + +integer :: e ! loop variable +real(r8) :: v_p(0:grid(id)%bt,ens_size) +real(r8) :: v_h(0:grid(id)%bt,ens_size) +logical :: lev0 + +fail = .false. + +select case (which_vert) + case(VERTISLEVEL) + zloc(:) = lon_lat_vert(3) + level_below(:) = nint(zloc(:)) + fail = .false. + case(VERTISPRESSURE) + call get_model_pressure_profile(id, ll, ul, lr, ur, dx, dy, dxm, dym, ens_size, state_handle, v_p) + do e = 1, ens_size + call pres_to_zk(lon_lat_vert(3), v_p(:,e), grid(id)%bt, zloc(e), level_below(e), lev0, fail) + if (fail) return + if (lev0) then ! pressure obs below lowest sigma + if (.not. allow_obs_below_vol) then + fail = .true. + return + endif + endif + enddo + case(VERTISHEIGHT) + call get_model_height_profile(ll, ul, lr, ur, dx, dy, dxm, dym, id, v_h, state_handle, ens_size) + do e = 1, ens_size + call height_to_zk(lon_lat_vert(3), v_h(:, e), grid(id)%bt, zloc(e), level_below(e), lev0, fail) + if (fail) return + if (lev0) then ! height obs below lowest sigma + if (.not. allow_obs_below_vol) then + fail = .true. + return + endif + endif + enddo + case(VERTISSURFACE) + zloc(:) = 1.0_r8 + level_below(:) = 1 + case(VERTISUNDEF) + zloc = 0.0_r8 + level_below(:) = 1 + case default + fail = .true. +end select + +end subroutine get_level_below_obs + +!------------------------------------------------------------------ +! Check that the surface elevation difference between the model and +! the observation is within acceptable limits +pure function surface_elevation_difference_too_big(id, ll, ul, lr, ur, dx, dy, dxm, dym, elev, tol) result(too_big) + +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) mass grid corners +real(r8), intent(in) :: dx, dy, dxm, dym +real(r8), intent(in) :: elev ! observation surface elevation +real(r8), intent(in) :: tol ! tolerance for elevation difference +logical :: too_big + +real(r8) :: sfc_elevation + +if ( tol < 0.0_r8 ) then + too_big = .false. + return +endif + +sfc_elevation = dym*( dxm*stat_dat(id)%hgt(ll(1), ll(2)) + & + dx*stat_dat(id)%hgt(lr(1), lr(2)) ) + & + dy*( dxm*stat_dat(id)%hgt(ul(1), ul(2)) + & + dx*stat_dat(id)%hgt(ur(1), ur(2)) ) + +too_big = (abs(elev-sfc_elevation) > tol) + +end function surface_elevation_difference_too_big + +!------------------------------------------------------------------ +! Calculate the model pressure profile on half (mass) levels, +! horizontally interpolated at the observation location. +subroutine get_model_pressure_profile(id, ll, ul, lr, ur, dx, dy, dxm, dym, & + ens_size, state_handle, v_p) + +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) mass grid corners +real(r8), intent(in) :: dx, dy, dxm, dym +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +real(r8), intent(out) :: v_p(0:grid(id)%bt, ens_size) + +integer(i8) :: ill, ilr, iul, iur +real(r8), dimension(ens_size) :: x_ill, x_ilr, x_iul, x_iur +real(r8), dimension(ens_size) :: pres1, pres2, pres3, pres4 +real(r8), dimension(ens_size) :: lev2_pres1, lev2_pres2, lev2_pres3, lev2_pres4 + +integer :: var_id, levk +integer :: k(ens_size) + +do levk=1, grid(id)%bt ! number of mass levels + + k(:) = levk + pres1 = model_pressure_t(ll(1), ll(2), k, id, state_handle, ens_size) + pres2 = model_pressure_t(lr(1), lr(2), k, id, state_handle, ens_size) + pres3 = model_pressure_t(ul(1), ul(2), k, id, state_handle, ens_size) + pres4 = model_pressure_t(ur(1), ur(2), k, id, state_handle, ens_size) + + v_p(levk, :) = interp_4pressure(pres1, pres2, pres3, pres4, dx, dxm, dy, dym, ens_size) + + if (levk == 2) then ! store result for extrapolation + lev2_pres1(:) = pres1(:) + lev2_pres2(:) = pres2(:) + lev2_pres3(:) = pres3(:) + lev2_pres4(:) = pres4(:) + endif + +enddo + +var_id = get_varid_from_kind(wrf_dom(id), QTY_SURFACE_PRESSURE) + +if (var_id > 0) then ! surface pressure in domain so get v_p(0,:) from surface pressure + + ill = get_dart_vector_index(ll(1), ll(2), 1, wrf_dom(id), var_id) + ilr = get_dart_vector_index(lr(1), lr(2), 1, wrf_dom(id), var_id) + iul = get_dart_vector_index(ul(1), ul(2), 1, wrf_dom(id), var_id) + iur = get_dart_vector_index(ur(1), ur(2), 1, wrf_dom(id), var_id) + + x_ill = get_state(ill, state_handle) + x_ilr = get_state(ilr, state_handle) + x_iul = get_state(iul, state_handle) + x_iur = get_state(iur, state_handle) + + v_p(0,:) = interp_4pressure(x_ill, x_ilr, x_iul, x_iur, dx, dxm, dy, dym, ens_size) + +!!! Old code: has a check for 0.0 surface pressure +!!! https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L4600-L4606 + +else !extrapolate v_p(0:,) from pressure level 2 and v_p(1:,:) + + v_p(0,:) = extrap_4pressure(lev2_pres1(:), lev2_pres2(:), lev2_pres3(:), lev2_pres4(:), dx, dxm, dy, dym, ens_size, & + edgep=v_p(1,:)) + +endif + +end subroutine get_model_pressure_profile + +!------------------------------------------------------------------ +! returns pressure at a point on the mass grid +function model_pressure_t(i,j,k,id,state_handle, ens_size) + +integer, intent(in) :: ens_size +integer, intent(in) :: i,j,k(ens_size),id +type(ensemble_type), intent(in) :: state_handle +real(r8) :: model_pressure_t(ens_size) + +integer(i8), dimension(ens_size) :: iqv, it +real(r8), dimension(ens_size) :: qvf1, rho, x_iqv, x_it + +integer :: var_idv, var_idt, e + +! Adapted the code from WRF module_big_step_utilities_em.F ---- +! subroutine calc_p_rho_phi Y.-R. Guo (10/20/2004) + +! Simplification: alb*mub = (phb(i,j,k+1) - phb(i,j,k))/dnw(k) + +var_idv = get_varid_from_kind(wrf_dom(id), QTY_VAPOR_MIXING_RATIO) +var_idt = get_varid_from_kind(wrf_dom(id), QTY_POTENTIAL_TEMPERATURE) ! HK original code type_t is this always QTY_POTENTIAL_TEMPERATURE + +if (var_idv < 0 .or. var_idt < 0 ) then + call error_handler(E_ERR, 'model_pressure_t:', 'BOTH QVAPOR and T must be in state vector to compute total pressure', source) +endif + + +do e = 1, ens_size + iqv = get_dart_vector_index(i,j,k(e), wrf_dom(id), var_idv) + it = get_dart_vector_index(i,j,k(e), wrf_dom(id), var_idt) +enddo + +call get_state_array(x_iqv, iqv, state_handle) +call get_state_array(x_it, it, state_handle) + +qvf1(:) = 1.0_r8 + x_iqv(:) / rd_over_rv + +rho(:) = model_rho_t(i,j,k,id,state_handle, ens_size) + +model_pressure_t(:) = ps0 * ( (gas_constant*(ts0+x_it)*qvf1) / & + (ps0/rho(:)) )**cpovcv + +end function model_pressure_t + +!------------------------------------------------------------------ +! returns surface pressure at a point on the mass grid +function model_pressure_s(i, j, id, state_handle, ens_size) + +integer, intent(in) :: i,j,id +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size + +real(r8) :: model_pressure_s(ens_size) + +integer(i8) :: ips, imu +integer :: var_id_psfc, var_id_mu +real(r8) :: x_imu(ens_size), x_ips(ens_size) + +var_id_psfc = get_varid_from_varname(wrf_dom(id), 'PSFC') +var_id_mu = get_varid_from_varname(wrf_dom(id), 'MU') + +if ( var_id_PSFC > 0 ) then + ips = get_dart_vector_index(i,j,1, wrf_dom(id), var_id_psfc) + x_ips = get_state(ips, state_handle) + model_pressure_s = x_ips + +elseif (var_id_mu > 0) then + imu = get_dart_vector_index(i,j,1, wrf_dom(id), var_id_mu) + x_imu = get_state(imu, state_handle) + model_pressure_s = stat_dat(id)%p_top + stat_dat(id)%mub(i,j) + x_imu + +else + call error_handler(E_ERR, 'model_pressure_s:', & + 'One of MU (QTY_PRESSURE) or PSFC (QTY_SURFACE_PRESSURE) must be in state vector to compute surface pressure', & + source) +endif + +end function model_pressure_s + + +!------------------------------------------------------------------ +! Calculate the model level "zk" on half (mass) levels, +! corresponding to pressure "pres" +subroutine pres_to_zk(pres, mdl_v, n3, zk, level_below, lev0, fail) + +real(r8), intent(in) :: pres +real(r8), intent(in) :: mdl_v(0:n3) +integer, intent(in) :: n3 +real(r8), intent(out) :: zk +integer, intent(out) :: level_below +logical, intent(out) :: lev0 +logical, intent(out) :: fail + +integer :: k + +lev0 = .false. +zk = MISSING_R8 +fail = .false. + +! if out of range completely, return missing_r8 and lev0 false +if (pres > mdl_v(0) .or. pres < mdl_v(n3)) then + fail = .true. + return +endif + +! if above surface but below lowest sigma level, return the +! sigma value but set lev0 true +if(pres <= mdl_v(0) .and. pres > mdl_v(1)) then + lev0 = .true. + level_below = 1 + if (log_vert_interp) then + zk = (log(mdl_v(0)) - log(pres))/(log(mdl_v(0)) - log(mdl_v(1))) + else + zk = (mdl_v(0) - pres)/(mdl_v(0) - mdl_v(1)) + endif + return + endif + +! find the 2 sigma levels the value is between and return that +! as a real number, including the fraction between the levels. +do k = 1,n3-1 + if(pres <= mdl_v(k) .and. pres >= mdl_v(k+1)) then + level_below = k + if (log_vert_interp) then + zk = real(k) + (log(mdl_v(k)) - log(pres))/(log(mdl_v(k)) - log(mdl_v(k+1))) + else + zk = real(k) + (mdl_v(k) - pres)/(mdl_v(k) - mdl_v(k+1)) + endif + exit + endif +enddo + +end subroutine pres_to_zk + +!------------------------------------------------------------------ +! Calculate the model height profile on half (mass) levels, +! horizontally interpolated at the observation location. +subroutine get_model_height_profile(ll, ul, lr, ur, dx, dy, dxm, dym, id, v_h, state_handle, ens_size) + +integer, dimension(2), intent(in) :: ll, ul, lr, ur ! (x,y) mass grid corners +integer, intent(in) :: id +real(r8), intent(in) :: dx,dy,dxm,dym +integer, intent(in) :: ens_size +real(r8), intent(out) :: v_h(0:grid(id)%bt, ens_size) +type(ensemble_type), intent(in) :: state_handle +integer e !< for ensemble loop + +real(r8) :: fll(grid(id)%bts, ens_size), geop(ens_size), lat(ens_size) +integer(i8) :: ill, iul, ilr, iur +integer :: k, rc + +real(r8), dimension(ens_size) :: x_ill, x_ilr, x_iul, x_iur +integer :: var_id + +var_id = get_varid_from_kind(wrf_dom(id), QTY_GEOPOTENTIAL_HEIGHT) + +do k = 1, grid(id)%bts ! geopotential height (PH) is on bottom_top_stag + + ill = get_dart_vector_index(ll(1), ll(2), k, wrf_dom(id), var_id) + iul = get_dart_vector_index(ul(1), ul(2), k, wrf_dom(id), var_id) + ilr = get_dart_vector_index(lr(1), lr(2), k, wrf_dom(id), var_id) + iur = get_dart_vector_index(ur(1), ur(2), k, wrf_dom(id), var_id) + + x_ill = get_state(ill, state_handle) + x_ilr = get_state(ilr, state_handle) + x_iul = get_state(iul, state_handle) + x_iur = get_state(iur, state_handle) + + geop(:) = ( dym*( dxm*( stat_dat(id)%phb(ll(1),ll(2),k) + x_ill ) + & + dx*( stat_dat(id)%phb(lr(1),lr(2),k) + x_ilr ) ) + & + dy*( dxm*( stat_dat(id)%phb(ul(1),ul(2),k) + x_iul ) + & + dx*( stat_dat(id)%phb(ur(1),ur(2),k) + x_iur ) ) )/gravity + + lat(:) = ( grid(id)%latitude(ll(1),ll(2)) + & + grid(id)%latitude(lr(1),lr(2)) + & + grid(id)%latitude(ul(1),ul(2)) + & + grid(id)%latitude(ur(1),ur(2)) ) / 4.0_r8 + + do e = 1, ens_size + fll(k, e) = compute_geometric_height(geop(e), lat(e)) + enddo +end do + +do k = 1, grid(id)%bt + v_h(k, :) = 0.5_r8*(fll(k, :) + fll(k+1, :) ) +end do + +v_h(0, :) = dym*( dxm*stat_dat(id)%hgt(ll(1), ll(2)) + & + dx*stat_dat(id)%hgt(lr(1), lr(2)) ) + & + dy*( dxm*stat_dat(id)%hgt(ul(1), ul(2)) + & + dx*stat_dat(id)%hgt(ur(1), ur(2)) ) + + + +end subroutine get_model_height_profile + +!------------------------------------------------------------------ +! Calculate the model level zk on half (mass) levels, +! corresponding to height obs_v. +subroutine height_to_zk(obs_v, mdl_v, n3, zk, level_below, lev0, fail) + +real(r8), intent(in) :: obs_v +integer, intent(in) :: n3 +real(r8), intent(in) :: mdl_v(0:n3) +real(r8), intent(out) :: zk +integer, intent(out) :: level_below +logical, intent(out) :: lev0 +logical, intent(out) :: fail + +integer :: k + +zk = MISSING_R8 +lev0 = .false. +fail = .false. + +! if out of range completely, return missing_r8 and lev0 false +if (obs_v < mdl_v(0) .or. obs_v > mdl_v(n3)) then + fail = .true. + return +endif + +! if above surface but below lowest 3-d height level, return the +! height value but set lev0 true +if(obs_v >= mdl_v(0) .and. obs_v < mdl_v(1)) then + lev0 = .true. + level_below = 1 + zk = (mdl_v(0) - obs_v)/(mdl_v(0) - mdl_v(1)) + return +endif + +! find the 2 height levels the value is between and return that +! as a real number, including the fraction between the levels. +do k = 1,n3-1 + if(obs_v >= mdl_v(k) .and. obs_v <= mdl_v(k+1)) then + level_below = k + zk = real(k) + (mdl_v(k) - obs_v)/(mdl_v(k) - mdl_v(k+1)) + exit + endif +enddo + +end subroutine height_to_zk + +!------------------------------------------------------------------ +! Interpolate pressure inside quad +! Four corners of a quad: +! p1 lower left +! p2 lower right +! p3 upper left +! p4 upper right +! dx is the distance in x, dxm is 1.0-dx +! dy is distance in y, dym is 1.0-dy +function interp_4pressure(p1, p2, p3, p4, dx, dxm, dy, dym, ens_size) + +integer, intent(in) :: ens_size +real(r8), intent(in) :: p1(ens_size), p2(ens_size), p3(ens_size), p4(ens_size) +real(r8), intent(in) :: dx, dxm, dy, dym +real(r8) :: interp_4pressure(ens_size) + +real(r8) :: l1(ens_size), l2(ens_size), l3(ens_size), l4(ens_size) + +if (log_horz_interpQ) then + l1 = log(p1) + l2 = log(p2) + l3 = log(p3) + l4 = log(p4) +endif + +! once we like the results, remove the log_horz_interpQ test. +if (log_horz_interpQ) then + interp_4pressure = exp(dym*( dxm*l1 + dx*l2 ) + dy*( dxm*l3 + dx*l4 )) +else + interp_4pressure = dym*( dxm*p1 + dx*p2 ) + dy*( dxm*p3 + dx*p4 ) +endif + +end function interp_4pressure + +!------------------------------------------------------------------ +! extrapolate quad where edgep is the edge pressure +function extrap_4pressure(p1, p2, p3, p4, dx, dxm, dy, dym, ens_size, edgep) + +integer, intent(in) :: ens_size +real(r8), intent(in) :: p1(ens_size), p2(ens_size), p3(ens_size), p4(ens_size) +real(r8), intent(in) :: dx, dxm, dy, dym +real(r8), intent(in) :: edgep(ens_size) +real(r8) :: extrap_4pressure(ens_size) + +real(r8) :: intermediate(ens_size) +real(r8) :: l1(ens_size), l2(ens_size), l3(ens_size), l4(ens_size) + +if (log_horz_interpQ) then + l1 = log(p1) + l2 = log(p2) + l3 = log(p3) + l4 = log(p4) +endif + +! once we like the results, remove the log_horz_interpQ test. +if (log_horz_interpQ) then + intermediate = (3.0_r8*log(edgep) - & + dym*( dxm*l1 + dx*l2 ) - dy*( dxm*l3 + dx*l4 ))/2.0_r8 + + where (intermediate <= 0.0_r8) + extrap_4pressure = edgep + else where + extrap_4pressure = exp(intermediate) + end where +else + extrap_4pressure = (3.0_r8*edgep - & + dym*( dxm*p1 + dx*p2 ) - dy*( dxm*p3 + dx*p4 ))/2.0_r8 +endif + +end function extrap_4pressure + +!------------------------------------------------------------------ +! Calculate the total density on mass point (half (mass) levels, T-point). +function model_rho_t(i,j,k,id,state_handle, ens_size) + +integer, intent(in) :: ens_size +integer, intent(in) :: i,j,k(ens_size),id +type(ensemble_type), intent(in) :: state_handle +real(r8) :: model_rho_t(ens_size) + +integer(i8), dimension(ens_size) :: imu,iph,iphp1 +real(r8), dimension(ens_size) :: ph_e, x_imu, x_iph, x_iphp1 +integer :: var_id_mu, var_id_ph, e + +! Adapted the code from WRF module_big_step_utilities_em.F ---- +! subroutine calc_p_rho_phi Y.-R. Guo (10/20/2004) + +! Simplification: alb*mub = (phb(i,j,k+1) - phb(i,j,k))/dnw(k) + +var_id_mu = get_varid_from_varname(wrf_dom(id), 'MU') +var_id_ph = get_varid_from_kind(wrf_dom(id), QTY_GEOPOTENTIAL_HEIGHT) + +if (var_id_mu < 0 .or. var_id_ph < 0) then + call error_handler(E_ERR, 'model_rho_t:', 'BOTH MU and PH must be in state vector to compute total density', source) +endif + +do e = 1, ens_size + imu = get_dart_vector_index(i,j,1, wrf_dom(id), var_id_mu) + iph = get_dart_vector_index(i,j,k(e), wrf_dom(id), var_id_ph) + iphp1 = get_dart_vector_index(i,j,k(e)+1, wrf_dom(id), var_id_ph) +enddo + +call get_state_array(x_imu, imu, state_handle) +call get_state_array(x_iph, iph, state_handle) +call get_state_array(x_iphp1, iphp1, state_handle) + +do e = 1, ens_size + ph_e(e) = ( (x_iphp1(e) + stat_dat(id)%phb(i,j,k(e)+1)) & + - (x_iph(e) + stat_dat(id)%phb(i,j,k(e) )) ) / stat_dat(id)%dnw(k(e)) +enddo + +! rho = - mu / dphi/deta +model_rho_t(:) = - (stat_dat(id)%mub(i,j)+x_imu) / ph_e + +end function model_rho_t + + +!------------------------------------------------------------------ +function density_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: density_interpolate(ens_size) + +real(r8), dimension(ens_size) :: rho1, rho2, rho3, rho4 + +rho1 = model_rho_t(ll(1), ll(2), k, id, state_handle, ens_size) +rho2 = model_rho_t(lr(1), lr(2), k, id, state_handle, ens_size) +rho3 = model_rho_t(ul(1), ul(2), k, id, state_handle, ens_size) +rho4 = model_rho_t(ur(1), ur(2), k, id, state_handle, ens_size) + +density_interpolate = dym*( dxm*rho1(:) + dx*rho2(:) ) + dy*( dxm*rho3(:) + dx*rho4(:) ) + +end function density_interpolate + +!------------------------------------------------------------------ +! wrfinput land mask XLAND 1 = land, 2 = water +! obs_def_rttov_mod 0 = land, 1 = water, 2 = sea ice +function surface_type_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: surface_type_interpolate(ens_size) ! same across the ensemble + +surface_type_interpolate(:) = -1.d0 + landmask_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + +end function surface_type_interpolate + +!------------------------------------------------------------------ +! land mask is static data XLAND 1 = land, 0 = wate +function landmask_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: landmask_interpolate(ens_size) ! same across the ensemble + +landmask_interpolate = dym*( dxm*stat_dat(id)%land(ll(1), ll(2)) + & + dx*stat_dat(id)%land(lr(1), lr(2)) ) + & + dy*( dxm*stat_dat(id)%land(ul(1), ul(2)) + & + dx*stat_dat(id)%land(ur(1), ur(2)) ) + +end function landmask_interpolate + + +!------------------------------------------------------------------ +function surface_elevation_interpolate(ens_size, id, ll, ul, lr, ur, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: surface_elevation_interpolate(ens_size) + +surface_elevation_interpolate(:) = dym*( dxm*stat_dat(id)%hgt(ll(1), ll(2)) + & + dx*stat_dat(id)%hgt(lr(1), lr(2)) ) + & + dy*( dxm*stat_dat(id)%hgt(ul(1), ul(2)) + & + dx*stat_dat(id)%hgt(ur(1), ur(2)) ) + +end function surface_elevation_interpolate + +!------------------------------------------------------------------ +function geopotential_height_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: qty +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: geopotential_height_interpolate(ens_size) + +real(r8), dimension(ens_size) :: a1 + +a1 = simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +! phb is constant across the ensemble, so use k(1) +geopotential_height_interpolate = ( a1 + & + dym* ( dxm*stat_dat(id)%phb(ll(1), ll(2), k(1) ) + & + dx * stat_dat(id)%phb(lr(1), lr(2), k(1)) ) + & + dy * ( dxm*stat_dat(id)%phb(ul(1), ul(2), k(1) ) + & + dx * stat_dat(id)%phb(ur(1), ur(2), k(1)) ) ) / gravity + +end function geopotential_height_interpolate + +!------------------------------------------------------------------ +function temperature_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: temperature_interpolate(ens_size) + +real(r8), dimension(ens_size) :: a1, pres, pres1, pres2, pres3, pres4 + +! In terms of perturbation potential temperature +a1 = simple_interpolation(ens_size, state_handle, QTY_POTENTIAL_TEMPERATURE, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +pres1 = model_pressure_t(ll(1), ll(2), k, id, state_handle, ens_size) +pres2 = model_pressure_t(lr(1), lr(2), k, id, state_handle, ens_size) +pres3 = model_pressure_t(ul(1), ul(2), k, id, state_handle, ens_size) +pres4 = model_pressure_t(ur(1), ur(2), k, id, state_handle, ens_size) + +! Pressure at location +pres = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) + +! Full sensible temperature field +temperature_interpolate = (ts0 + a1(:))*(pres(:)/ps0)**kappa + +end function temperature_interpolate + +!------------------------------------------------------------------ +function pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: pressure_interpolate(ens_size) + +real(r8), dimension(ens_size) :: pres1, pres2, pres3, pres4 + +pres1 = model_pressure_t(ll(1), ll(2), k, id, state_handle, ens_size) +pres2 = model_pressure_t(lr(1), lr(2), k, id, state_handle, ens_size) +pres3 = model_pressure_t(ul(1), ul(2), k, id, state_handle, ens_size) +pres4 = model_pressure_t(ur(1), ur(2), k, id, state_handle, ens_size) + +! Pressure at location +pressure_interpolate = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) + +end function pressure_interpolate + +!------------------------------------------------------------------ + +function pressure_interpolate_vertically(ens_size, state_handle, id, xy_point, k, zloc) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: xy_point(2) ! (x,y) at one corner +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: zloc(ens_size) + +real(r8) :: pressure_interpolate_vertically(ens_size) + +real(r8), dimension(ens_size) :: pres_k, pres_kp1 + +pres_k = model_pressure_t(xy_point(1), xy_point(2), k, id, state_handle, ens_size) +pres_kp1 = model_pressure_t(xy_point(1), xy_point(2), k+1, id, state_handle, ens_size) + +pressure_interpolate_vertically = vertical_interpolation(ens_size, zloc, pres_k, pres_kp1) + +end function pressure_interpolate_vertically + + + +!------------------------------------------------------------------ +function surface_pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: surface_pressure_interpolate(ens_size) + +real(r8), dimension(ens_size) :: x_ill, x_iul, x_ilr, x_iur +integer :: e + +x_ill = model_pressure_s(ll(1), ll(2), id, state_handle, ens_size) +x_iul = model_pressure_s(ul(1), ul(2), id, state_handle, ens_size) +x_ilr = model_pressure_s(lr(1), lr(2), id, state_handle, ens_size) +x_iur = model_pressure_s(ur(1), ur(2), id, state_handle, ens_size) + +! Pressure at location + +do e = 1, ens_size + ! HK todo original code comment: + ! I'm not quite sure where this comes from, but I will trust them on it.... + if ( x_ill(e) /= 0.0_r8 .and. x_ilr(e) /= 0.0_r8 .and. x_iul(e) /= 0.0_r8 .and. & + x_iur(e) /= 0.0_r8 ) then + + surface_pressure_interpolate(e) = dym*( dxm*x_ill(e) + dx*x_ilr(e) ) + dy*( dxm*x_iul(e) + dx*x_iur(e) ) + endif !HK todo initialize to missing_r8? +enddo + + +end function surface_pressure_interpolate + +!------------------------------------------------------------------ +function specific_humidity_interpolate(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: qty +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: specific_humidity_interpolate(ens_size) + +real(r8), dimension(ens_size) :: a1 +a1 = simple_interpolation(ens_size, state_handle, QTY_VAPOR_MIXING_RATIO, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +specific_humidity_interpolate = a1(:) /(1.0_r8 + a1(:)) + +end function specific_humidity_interpolate + +!------------------------------------------------------------------ +function wind_interpolate(ens_size, state_handle, qty, id, k, xloc, yloc, i, j, dxm, dx, dy, dym, lon) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: qty +integer, intent(in) :: id +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: xloc, yloc ! location on mass grid +integer, intent(in) :: i,j ! ll corners of mass grid +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8), intent(in) :: lon ! Longitude of point in degrees +real(r8) :: wind_interpolate(ens_size) + +real(r8), dimension(ens_size) :: u_wind_grid, v_wind_grid, u_wind, v_wind +real(r8) :: xloc_u, yloc_v ! x ugrid, y vgrid +integer :: u_qty, v_qty ! quantity ids for u and v components +real(r8) :: dx_u, dxm_u, dy_v, dym_v +integer :: i_u, j_v +integer :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer :: e, rc + +! HK TODO relationship between mass grid and u grid and v grid +! Original code adds 0.5 to xloc, yloc. But what if you are on the edge of a domain? +! https://github.com/NCAR/DART/blob/70e6af803a52d14b9f77f872c94b1fe11d5dc2d9/models/wrf/model_mod.f90#L1425-L1432 + +! xloc and yloc are indices on mass-grid. If we are on a periodic longitude domain, +! then xloc can range from [1 wes). This means that simply adding 0.5 to xloc has +! the potential to render xloc_u out of the valid mass-grid index bounds (>wes). +! Only add 0.5 offset for 3D winds, not surface winds +if (qty == QTY_U_WIND_COMPONENT .or. qty == QTY_V_WIND_COMPONENT) then + xloc_u = xloc + 0.5_r8 + yloc_v = yloc + 0.5_r8 + u_qty = QTY_U_WIND_COMPONENT + v_qty = QTY_V_WIND_COMPONENT +else + ! For surface winds (QTY_10M_U_WIND_COMPONENT, QTY_10M_V_WIND_COMPONENT), no offset + xloc_u = xloc + yloc_v = yloc + u_qty = QTY_10M_U_WIND_COMPONENT + v_qty = QTY_10M_V_WIND_COMPONENT +endif + +! HK TODO what about periodic_y? +if ( grid(id)%periodic_x .and. xloc_u > real(grid(id)%wes,r8) ) xloc_u = xloc_u - real(grid(id)%we,r8) + +call toGrid(xloc_u,i_u,dx_u,dxm_u) +call toGrid(yloc_v,j_v,dy_v,dym_v) + +call getCorners(i_u, j, id, u_qty, ll, ul, lr, ur, rc) +u_wind_grid = simple_interpolation(ens_size, state_handle, u_qty, id, ll, ul, lr, ur, k, dxm_u, dx_u, dy, dym) +call getCorners(i, j_v, id, v_qty, ll, ul, lr, ur, rc) +v_wind_grid = simple_interpolation(ens_size, state_handle, v_qty, id, ll, ul, lr, ur, k, dxm, dx, dy_v, dym_v) + +do e = 1, ens_size + call gridwind_to_truewind(lon, grid(id)%proj, u_wind_grid(e), v_wind_grid(e), u_wind(e), v_wind(e)) +enddo + +if ( qty == u_qty ) then + wind_interpolate = u_wind +else + wind_interpolate = v_wind +endif + +end function wind_interpolate + +!------------------------------------------------------------------ +! If there are other domains in the state: +! wrf domain id =/ state domain id +function get_wrf_domain(state_id) + +integer, intent(in) :: state_id +integer :: get_wrf_domain + +integer :: i + +do i = 1, num_state_domains() + if (wrf_dom(i) == state_id) then + get_wrf_domain = mod(i-1, num_domains) + 1 + return + endif +enddo + +get_wrf_domain = -1 ! not found + +end function get_wrf_domain + +!------------------------------------------------------------------ +subroutine convert_vertical_state(state_handle, num, locs, loc_qtys, loc_indx, & + which_vert, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: num +type(location_type), intent(inout) :: locs(num) ! location for each state element +integer, intent(in) :: loc_qtys(num) ! qty for each state element +integer(i8), intent(in) :: loc_indx(num) ! index into the state vector +integer, intent(in) :: which_vert ! vertical coordinate to be converted to +integer, intent(out) :: istatus + +integer :: ip, jp, kp, id ! x,y,z of state index +integer :: var_id, state_id +real(r8) :: vert ! vertical after conversion +integer :: i +integer :: vert_coord_in +real(r8) :: lon_lat_vert(3) + +istatus = 0 ! can not fail + +do i = 1, num + + lon_lat_vert = get_location(locs(i)) + vert_coord_in = nint(query_location(locs(i))) + + if (vert_coord_in == which_vert) then ! conversion is already done + cycle + endif + + call get_model_variable_indices(loc_indx(i), ip, jp, kp, var_id=var_id, dom_id=state_id) + id = get_wrf_domain(state_id) + + if (which_vert == VERTISLEVEL) then + + if (on_w_grid(state_id, var_id)) then + vert = real(kp) - 0.5_r8 + else + vert = real(kp) + endif + + elseif (which_vert == VERTISPRESSURE) then + + vert = model_pressure(ip, jp, kp, id, var_id, state_id, state_handle) + + elseif (which_vert == VERTISHEIGHT) then + + vert = model_height(ip, jp, kp, id, loc_qtys(i), var_id, state_id, state_handle) + + elseif (which_vert == VERTISSCALEHEIGHT) then + + vert = -log(model_pressure(ip, jp, kp, id, var_id, state_id, state_handle) / & + model_surface_pressure(ip, jp, id, var_id, state_id, state_handle)) + + endif + + call set_vertical(locs(i), vert, which_vert) + +enddo + + + +end subroutine convert_vertical_state + +!------------------------------------------------------------------ +subroutine convert_vertical_obs(state_handle, num, locs, loc_qtys, loc_types, & + which_vert, istatus) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: num +type(location_type), intent(inout) :: locs(num) +integer, intent(in) :: loc_qtys(num) +integer, intent(in) :: loc_types(num) +integer, intent(in) :: which_vert +integer, intent(out) :: istatus(num) + +integer :: vert_coord_in, id +real(r8) :: lon_lat_vert(3) +integer :: ll(2), ul(2), lr(2), ur(2) !(x,y) of four corners +integer :: i, j ! grid +real(r8) :: xloc, yloc ! WRF i,j in the grid +real(r8) :: dx, dxm, dy, dym ! grid fractions +real(r8) :: zout(1), zk(1), zk1(1), geop(1), zloc(1) +real(r8) :: h1(1), h2(1), h3(1), h4(1) +real(r8) :: pres1(1), pres2(1), pres3(1), pres4(1) +integer :: ens_size, rc, k(1) +logical :: fail +integer :: ob ! loop variable +real(r8) :: dz, dzm ! dummys for toGrid call + +integer, parameter :: FAILED_BOUNDS_CHECK = 144 +integer, parameter :: CANNOT_INTERPOLATE_QTY = 155 +integer, parameter :: NOT_IN_ANY_DOMAIN = 111 +integer, parameter :: VERTICAL_LOCATION_FAIL = 166 + +ens_size = 1 ! working with the mean state + +do ob = 1, num + + lon_lat_vert = get_location(locs(ob)) + vert_coord_in = nint(query_location(locs(ob))) + + if (vert_coord_in == which_vert) then ! conversion is already done + istatus(ob) = 0 + cycle + endif + + if (lon_lat_vert(3) == VERTISUNDEF) then ! no vertical, no conversion + istatus(ob) = 0 + cycle + endif + + if (lon_lat_vert(3) == MISSING_R8) then ! vertical is missing, no conversion + call set_vertical(locs(ob), MISSING_R8, which_vert) + istatus(ob) = 1 ! HK todo original code does not set success for this + cycle + endif + + ! convert to which_vert + call get_domain_info(lon_lat_vert(1),lon_lat_vert(2),id,xloc,yloc) + if (id == 0) then + call set_vertical(locs(ob), MISSING_R8, which_vert) !HK original code - why set to missing? + istatus(ob) = NOT_IN_ANY_DOMAIN + cycle + endif + + ! horizontal location mass point + call toGrid(xloc,i,dx,dxm) + call toGrid(yloc,j,dy,dym) + + if ( .not. within_bounds_horizontal(i, j, id, QTY_POTENTIAL_TEMPERATURE) ) then ! HK origianal code mass grid qty + call set_vertical(locs(ob), MISSING_R8, which_vert) + istatus(ob) = FAILED_BOUNDS_CHECK + cycle + endif + + call getCorners(i, j, id, QTY_POTENTIAL_TEMPERATURE, ll, ul, lr, ur, rc) !HK todo what qty to pick? + + ! vertical location mass level + call get_level_below_obs(vert_coord_in, id, lon_lat_vert, ens_size, state_handle, ll, ul, lr, ur, dx, dy, dxm, dym, k, zloc, fail) + if (fail) then + istatus(ob) = VERTICAL_LOCATION_FAIL + ! set vertical? + cycle + endif + + select case (which_vert) + + case (VERTISLEVEL) + + if (vert_coord_in == VERTISSURFACE) then + zout = 1.0_r8 + else + zout = zloc + endif + + case (VERTISPRESSURE) + + if (vert_coord_in == VERTISSURFACE) then + ! compute surface pressure at all neighboring mass points + pres1 = model_pressure_s(ll(1), ll(2), id, state_handle, ens_size) + pres2 = model_pressure_s(lr(1), lr(2), id, state_handle, ens_size) + pres3 = model_pressure_s(ul(1), ul(2), id, state_handle, ens_size) + pres4 = model_pressure_s(ur(1), ur(2), id, state_handle, ens_size) + zout = dym*( dxm*pres1(1) + dx*pres2(1) ) + dy*( dxm*pres3(1) + dx*pres4(1) ) + + else + zk = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + zk1 = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + zout = vertical_interpolation(ens_size, zloc, zk, zk1) + endif + + case (VERTISHEIGHT) + + if (vert_coord_in == VERTISSURFACE) then + ! a surface ob is assumed to have height as vertical coordinate. + ! this code needs to be revised if this is not true + ! (in that case uncomment lines below to get terrain height + ! from model) + zout = lon_lat_vert(3) + !! or: directly interpolate terrain height at neighboring mass points + !zout = dym*( dxm*stat_dat(id)%hgt(i, j) + & + ! dx*stat_dat(id)%hgt(i+1,j) ) + & + ! dy*( dxm*stat_dat(id)%hgt(i, j+1) + & + ! dx*stat_dat(id)%hgt(i+1,j+1) ) + else + ! adding 0.5 to get to the staggered vertical grid for height + zloc = zloc + 0.5_r8 ! Adjust zloc for staggered + call toGrid(zloc(1),k(1),dz,dzm) + + ! This method does not give bitwise with main + !zk = geopotential_height_interpolate(ens_size, state_handle, QTY_GEOPOTENTIAL_HEIGHT, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + !zk1 = geopotential_height_interpolate(ens_size, state_handle, QTY_GEOPOTENTIAL_HEIGHT, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + !geop = vertical_interpolation(ens_size, zloc, zk, zk1) + !zout = compute_geometric_height(geop(1), lon_lat_vert(2)) + + ! horizontal then vertical interpolation, not what main code does + zk = interpolate_geometric_height(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + zk1 = interpolate_geometric_height(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + zout = vertical_interpolation(ens_size, zloc, zk, zk1) + + ! vertical then horizontal interpolation + h1 = interpolate_geometric_height_vertically(ens_size, state_handle, id, ll, k, zloc) + h2 = interpolate_geometric_height_vertically(ens_size, state_handle, id, lr, k, zloc) + h3 = interpolate_geometric_height_vertically(ens_size, state_handle, id, ul, k, zloc) + h4 = interpolate_geometric_height_vertically(ens_size, state_handle, id, ur, k, zloc) + + zout = dym*( dxm*h1(1) + dx*h2(1) ) + dy*( dxm*h3(1) + dx*h4(1) ) + + endif + + case (VERTISSCALEHEIGHT) + if (vert_coord_in == VERTISSURFACE) then + zout = -log(1.0_r8) + elseif (vert_coord_in == VERTISPRESSURE) then + ! surface pressure + pres1 = model_pressure_s(ll(1), ll(2), id, state_handle, ens_size) + pres2 = model_pressure_s(lr(1), lr(2), id, state_handle, ens_size) + pres3 = model_pressure_s(ul(1), ul(2), id, state_handle, ens_size) + pres4 = model_pressure_s(ur(1), ur(2), id, state_handle, ens_size) + + zout = -log(lon_lat_vert(3) / (dym*( dxm*pres1(1) + dx*pres2(1) ) + dy*( dxm*pres3(1) + dx*pres4(1) ))) + + else ! vert_cood_in == VERTISHEIGHT + + ! horizontal interpolation, then vertical interpolation + !zk = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + !zk1 = pressure_interpolate(ens_size, state_handle, id, ll, ul, lr, ur, k+1, dxm, dx, dy, dym) + !zout = vertical_interpolation(ens_size, zloc, zk, zk1) + + ! vertical interpolation, then horizontal interpolation + pres1 = pressure_interpolate_vertically(ens_size, state_handle, id, ll, k, zloc) + pres2 = pressure_interpolate_vertically(ens_size, state_handle, id, lr, k, zloc) + pres3 = pressure_interpolate_vertically(ens_size, state_handle, id, ul, k, zloc) + pres4 = pressure_interpolate_vertically(ens_size, state_handle, id, ur, k, zloc) + + zout = -log(zout / (dym*( dxm*pres1(1) + dx*pres2(1) ) + dy*( dxm*pres3(1) + dx*pres4(1) ))) + + endif + + + end select + + !HK original code uses set_location, you could use set_vertical here see #621 + locs(ob) = set_location(lon_lat_vert(1), lon_lat_vert(2), zout(1), which_vert) + istatus(ob) = 0 + +enddo + +end subroutine convert_vertical_obs + + +!------------------------------------------------------------------ +function interpolate_geometric_height(ens_size, state_handle, id, ll, ul, lr, ur, k, dxm, dx, dy, dym) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: dxm, dx, dy, dym +real(r8) :: interpolate_geometric_height(ens_size) + +integer :: e ! loop variable +! lower left, upper left, lower right, upper right +integer(i8), dimension(ens_size) :: ill, iul, ilr, iur ! dart index at four corners +real(r8), dimension(ens_size) :: x_ill, x_iul, x_ilr, x_iur ! state value at four corners +real(r8), dimension(ens_size) :: geop_ll, geop_ul, geop_lr, geop_ur ! geopotential height at four corners +real(r8), dimension(ens_size) :: geomet_ll, geomet_ul, geomet_lr, geomet_ur ! geometric height at four corners +integer :: var_id + +var_id = get_varid_from_kind(wrf_dom(id), QTY_GEOPOTENTIAL_HEIGHT) + +do e = 1, ens_size + ! x, y, z, domain, variable + ill(e) = get_dart_vector_index(ll(1), ll(2), k(e), wrf_dom(id), var_id) + iul(e) = get_dart_vector_index(ul(1), ul(2), k(e), wrf_dom(id), var_id) + ilr(e) = get_dart_vector_index(lr(1), lr(2), k(e), wrf_dom(id), var_id) + iur(e) = get_dart_vector_index(ur(1), ur(2), k(e), wrf_dom(id), var_id) + +enddo + +call get_state_array(x_ill, ill, state_handle) +call get_state_array(x_iul, iul, state_handle) +call get_state_array(x_ilr, ilr, state_handle) +call get_state_array(x_iur, iur, state_handle) + +do e = 1, ens_size + geop_ll(e) = (stat_dat(id)%phb(ll(1), ll(2), k(e)) + x_ill(e)) / gravity + geop_ul(e) = (stat_dat(id)%phb(ul(1), ul(2), k(e)) + x_iul(e)) / gravity + geop_lr(e) = (stat_dat(id)%phb(lr(1), lr(2), k(e)) + x_ilr(e)) / gravity + geop_ur(e) = (stat_dat(id)%phb(ur(1), ur(2), k(e)) + x_iur(e)) / gravity + + geomet_ll(e) = compute_geometric_height(geop_ll(e), grid(id)%latitude(ll(1), ll(2))) + geomet_ul(e) = compute_geometric_height(geop_ul(e), grid(id)%latitude(ul(1), ul(2))) + geomet_lr(e) = compute_geometric_height(geop_lr(e), grid(id)%latitude(lr(1), lr(2))) + geomet_ur(e) = compute_geometric_height(geop_ur(e), grid(id)%latitude(ur(1), ur(2))) + +enddo + +interpolate_geometric_height = dym*( dxm*geomet_ll(:) + dx*geomet_lr(:) ) + dy*( dxm*geomet_ul(:) + dx*geomet_ur(:) ) + +end function interpolate_geometric_height + +!------------------------------------------------------------------ +function interpolate_geometric_height_vertically(ens_size, state_handle, id, xy_point, k, zloc) + +integer, intent(in) :: ens_size +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: id +integer, intent(in) :: xy_point(2) ! (x,y) at one corner +integer, intent(in) :: k(ens_size) ! k may be different across the ensemble +real(r8), intent(in) :: zloc(ens_size) + +real(r8) :: interpolate_geometric_height_vertically(ens_size) + +integer :: e ! loop variable +! lower left, upper left, lower right, upper right +integer(i8), dimension(ens_size) :: idx_k, idx_kp1 ! dart index at point +real(r8), dimension(ens_size) :: x_k, x_kp1 ! state value at point +real(r8), dimension(ens_size) :: geop_k, geop_kp1 ! geopotential height at point +real(r8), dimension(ens_size) :: geomet_k, geomet_kp1 ! geometric height at point +integer :: var_id + +var_id = get_varid_from_kind(wrf_dom(id), QTY_GEOPOTENTIAL_HEIGHT) + +do e = 1, ens_size + ! x, y, z, domain, variable + idx_k(e) = get_dart_vector_index(xy_point(1), xy_point(2), k(e), wrf_dom(id), var_id) + idx_kp1(e) = get_dart_vector_index(xy_point(1), xy_point(2), k(e)+1, wrf_dom(id), var_id) + +enddo + +call get_state_array(x_k, idx_k, state_handle) +call get_state_array(x_kp1, idx_kp1, state_handle) + +do e = 1, ens_size + geop_k(e) = (stat_dat(id)%phb(xy_point(1), xy_point(2), k(e)) + x_k(e)) / gravity + geop_kp1(e) = (stat_dat(id)%phb(xy_point(1), xy_point(2), k(e)+1) + x_kp1(e)) / gravity + geomet_k(e) = compute_geometric_height(geop_k(e), grid(id)%latitude(xy_point(1), xy_point(2))) + geomet_kp1(e) = compute_geometric_height(geop_kp1(e), grid(id)%latitude(xy_point(1), xy_point(2))) +enddo + +interpolate_geometric_height_vertically = vertical_interpolation(ens_size, zloc, geomet_k, geomet_kp1) + +end function interpolate_geometric_height_vertically + +!------------------------------------------------------------------ +! model height any grid u,v,w,t +! used in convert_vertical_state so ens_size = 1 +function model_height(i, j, k, id, qty, var_id, state_id, state_handle) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: i, j, k, id, qty +integer, intent(in) :: var_id,state_id +real(r8) :: model_height + +integer(i8) :: i1, i2, i3, i4 +integer :: off +real(r8) :: x_i1(1), x_i2(1), x_i3(1), x_i4(1) +real(r8) :: geop, lat +integer :: gz_id + +gz_id = get_varid_from_kind(state_id, QTY_GEOPOTENTIAL_HEIGHT) + +!HK todo for these special cases would it be better to check by variable name +! instead of QTY? +if( qty == QTY_PRESSURE .or. & ! MU + qty == QTY_SURFACE_PRESSURE .or. & ! PSFC SFC PRESSUR + qty == QTY_SKIN_TEMPERATURE) then ! TSK SURFACE SKIN TEMPERATURE + + model_height = stat_dat(id)%hgt(i,j) + +elseif( qty == QTY_SOIL_TEMPERATURE .or. & ! TSLB SOIL TEMPERATURE + qty == QTY_SOIL_MOISTURE .or. & ! SMOIS SOIL MOISTURE + qty == QTY_SOIL_LIQUID_WATER ) then ! SH2O SOIL LIQUID WATER + + model_height = stat_dat(id)%hgt(i,j) - stat_dat(id)%zs(k) + +elseif( qty == QTY_10M_U_WIND_COMPONENT .or. & + qty == QTY_10M_V_WIND_COMPONENT ) then + + model_height = stat_dat(id)%hgt(i,j) + 10.0_r8 + +elseif( qty == QTY_2M_TEMPERATURE .or. & + qty == QTY_2M_POTENTIAL_TEMPERATURE .or. & ! TH2 POT TEMP at 2 M + qty == QTY_2M_SPECIFIC_HUMIDITY ) then ! Q2 QV at 2 M + + model_height = stat_dat(id)%hgt(i,j) + 2.0_r8 + +! If W-grid (on ZNW levels), native to GZ +elseif( on_w_grid(state_id, var_id) ) then + + i1 = get_dart_vector_index(i,j,k, state_id, gz_id) + x_i1 = get_state(i1, state_handle) + + geop = minval(stat_dat(id)%phb(i,j,k)+x_i1)/gravity + model_height = compute_geometric_height(geop, grid(id)%latitude(i, j)) + +! If U-grid, then height is defined between U points, both in horizontal +! and in vertical, so average -- averaging depends on longitude periodicity +elseif( on_u_grid(state_id, var_id) ) then + + if( i == grid(id)%wes ) then + + ! Check to see if periodic in longitude + if ( grid(id)%periodic_x ) then + + ! We are at the seam in longitude, so take first and last mass points + i1 = get_dart_vector_index(i-1,j,k , state_id, gz_id) + i2 = get_dart_vector_index(i-1,j,k+1, state_id, gz_id) + i3 = get_dart_vector_index(1, j,k , state_id, gz_id) + i4 = get_dart_vector_index(1, j,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + +! HK todo what is minval for? Is it just for converting an array to a scalar? + geop = minval(( (stat_dat(id)%phb(i-1,j,k ) + x_i1) & + +(stat_dat(id)%phb(i-1,j,k+1) + x_i2) & + +(stat_dat(id)%phb(1 ,j,k ) + x_i3) & + +(stat_dat(id)%phb(1 ,j,k+1) + x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(i-1,j) & + +grid(id)%latitude(i-1,j) & + +grid(id)%latitude(1 ,j) & + +grid(id)%latitude(1 ,j) ) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + else + + ! If not periodic, then try extrapolating + i1 = get_dart_vector_index(i-1,j,k , state_id, gz_id) + i2 = get_dart_vector_index(i-1,j,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i1 -1, state_handle) + x_i4 = get_state(i2 -1, state_handle) + + + geop = minval(( 3.0_r8*(stat_dat(id)%phb(i-1,j,k )+x_i1) & + +3.0_r8*(stat_dat(id)%phb(i-1,j,k+1)+x_i2) & + -(stat_dat(id)%phb(i-2,j,k )+x_i3) & + -(stat_dat(id)%phb(i-2,j,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( 3.0_r8*grid(id)%latitude(i-1,j) & + +3.0_r8*grid(id)%latitude(i-1,j) & + -grid(id)%latitude(i-2,j) & + -grid(id)%latitude(i-2,j)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + + elseif( i == 1 ) then + + ! Check to see if periodic in longitude + if ( grid(id)%periodic_x ) then + + ! We are at the seam in longitude, so take first and last mass points + off = grid(id)%we + i1 = get_dart_vector_index(i ,j,k ,state_id, gz_id) + i2 = get_dart_vector_index(i ,j,k+1,state_id, gz_id) + i3 = get_dart_vector_index(off,j,k ,state_id, gz_id) + i4 = get_dart_vector_index(off,j,k+1,state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( (stat_dat(id)%phb(i ,j,k ) + x_i1) & + +(stat_dat(id)%phb(i ,j,k+1) + x_i2) & + +(stat_dat(id)%phb(off,j,k ) + x_i3) & + +(stat_dat(id)%phb(off,j,k+1) + x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(i ,j) & + +grid(id)%latitude(i ,j) & + +grid(id)%latitude(off,j) & + +grid(id)%latitude(off,j)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + else + + ! If not periodic, then try extrapolating + i1 = get_dart_vector_index(i,j,k ,state_id, gz_id) + i2 = get_dart_vector_index(i,j,k+1,state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i1 +1, state_handle) + x_i4 = get_state(i2 +1, state_handle) + + + geop = minval(( 3.0_r8*(stat_dat(id)%phb(i ,j,k )+x_i1) & + +3.0_r8*(stat_dat(id)%phb(i ,j,k+1)+x_i2) & + -(stat_dat(id)%phb(i+1,j,k )+x_i3) & + -(stat_dat(id)%phb(i+1,j,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( 3.0_r8*grid(id)%latitude(i ,j) & + +3.0_r8*grid(id)%latitude(i ,j) & + -grid(id)%latitude(i+1,j) & + -grid(id)%latitude(i+1,j)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + + else + + i1 = get_dart_vector_index(i,j,k ,state_id, gz_id) + i2 = get_dart_vector_index(i,j,k+1,state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i1 -1, state_handle) + x_i4 = get_state(i2 -1, state_handle) + + + geop = minval(( (stat_dat(id)%phb(i ,j,k )+x_i1) & + +(stat_dat(id)%phb(i ,j,k+1)+x_i2) & + +(stat_dat(id)%phb(i-1,j,k )+x_i3) & + +(stat_dat(id)%phb(i-1,j,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(i ,j) & + +grid(id)%latitude(i ,j) & + +grid(id)%latitude(i-1,j) & + +grid(id)%latitude(i-1,j)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + +! If V-grid, then pressure is defined between V points, both in horizontal +! and in vertical, so average -- averaging depends on polar periodicity +elseif( on_v_grid(state_id, var_id) ) then + + if( j == grid(id)%sns ) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The upper corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + i1 = get_dart_vector_index(off,j-1,k , state_id, gz_id) + i2 = get_dart_vector_index(off,j-1,k+1, state_id, gz_id) + i3 = get_dart_vector_index(i ,j-1,k , state_id, gz_id) + i4 = get_dart_vector_index(i ,j-1,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( (stat_dat(id)%phb(off,j-1,k )+x_i1) & + +(stat_dat(id)%phb(off,j-1,k+1)+x_i2) & + +(stat_dat(id)%phb(i ,j-1,k )+x_i3) & + +(stat_dat(id)%phb(i ,j-1,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(off,j-1) & + +grid(id)%latitude(off,j-1) & + +grid(id)%latitude(i ,j-1) & + +grid(id)%latitude(i ,j-1)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + else + + ! If not periodic, then try extrapolating + i1 = get_dart_vector_index(i,j-1,k , state_id, gz_id) + i2 = get_dart_vector_index(i,j-1,k+1, state_id, gz_id) + i3 = get_dart_vector_index(i,j-2,k , state_id, gz_id) + i4 = get_dart_vector_index(i,j-2,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( 3.0_r8*(stat_dat(id)%phb(i,j-1,k )+x_i1) & + +3.0_r8*(stat_dat(id)%phb(i,j-1,k+1)+x_i2) & + -(stat_dat(id)%phb(i,j-2,k )+x_i3) & + -(stat_dat(id)%phb(i,j-2,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( 3.0_r8*grid(id)%latitude(i,j-1) & + +3.0_r8*grid(id)%latitude(i,j-1) & + -grid(id)%latitude(i,j-2) & + -grid(id)%latitude(i,j-2)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + + elseif( j == 1 ) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The lower corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + i1 = get_dart_vector_index(off,j,k , state_id, gz_id) + i2 = get_dart_vector_index(off,j,k+1, state_id, gz_id) + i3 = get_dart_vector_index(i ,j,k , state_id, gz_id) + i4 = get_dart_vector_index(i ,j,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( (stat_dat(id)%phb(off,j,k )+x_i1) & + +(stat_dat(id)%phb(off,j,k+1)+x_i2) & + +(stat_dat(id)%phb(i ,j,k )+x_i3) & + +(stat_dat(id)%phb(i ,j,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(off,j) & + +grid(id)%latitude(off,j) & + +grid(id)%latitude(i ,j) & + +grid(id)%latitude(i ,j)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + else + + ! If not periodic, then try extrapolating + i1 = get_dart_vector_index(i,j ,k , state_id, gz_id) + i2 = get_dart_vector_index(i,j ,k+1, state_id, gz_id) + i3 = get_dart_vector_index(i,j+1,k , state_id, gz_id) + i4 = get_dart_vector_index(i,j+1,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( 3.0_r8*(stat_dat(id)%phb(i,j ,k )+x_i1) & + +3.0_r8*(stat_dat(id)%phb(i,j ,k+1)+x_i2) & + -(stat_dat(id)%phb(i,j+1,k )+x_i3) & + -(stat_dat(id)%phb(i,j+1,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( 3.0_r8*grid(id)%latitude(i,j ) & + +3.0_r8*grid(id)%latitude(i,j ) & + -grid(id)%latitude(i,j+1) & + -grid(id)%latitude(i,j+1)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + + else + + i1 = get_dart_vector_index(i,j ,k , state_id, gz_id) + i2 = get_dart_vector_index(i,j ,k+1, state_id, gz_id) + i3 = get_dart_vector_index(i,j-1,k , state_id, gz_id) + i4 = get_dart_vector_index(i,j-1,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + x_i3 = get_state(i3, state_handle) + x_i4 = get_state(i4, state_handle) + + geop = minval(( (stat_dat(id)%phb(i,j ,k )+x_i1) & + +(stat_dat(id)%phb(i,j ,k+1)+x_i2) & + +(stat_dat(id)%phb(i,j-1,k )+x_i3) & + +(stat_dat(id)%phb(i,j-1,k+1)+x_i4) )/(4.0_r8*gravity)) + + lat = ( grid(id)%latitude(i,j ) & + +grid(id)%latitude(i,j ) & + +grid(id)%latitude(i,j-1) & + +grid(id)%latitude(i,j-1)) / 4.0_r8 + + model_height = compute_geometric_height(geop, lat) + + endif + +else + + i1 = get_dart_vector_index(i,j,k , state_id, gz_id) + i2 = get_dart_vector_index(i,j,k+1, state_id, gz_id) + + x_i1 = get_state(i1, state_handle) + x_i2 = get_state(i2, state_handle) + + geop = minval(( (stat_dat(id)%phb(i,j,k )+x_i1) & + +(stat_dat(id)%phb(i,j,k+1)+x_i2) )/(2.0_r8*gravity)) + + lat = grid(id)%latitude(i,j) + + model_height = compute_geometric_height(geop, lat) + +endif + +end function model_height + +!------------------------------------------------------------------ +pure function interp_pressure(p1, p2, use_log) + +real(r8), intent(in) :: p1(1), p2(1) +logical, intent(in) :: use_log +real(r8) :: interp_pressure + +if (use_log) then + interp_pressure = exp((log(p1(1)) + log(p2(1)))/2.0_r8) +else + interp_pressure = (p1(1) + p2(1))/2.0_r8 +endif + +end function interp_pressure + +!------------------------------------------------------------------ +pure function extrap_pressure(p1, p2, use_log) + +real(r8), intent(in) :: p1(1), p2(1) +logical, intent(in) :: use_log + +real(r8) :: extrap_pressure +real(r8) :: intermediate + +if (use_log) then + intermediate = (3.0_r8*log(p1(1)) - log(p2(1)))/2.0_r8 + if (intermediate <= 0.0_r8) then + extrap_pressure = p1(1) + else + extrap_pressure = exp(intermediate) + endif +else + extrap_pressure = (3.0_r8*p1(1) - p2(1))/2.0_r8 +endif + +end function extrap_pressure + +!------------------------------------------------------------------ +! model pressure any grid (w,u,v,t) +! used in convert_vertical_state so ens_size = 1 +function model_pressure(i, j, kp, id, var_id, state_id, state_handle) + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: i,j,kp,id +integer, intent(in) :: var_id +integer, intent(in) :: state_id +real(r8) :: model_pressure + +real(r8) :: p(1), p_one(1), p_two(1) ! only using the mean, so calling model_pressure_t with ens_size=1 +integer :: k(1), off, n + +k(1) = kp ! array version + +! HK this is for soil variables +do n = 1, get_num_dims(state_id, var_id) + if ( get_dim_name(state_id, var_id, n) == 'soil_layers_stag' ) then + p = model_pressure_s(i, j, id, state_handle, 1) + model_pressure = p(1) + return + endif +enddo + +if (get_num_dims(state_id, var_id) == 2) then ! surface qty + p = model_pressure_s(i, j, id, state_handle, 1) + model_pressure = p(1) + return +endif + +if (on_w_grid(state_id, var_id)) then ! average in the vertical + + if (kp==1) then ! on boundary, extrapolate in the vertical + + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(i, j, k+1, id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_vert_interp) + + elseif ( kp==grid(id)%bts ) then ! on boundary, extrapolate in the vertical + + p_one = model_pressure_t(i, j, k-1, id, state_handle, 1) + p_two = model_pressure_t(i, j, k-2, id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_vert_interp) + + else ! interpolate + + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(i, j, k-1,id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_vert_interp) + + endif + +elseif (on_u_grid(state_id, var_id)) then ! average in the horizontal u direction + + if (i==grid(id)%wes) then + + if ( grid(id)%periodic_x ) then + + ! We are at seam in longitude, take first and last M-grid points + p_one = model_pressure_t(i-1,j, k, id, state_handle, 1) + p_two = model_pressure_t(1, j, k, id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + + else + + ! If not periodic, then try extrapolating + p_one = model_pressure_t(i-1, j, k, id, state_handle, 1) + p_two = model_pressure_t(i-2, j, k, id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_horz_interpM) + + endif + + elseif (i==1) then + + if ( grid(id)%periodic_x ) then + + ! We are at seam in longitude, take first and last M-grid points + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(grid(id)%we, j, k, id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + else + + ! If not periodic, then try extrapolating + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(i+1, j, k, id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_horz_interpM) + + endif + + else + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(i-1, j, k, id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + endif + +elseif (on_v_grid(state_id, var_id)) then ! average in the horizontal v direction + + if (j==grid(id)%sns) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The upper corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + p_one = model_pressure_t(off, j-1, k, id, state_handle, 1) + p_two = model_pressure_t(i, j-1, k, id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + + ! If not periodic, then try extrapolating + else + + p_one = model_pressure_t(i,j-1,k,id, state_handle, 1) + p_two = model_pressure_t(i,j-2,k,id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_horz_interpM) + + endif + + elseif (j==1) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The lower corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + p_one = model_pressure_t(off,j,k,id, state_handle, 1) + p_two = model_pressure_t(i, j,k,id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + + ! If not periodic, then try extrapolating + else + p_one = model_pressure_t(i, j, k, id, state_handle, 1) + p_two = model_pressure_t(i, j+1, k, id, state_handle, 1) + model_pressure = extrap_pressure(p_one, p_two, log_horz_interpM) + endif + + else + p_one = model_pressure_t(i, j ,k, id, state_handle, 1) + p_two = model_pressure_t(i, j-1, k, id, state_handle, 1) + model_pressure = interp_pressure(p_one, p_two, log_horz_interpM) + + endif + +else ! on t_grid + + p = model_pressure_t(i, j, k, id, state_handle, 1) + model_pressure = p(1) + +endif + +end function model_pressure + +!------------------------------------------------------------------ +! model surface pressure any grid u,v,w,t +function model_surface_pressure(i, j, id, var_id, state_id, state_handle) + +! Calculate the surface pressure at grid point (i,j), domain id. +! The grid is defined according to var_type. + +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: i, j, id +integer, intent(in) :: var_id +integer, intent(in) :: state_id +real(r8) :: model_surface_pressure + +integer :: off +integer :: ens_size ! only using the mean, so calling model_pressure_s with ens_size=1 +real(r8) :: pres1(1), pres2(1) + +ens_size = 1 + +model_surface_pressure = missing_r8 + +! If U-grid, then pressure is defined between U points, so average -- +! averaging depends on longitude periodicity +if( on_u_grid(state_id, var_id) ) then + + if (i==grid(id)%wes) then + + ! Check to see if periodic in longitude + if ( grid(id)%periodic_x ) then + + ! We are at seam in longitude, take first and last M-grid points + pres1 = model_pressure_s(i-1,j,id, state_handle, ens_size) + pres2 = model_pressure_s(1, j,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + else + + ! If not periodic, then try extrapolating + pres1 = model_pressure_s(i-1,j,id, state_handle, ens_size) + pres2 = model_pressure_s(i-2,j,id, state_handle, ens_size) + model_surface_pressure = extrap_pressure(pres1, pres2, log_horz_interpM) + + endif + + elseif( i == 1 ) then + + ! Check to see if periodic in longitude + if ( grid(id)%periodic_x ) then + + ! We are at seam in longitude, take first and last M-grid points + pres1 = model_pressure_s(i, j,id, state_handle, ens_size) + pres2 = model_pressure_s(grid(id)%we,j,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + else + + ! If not periodic, then try extrapolating + pres1 = model_pressure_s(i, j,id, state_handle, ens_size) + pres2 = model_pressure_s(i+1,j,id, state_handle, ens_size) + model_surface_pressure = extrap_pressure(pres1, pres2, log_horz_interpM) + + endif + + else + + pres1 = model_pressure_s(i, j,id, state_handle, ens_size) + pres2 = model_pressure_s(i-1,j,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + endif + +! If V-grid, then pressure is defined between V points, so average -- +! averaging depends on polar periodicity +elseif( on_v_grid(state_id, var_id) ) then + + if (j==grid(id)%sns) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The upper corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + pres1 = model_pressure_s(off,j-1,id, state_handle, ens_size) + pres2 = model_pressure_s(i ,j-1,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + ! If not periodic, then try extrapolating + else + + pres1 = model_pressure_s(i,j-1,id, state_handle, ens_size) + pres2 = model_pressure_s(i,j-2,id, state_handle, ens_size) + model_surface_pressure = extrap_pressure(pres1, pres2, log_horz_interpM) + + endif + + elseif( j == 1 ) then + + ! Check to see if periodic in latitude (polar) + if ( grid(id)%polar ) then + + ! The lower corner is 180 degrees of longitude away + off = i + grid(id)%we/2 + if ( off > grid(id)%we ) off = off - grid(id)%we + + pres1 = model_pressure_s(off,j,id, state_handle, ens_size) + pres2 = model_pressure_s(i, j,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + ! If not periodic, then try extrapolating + else + + pres1 = model_pressure_s(i,j, id, state_handle, ens_size) + pres2 = model_pressure_s(i,j+1,id, state_handle, ens_size) + model_surface_pressure = extrap_pressure(pres1, pres2, log_horz_interpM) + + endif + + else + + pres1 = model_pressure_s(i,j, id, state_handle, ens_size) + pres2 = model_pressure_s(i,j-1,id, state_handle, ens_size) + model_surface_pressure = interp_pressure(pres1, pres2, log_horz_interpM) + + endif + +else + + pres1 = model_pressure_s(i,j,id, state_handle, ens_size) + model_surface_pressure = pres1(1) + +endif + +end function model_surface_pressure + +!------------------------------------------------------------------ +function convert_indices_to_lon_lat_lev(i, j, k, var_id, state_id) + +integer, intent(in) :: i, j, k, var_id, state_id +type(location_type) :: convert_indices_to_lon_lat_lev + +real(r8) :: long, lat, lev +integer :: dom_id + +dom_id = get_wrf_domain(state_id) + +if ( on_u_grid(dom_id, var_id) ) then + long = grid(dom_id)%longitude_u(i,j) + lat = grid(dom_id)%latitude_u(i,j) +elseif ( on_v_grid(dom_id, var_id) ) then + long = grid(dom_id)%longitude_v(i,j) + lat = grid(dom_id)%latitude_v(i,j) +else ! on mass grid + long = grid(dom_id)%longitude(i,j) + lat = grid(dom_id)%latitude(i,j) +endif + +! dart expects longitude [0,360] +do while (long < 0.0_r8) + long = long + 360.0_r8 +end do +do while (long > 360.0_r8) + long = long - 360.0_r8 +end do + + +if ( on_w_grid(state_id, var_id) ) then + lev = real(k) - 0.5_r8 +else + lev = real(k) +endif + +convert_indices_to_lon_lat_lev = set_location(long,lat,lev, VERTISLEVEL) + +end function convert_indices_to_lon_lat_lev + +!------------------------------------------------------------------ +! which grid a variable is on. +! querying dimension here, could do by qty? +!------------------------------------------------------------------ +function on_u_grid(state_id, ivar) +integer, intent(in) :: state_id, ivar +logical :: on_u_grid + +on_u_grid = (get_dim_name(state_id, ivar, 1) == 'west_east_stag') + +end function + +!------------------------------------------------------------------ +function on_v_grid(state_id, ivar) +integer, intent(in) :: state_id, ivar +logical :: on_v_grid + +on_v_grid = (get_dim_name(state_id, ivar, 2) == 'south_north_stag') + +end function + +!------------------------------------------------------------------ +function on_w_grid(state_id, ivar) +integer, intent(in) :: state_id, ivar +logical :: on_w_grid + +if (get_num_dims(state_id, ivar) > 2) then + on_w_grid = (get_dim_name(state_id, ivar, 3) == 'bottom_top_stag') +else + on_w_grid = .false. +endif + +end function on_w_grid + +!------------------------------------------------------------------ +function on_t_grid(state_id, ivar) +integer, intent(in) :: state_id, ivar +logical :: on_t_grid + +on_t_grid = (get_dim_name(state_id, ivar, 1) == 'west_east') .and. & + (get_dim_name(state_id, ivar, 2) == 'south_north') + +end function on_t_grid + + +!------------------------------------------------------------------ +!------------------------------------------------------------------ + +function within_bounds_horizontal(i, j, id, qty_in) + +integer, intent(in) :: i, j +integer, intent(in) :: id +integer, intent(in) :: qty_in +logical :: within_bounds_horizontal + +integer :: var_id, qty + +! Some qtys we can interpolate, but the qty is not in the state. +! using QTY_POTENTIAL_TEMPERATURE because this is on the mass grid. +! internal to model_mod, QTY_TEMPERATURE is QTY_POTENTIAL_TEMPERATURE +select case (qty_in) + + case (QTY_SURFACE_TYPE); qty = QTY_POTENTIAL_TEMPERATURE ! uses land mask XLAND which is static data on mass grid + case (QTY_LANDMASK); qty = QTY_POTENTIAL_TEMPERATURE ! land mask XLAND is static data on mass grid + case (QTY_SURFACE_ELEVATION); qty = QTY_POTENTIAL_TEMPERATURE ! terrain height HGT is static data on mass grid + case (QTY_TEMPERATURE); qty = QTY_POTENTIAL_TEMPERATURE ! QTY_TEMPERATURE is not in the state, use QTY_POTENTIAL_TEMPERATURE + case (QTY_SPECIFIC_HUMIDITY); qty = QTY_VAPOR_MIXING_RATIO ! we use vapor mixing ratio to compute specific humidity + case (QTY_DENSITY); qty = QTY_GEOPOTENTIAL_HEIGHT ! density interpolated from geopotential height and pressure + case default + qty = qty_in +end select + +var_id = get_varid_from_kind(wrf_dom(id), qty) +if (var_id <= 0) then + within_bounds_horizontal = .false. + call error_handler(E_ERR, 'within_bounds_horizontal', & + 'Quantity ' // trim(adjustl(get_name_for_quantity(qty))) // & + ' not in state') + return +endif + +within_bounds_horizontal = (bounds_check_lon(i, id, var_id) .and. bounds_check_lat(j, id, var_id)) + +end function within_bounds_horizontal + +!------------------------------------------------------------------ +function able_to_interpolate_qty(id, qty) + +integer, intent(in) :: id ! grid domain id +integer, intent(in) :: qty + +logical :: able_to_interpolate_qty + +select case (qty) + case (QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT) + able_to_interpolate_qty = qty_in_domain(id, QTY_U_WIND_COMPONENT) .and. & + qty_in_domain(id, QTY_V_WIND_COMPONENT) + + case (QTY_10M_U_WIND_COMPONENT, QTY_10M_V_WIND_COMPONENT) + able_to_interpolate_qty = qty_in_domain(id, QTY_10M_U_WIND_COMPONENT) .and. & + qty_in_domain(id, QTY_10M_V_WIND_COMPONENT) + + case (QTY_DENSITY) + able_to_interpolate_qty = qty_in_domain(id, QTY_GEOPOTENTIAL_HEIGHT) .and. & + qty_in_domain(id, QTY_PRESSURE) + + case (QTY_SURFACE_TYPE) + able_to_interpolate_qty = .true. ! land mask XLAND is static data + + case (QTY_LANDMASK) + able_to_interpolate_qty = .true. ! land mask XLAND is static data + + case (QTY_SURFACE_ELEVATION) + able_to_interpolate_qty = .true. ! terrain height HGT is static data + + case (QTY_TEMPERATURE, QTY_POTENTIAL_TEMPERATURE) + able_to_interpolate_qty = qty_in_domain(id, QTY_POTENTIAL_TEMPERATURE) + + case (QTY_SPECIFIC_HUMIDITY) ! we use vapor mixing ratio to compute specific humidity + able_to_interpolate_qty = qty_in_domain(id, QTY_VAPOR_MIXING_RATIO) + + case default + if (chemistry_separate_file) then + able_to_interpolate_qty = qty_in_domain(id, qty) .or. & + qty_in_domain(id+num_domains, qty) ! chemistry variables stored in separate file + else + able_to_interpolate_qty = qty_in_domain(id, qty) + endif + +end select + + +end function able_to_interpolate_qty + +!------------------------------------------------------------------ +function qty_in_domain(id, qty) + +integer, intent(in) :: id +integer, intent(in) :: qty + +logical :: qty_in_domain + +integer :: varid + +varid = get_varid_from_kind(wrf_dom(id), qty) + +if (varid > 0) then + qty_in_domain = .true. +else + qty_in_domain = .false. +endif + +end function qty_in_domain + +!------------------------------------------------------------------ +! returns closest domain id and horizontal mass point grid points (iloc,jloc) +subroutine get_domain_info(obslon,obslat,id,iloc,jloc,domain_id_start) + +real(r8), intent(in) :: obslon, obslat +integer, intent(out) :: id +real(r8), intent(out) :: iloc, jloc +integer, intent(in), optional :: domain_id_start ! HK this is used in wrf_dart_obs_preprocess.f90 + +integer :: n ! domain to start from + +if (present(domain_id_start)) then + n = domain_id_start +else + n = num_domains +endif + +do id = n, 1, -1 + +! From module_map_utils.f90 +! latlon_to_ij(proj, lat, lon, i, j) +! ij_to_latlon(proj, i, j, lat, lon) +! +! It is incumbent upon the calling routine to determine whether or +! not the values returned are within your domain's bounds. All values +! of i, j, lat, and lon are REAL values. + + call latlon_to_ij(grid(id)%proj,min(max(obslat,-89.9999999_r8),89.9999999_r8),obslon,iloc,jloc) + + if (found_in_domain(id, iloc,jloc)) return + +enddo + +! domain not found +id=0 + +end subroutine get_domain_info + +!------------------------------------------------------------------ +function found_in_domain(id, i,j) +integer, intent(in) :: id +real(r8), intent(in) :: i, j +logical :: found_in_domain + +found_in_domain = .false. + +if (id > 1) then + + found_in_domain = ( i >= 1.0_r8 .and. i <= real(grid(id)%we,r8) .and. & + j >= 1.0_r8 .and. j <= real(grid(id)%sn,r8) ) + +else ! have to check periodic + +! Array bound checking depends on whether periodic or not -- these are +! real-valued indices here, so we cannot use boundsCheck :( + + if ( grid(id)%periodic_x .and. .not. grid(id)%periodic_y ) then + if ( grid(id)%polar ) then + ! Periodic X & M_grid ==> [1 we+1) + ! Periodic Y & M_grid ==> [0.5 sn+0.5] + found_in_domain = ( i >= 1.0_r8 .and. i < real(grid(id)%we,r8)+1.0_r8 .and. & + j >= 0.5_r8 .and. j <= real(grid(id)%sn,r8)+0.5_r8 ) + else + ! Periodic X & M_grid ==> [1 we+1) + ! NOT Periodic Y & M_grid ==> [1 sn] + found_in_domain = ( i >= 1.0_r8 .and. i < real(grid(id)%we,r8)+1.0_r8 .and. & + j >= 1.0_r8 .and. j <= real(grid(id)%sn,r8) ) + endif + + elseif ( grid(id)%periodic_x .and. grid(id)%periodic_y ) then + ! Periodic X & M_grid ==> [1 we+1) + ! Periodic Y & M_grid ==> [1 sn+1] + found_in_domain = ( i >= 1.0_r8 .and. i < real(grid(id)%we,r8)+1.0_r8 .and. & + j >= 1.0_r8 .and. j <= real(grid(id)%sn,r8)+1.0_r8 ) + + else + if ( grid(id)%polar ) then + ! NOT Periodic X & M_grid ==> [1 we] + ! Periodic Y & M_grid ==> [0.5 sn+0.5] + found_in_domain = ( i >= 1.0_r8 .and. i <= real(grid(id)%we,r8) .and. & + j >= 0.5_r8 .and. j <= real(grid(id)%sn,r8)+0.5_r8 ) + else + ! NOT Periodic X & M_grid ==> [1 we] + ! NOT Periodic Y & M_grid ==> [1 sn] + found_in_domain = ( i >= 1.0_r8 .and. i <= real(grid(id)%we,r8) .and. & + j >= 1.0_r8 .and. j <= real(grid(id)%sn,r8) ) + endif + endif +endif + + +end function found_in_domain + +!------------------------------------------------------------------ +subroutine getCorners(i, j, id, qty, ll, ul, lr, ur, rc) + +integer, intent(in) :: i, j, id, qty +integer, dimension(2), intent(out) :: ll, ul, lr, ur ! (x,y) of each corner +integer, intent(out) :: rc + +integer :: var_id + +! set return code to 0, and change this if necessary +rc = 0 + +var_id = get_varid_from_kind(wrf_dom(id), qty) + + +!---------------- +! LOWER LEFT ll +!---------------- +! i and j are the lower left (ll) corner already +! +! NOTE :: once we allow for polar periodicity, the incoming j index could actually +! be 0, which would imply a ll(2) value of 1, with a ll(1) value 180 degrees +! of longitude away from the incoming i index! But we have not included +! this possibility yet. + +! As of 22 Oct 2007, this option is not allowed! +! Note that j = 0 can only happen if we are on the M (or U) wrt to latitude +if ( grid(id)%polar .and. j == 0 .and. .not. restrict_polar ) then + + ! j = 0 should be mapped to j = 1 (ll is on other side of globe) + ll(2) = 1 + + ! Need to map i index 180 degrees away + ll(1) = i + grid(id)%we/2 + + ! Check validity of bounds & adjust by periodicity if necessary + if ( ll(1) > grid(id)%we ) ll(1) = ll(1) - grid(id)%we + + ! We shouldn't be able to get this return code if restrict_polar = .true. + rc = 1 + print*, 'model_mod.f90 :: getCorners :: Tried to do polar bc -- rc = ', rc + +else + + ll(1) = i + ll(2) = j + +endif + + +!---------------- +! LOWER RIGHT lr +!---------------- + +! Most of the time, the lower right (lr) corner will simply be (i+1,j), but we need to check +! Summary of x-direction corners: +! Periodic & M_grid has ind = [1 wes) +! ind = [1 we) ==> ind_p_1 = ind + 1 +! ind = [we wes) ==> ind_p_1 = 1 +! Periodic & U_grid has ind = [1 wes) +! ind = [1 we) ==> ind_p_1 = ind + 1 +! ind = [we wes) ==> ind_p_1 = wes ( keep in mind that U(1) = U(wes) if periodic ) +! NOT Periodic & M_grid has ind = [1 we) +! ind = [1 we-1) ==> ind_p_1 = ind + 1 +! ind = [we-1 we) ==> ind_p_1 = we +! NOT Periodic & U_grid has ind = [1 wes) +! ind = [1 we) ==> ind_p_1 = ind + 1 +! ind = [we wes) ==> ind_p_1 = wes + +if ( grid(id)%periodic_x ) then + + ! Check to see what grid we have, M vs. U + if (on_u_grid(wrf_dom(id), var_id) ) then + ! U-grid is always i+1 -- do this in reference to already adjusted ll points + lr(1) = ll(1) + 1 + lr(2) = ll(2) + else + ! M-grid is i+1 except if we <= ind < wes, in which case it's 1 + if ( i < grid(id)%we ) then + lr(1) = ll(1) + 1 + else + lr(1) = 1 + endif + lr(2) = ll(2) + endif + +else + + ! Regardless of grid, NOT Periodic always has i+1 + lr(1) = ll(1) + 1 + lr(2) = ll(2) + +endif + + +!---------------- +! UPPER LEFT ul +!---------------- + +!** NOTE: For now are disallowing observation locations that occur poleward of the +! first and last M-grid gridpoints. This need not be the case because +! the information should be available for proper interpolation across the +! poles, but it will require more clever thinking. Hopefully this can +! be added in later. + +! Most of the time, the upper left (ul) corner will simply be (i,j+1), but we need to check +! Summary of y-direction corners: +! Periodic & M_grid has ind = [0 sns) *though in practice, [1 sn)* +! ind = [1 sn-1) ==> ind_p_1 = ind + 1 +! ind = [sn-1 sn) ==> ind_p_1 = sn +! Periodic & V_grid has ind = [1 sns) +! ind = [1 sn) ==> ind_p_1 = ind + 1 +! ind = [sn sns) ==> ind_p_1 = sns +! NOT Periodic & M_grid has ind = [1 sn) +! ind = [1 sn-1) ==> ind_p_1 = ind + 1 +! ind = [sn-1 sn) ==> ind_p_1 = sn +! NOT Periodic & V_grid has ind = [1 sns) +! ind = [1 sn) ==> ind_p_1 = ind + 1 +! ind = [sn sns) ==> ind_p_1 = sns +! +! Hence, with our current polar obs restrictions, all four possible cases DO map into +! ul = (i,j+1). But this will not always be the case. + +if ( grid(id)%polar ) then + + ! Check to see what grid we have, M vs. V + if ( on_v_grid(wrf_dom(id), var_id) ) then + ! V-grid is always j+1, even if we allow for full [1 sns) range + ul(1) = ll(1) + ul(2) = ll(2) + 1 + else + ! M-grid changes depending on polar restriction + if ( restrict_polar ) then + ! If restricted, then we can simply add 1 + ul(1) = ll(1) + ul(2) = ll(2) + 1 + else + ! If not restricted, then we can potentially wrap over the north pole, which + ! means that ul(2) is set to sn and ul(1) is shifted by 180 deg. + + if ( j == grid(id)%sn ) then + ! j > sn should be mapped to j = sn (ul is on other side of globe) + ul(2) = grid(id)%sn + + ! Need to map i index 180 degrees away + ul(1) = ll(1) + grid(id)%we/2 + + ! Check validity of bounds & adjust by periodicity if necessary + if ( ul(1) > grid(id)%we ) ul(1) = ul(1) - grid(id)%we + + ! We shouldn't be able to get this return code if restrict_polar = .true. + rc = 1 + print*, 'model_mod.f90 :: getCorners :: Tried to do polar bc -- rc = ', rc + + elseif ( j == 0 ) then + ! In this case, we have place ll on the other side of the globe, so we + ! cannot reference ul to ll + ul(1) = i + ul(2) = 1 + + else + ! We can confidently set to j+1 + ul(1) = ll(1) + ul(2) = ll(2) + 1 + endif + + endif + endif + +elseif ( grid(id)%periodic_y ) then + + ! Check to see what grid we have, M vs. V + if ( on_v_grid(wrf_dom(id), var_id) ) then + ! V-grid is always j+1 -- do this in reference to already adjusted ll points + ul(1) = ll(1) + ul(2) = ll(2)+1 + else + ! M-grid is j+1 except if we <= ind < wes, in which case it's 1 + if ( j < grid(id)%sn ) then + ul(2) = ll(2) + 1 + else + ul(2) = 1 + endif + ul(1) = ll(1) + endif + +else + + ! Regardless of grid, NOT Periodic always has j+1 + ul(1) = ll(1) + ul(2) = ll(2) + 1 + +endif + + +!---------------- +! UPPER RIGHT ur +!---------------- + +!*** NOTE: For now are disallowing observation locations that occur poleward of the +! first and last M-grid gridpoints. This need not be the case because +! the information should be available for proper interpolation across the +! poles, but it will require more clever thinking. Hopefully this can +! be added in later. + +! Most of the time, the upper right (ur) corner will simply be (i+1,j+1), but we need to check +! In fact, we can largely get away with ur = (lr(1),ul(2)). Where this will NOT work is +! where we have had to re-map the i index to the other side of the globe (180 deg) due to +! the polar boundary condition. There are no situations where ur(2) will not be equal to +! ul(2). + +ur(2) = ul(2) + +! Need to check if ur(1) .ne. lr(1) +if ( grid(id)%polar .and. .not. restrict_polar ) then + + ! Only if j == 0 or j == sn + if ( j == 0 .or. j == grid(id)%sn) then + ! j == 0 means that ll(1) = i + 180 deg, so we cannot use lr(1) -- hence, we will + ! add 1 to ul(1), unless doing so spans the longitude seam point. + ! j == sn means that ul(1) = i + 180 deg. Here we cannot use lr(1) either because + ! it will be half a domain away from ul(1)+1. Be careful of longitude seam point. + + ! Here we need to check longitude periodicity and the type of grid + if ( grid(id)%periodic_x ) then + + ! Check to see what grid we have, M vs. U + if ( on_u_grid(wrf_dom(id), var_id) ) then + ! U-grid is always i+1 -- do this in reference to already adjusted ll points + ur(1) = ul(1) + 1 + else + ! M-grid is i+1 except if we <= ind < wes, in which case it's 1 + if ( ul(1) < grid(id)%we ) then + ur(1) = ul(1) + 1 + else + ur(1) = 1 + endif + endif + + else + + ! Regardless of grid, NOT Periodic always has i+1 + ur(1) = ul(1) + 1 + + endif + + ! If not a special j value, then we are set for the ur(1) = lr(1) + else + + ur(1) = lr(1) + + endif + +! If not an unrestricted polar periodic domain, then we have nothing to worry about +else + + ur(1) = lr(1) + +endif + +end subroutine getCorners + + +!------------------------------------------------------------------ +! Transfer obs. x to grid j and calculate its +! distance to grid j and j+1 +subroutine toGrid (x, j, dx, dxm) + +real(r8), intent(in) :: x +real(r8), intent(out) :: dx, dxm +integer, intent(out) :: j + +j = int(x) +dx = x - real(j) ! HK todo this should be real(j, KIND=r8)? +dxm= 1.0_r8 - dx + +end subroutine toGrid + +!------------------------------------------------------------------ +subroutine setup_map_projection(id) + +integer, intent(in) :: id +logical, parameter :: debug = .false. + +real(r8) :: latinc,loninc,stdlon +real(r8) :: truelat1, truelat2 + +call map_init(grid(id)%proj) + +! Populate the map projection structure + +!nc -- added in case structures for CASSINI and CYL +!nc -- global wrfinput_d01 has truelat1 = 1.e20, so we need to change it where needed +!nc -- for PROJ_LATLON stdlon and truelat1 have different meanings -- +!nc -- stdlon --> loninc and truelat1 --> latinc +!JPH --- this latinc/loninc calculations are only valid for global domains + +!if ( wrf%dom(id)%scm ) then +!! JPH -- set to zero which should cause the map utils to return NaNs if called +! latinc = 0.0_r8 +! loninc = 0.0_r8 +!else +! latinc = 180.0_r8/wrf%dom(id)%sn +! loninc = 360.0_r8/wrf%dom(id)%we +!endif + +latinc = 180.0_r8/size(grid(id)%longitude(:,1)) ! west_east +loninc = 360.0_r8/size(grid(id)%longitude(1,:)) ! north_south + +if(grid(id)%map_proj == PROJ_LATLON) then !HK why are these different to the wrfinput file? + truelat1 = latinc + stdlon = loninc +else + truelat1 = grid(id)%truelat1 + truelat2 = grid(id)%truelat2 + stdlon = grid(id)%stand_lon +endif + +!nc -- specified inputs to hopefully handle ALL map projections -- hopefully map_set will +! just ignore the inputs it doesn't need for its map projection of interest (?) +! +! NOTE:: We are NOT yet supporting the Gaussian grid or the Rotated Lat/Lon, so we +! are going to skip the entries: nlon, nlat, ixdim, jydim, stagger, phi, lambda +! +! + Gaussian grid uses nlat & nlon +! + Rotated Lat/Lon uses ixdim, jydim, stagger, phi, & lambda +! +call map_set( proj_code=grid(id)%map_proj, & + proj=grid(id)%proj, & + lat1=grid(id)%latitude(1,1), & + lon1=grid(id)%longitude(1,1), & + lat0=90.0_r8, & + lon0=0.0_r8, & + knowni=1.0_r8, & + knownj=1.0_r8, & + dx=grid(id)%dx, & + latinc=latinc, & + loninc=loninc, & + stdlon=stdlon, & + truelat1=truelat1, & + truelat2=truelat2 ) + +end subroutine setup_map_projection + + +!------------------------------------------------------------------ +!------------------------------------------------------------------ +! Bounds checking routines: + ! Summary of Allowable REAL-VALUED Index Values ==> INTEGER Index Values + ! + ! In longitude (x) direction + ! Periodic & M_grid ==> [1 we+1) ==> [1 wes) + ! Periodic & U_grid ==> [1 wes) ==> [1 wes) + ! NOT Periodic & M_grid ==> [1 we] ==> [1 we) + ! NOT Periodic & U_grid ==> [1.5 we+0.5] ==> [1 wes) + ! In latitude (y) direction + ! Periodic & M_grid ==> [0.5 sn+0.5] ==> [0 sns) *though in practice, [1 sn)* + ! Periodic & V_grid ==> [1 sns] ==> [1 sns) *though allowable range, [1.5 sn+.5]* + ! NOT Periodic & M_grid ==> [1 sn] ==> [1 sn) + ! NOT Periodic & V_grid ==> [1.5 sn+0.5] ==> [1 sns) + ! In vertical (z) direction + ! M_grid ==> [1 bt] ==> [1 bt) + ! W_grid ==> [1.5 bt+0.5] ==> [1 bts) + +!------------------------------------------------------------------ + +!------------------------------------------------------------------ +! Determines whether real-valued location indices are +! within a sensible range based on the assumed (un)staggered grid and based on +! whether the domain is assumed to be periodic in a given direction. +function bounds_check_lon(ind, id, var_id) ! or variable? + +integer, intent(in) :: ind ! index into the grid +integer, intent(in) :: id ! domain id +integer, intent(in) :: var_id + +logical :: bounds_check_lon + +! Consider cases in REAL-VALUED indexing: +! +! I. Longitude -- x-direction +! A. PERIODIC (period_x = .true.) +! +! Consider Mass-grid (& V-grid) longitude grid with 4 west-east gridpoints +! Values :: [ -135 -45 45 135 ] .. {225} +! Indices :: [ 1 2 3 4 ] .. {1,5} +! Complementary U-grid +! Values :: [ -180 -90 0 90 180 ] +! Indices :: [ 1 2 3 4 5 ] +! +! What are the allowable values for a real-valued index on each of these grids? +! 1. M-grid ---> [1 5) ---> [1 we+1) +! ---> [-135 225) +! 2. U-grid ---> [1 5) ---> [1 wes) +! ---> [-180 180) +! [Note that above "allowable values" reflect that one should be able to have +! an observation anywhere on a given longitude circle -- the information +! exists in order to successfully interpolate to anywhere over [0 360).] +! +! It is up to the routine calling "boundsCheck" to have handled the 0.5 offset +! in indices between the M-grid & U-grid. Hence, two examples: +! a. If there is an observation location at -165 longitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 4.667 +! * An observation of TYPE_U (on the U-grid) would have ind = 1.167 +! b. If there is an observation location at 0 longitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 2.5 +! * An observation of TYPE_U (on the U-grid) would have ind = 3.0 +! +! B. NOT periodic (period_x = .false.) +! +! Consider Mass-grid (& V-grid) longitude grid with 4 west-east gridpoints +! Values :: [ 95 105 115 125 ] +! Indices :: [ 1 2 3 4 ] +! Complementary U-grid +! Values :: [ 90 100 110 120 130 ] +! Indices :: [ 1 2 3 4 5 ] +! +! What are the allowable values for a real-valued index on each of these grids? +! 1. M-grid ---> [1 4] ---> [1 we] +! ---> [95 125] +! 2. U-grid ---> [1.5 4.5] ---> [1.5 we+0.5] +! ---> [95 125] +! [Note that above "allowable values" reflect that one should only be able to +! have an observation within the M-grid, since that is the only way to +! guarantee that the necessary information exists in order to successfully +! interpolate to a specified location.] +! +! It is up to the routine calling "boundsCheck" to have handled the 0.5 offset +! in indices between the M-grid & U-grid. Hence, two examples: +! a. If there is an observation location at 96 longitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 1.1 +! * An observation of TYPE_U (on the U-grid) would have ind = 1.6 +! b. If there is an observation location at 124 longitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 3.9 +! * An observation of TYPE_U (on the U-grid) would have ind = 4.4 +! + +bounds_check_lon = .false. + +! Next check periodicity +if ( grid(id)%periodic_x ) then + + ! If periodic in longitude, then no need to check staggering because both + ! M and U grids allow integer indices from [1 wes) + if ( ind >= 1 .and. ind < grid(id)%wes ) bounds_check_lon = .true. + +else + + ! If NOT periodic in longitude, then we need to check staggering because + ! M and U grids allow different index ranges + + ! Check staggering by comparing var_size(dim,type) to the staggered dimension + if ( on_u_grid(wrf_dom(id),var_id) ) then + ! U-grid allows integer range of [1 wes) + if ( ind >= 1 .and. ind < grid(id)%wes ) bounds_check_lon = .true. + else + ! M & V-grid allow [1 we) + if ( ind >= 1 .and. ind < grid(id)%we ) bounds_check_lon = .true. + endif + +endif + +end function bounds_check_lon + +!------------------------------------------------------------------ +function bounds_check_lat(ind, id, var_id) + +integer, intent(in) :: ind ! index into the grid +integer, intent(in) :: id ! domain id +integer, intent(in) :: var_id + +logical :: bounds_check_lat +! II. Latitude -- y-direction +! A. PERIODIC (polar = .true.) +! +! Consider Mass-grid (& U-Grid) latitude grid with 4 south-north gridpoints +! Values :: [ -67.5 -22.5 22.5 67.5 ] +! Indices :: [ 1 2 3 4 ] +! Complementary V-grid +! Values :: [ -90 -45 0 45 90 ] +! Indices :: [ 1 2 3 4 5 ] +! +! What are the allowable values for a real-valued index on each of these grids? +! 1. M-grid ---> [0.5 4.5] ---> [0.5 sn+0.5] +! ---> [-90 90] +! 2. U-grid ---> [1 5] ---> [1 sns] +! ---> [-90 90] +! [Note that above "allowable values" reflect that one should be able to have +! an observation anywhere along a give latitude circle -- the information +! exists in order to successfully interpolate to anywhere over [-90 90]; +! however, in latitude this poses a special challenge since the seams join +! two separate columns of data over the pole, as opposed to in longitude +! where the seam wraps back on a single row of data.] +! +! It is up to the routine calling "boundsCheck" to have handled the 0.5 offset +! in indices between the M-grid & V-grid. Hence, two examples: +! a. If there is an observation location at -75 latitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 0.833 +! * An observation of TYPE_V (on the V-grid) would have ind = 1.333 +! b. If there is an observation location at 0 latitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 2.5 +! * An observation of TYPE_V (on the V-grid) would have ind = 3.0 +! +! B. NOT periodic (polar = .false.) +! +! Consider Mass-grid (& U-Grid) latitude grid with 4 south-north gridpoints +! Values :: [ 10 20 30 40 ] +! Indices :: [ 1 2 3 4 ] +! Complementary V-grid +! Values :: [ 5 15 25 35 45 ] +! Indices :: [ 1 2 3 4 5 ] +! +! What are the allowable values for a real-valued index on each of these grids? +! 1. M-grid ---> [1 4] ---> [1 sn] +! ---> [10 40] +! 2. U-grid ---> [1.5 4.5] ---> [1.5 sn+0.5] +! ---> [10 40] +! [Note that above "allowable values" reflect that one should only be able to +! have an observation within the M-grid, since that is the only way to +! guarantee that the necessary information exists in order to successfully +! interpolate to a specified location.] +! +! It is up to the routine calling "boundsCheck" to have handled the 0.5 offset +! in indices between the M-grid & V-grid. Hence, two examples: +! a. If there is an observation location at 11 latitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 1.1 +! * An observation of TYPE_V (on the V-grid) would have ind = 1.6 +! b. If there is an observation location at 25 latitude, then: +! * An observation of TYPE_T (on the M-grid) would have ind = 2.5 +! * An observation of TYPE_V (on the V-grid) would have ind = 3.0 +! + +if ( grid(id)%periodic_y ) then + + ! We need to check staggering because M and V grids allow different indices + +!*** NOTE: For now are disallowing observation locations that occur poleward of the +! first and last M-grid gridpoints. This means that this function will +! return false for polar observations. This need not be the case because +! the information should be available for proper interpolation across the +! poles, but it will require more clever thinking. Hopefully this can +! be added in later. + + ! Check staggering by comparing var_size(dim,type) to the staggered dimension + if ( on_v_grid(wrf_dom(id), var_id) ) then + ! V-grid allows integer range [1 sns) + if ( ind >= 1 .and. ind < grid(id)%sns ) bounds_check_lat = .true. + else + ! For now we will set a logical flag to more restrictively check the array + ! bounds under our no-polar-obs assumptions + if ( restrict_polar ) then + ! M & U-grid allow integer range [1 sn) in practice (though properly, [0 sns) ) + if ( ind >= 1 .and. ind < grid(id)%sn ) bounds_check_lat = .true. + else + ! M & U-grid allow integer range [0 sns) in unrestricted circumstances + if ( ind >= 0 .and. ind < grid(id)%sns ) bounds_check_lat = .true. + endif + endif + + else + + ! We need to check staggering because M and V grids allow different indices + if ( on_v_grid(wrf_dom(id), var_id) ) then + ! V-grid allows [1 sns) + if ( ind >= 1 .and. ind < grid(id)%sns ) bounds_check_lat = .true. + else + ! M & U-grid allow [1 sn) + if ( ind >= 1 .and. ind < grid(id)%sn ) bounds_check_lat = .true. + endif + + endif + +end function bounds_check_lat + + +!------------------------------------------------------------------ +function bounds_check_vertical(ind, id, var_id) + +integer, intent(in) :: ind ! index into the grid +integer, intent(in) :: id ! domain id +integer, intent(in) :: var_id + +logical :: bounds_check_vertical +! III. Vertical -- z-direction (periodicity not an issue) +! +! Consider Mass vertical grid with 4 bottom-top gridpoints +! Values :: [ 0.875 0.625 0.375 0.125 ] +! Indices :: [ 1 2 3 4 ] +! Complementary W-grid +! Values :: [ 1 0.75 0.50 0.25 0 ] +! Indices :: [ 1 2 3 4 5 ] +! +! What are the allowable values for a real-valued index on each of these grids? +! 1. M-grid ---> [1 4] ---> [1 bt] +! ---> [0.875 0.125] +! 2. W-grid ---> [1.5 4.5] ---> [1.5 bt+0.5] +! ---> [0.875 0.125] +! +! [Note that above "allowable values" reflect that one should only be able to +! have an observation within the M-grid, since that is the only way to +! guarantee that the necessary information exists in order to successfully +! interpolate to a specified location.] + +bounds_check_vertical = .false. + +if ( on_w_grid(wrf_dom(id), var_id) ) then + ! W vertical grid allows [1 bts) + if ( ind >= 1 .and. ind < grid(id)%bts ) bounds_check_vertical = .true. +else + ! M vertical grid allows [1 bt) + if ( ind >= 1 .and. ind < grid(id)%bt ) bounds_check_vertical = .true. +endif + +end function bounds_check_vertical + +!------------------------------------------------------------------ +! Summary of Allowable REAL-VALUED Index Values ==> INTEGER Index Values +! +! In longitude (x) direction +! Periodic & M_grid ==> [1 we+1) ==> [1 wes) +! Periodic & U_grid ==> [1 wes) ==> [1 wes) +! NOT Periodic & M_grid ==> [1 we] ==> [1 we) +! NOT Periodic & U_grid ==> [1.5 we+0.5] ==> [1 wes) +! In latitude (y) direction +! Periodic & M_grid ==> [0.5 sn+0.5] ==> [0 sns) *though in practice, [1 sn)* +! Periodic & V_grid ==> [1 sns] ==> [1 sns) *though allowable range, [1.5 sn+.5]* +! NOT Periodic & M_grid ==> [1 sn] ==> [1 sn) +! NOT Periodic & V_grid ==> [1.5 sn+0.5] ==> [1 sns) +! In vertical (z) direction +! M_grid ==> [1 bt] ==> [1 bt) +! W_grid ==> [1.5 bt+0.5] ==> [1 bts) +!------------------------------------------------------------------ + +!------------------------------------------------------------------ +! Vortex + +subroutine vortex() + +print*, 'Do vortex' + +end subroutine vortex + +end module model_mod diff --git a/models/wrf_unified/model_mod.nml b/models/wrf_unified/model_mod.nml new file mode 100644 index 000000000..c41e23f81 --- /dev/null +++ b/models/wrf_unified/model_mod.nml @@ -0,0 +1,14 @@ +&model_nml + wrf_state_variables = 'NULL' + wrf_state_bounds = 'NULL' + num_domains = 1 + calendar_type = 3 + assimilation_period_seconds = 21600 + allow_obs_below_vol = .false. + vert_localization_coord = 3 + sfc_elev_max_diff = -1.0 + polar = .false. + periodic_x = .false. + periodic_y = .false. + allow_perturbed_ics = .false. + / diff --git a/models/wrf_unified/module_map_utils.f90 b/models/wrf_unified/module_map_utils.f90 new file mode 100644 index 000000000..f874e7a9e --- /dev/null +++ b/models/wrf_unified/module_map_utils.f90 @@ -0,0 +1,2348 @@ +! This code is not protected by the DART copyright agreement. + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! MODULE CONSTANTS_MODULE +! +! This module defines constants that are used by other modules +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +module constants_module + + use types_mod, only : r8, digits12 + + real (kind=r8), parameter :: PI = 3.141592653589793_r8 + real (kind=r8), parameter :: OMEGA_E = 0.00007292_r8 ! Angular rotation rate of the earth + + real (kind=r8), parameter :: DEG_PER_RAD = 180.0_r8 / PI + real (kind=r8), parameter :: RAD_PER_DEG = PI / 180.0_r8 + + ! Mean Earth Radius in m. The value below is consistent + ! with NCEP's routines and grids. + real (kind=r8), parameter :: A_WGS84 = 6378137.0_r8 + real (kind=r8), parameter :: B_WGS84 = 6356752.314_r8 + real (kind=r8), parameter :: RE_WGS84 = A_WGS84 + real (kind=r8), parameter :: E_WGS84 = 0.081819192_r8 + + real (kind=r8), parameter :: A_NAD83 = 6378137.0_r8 + real (kind=r8), parameter :: RE_NAD83 = A_NAD83 + real (kind=r8), parameter :: E_NAD83 = 0.0818187034_r8 + + real (kind=r8), parameter :: EARTH_RADIUS_M = 6370000.0_r8 ! same as MM5 system + real (kind=r8), parameter :: EARTH_CIRC_M = 2.0_r8 * PI * EARTH_RADIUS_M + +end module constants_module + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! MODULE MISC_DEFINITIONS_MODULE +! +! This module defines various non-meteorological constants that are used +! by other modules for readability. +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +module misc_definitions_module + + use types_mod, only : r8 + + real (kind=r8), parameter :: NAN=1.E20_r8 + + real (kind=r8), parameter :: NOT_MASKED = -2.0_r8, & + MASKED_BOTH = -1.0_r8, & + MASKED_WATER = 0.0_r8, & + MASKED_LAND = 1.0_r8 + + integer, parameter :: OUTSIDE_DOMAIN=1E8, NOT_PROCESSED=1E9, INVALID=1E9 + + integer, parameter :: SIXTEEN_POINT=1, FOUR_POINT=2, N_NEIGHBOR=3, & + AVERAGE4=4, AVERAGE16=5, W_AVERAGE4=6, W_AVERAGE16=7, & + SEARCH=8 + + integer, parameter :: BOTTOM_TOP=1, TOP_BOTTOM=2 + + integer, parameter :: CONTINUOUS=0, CATEGORICAL=1, SP_CONTINUOUS=2 + + integer, parameter :: M=1, U=2, V=3, HH=4, VV=5 + + integer, parameter :: ONETWOONE=1, SMTHDESMTH=2, SMTHDESMTH_SPECIAL=3 + + integer, parameter :: BINARY=1, NETCDF=2, GRIB1=3, HDF=4 + + integer, parameter :: BIG_ENDIAN=0, LITTLE_ENDIAN=1 + + ! Projection codes for proj_info structure: + INTEGER, PUBLIC, PARAMETER :: PROJ_LATLON = 0 + INTEGER, PUBLIC, PARAMETER :: PROJ_LC = 1 + INTEGER, PUBLIC, PARAMETER :: PROJ_PS = 2 + INTEGER, PUBLIC, PARAMETER :: PROJ_PS_WGS84 = 102 + INTEGER, PUBLIC, PARAMETER :: PROJ_MERC = 3 + INTEGER, PUBLIC, PARAMETER :: PROJ_GAUSS = 4 + INTEGER, PUBLIC, PARAMETER :: PROJ_CYL = 5 + INTEGER, PUBLIC, PARAMETER :: PROJ_CASSINI = 6 + INTEGER, PUBLIC, PARAMETER :: PROJ_ALBERS_NAD83 = 105 + INTEGER, PUBLIC, PARAMETER :: PROJ_ROTLL = 203 + +end module misc_definitions_module + + +MODULE map_utils + +! Module that defines constants, data structures, and +! subroutines used to convert grid indices to lat/lon +! and vice versa. +! +! SUPPORTED PROJECTIONS +! --------------------- +! Cylindrical Lat/Lon (code = PROJ_LATLON) +! Mercator (code = PROJ_MERC) +! Lambert Conformal (code = PROJ_LC) +! Gaussian (code = PROJ_GAUSS) +! Polar Stereographic (code = PROJ_PS) +! Rotated Lat/Lon (code = PROJ_ROTLL) +! +! REMARKS +! ------- +! The routines contained within were adapted from routines +! obtained from NCEP's w3 library. The original NCEP routines were less +! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S) +! than what we needed, so modifications based on equations in Hoke, Hayes, and +! Renninger (AFGWC/TN/79-003) were added to improve the flexibility. +! Additionally, coding was improved to F90 standards and the routines were +! combined into this module. +! +! ASSUMPTIONS +! ----------- +! Grid Definition: +! For mercator, lambert conformal, and polar-stereographic projections, +! the routines within assume the following: +! +! 1. Grid is dimensioned (i,j) where i is the East-West direction, +! positive toward the east, and j is the north-south direction, +! positive toward the north. +! 2. Origin is at (1,1) and is located at the southwest corner, +! regardless of hemispere. +! 3. Grid spacing (dx) is always positive. +! 4. Values of true latitudes must be positive for NH domains +! and negative for SH domains. +! +! For the latlon and Gaussian projection, the grid origin may be at any +! of the corners, and the deltalat and deltalon values can be signed to +! account for this using the following convention: +! Origin Location Deltalat Sign Deltalon Sign +! --------------- ------------- ------------- +! SW Corner + + +! NE Corner - - +! NW Corner - + +! SE Corner + - +! +! Data Definitions: +! 1. Any arguments that are a latitude value are expressed in +! degrees north with a valid range of -90 -> 90 +! 2. Any arguments that are a longitude value are expressed in +! degrees east with a valid range of -180 -> 180. +! 3. Distances are in meters and are always positive. +! 4. The standard longitude (stdlon) is defined as the longitude +! line which is parallel to the grid's y-axis (j-direction), along +! which latitude increases (NOT the absolute value of latitude, but +! the actual latitude, such that latitude increases continuously +! from the south pole to the north pole) as j increases. +! 5. One true latitude value is required for polar-stereographic and +! mercator projections, and defines at which latitude the +! grid spacing is true. For lambert conformal, two true latitude +! values must be specified, but may be set equal to each other to +! specify a tangent projection instead of a secant projection. +! +! USAGE +! ----- +! To use the routines in this module, the calling routines must have the +! following statement at the beginning of its declaration block: +! USE map_utils +! +! The use of the module not only provides access to the necessary routines, +! but also defines a structure of TYPE (proj_info) that can be used +! to declare a variable of the same type to hold your map projection +! information. It also defines some integer parameters that contain +! the projection codes so one only has to use those variable names rather +! than remembering the acutal code when using them. The basic steps are +! as follows: +! +! 1. Ensure the "USE map_utils" is in your declarations. +! 2. Declare the projection information structure as type(proj_info): +! TYPE(proj_info) :: proj +! 3. Populate your structure by calling the map_set routine: +! CALL map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj) +! where: +! code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS, +! PROJ_GAUSS, or PROJ_ROTLL +! lat1 (input) = Latitude of grid origin point (i,j)=(1,1) +! (see assumptions!) +! lon1 (input) = Longitude of grid origin +! knowni (input) = origin point, x-location +! knownj (input) = origin point, y-location +! dx (input) = grid spacing in meters (ignored for LATLON projections) +! stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC, +! deltalon (see assumptions) for PROJ_LATLON, +! ignored for PROJ_MERC +! truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and +! PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON +! truelat2 (input) = 2nd true latitude for PROJ_LC, +! ignored for all others. +! proj (output) = The structure of type (proj_info) that will be fully +! populated after this call +! +! 4. Now that the proj structure is populated, you may call either +! of the following routines: +! +! latlon_to_ij(proj, lat, lon, i, j) +! ij_to_latlon(proj, i, j, lat, lon) +! +! It is incumbent upon the calling routine to determine whether or +! not the values returned are within your domain's bounds. All values +! of i, j, lat, and lon are REAL values. +! +! +! REFERENCES +! ---------- +! Hoke, Hayes, and Renninger, "Map Projections and Grid Systems for +! Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather +! Service, 1985. +! +! NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F +! NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12 +! +! HISTORY +! ------- +! 27 Mar 2001 - Original Version +! Brent L. Shaw, NOAA/FSL (CSU/CIRA) +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! other modules included in this .f90 file + use constants_module + use misc_definitions_module + + ! external modules + use types_mod, only : r8 +! use utilities_mod, only : register_module + + ! Define some private constants + INTEGER, PRIVATE, PARAMETER :: HIGH = digits12 + + TYPE proj_info + + INTEGER :: code ! Integer code for projection TYPE + INTEGER :: nlat ! For Gaussian -- number of latitude points + ! north of the equator + INTEGER :: nlon ! + ! + INTEGER :: ixdim ! For Rotated Lat/Lon -- number of mass points + ! in an odd row + INTEGER :: jydim ! For Rotated Lat/Lon -- number of rows + INTEGER :: stagger ! For Rotated Lat/Lon -- mass or velocity grid + REAL(r8) :: phi ! For Rotated Lat/Lon -- domain half-extent in + ! degrees latitude + REAL(r8) :: lambda ! For Rotated Lat/Lon -- domain half-extend in + ! degrees longitude + REAL(r8) :: lat1 ! SW latitude (1,1) in degrees (-90->90N) + REAL(r8) :: lon1 ! SW longitude (1,1) in degrees (-180->180E) + REAL(r8) :: lat0 ! For Cassini, latitude of projection pole + REAL(r8) :: lon0 ! For Cassini, longitude of projection pole + REAL(r8) :: dx ! Grid spacing in meters at truelats, used + ! only for ps, lc, and merc projections + REAL(r8) :: dy ! Grid spacing in meters at truelats, used + ! only for ps, lc, and merc projections + REAL(r8) :: latinc ! Latitude increment for cylindrical lat/lon + REAL(r8) :: loninc ! Longitude increment for cylindrical lat/lon + ! also the lon increment for Gaussian grid + REAL(r8) :: dlat ! Lat increment for lat/lon grids + REAL(r8) :: dlon ! Lon increment for lat/lon grids + REAL(r8) :: stdlon ! Longitude parallel to y-axis (-180->180E) + REAL(r8) :: truelat1 ! First true latitude (all projections) + REAL(r8) :: truelat2 ! Second true lat (LC only) + REAL(r8) :: hemi ! 1 for NH, -1 for SH + REAL(r8) :: cone ! Cone factor for LC projections + REAL(r8) :: polei ! Computed i-location of pole point + REAL(r8) :: polej ! Computed j-location of pole point + REAL(r8) :: rsw ! Computed radius to SW corner + REAL(r8) :: rebydx ! Earth radius divided by dx + REAL(r8) :: knowni ! X-location of known lat/lon + REAL(r8) :: knownj ! Y-location of known lat/lon + REAL(r8) :: re_m ! Radius of spherical earth, meters + REAL(r8) :: rho0 ! For Albers equal area + REAL(r8) :: nc ! For Albers equal area + REAL(r8) :: bigc ! For Albers equal area + LOGICAL :: init ! Flag to indicate if this struct is + ! ready for use + LOGICAL :: wrap ! For Gaussian -- flag to indicate wrapping + ! around globe? + REAL(r8), POINTER, DIMENSION(:) :: gauss_lat ! Latitude array for Gaussian grid + + END TYPE proj_info + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + CONTAINS + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE map_init(proj) + ! Initializes the map projection structure to missing values + + IMPLICIT NONE + TYPE(proj_info), INTENT(INOUT) :: proj + + proj%lat1 = -999.9_r8 + proj%lon1 = -999.9_r8 + proj%lat0 = -999.9_r8 + proj%lon0 = -999.9_r8 + proj%dx = -999.9_r8 + proj%dy = -999.9_r8 + proj%latinc = -999.9_r8 + proj%loninc = -999.9_r8 + proj%stdlon = -999.9_r8 + proj%truelat1 = -999.9_r8 + proj%truelat2 = -999.9_r8 + proj%phi = -999.9_r8 + proj%lambda = -999.9_r8 + proj%ixdim = -999 + proj%jydim = -999 + proj%stagger = HH + proj%nlat = 0 + proj%nlon = 0 + proj%hemi = 0.0_r8 + proj%cone = -999.9_r8 + proj%polei = -999.9_r8 + proj%polej = -999.9_r8 + proj%rsw = -999.9_r8 + proj%knowni = -999.9_r8 + proj%knownj = -999.9_r8 + proj%re_m = EARTH_RADIUS_M + proj%init = .FALSE. + proj%wrap = .FALSE. + proj%rho0 = 0.0_r8 + proj%nc = 0.0_r8 + proj%bigc = 0.0_r8 + nullify(proj%gauss_lat) + + END SUBROUTINE map_init + + +!nc -- Global WRF assumes proj_code = PROJ_CASSINI (=6 in misc_definitions_module) + + SUBROUTINE map_set(proj_code, proj, lat1, lon1, lat0, lon0, knowni, knownj, dx, latinc, & + loninc, stdlon, truelat1, truelat2, nlat, nlon, ixdim, jydim, & + stagger, phi, lambda, r_earth) + ! Given a partially filled proj_info structure, this routine computes + ! polei, polej, rsw, and cone (if LC projection) to complete the + ! structure. This allows us to eliminate redundant calculations when + ! calling the coordinate conversion routines multiple times for the + ! same map. + ! This will generally be the first routine called when a user wants + ! to be able to use the coordinate conversion routines, and it + ! will call the appropriate subroutines based on the + ! proj%code which indicates which projection type this is. + + IMPLICIT NONE + + ! Declare arguments + INTEGER, INTENT(IN) :: proj_code + INTEGER, INTENT(IN), OPTIONAL :: nlat + INTEGER, INTENT(IN), OPTIONAL :: nlon + INTEGER, INTENT(IN), OPTIONAL :: ixdim + INTEGER, INTENT(IN), OPTIONAL :: jydim + INTEGER, INTENT(IN), OPTIONAL :: stagger + REAL(r8), INTENT(IN), OPTIONAL :: latinc + REAL(r8), INTENT(IN), OPTIONAL :: loninc + REAL(r8), INTENT(IN), OPTIONAL :: lat1 + REAL(r8), INTENT(IN), OPTIONAL :: lon1 + REAL(r8), INTENT(IN), OPTIONAL :: lat0 + REAL(r8), INTENT(IN), OPTIONAL :: lon0 + REAL(r8), INTENT(IN), OPTIONAL :: dx + REAL(r8), INTENT(IN), OPTIONAL :: stdlon + REAL(r8), INTENT(IN), OPTIONAL :: truelat1 + REAL(r8), INTENT(IN), OPTIONAL :: truelat2 + REAL(r8), INTENT(IN), OPTIONAL :: knowni + REAL(r8), INTENT(IN), OPTIONAL :: knownj + REAL(r8), INTENT(IN), OPTIONAL :: phi + REAL(r8), INTENT(IN), OPTIONAL :: lambda + REAL(r8), INTENT(IN), OPTIONAL :: r_earth + TYPE(proj_info), INTENT(OUT) :: proj + + INTEGER :: iter + REAL(r8) :: dummy_lon1 + REAL(r8) :: dummy_lon0 + REAL(r8) :: dummy_stdlon + + ! First, verify that mandatory parameters are present for the specified proj_code + IF ( proj_code == PROJ_LC ) THEN + IF ( .NOT.PRESENT(truelat1) .OR. & + .NOT.PRESENT(truelat2) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(stdlon) .OR. & + .NOT.PRESENT(dx) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knowni, knownj, stdlon, dx' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_PS ) THEN + IF ( .NOT.PRESENT(truelat1) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(stdlon) .OR. & + .NOT.PRESENT(dx) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_PS_WGS84 ) THEN + IF ( .NOT.PRESENT(truelat1) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(stdlon) .OR. & + .NOT.PRESENT(dx) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_ALBERS_NAD83 ) THEN + IF ( .NOT.PRESENT(truelat1) .OR. & + .NOT.PRESENT(truelat2) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(stdlon) .OR. & + .NOT.PRESENT(dx) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knonwi, knownj, stdlon, dx' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_MERC ) THEN + IF ( .NOT.PRESENT(truelat1) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(dx) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' truelat1, lat1, lon1, knowni, knownj, dx' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_LATLON ) THEN + IF ( .NOT.PRESENT(latinc) .OR. & + .NOT.PRESENT(loninc) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' latinc, loninc, knowni, knownj, lat1, lon1' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_CYL ) THEN + IF ( .NOT.PRESENT(latinc) .OR. & + .NOT.PRESENT(loninc) .OR. & + .NOT.PRESENT(stdlon) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' latinc, loninc, stdlon' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_CASSINI ) THEN + IF ( .NOT.PRESENT(latinc) .OR. & + .NOT.PRESENT(loninc) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(lat0) .OR. & + .NOT.PRESENT(lon0) .OR. & + .NOT.PRESENT(knowni) .OR. & + .NOT.PRESENT(knownj) .OR. & + .NOT.PRESENT(stdlon) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' latinc, loninc, lat1, lon1, knowni, knownj, lat0, lon0, stdlon' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_GAUSS ) THEN + IF ( .NOT.PRESENT(nlat) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(loninc) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' nlat, lat1, lon1, loninc' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE IF ( proj_code == PROJ_ROTLL ) THEN + IF ( .NOT.PRESENT(ixdim) .OR. & + .NOT.PRESENT(jydim) .OR. & + .NOT.PRESENT(phi) .OR. & + .NOT.PRESENT(lambda) .OR. & + .NOT.PRESENT(lat1) .OR. & + .NOT.PRESENT(lon1) .OR. & + .NOT.PRESENT(stagger) ) THEN + PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code + PRINT '(A)', ' ixdim, jydim, phi, lambda, lat1, lon1, stagger' + write(6,*) 'ERROR: MAP_INIT' + END IF + ELSE + PRINT '(A,I2)', 'Unknown projection code: ', proj_code + write(6,*) 'ERROR: MAP_INIT' + END IF + + ! Check for validity of mandatory variables in proj + IF ( PRESENT(lat1) ) THEN + IF ( ABS(lat1) .GT. 90.0_r8 ) THEN + PRINT '(A)', 'Latitude of origin corner required as follows:' + PRINT '(A)', ' -90N <= lat1 < = 90.N' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + + IF ( PRESENT(lon1) ) THEN + dummy_lon1 = lon1 + IF ( ABS(dummy_lon1) .GT. 180.0_r8 ) THEN + iter = 0 + DO WHILE (ABS(dummy_lon1) > 180.0_r8 .AND. iter < 10) + IF (dummy_lon1 < -180.0_r8) dummy_lon1 = dummy_lon1 + 360.0_r8 + IF (dummy_lon1 > 180.0_r8) dummy_lon1 = dummy_lon1 - 360.0_r8 + iter = iter + 1 + END DO + IF (abs(dummy_lon1) > 180.0_r8) THEN + PRINT '(A)', 'Longitude of origin required as follows:' + PRINT '(A)', ' -180E <= lon1 <= 180W' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + ENDIF + + IF ( PRESENT(lon0) ) THEN + dummy_lon0 = lon0 + IF ( ABS(dummy_lon0) .GT. 180.0_r8) THEN + iter = 0 + DO WHILE (ABS(dummy_lon0) > 180.0_r8 .AND. iter < 10) + IF (dummy_lon0 < -180.0_r8) dummy_lon0 = dummy_lon0 + 360.0_r8 + IF (dummy_lon0 > 180.0_r8) dummy_lon0 = dummy_lon0 - 360.0_r8 + iter = iter + 1 + END DO + IF (abs(dummy_lon0) > 180.0_r8) THEN + PRINT '(A)', 'Longitude of pole required as follows:' + PRINT '(A)', ' -180E <= lon0 <= 180W' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + ENDIF + + IF ( PRESENT(dx) ) THEN + IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN + PRINT '(A)', 'Require grid spacing (dx) in meters be positive!' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + + IF ( PRESENT(stdlon) ) THEN + dummy_stdlon = stdlon + IF ((ABS(dummy_stdlon) > 180.0_r8).AND.(proj_code /= PROJ_MERC)) THEN + iter = 0 + DO WHILE (ABS(dummy_stdlon) > 180.0_r8 .AND. iter < 10) + IF (dummy_stdlon < -180.0_r8) dummy_stdlon = dummy_stdlon + 360.0_r8 + IF (dummy_stdlon > 180.0_r8) dummy_stdlon = dummy_stdlon - 360.0_r8 + iter = iter + 1 + END DO + IF (abs(dummy_stdlon) > 180.0_r8) THEN + PRINT '(A)', 'Need orientation longitude (stdlon) as: ' + PRINT '(A)', ' -180E <= stdlon <= 180W' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + ENDIF + + IF ( PRESENT(truelat1) ) THEN + IF (ABS(truelat1).GT.90.0_r8) THEN + PRINT '(A)', 'Set true latitude 1 for all projections!' + write(6,*) 'ERROR: MAP_INIT' + ENDIF + ENDIF + + CALL map_init(proj) + proj%code = proj_code + IF ( PRESENT(lat1) ) proj%lat1 = lat1 + IF ( PRESENT(lon1) ) proj%lon1 = dummy_lon1 + IF ( PRESENT(lat0) ) proj%lat0 = lat0 + IF ( PRESENT(lon0) ) proj%lon0 = dummy_lon0 + IF ( PRESENT(latinc) ) proj%latinc = latinc + IF ( PRESENT(loninc) ) proj%loninc = loninc + IF ( PRESENT(knowni) ) proj%knowni = knowni + IF ( PRESENT(knownj) ) proj%knownj = knownj + IF ( PRESENT(dx) ) proj%dx = dx + IF ( PRESENT(stdlon) ) proj%stdlon = dummy_stdlon + IF ( PRESENT(truelat1) ) proj%truelat1 = truelat1 + IF ( PRESENT(truelat2) ) proj%truelat2 = truelat2 + IF ( PRESENT(nlat) ) proj%nlat = nlat + IF ( PRESENT(nlon) ) proj%nlon = nlon + IF ( PRESENT(ixdim) ) proj%ixdim = ixdim + IF ( PRESENT(jydim) ) proj%jydim = jydim + IF ( PRESENT(stagger) ) proj%stagger = stagger + IF ( PRESENT(phi) ) proj%phi = phi + IF ( PRESENT(lambda) ) proj%lambda = lambda + IF ( PRESENT(r_earth) ) proj%re_m = r_earth + + IF ( PRESENT(dx) ) THEN + IF ( (proj_code == PROJ_LC) .OR. (proj_code == PROJ_PS) .OR. & + (proj_code == PROJ_PS_WGS84) .OR. (proj_code == PROJ_ALBERS_NAD83) .OR. & + (proj_code == PROJ_MERC) ) THEN + proj%dx = dx + IF (truelat1 .LT. 0.0_r8) THEN + proj%hemi = -1.0_r8 + ELSE + proj%hemi = 1.0_r8 + ENDIF + proj%rebydx = proj%re_m / dx + ENDIF + ENDIF + + pick_proj: SELECT CASE(proj%code) + + CASE(PROJ_PS) + CALL set_ps(proj) + + CASE(PROJ_PS_WGS84) + CALL set_ps_wgs84(proj) + + CASE(PROJ_ALBERS_NAD83) + CALL set_albers_nad83(proj) + + CASE(PROJ_LC) + IF (ABS(proj%truelat2) .GT. 90.0_r8) THEN + proj%truelat2=proj%truelat1 + ENDIF + CALL set_lc(proj) + + CASE (PROJ_MERC) + CALL set_merc(proj) + + CASE (PROJ_LATLON) + + CASE (PROJ_GAUSS) + CALL set_gauss(proj) + + CASE (PROJ_CYL) + CALL set_cyl(proj) + + CASE (PROJ_CASSINI) + CALL set_cassini(proj) + + CASE (PROJ_ROTLL) + + END SELECT pick_proj + proj%init = .TRUE. + + RETURN + + END SUBROUTINE map_set + + + SUBROUTINE latlon_to_ij(proj, lat, lon, i, j) + ! Converts input lat/lon values to the cartesian (i,j) value + ! for the given projection. + + IMPLICIT NONE + TYPE(proj_info), INTENT(IN) :: proj + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + REAL(r8), INTENT(OUT) :: i + REAL(r8), INTENT(OUT) :: j + + IF (.NOT.proj%init) THEN + PRINT '(A)', 'You have not called map_set for this projection!' + write(6,*) 'ERROR: LATLON_TO_IJ' + ENDIF + + SELECT CASE(proj%code) + + CASE(PROJ_LATLON) + CALL llij_latlon(lat,lon,proj,i,j) + + CASE(PROJ_MERC) + CALL llij_merc(lat,lon,proj,i,j) + + CASE(PROJ_PS) + CALL llij_ps(lat,lon,proj,i,j) + + CASE(PROJ_PS_WGS84) + CALL llij_ps_wgs84(lat,lon,proj,i,j) + + CASE(PROJ_ALBERS_NAD83) + CALL llij_albers_nad83(lat,lon,proj,i,j) + + CASE(PROJ_LC) + CALL llij_lc(lat,lon,proj,i,j) + + CASE(PROJ_GAUSS) + CALL llij_gauss(lat,lon,proj,i,j) + + CASE(PROJ_CYL) + CALL llij_cyl(lat,lon,proj,i,j) + + CASE(PROJ_CASSINI) + CALL llij_cassini(lat,lon,proj,i,j) + + CASE(PROJ_ROTLL) + CALL llij_rotlatlon(lat,lon,proj,i,j) + + CASE DEFAULT + PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code + write(6,*) 'ERROR: LATLON_TO_IJ' + + END SELECT + + RETURN + + END SUBROUTINE latlon_to_ij + + + SUBROUTINE ij_to_latlon(proj, i, j, lat, lon) + ! Computes geographical latitude and longitude for a given (i,j) point + ! in a grid with a projection of proj + + IMPLICIT NONE + TYPE(proj_info),INTENT(IN) :: proj + REAL(r8), INTENT(IN) :: i + REAL(r8), INTENT(IN) :: j + REAL(r8), INTENT(OUT) :: lat + REAL(r8), INTENT(OUT) :: lon + + IF (.NOT.proj%init) THEN + PRINT '(A)', 'You have not called map_set for this projection!' + write(6,*) 'ERROR: IJ_TO_LATLON' + ENDIF + SELECT CASE (proj%code) + + CASE (PROJ_LATLON) + CALL ijll_latlon(i, j, proj, lat, lon) + + CASE (PROJ_MERC) + CALL ijll_merc(i, j, proj, lat, lon) + + CASE (PROJ_PS) + CALL ijll_ps(i, j, proj, lat, lon) + + CASE (PROJ_PS_WGS84) + CALL ijll_ps_wgs84(i, j, proj, lat, lon) + + CASE (PROJ_ALBERS_NAD83) + CALL ijll_albers_nad83(i, j, proj, lat, lon) + + CASE (PROJ_LC) + CALL ijll_lc(i, j, proj, lat, lon) + + CASE (PROJ_CYL) + CALL ijll_cyl(i, j, proj, lat, lon) + + CASE (PROJ_CASSINI) + CALL ijll_cassini(i, j, proj, lat, lon) + + CASE (PROJ_ROTLL) + CALL ijll_rotlatlon(i, j, proj, lat, lon) + + CASE DEFAULT + PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code + write(6,*) 'ERROR: IJ_TO_LATLON' + + END SELECT + RETURN + END SUBROUTINE ij_to_latlon + + + SUBROUTINE set_ps(proj) + ! Initializes a polar-stereographic map projection from the partially + ! filled proj structure. This routine computes the radius to the + ! southwest corner and computes the i/j location of the pole for use + ! in llij_ps and ijll_ps. + IMPLICIT NONE + + ! Declare args + TYPE(proj_info), INTENT(INOUT) :: proj + + ! Local vars + REAL(r8) :: ala1 + REAL(r8) :: alo1 + REAL(r8) :: reflon + REAL(r8) :: scale_top + + ! Executable code + + ! Thanks to Kevin Manning for the 'cone' fix. It must be set to 1.0 + ! to fix wind rotation for polar stereographic projection + reflon = proj%stdlon + 90.0_r8 + proj%cone = 1.0_r8 + + ! Compute numerator term of map scale factor + scale_top = 1.0_r8 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) + + ! Compute radius to lower-left (SW) corner + ala1 = proj%lat1 * rad_per_deg + proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.0_r8+proj%hemi*SIN(ala1)) + + ! Find the pole point + alo1 = (proj%lon1 - reflon) * rad_per_deg + proj%polei = proj%knowni - proj%rsw * COS(alo1) + proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1) + + RETURN + + END SUBROUTINE set_ps + + + SUBROUTINE llij_ps(lat,lon,proj,i,j) + ! Given latitude (-90 to 90), longitude (-180 to 180), and the + ! standard polar-stereographic projection information via the + ! public proj structure, this routine returns the i/j indices which + ! if within the domain range from 1->nx and 1->ny, respectively. + + IMPLICIT NONE + + ! Delcare input arguments + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + TYPE(proj_info),INTENT(IN) :: proj + + ! Declare output arguments + REAL(r8), INTENT(OUT) :: i !(x-index) + REAL(r8), INTENT(OUT) :: j !(y-index) + + ! Declare local variables + + REAL(r8) :: reflon + REAL(r8) :: scale_top + REAL(r8) :: ala + REAL(r8) :: alo + REAL(r8) :: rm + + ! BEGIN CODE + + reflon = proj%stdlon + 90.0_r8 + + ! Compute numerator term of map scale factor + + scale_top = 1.0_r8 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) + + ! Find radius to desired point + ala = lat * rad_per_deg + rm = proj%rebydx * COS(ala) * scale_top/(1.0_r8 + proj%hemi *SIN(ala)) + alo = (lon - reflon) * rad_per_deg + i = proj%polei + rm * COS(alo) + j = proj%polej + proj%hemi * rm * SIN(alo) + + RETURN + + END SUBROUTINE llij_ps + + + SUBROUTINE ijll_ps(i, j, proj, lat, lon) + + ! This is the inverse subroutine of llij_ps. It returns the + ! latitude and longitude of an i/j point given the projection info + ! structure. + + IMPLICIT NONE + + ! Declare input arguments + REAL(r8), INTENT(IN) :: i ! Column + REAL(r8), INTENT(IN) :: j ! Row + TYPE (proj_info), INTENT(IN) :: proj + + ! Declare output arguments + REAL(r8), INTENT(OUT) :: lat ! -90 -> 90 north + REAL(r8), INTENT(OUT) :: lon ! -180 -> 180 East + + ! Local variables + REAL(r8) :: reflon + REAL(r8) :: scale_top + REAL(r8) :: xx,yy + REAL(r8) :: gi2, r2 + REAL(r8) :: arccos + + ! Begin Code + + ! Compute the reference longitude by rotating 90 degrees to the east + ! to find the longitude line parallel to the positive x-axis. + reflon = proj%stdlon + 90.0_r8 + + ! Compute numerator term of map scale factor + scale_top = 1.0_r8 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) + + ! Compute radius to point of interest + xx = i - proj%polei + yy = (j - proj%polej) * proj%hemi + r2 = xx**2 + yy**2 + + ! Now the magic code + IF (r2 .EQ. 0.0_r8) THEN + lat = proj%hemi * 90.0_r8 + lon = reflon + ELSE + gi2 = (proj%rebydx * scale_top)**2. + lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2)) + arccos = ACOS(xx/SQRT(r2)) + IF (yy .GT. 0.0_r8) THEN + lon = reflon + deg_per_rad * arccos + ELSE + lon = reflon - deg_per_rad * arccos + ENDIF + ENDIF + + ! Convert to a -180 -> 180 East convention + IF (lon .GT. 180.0_r8) lon = lon - 360.0_r8 + IF (lon .LT. -180.0_r8) lon = lon + 360.0_r8 + + RETURN + + END SUBROUTINE ijll_ps + + + SUBROUTINE set_ps_wgs84(proj) + ! Initializes a polar-stereographic map projection (WGS84 ellipsoid) + ! from the partially filled proj structure. This routine computes the + ! radius to the southwest corner and computes the i/j location of the + ! pole for use in llij_ps and ijll_ps. + + IMPLICIT NONE + + ! Arguments + TYPE(proj_info), INTENT(INOUT) :: proj + + ! Local variables + real(r8) :: h, mc, tc, t, rho + + h = proj%hemi + + mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0_r8-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0) + tc = sqrt(((1.0_r8-sin(h*proj%truelat1*rad_per_deg))/(1.0_r8+sin(h*proj%truelat1*rad_per_deg)))* & + (((1.0_r8+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0_r8 - & + E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 )) + + ! Find the i/j location of reference lat/lon with respect to the pole of the projection + t = sqrt(((1.0_r8-sin(h*proj%lat1*rad_per_deg))/(1.0_r8+sin(h*proj%lat1*rad_per_deg)))* & + (((1.0_r8+E_WGS84*sin(h*proj%lat1*rad_per_deg))/(1.0_r8 - & + E_WGS84*sin(h*proj%lat1*rad_per_deg)) )**E_WGS84 ) ) + rho = h * (A_WGS84 / proj%dx) * mc * t / tc + proj%polei = rho * sin((h*proj%lon1 - h*proj%stdlon)*rad_per_deg) + proj%polej = -rho * cos((h*proj%lon1 - h*proj%stdlon)*rad_per_deg) + + RETURN + + END SUBROUTINE set_ps_wgs84 + + + SUBROUTINE llij_ps_wgs84(lat,lon,proj,i,j) + ! Given latitude (-90 to 90), longitude (-180 to 180), and the + ! standard polar-stereographic projection information via the + ! public proj structure, this routine returns the i/j indices which + ! if within the domain range from 1->nx and 1->ny, respectively. + + IMPLICIT NONE + + ! Arguments + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + REAL(r8), INTENT(OUT) :: i !(x-index) + REAL(r8), INTENT(OUT) :: j !(y-index) + TYPE(proj_info),INTENT(IN) :: proj + + ! Local variables + real(r8) :: h, mc, tc, t, rho + + h = proj%hemi + + mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0_r8-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0) + tc = sqrt(((1.0_r8-sin(h*proj%truelat1*rad_per_deg))/(1.0_r8+sin(h*proj%truelat1*rad_per_deg)))* & + (((1.0_r8+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0_r8 - & + E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 )) + + t = sqrt(((1.0_r8-sin(h*lat*rad_per_deg))/(1.0_r8+sin(h*lat*rad_per_deg))) * & + (((1.0_r8+E_WGS84*sin(h*lat*rad_per_deg))/(1.0_r8 - & + E_WGS84*sin(h*lat*rad_per_deg)))**E_WGS84)) + + ! Find the x/y location of the requested lat/lon with respect to the pole of the projection + rho = (A_WGS84 / proj%dx) * mc * t / tc + i = h * rho * sin((h*lon - h*proj%stdlon)*rad_per_deg) + j = h *(-rho)* cos((h*lon - h*proj%stdlon)*rad_per_deg) + + ! Get i/j relative to reference i/j + i = proj%knowni + (i - proj%polei) + j = proj%knownj + (j - proj%polej) + + RETURN + + END SUBROUTINE llij_ps_wgs84 + + + SUBROUTINE ijll_ps_wgs84(i, j, proj, lat, lon) + + ! This is the inverse subroutine of llij_ps. It returns the + ! latitude and longitude of an i/j point given the projection info + ! structure. + + implicit none + + ! Arguments + REAL(r8), INTENT(IN) :: i ! Column + REAL(r8), INTENT(IN) :: j ! Row + REAL(r8), INTENT(OUT) :: lat ! -90 -> 90 north + REAL(r8), INTENT(OUT) :: lon ! -180 -> 180 East + TYPE (proj_info), INTENT(IN) :: proj + + ! Local variables + real(r8) :: h, mc, tc, t, rho, x, y + real(r8) :: chi, a, b, c, d + + h = proj%hemi + x = (i - proj%knowni + proj%polei) + y = (j - proj%knownj + proj%polej) + + mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0_r8-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0) + tc = sqrt(((1.0_r8-sin(h*proj%truelat1*rad_per_deg))/(1.0_r8+sin(h*proj%truelat1*rad_per_deg))) * & + (((1.0_r8+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0_r8 - & + E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 )) + + rho = sqrt((x*proj%dx)**2.0 + (y*proj%dx)**2.0) + t = rho * tc / (A_WGS84 * mc) + + lon = h*proj%stdlon + h*atan2(h*x,h*(-y)) + + chi = PI/2.0_r8-2.0_r8*atan(t) + a = (1.0_r8/2.0_r8)*E_WGS84**2. + (5.0_r8/24.0_r8)*E_WGS84**4. + (1.0_r8/40.0_r8)*E_WGS84**6. + & + (73.0_r8/2016.0_r8)*E_WGS84**8. + b = (7.0_r8/24.0_r8)*E_WGS84**4. + (29.0_r8/120.0_r8)*E_WGS84**6. + & + (54113.0_r8/40320.0_r8)*E_WGS84**8. + c = (7.0_r8/30.0_r8)*E_WGS84**6. + (81.0_r8/280.0_r8)*E_WGS84**8. + d = (4279.0_r8/20160.0_r8)*E_WGS84**8. + + lat = chi + sin(2.0_r8*chi)*(a + cos(2.0_r8*chi)*(b + cos(2.0_r8*chi)*(c + d*cos(2.0_r8*chi)))) + lat = h * lat + + lat = lat*deg_per_rad + lon = lon*deg_per_rad + + RETURN + + END SUBROUTINE ijll_ps_wgs84 + + + SUBROUTINE set_albers_nad83(proj) + ! Initializes an Albers equal area map projection (NAD83 ellipsoid) + ! from the partially filled proj structure. This routine computes the + ! radius to the southwest corner and computes the i/j location of the + ! pole for use in llij_albers_nad83 and ijll_albers_nad83. + + IMPLICIT NONE + + ! Arguments + TYPE(proj_info), INTENT(INOUT) :: proj + + ! Local variables + real(r8) :: h, m1, m2, q1, q2, theta, q, sinphi + + h = proj%hemi + + m1 = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0_r8-(E_NAD83*sin(h*proj%truelat1*rad_per_deg))**2.0) + m2 = cos(h*proj%truelat2*rad_per_deg)/sqrt(1.0_r8-(E_NAD83*sin(h*proj%truelat2*rad_per_deg))**2.0) + + sinphi = sin(proj%truelat1*rad_per_deg) + q1 = (1.0_r8-E_NAD83**2.0) * & + ((sinphi/(1.0_r8-(E_NAD83*sinphi)**2.0)) - 1.0_r8/(2.0_r8*E_NAD83) * & + log((1.0_r8-E_NAD83*sinphi)/(1.0_r8+E_NAD83*sinphi))) + + sinphi = sin(proj%truelat2*rad_per_deg) + q2 = (1.0_r8-E_NAD83**2.0) * & + ((sinphi/(1.0_r8-(E_NAD83*sinphi)**2.0)) - 1.0_r8/(2.0_r8*E_NAD83) * & + log((1.0_r8-E_NAD83*sinphi)/(1.0_r8+E_NAD83*sinphi))) + + if (proj%truelat1 == proj%truelat2) then + proj%nc = sin(proj%truelat1*rad_per_deg) + else + proj%nc = (m1**2.0 - m2**2.0) / (q2 - q1) + end if + + proj%bigc = m1**2.0 + proj%nc*q1 + + ! Find the i/j location of reference lat/lon with respect to the pole of the projection + sinphi = sin(proj%lat1*rad_per_deg) + q = (1.0_r8-E_NAD83**2.0) * & + ((sinphi/(1.0_r8-(E_NAD83*sinphi)**2.0)) - 1.0_r8/(2.0_r8*E_NAD83) * & + log((1.0_r8-E_NAD83*sinphi)/(1.0_r8+E_NAD83*sinphi))) + + proj%rho0 = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc + theta = proj%nc*(proj%lon1 - proj%stdlon)*rad_per_deg + + proj%polei = proj%rho0 * sin(h*theta) + proj%polej = proj%rho0 - proj%rho0 * cos(h*theta) + + RETURN + + END SUBROUTINE set_albers_nad83 + + + SUBROUTINE llij_albers_nad83(lat,lon,proj,i,j) + ! Given latitude (-90 to 90), longitude (-180 to 180), and the + ! standard projection information via the + ! public proj structure, this routine returns the i/j indices which + ! if within the domain range from 1->nx and 1->ny, respectively. + + IMPLICIT NONE + + ! Arguments + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + REAL(r8), INTENT(OUT) :: i !(x-index) + REAL(r8), INTENT(OUT) :: j !(y-index) + TYPE(proj_info),INTENT(IN) :: proj + + ! Local variables + real(r8) :: h, q, rho, theta, sinphi + + h = proj%hemi + + sinphi = sin(h*lat*rad_per_deg) + + ! Find the x/y location of the requested lat/lon with respect to the pole of the projection + q = (1.0_r8-E_NAD83**2.0) * & + ((sinphi/(1.0_r8-(E_NAD83*sinphi)**2.0)) - 1.0_r8/(2.0_r8*E_NAD83) * & + log((1.0_r8-E_NAD83*sinphi)/(1.0_r8+E_NAD83*sinphi))) + + rho = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc + theta = proj%nc * (h*lon - h*proj%stdlon)*rad_per_deg + + i = h*rho*sin(theta) + j = h*proj%rho0 - h*rho*cos(theta) + + ! Get i/j relative to reference i/j + i = proj%knowni + (i - proj%polei) + j = proj%knownj + (j - proj%polej) + + RETURN + + END SUBROUTINE llij_albers_nad83 + + + SUBROUTINE ijll_albers_nad83(i, j, proj, lat, lon) + + ! This is the inverse subroutine of llij_albers_nad83. It returns the + ! latitude and longitude of an i/j point given the projection info + ! structure. + + implicit none + + ! Arguments + REAL(r8), INTENT(IN) :: i ! Column + REAL(r8), INTENT(IN) :: j ! Row + REAL(r8), INTENT(OUT) :: lat ! -90 -> 90 north + REAL(r8), INTENT(OUT) :: lon ! -180 -> 180 East + TYPE (proj_info), INTENT(IN) :: proj + + ! Local variables + real(r8) :: h, q, rho, theta, beta, x, y + real(r8) :: a, b, c + + h = proj%hemi + + x = (i - proj%knowni + proj%polei) + y = (j - proj%knownj + proj%polej) + + rho = sqrt(x**2.0 + (proj%rho0 - y)**2.0) + theta = atan2(x, proj%rho0-y) + + q = (proj%bigc - (rho*proj%nc*proj%dx/A_NAD83)**2.0) / proj%nc + + beta = asin(q/(1.0_r8 - & + log((1.0_r8-E_NAD83)/(1.0_r8+E_NAD83))*(1.0_r8-E_NAD83**2.0)/(2.0_r8*E_NAD83))) + a = (1.0_r8/3.0_r8)*E_NAD83**2 + (31.0_r8/180.0_r8)*E_NAD83**4 + (517.0_r8/5040.0_r8)*E_NAD83**6 + b = (23.0_r8/360.0_r8)*E_NAD83**4 + (251.0_r8/3780.0_r8)*E_NAD83**6 + c = (761.0_r8/45360.0_r8)*E_NAD83**6 + + lat = beta + a*sin(2.0_r8*beta) + b*sin(4.0_r8*beta) + c*sin(6.0_r8*beta) + + lat = h*lat*deg_per_rad + lon = proj%stdlon + theta*deg_per_rad/proj%nc + + RETURN + + END SUBROUTINE ijll_albers_nad83 + + + SUBROUTINE set_lc(proj) + ! Initialize the remaining items in the proj structure for a + ! lambert conformal grid. + + IMPLICIT NONE + + TYPE(proj_info), INTENT(INOUT) :: proj + + REAL(r8) :: arg + REAL(r8) :: deltalon1 + REAL(r8) :: tl1r + REAL(r8) :: ctl1r + + ! Compute cone factor + CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone) + + ! Compute longitude differences and ensure we stay out of the + ! forbidden "cut zone" + deltalon1 = proj%lon1 - proj%stdlon + IF (deltalon1 .GT. +180.0_r8) deltalon1 = deltalon1 - 360.0_r8 + IF (deltalon1 .LT. -180.0_r8) deltalon1 = deltalon1 + 360.0_r8 + + ! Convert truelat1 to radian and compute COS for later use + tl1r = proj%truelat1 * rad_per_deg + ctl1r = COS(tl1r) + + ! Compute the radius to our known lower-left (SW) corner + proj%rsw = proj%rebydx * ctl1r/proj%cone * & + (TAN((90.0_r8*proj%hemi-proj%lat1)*rad_per_deg/2.0_r8) / & + TAN((90.0_r8*proj%hemi-proj%truelat1)*rad_per_deg/2.0_r8))**proj%cone + + ! Find pole point + arg = proj%cone*(deltalon1*rad_per_deg) + proj%polei = proj%hemi*proj%knowni - proj%hemi * proj%rsw * SIN(arg) + proj%polej = proj%hemi*proj%knownj + proj%rsw * COS(arg) + + RETURN + + END SUBROUTINE set_lc + + + SUBROUTINE lc_cone(truelat1, truelat2, cone) + + ! Subroutine to compute the cone factor of a Lambert Conformal projection + + IMPLICIT NONE + + ! Input Args + REAL(r8), INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N) + REAL(r8), INTENT(IN) :: truelat2 ! " " " " " + + ! Output Args + REAL(r8), INTENT(OUT) :: cone + + ! Locals + + ! BEGIN CODE + + ! First, see if this is a secant or tangent projection. For tangent + ! projections, truelat1 = truelat2 and the cone is tangent to the + ! Earth's surface at this latitude. For secant projections, the cone + ! intersects the Earth's surface at each of the distinctly different + ! latitudes + IF (ABS(truelat1-truelat2) .GT. 0.10_r8) THEN + cone = LOG10(COS(truelat1*rad_per_deg)) - & + LOG10(COS(truelat2*rad_per_deg)) + cone = cone /(LOG10(TAN((45.0_r8 - ABS(truelat1)/2.0_r8) * rad_per_deg)) - & + LOG10(TAN((45.0_r8 - ABS(truelat2)/2.0_r8) * rad_per_deg))) + ELSE + cone = SIN(ABS(truelat1)*rad_per_deg ) + ENDIF + + RETURN + + END SUBROUTINE lc_cone + + + SUBROUTINE ijll_lc( i, j, proj, lat, lon) + + ! Subroutine to convert from the (i,j) cartesian coordinate to the + ! geographical latitude and longitude for a Lambert Conformal projection. + + ! History: + ! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL + ! + IMPLICIT NONE + + ! Input Args + REAL(r8), INTENT(IN) :: i ! Cartesian X coordinate + REAL(r8), INTENT(IN) :: j ! Cartesian Y coordinate + TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure + + ! Output Args + REAL(r8), INTENT(OUT) :: lat ! Latitude (-90->90 deg N) + REAL(r8), INTENT(OUT) :: lon ! Longitude (-180->180 E) + + ! Locals + REAL(r8) :: inew + REAL(r8) :: jnew + REAL(r8) :: r + REAL(r8) :: chi,chi1,chi2 + REAL(r8) :: r2 + REAL(r8) :: xx + REAL(r8) :: yy + + ! BEGIN CODE + + chi1 = (90.0_r8 - proj%hemi*proj%truelat1)*rad_per_deg + chi2 = (90.0_r8 - proj%hemi*proj%truelat2)*rad_per_deg + + ! See if we are in the southern hemispere and flip the indices + ! if we are. + inew = proj%hemi * i + jnew = proj%hemi * j + + ! Compute radius**2 to i/j location + xx = inew - proj%polei + yy = proj%polej - jnew + r2 = (xx*xx + yy*yy) + r = SQRT(r2)/proj%rebydx + + ! Convert to lat/lon + IF (r2 .EQ. 0.0_r8) THEN + lat = proj%hemi * 90.0_r8 + lon = proj%stdlon + ELSE + + ! Longitude + lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone + lon = MOD(lon+360.0_r8, 360.0_r8) + + ! Latitude. Latitude determined by solving an equation adapted + ! from: + ! Maling, D.H., 1973: Coordinate Systems and Map Projections + ! Equations #20 in Appendix I. + + IF (chi1 .EQ. chi2) THEN + chi = 2.0_r8*ATAN( ( r/TAN(chi1) )**(1.0_r8/proj%cone) * TAN(chi1*0.5_r8) ) + ELSE + chi = 2.0_r8*ATAN( (r*proj%cone/SIN(chi1))**(1.0_r8/proj%cone) * TAN(chi1*0.5_r8)) + ENDIF + lat = (90.0_r8-chi*deg_per_rad)*proj%hemi + + ENDIF + + IF (lon .GT. +180.0_r8) lon = lon - 360.0_r8 + IF (lon .LT. -180.0_r8) lon = lon + 360.0_r8 + + RETURN + + END SUBROUTINE ijll_lc + + + SUBROUTINE llij_lc( lat, lon, proj, i, j) + + ! Subroutine to compute the geographical latitude and longitude values + ! to the cartesian x/y on a Lambert Conformal projection. + + IMPLICIT NONE + + ! Input Args + REAL(r8), INTENT(IN) :: lat ! Latitude (-90->90 deg N) + REAL(r8), INTENT(IN) :: lon ! Longitude (-180->180 E) + TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure + + ! Output Args + REAL(r8), INTENT(OUT) :: i ! Cartesian X coordinate + REAL(r8), INTENT(OUT) :: j ! Cartesian Y coordinate + + ! Locals + REAL(r8) :: arg + REAL(r8) :: deltalon + REAL(r8) :: tl1r + REAL(r8) :: rm + REAL(r8) :: ctl1r + + + ! BEGIN CODE + + ! Compute deltalon between known longitude and standard lon and ensure + ! it is not in the cut zone + deltalon = lon - proj%stdlon + IF (deltalon .GT. +180.0_r8) deltalon = deltalon - 360.0_r8 + IF (deltalon .LT. -180.0_r8) deltalon = deltalon + 360.0_r8 + + ! Convert truelat1 to radian and compute COS for later use + tl1r = proj%truelat1 * rad_per_deg + ctl1r = COS(tl1r) + + ! Radius to desired point + rm = proj%rebydx * ctl1r/proj%cone * & + (TAN((90.0_r8*proj%hemi-lat)*rad_per_deg/2.0_r8) / & + TAN((90.0_r8*proj%hemi-proj%truelat1)*rad_per_deg/2.0_r8))**proj%cone + + arg = proj%cone*(deltalon*rad_per_deg) + i = proj%polei + proj%hemi * rm * SIN(arg) + j = proj%polej - rm * COS(arg) + + ! Finally, if we are in the southern hemisphere, flip the i/j + ! values to a coordinate system where (1,1) is the SW corner + ! (what we assume) which is different than the original NCEP + ! algorithms which used the NE corner as the origin in the + ! southern hemisphere (left-hand vs. right-hand coordinate?) + i = proj%hemi * i + j = proj%hemi * j + + RETURN + END SUBROUTINE llij_lc + + + SUBROUTINE set_merc(proj) + + ! Sets up the remaining basic elements for the mercator projection + + IMPLICIT NONE + TYPE(proj_info), INTENT(INOUT) :: proj + REAL(r8) :: clain + + + ! Preliminary variables + + clain = COS(rad_per_deg*proj%truelat1) + proj%dlon = proj%dx / (proj%re_m * clain) + + ! Compute distance from equator to origin, and store in the + ! proj%rsw tag. + + proj%rsw = 0.0_r8 + IF (proj%lat1 .NE. 0.0_r8) THEN + proj%rsw = (LOG(TAN(0.50_r8*((proj%lat1+90.0_r8)*rad_per_deg))))/proj%dlon + ENDIF + + RETURN + + END SUBROUTINE set_merc + + + SUBROUTINE llij_merc(lat, lon, proj, i, j) + + ! Compute i/j coordinate from lat lon for mercator projection + + IMPLICIT NONE + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + TYPE(proj_info),INTENT(IN) :: proj + REAL(r8),INTENT(OUT) :: i + REAL(r8),INTENT(OUT) :: j + REAL(r8) :: deltalon, i2 + + deltalon = lon - proj%lon1 + IF (deltalon .LT. -180.0_r8) deltalon = deltalon + 360.0_r8 + IF (deltalon .GT. 180.0_r8) deltalon = deltalon - 360.0_r8 + i = proj%knowni + (deltalon/(proj%dlon*deg_per_rad)) + i2 = proj%knowni + ((deltalon+360.0_r8)/(proj%dlon*deg_per_rad)) + if ( i < 0.0_r8 ) i = i2 + j = proj%knownj + (LOG(TAN(0.50_r8*((lat + 90.0_r8) * rad_per_deg)))) / & + proj%dlon - proj%rsw + + RETURN + + END SUBROUTINE llij_merc + + + SUBROUTINE ijll_merc(i, j, proj, lat, lon) + + ! Compute the lat/lon from i/j for mercator projection + + IMPLICIT NONE + REAL(r8),INTENT(IN) :: i + REAL(r8),INTENT(IN) :: j + TYPE(proj_info),INTENT(IN) :: proj + REAL(r8), INTENT(OUT) :: lat + REAL(r8), INTENT(OUT) :: lon + + + lat = 2.00_r8*ATAN(EXP(proj%dlon*(proj%rsw + j-proj%knownj)))*deg_per_rad - 90.0_r8 + lon = (i-proj%knowni)*proj%dlon*deg_per_rad + proj%lon1 + IF (lon.GT.180.0_r8) lon = lon - 360.0_r8 + IF (lon.LT.-180.0_r8) lon = lon + 360.0_r8 + RETURN + + END SUBROUTINE ijll_merc + + + SUBROUTINE llij_latlon(lat, lon, proj, i, j) + + ! Compute the i/j location of a lat/lon on a LATLON grid. + IMPLICIT NONE + REAL(r8), INTENT(IN) :: lat + REAL(r8), INTENT(IN) :: lon + TYPE(proj_info), INTENT(IN) :: proj + REAL(r8), INTENT(OUT) :: i + REAL(r8), INTENT(OUT) :: j + + REAL(r8) :: deltalat + REAL(r8) :: deltalon + + ! Compute deltalat and deltalon as the difference between the input + ! lat/lon and the origin lat/lon + deltalat = lat - proj%lat1 + deltalon = lon - proj%lon1 + + ! Compute i/j + i = deltalon/proj%loninc + j = deltalat/proj%latinc + + i = i + proj%knowni + j = j + proj%knownj + + RETURN + + END SUBROUTINE llij_latlon + + + SUBROUTINE ijll_latlon(i, j, proj, lat, lon) + + ! Compute the lat/lon location of an i/j on a LATLON grid. + IMPLICIT NONE + REAL(r8), INTENT(IN) :: i + REAL(r8), INTENT(IN) :: j + TYPE(proj_info), INTENT(IN) :: proj + REAL(r8), INTENT(OUT) :: lat + REAL(r8), INTENT(OUT) :: lon + + REAL(r8) :: i_work, j_work + REAL(r8) :: deltalat + REAL(r8) :: deltalon + + i_work = i - proj%knowni + j_work = j - proj%knownj + + ! Compute deltalat and deltalon + deltalat = j_work*proj%latinc + deltalon = i_work*proj%loninc + + lat = proj%lat1 + deltalat + lon = proj%lon1 + deltalon + + RETURN + + END SUBROUTINE ijll_latlon + + + SUBROUTINE set_cyl(proj) + + implicit none + + ! Arguments + type(proj_info), intent(inout) :: proj + + proj%hemi = 1.00_r8 + + END SUBROUTINE set_cyl + + + SUBROUTINE llij_cyl(lat, lon, proj, i, j) + + implicit none + + ! Arguments + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: i, j + type(proj_info), intent(in) :: proj + + ! Local variables + real(r8) :: deltalat + real(r8) :: deltalon + + ! Compute deltalat and deltalon as the difference between the input + ! lat/lon and the origin lat/lon + deltalat = lat - proj%lat1 +! deltalon = lon - proj%stdlon + deltalon = lon - proj%lon1 + + if (deltalon < 0.0_r8) deltalon = deltalon + 360.0_r8 + if (deltalon > 360.0_r8) deltalon = deltalon - 360.0_r8 + + ! Compute i/j + i = deltalon/proj%loninc + j = deltalat/proj%latinc + + if (i <= 0.0_r8) i = i + 360.0_r8/proj%loninc + if (i > 360.0_r8/proj%loninc) i = i - 360.0_r8/proj%loninc + + i = i + proj%knowni +!nc -- typo in original I am assuming +!nc -- orig --> j = j + proj%knowni + j = j + proj%knownj + + END SUBROUTINE llij_cyl + + + SUBROUTINE ijll_cyl(i, j, proj, lat, lon) + + implicit none + + ! Arguments + real(r8), intent(in) :: i, j + real(r8), intent(out) :: lat, lon + type(proj_info), intent(in) :: proj + + ! Local variables + real(r8) :: deltalat + real(r8) :: deltalon + real(r8) :: i_work, j_work + + i_work = i - proj%knowni + j_work = j - proj%knownj + + if (i_work < 0.0_r8) i_work = i_work + 360.0_r8/proj%loninc + if (i_work >= 360.0_r8/proj%loninc) i_work = i_work - 360.0_r8/proj%loninc + + ! Compute deltalat and deltalon + deltalat = j_work*proj%latinc + deltalon = i_work*proj%loninc + + lat = deltalat + proj%lat1 +! lon = deltalon + proj%stdlon + lon = deltalon + proj%lon1 + + if (lon < -180.0_r8) lon = lon + 360.0_r8 + if (lon > 180.0_r8) lon = lon - 360.0_r8 + + END SUBROUTINE ijll_cyl + + + SUBROUTINE set_cassini(proj) + + implicit none + + ! Arguments + type(proj_info), intent(inout) :: proj + + ! Local variables + real(r8) :: comp_lat, comp_lon + logical :: global_domain + + proj%hemi = 1.00_r8 + + ! Try to determine whether this domain has global coverage + if (abs(proj%lat1 - proj%latinc/2.0_r8 + 90.0_r8) < 0.001_r8 .and. & + abs(mod(proj%lon1 - proj%loninc/2.0_r8 - proj%stdlon,360.0_r8)) < 0.001_r8) then + global_domain = .true. + else + global_domain = .false. + end if + + if (abs(proj%lat0) /= 90.0_r8 .and. .not.global_domain) then + call rotate_coords(proj%lat1,proj%lon1,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1) + proj%lat1 = comp_lat + proj%lon1 = comp_lon + end if + + END SUBROUTINE set_cassini + + + SUBROUTINE llij_cassini(lat, lon, proj, i, j) + + implicit none + + ! Arguments + real(r8), intent(in) :: lat, lon + real(r8), intent(out) :: i, j + type(proj_info), intent(in) :: proj + + ! Local variables + real(r8) :: comp_lat, comp_lon + + ! Convert geographic to computational lat/lon + if (abs(proj%lat0) /= 90.0_r8) then + call rotate_coords(lat,lon,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1) + else + comp_lat = lat + comp_lon = lon + end if + + ! Convert computational lat/lon to i/j + call llij_cyl(comp_lat, comp_lon, proj, i, j) + + END SUBROUTINE llij_cassini + + + SUBROUTINE ijll_cassini(i, j, proj, lat, lon) + + implicit none + + ! Arguments + real(r8), intent(in) :: i, j + real(r8), intent(out) :: lat, lon + type(proj_info), intent(in) :: proj + + ! Local variables + real(r8) :: comp_lat, comp_lon + + ! Convert i/j to computational lat/lon + call ijll_cyl(i, j, proj, comp_lat, comp_lon) + + ! Convert computational to geographic lat/lon + if (abs(proj%lat0) /= 90.0_r8) then + call rotate_coords(comp_lat,comp_lon,lat,lon,proj%lat0,proj%lon0,proj%stdlon,1) + else + lat = comp_lat + lon = comp_lon + end if + + END SUBROUTINE ijll_cassini + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Purpose: Converts between computational and geographic lat/lon for Cassini + ! + ! Notes: This routine was provided by Bill Skamarock, 2007-03-27 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE rotate_coords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction) + + IMPLICIT NONE + + REAL(r8), INTENT(IN ) :: ilat, ilon + REAL(r8), INTENT( OUT) :: olat, olon + REAL(r8), INTENT(IN ) :: lat_np, lon_np, lon_0 + INTEGER, INTENT(IN ), OPTIONAL :: direction + ! >=0, default : computational -> geographical + ! < 0 : geographical -> computational + + REAL(r8) :: rlat, rlon + REAL(r8) :: phi_np, lam_np, lam_0, dlam + REAL(r8) :: sinphi, cosphi, coslam, sinlam + + ! Convert all angles to radians + phi_np = lat_np * rad_per_deg + lam_np = lon_np * rad_per_deg + lam_0 = lon_0 * rad_per_deg + rlat = ilat * rad_per_deg + rlon = ilon * rad_per_deg + + IF (PRESENT(direction) .AND. (direction < 0)) THEN + ! The equations are exactly the same except for one small difference + ! with respect to longitude ... + dlam = PI - lam_0 + ELSE + dlam = lam_np + END IF + sinphi = COS(phi_np)*COS(rlat)*COS(rlon-dlam) + SIN(phi_np)*SIN(rlat) + cosphi = SQRT(1.0_r8-sinphi*sinphi) + coslam = SIN(phi_np)*COS(rlat)*COS(rlon-dlam) - COS(phi_np)*SIN(rlat) + sinlam = COS(rlat)*SIN(rlon-dlam) + IF ( cosphi /= 0.0_r8 ) THEN + coslam = coslam/cosphi + sinlam = sinlam/cosphi + END IF + olat = deg_per_rad*ASIN(sinphi) + olon = deg_per_rad*(ATAN2(sinlam,coslam)-dlam-lam_0+lam_np) + ! Both of my F90 text books prefer the DO-EXIT form, and claim it is faster + ! when optimization is turned on (as we will always do...) + DO + IF (olon >= -180.0_r8) EXIT + olon = olon + 360.0_r8 + END DO + DO + IF (olon <= 180.0_r8) EXIT + olon = olon - 360.0_r8 + END DO + + END SUBROUTINE rotate_coords + + + SUBROUTINE llij_rotlatlon(lat, lon, proj, i, j) + + IMPLICIT NONE + + ! Arguments + REAL(r8), INTENT(IN) :: lat, lon + REAL(r8), INTENT(OUT) :: i, j + TYPE (proj_info), INTENT(IN) :: proj + + ! Local variables -- pi is a parameter in constants_module! -- local value will + ! be "ppi" + REAL(KIND=HIGH) :: dphd,dlmd !Grid increments, degrees + INTEGER :: ii,jj,jmt,ncol,nrow + REAL(KIND=HIGH) :: glatd !Geographic latitude, positive north + REAL(KIND=HIGH) :: glond !Geographic longitude, positive west + REAL(KIND=HIGH) :: col,d1,d2,d2r,dlm,dlm1,dlm2,dph,glat,glon, & + ppi,r2d,row,tlat,tlat1,tlat2, & + tlon,tlon1,tlon2,tph0,tlm0,x,y,z + + glatd = REAL(lat,HIGH) ! HIGH + glond = -REAL(lon,HIGH) ! HIGH + + dphd = REAL(proj%phi,HIGH)/REAL((proj%jydim-1)/2,HIGH) ! HIGH + dlmd = REAL(proj%lambda,HIGH)/REAL(proj%ixdim-1,HIGH) ! HIGH + + ppi = ACOS(-1.0_HIGH) ! HIGH (same name as module variable) + d2r = ppi/180.0_HIGH ! HIGH + r2d = 1.0_HIGH/d2r ! HIGH + + jmt = proj%jydim/2+1 ! INTEGER + + glat = glatd*d2r ! HIGH + glon = glond*d2r ! HIGH + dph = dphd*d2r ! HIGH + dlm = dlmd*d2r ! HIGH + tph0 = REAL(proj%lat1,HIGH)*d2r ! HIGH + tlm0 = -REAL(proj%lon1,HIGH)*d2r ! HIGH + + x = COS(tph0)*COS(glat)*COS(glon-tlm0)+SIN(tph0)*SIN(glat) ! HIGH + y = -COS(glat)*SIN(glon-tlm0) ! HIGH + z = COS(tph0)*SIN(glat)-SIN(tph0)*COS(glat)*COS(glon-tlm0) ! HIGH + tlat = r2d*ATAN(z/SQRT(x*x+y*y)) ! HIGH + tlon = r2d*ATAN(y/x) ! HIGH + + row = tlat/dphd+REAL(jmt,HIGH) ! HIGH + col = tlon/dlmd+REAL(proj%ixdim,HIGH) ! HIGH + + nrow = NINT(row) ! INTEGER + ncol = NINT(col) ! INTEGER + + tlat = tlat*d2r ! HIGH + tlon = tlon*d2r ! HIGH + + IF (proj%stagger == HH) THEN + + IF ((abs(MOD(nrow,2)) == 1 .AND. abs(MOD(ncol,2)) == 1) .OR. & + (MOD(nrow,2) == 0 .AND. MOD(ncol,2) == 0)) THEN + + tlat1 = REAL((nrow-jmt),HIGH)*dph + tlat2 = tlat1+dph + tlon1 = REAL((ncol-proj%ixdim),HIGH)*dlm + tlon2 = tlon1+dlm + + dlm1 = tlon-tlon1 + dlm2 = tlon-tlon2 + d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1)) + d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2)) + + IF (d1 > d2) THEN + nrow = nrow+1 + ncol = ncol+1 + END IF + + ELSE + tlat1 = REAL((nrow+1-jmt),HIGH)*dph + tlat2 = tlat1-dph + tlon1 = REAL((ncol-proj%ixdim),HIGH)*dlm + tlon2 = tlon1+dlm + dlm1 = tlon-tlon1 + dlm2 = tlon-tlon2 + d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1)) + d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2)) + + IF (d1 < d2) THEN + nrow = nrow+1 + ELSE + ncol = ncol+1 + END IF + END IF + + ELSE IF (proj%stagger == VV) THEN + + IF ((MOD(nrow,2) == 0 .AND. abs(MOD(ncol,2)) == 1) .OR. & + (abs(MOD(nrow,2)) == 1 .AND. MOD(ncol,2) == 0)) THEN + tlat1 = REAL((nrow-jmt),HIGH)*dph + tlat2 = tlat1+dph + tlon1 = REAL((ncol-proj%ixdim),HIGH)*dlm + tlon2 = tlon1+dlm + dlm1 = tlon-tlon1 + dlm2 = tlon-tlon2 + d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1)) + d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2)) + + IF (d1 > d2) THEN + nrow = nrow+1 + ncol = ncol+1 + END IF + + ELSE + tlat1 = REAL((nrow+1-jmt),HIGH)*dph + tlat2 = tlat1-dph + tlon1 = REAL((ncol-proj%ixdim),HIGH)*dlm + tlon2 = tlon1+dlm + dlm1 = tlon-tlon1 + dlm2 = tlon-tlon2 + d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1)) + d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2)) + + IF (d1 < d2) THEN + nrow = nrow+1 + ELSE + ncol = ncol+1 + END IF + END IF + END IF + + +!!! Added next line as a Kludge - not yet understood why needed + if (ncol .le. 0) ncol=ncol-1 + + jj = nrow + ii = ncol/2 + + IF (proj%stagger == HH) THEN + IF (abs(MOD(jj,2)) == 1) ii = ii+1 + ELSE IF (proj%stagger == VV) THEN + IF (MOD(jj,2) == 0) ii=ii+1 + END IF + + i = REAL(ii,r8) + j = REAL(jj,r8) + + END SUBROUTINE llij_rotlatlon + + + SUBROUTINE ijll_rotlatlon(i, j, proj, lat,lon) + + IMPLICIT NONE + + ! Arguments + REAL(R8), INTENT(IN) :: i, j + REAL(R8), INTENT(OUT) :: lat, lon + TYPE (proj_info), INTENT(IN) :: proj + + ! Local variables + INTEGER :: ih,jh + INTEGER :: midcol,midrow,ncol + REAL(KIND=HIGH) :: dphd,dlmd !Grid increments, degrees + REAL(KIND=HIGH) :: arg1,arg2,d2r,fctr,glatr,glatd,glond,ppi, & + r2d,tlatd,tlond,tlatr,tlonr,tlm0,tph0 + + ih = NINT(i) + jh = NINT(j) + + dphd = REAL(proj%phi,HIGH)/REAL((proj%jydim-1)/2,HIGH) + dlmd = REAL(proj%lambda,HIGH)/REAL(proj%ixdim-1,HIGH) + + ppi = ACOS(-1.0_HIGH) + d2r = ppi/180.0_HIGH + r2d = 1.0_HIGH/d2r + tph0 = REAL(proj%lat1,HIGH)*d2r + tlm0 = -REAL(proj%lon1,HIGH)*d2r + + midrow = (proj%jydim+1)/2 + midcol = proj%ixdim + +! IF (proj%stagger == HH) THEN + +!!! ncol = 2*ih-1+MOD(jh+1,2) + ncol = 2*ih-1+abs(MOD(jh+1,2)) + tlatd = REAL((jh-midrow),HIGH)*dphd + tlond = REAL((ncol-midcol),HIGH)*dlmd + IF (proj%stagger == VV) THEN + if (mod(jh,2) .eq. 0) then + tlond = tlond - DLMD + else + tlond = tlond + DLMD + end if + END IF + + tlatr = tlatd*d2r + tlonr = tlond*d2r + arg1 = SIN(tlatr)*COS(tph0)+COS(tlatr)*SIN(tph0)*COS(tlonr) + glatr = ASIN(arg1) + + glatd = glatr*r2d + + arg2 = COS(tlatr)*COS(tlonr)/(COS(glatr)*COS(tph0))-TAN(glatr)*TAN(tph0) + IF (ABS(arg2) > 1.0_HIGH) arg2 = ABS(arg2)/arg2 + fctr = 1.0_HIGH + IF (tlond > 0.0_HIGH) fctr = -1.0_HIGH + + glond = tlm0*r2d+fctr*ACOS(arg2)*r2d + + lat = REAL(glatd,r8) + lon = -REAL(glond,r8) + + IF (lon .GT. +180.0_r8) lon = lon - 360.0_r8 + IF (lon .LT. -180.0_r8) lon = lon + 360.0_r8 + + END SUBROUTINE ijll_rotlatlon + + + SUBROUTINE set_gauss(proj) + + IMPLICIT NONE + + ! Argument + type (proj_info), intent(inout) :: proj + + ! Initialize the array that will hold the Gaussian latitudes. + + IF ( ASSOCIATED( proj%gauss_lat ) ) THEN + DEALLOCATE ( proj%gauss_lat ) + END IF + + ! Get the needed space for our array. + + ALLOCATE ( proj%gauss_lat(proj%nlat*2) ) + + ! Compute the Gaussian latitudes. + + CALL gausll( proj%nlat*2 , proj%gauss_lat ) + + ! Now, these could be upside down from what we want, so let's check. + ! We take advantage of the equatorial symmetry to remove any sort of + ! array re-ordering. + + IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01_r8 ) THEN + proj%gauss_lat = -1.0_r8 * proj%gauss_lat + END IF + + ! Just a sanity check. + + IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01_r8 ) THEN + PRINT '(A)','Oops, something is not right with the Gaussian latitude computation.' + PRINT '(A,F8.3,A)','The input data gave the starting latitude as ',proj%lat1,'.' + PRINT '(A,F8.3,A)','This routine computed the starting latitude as +-',ABS(proj%gauss_lat(1)),'.' + PRINT '(A,F8.3,A)','The difference is larger than 0.01 degrees, which is not expected.' + write(6,*) 'ERROR: Gaussian_latitude_computation' + END IF + + END SUBROUTINE set_gauss + + + SUBROUTINE gausll ( nlat , lat_sp ) + + IMPLICIT NONE + + INTEGER :: nlat , i + REAL (KIND=HIGH) , PARAMETER :: ppi = 3.141592653589793_HIGH + REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2 , lat + REAL(R8) , DIMENSION(nlat) :: lat_sp + + CALL lggaus(nlat, cosc, gwt, sinc, colat, wos2) + + DO i = 1, nlat + lat(i) = ACOS(sinc(i)) * 180.0_HIGH / ppi + IF (i.gt.nlat/2) lat(i) = -lat(i) + END DO + + lat_sp = REAL(lat,r8) + + END SUBROUTINE gausll + + + SUBROUTINE lggaus( nlat, cosc, gwt, sinc, colat, wos2 ) + + IMPLICIT NONE + + ! LGGAUS finds the Gaussian latitudes by finding the roots of the + ! ordinary Legendre polynomial of degree NLAT using Newton's + ! iteration method. + + ! On entry: + integer NLAT ! the number of latitudes (degree of the polynomial) + + ! On exit: for each Gaussian latitude + ! COSC - cos(colatitude) or sin(latitude) + ! GWT - the Gaussian weights + ! SINC - sin(colatitude) or cos(latitude) + ! COLAT - the colatitudes in radians + ! WOS2 - Gaussian weight over sin**2(colatitude) + + REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2 + REAL (KIND=HIGH) , PARAMETER :: ppi = 3.141592653589793_HIGH + + ! Convergence criterion for iteration of cos latitude + + REAL(KIND=HIGH) , PARAMETER :: xlim = 1.0E-14_HIGH + + INTEGER :: nzero, i, j + REAL (KIND=HIGH) :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d + + ! The number of zeros between pole and equator + + nzero = nlat/2 + + ! Set first guess for cos(colat) + + DO i=1,nzero + cosc(i) = SIN( (REAL(i,HIGH)-0.5_HIGH)*ppi/REAL(nlat,HIGH) + ppi*0.5_HIGH ) + END DO + + ! Constants for determining the derivative of the polynomial + fi = nlat + fi1 = fi + 1.0_HIGH + a = fi*fi1 / SQRT(4.0_HIGH * fi1*fi1 - 1.0_HIGH) + b = fi1*fi / SQRT(4.0_HIGH * fi*fi - 1.0_HIGH) + + ! Loop over latitudes, iterating the search for each root + + DO i=1,nzero + j=0 + + ! Determine the value of the ordinary Legendre polynomial for + ! the current guess root + + DO + CALL lgord( g, cosc(i), nlat ) + + ! Determine the derivative of the polynomial at this point + + CALL lgord( gm, cosc(i), nlat-1 ) + CALL lgord( gp, cosc(i), nlat+1 ) + gt = (cosc(i)*cosc(i)-1.0_HIGH) / (a*gp-b*gm) + + ! Update the estimate of the root + + delta = g*gt + cosc(i) = cosc(i) - delta + + ! If convergence criterion has not been met, keep trying + + j = j+1 + IF( ABS(delta).GT.xlim ) CYCLE + + ! Determine the Gaussian weights + + c = 2.0_HIGH *( 1.0_HIGH - cosc(i)*cosc(i) ) + CALL lgord( d, cosc(i), nlat-1 ) + d = d*d*fi*fi + gwt(i) = c *( fi-0.5_HIGH ) / d + EXIT + + END DO + + END DO + + ! Determine the colatitudes and sin(colat) and weights over sin**2 + + DO i=1,nzero + colat(i)= ACOS(cosc(i)) + sinc(i) = SIN(colat(i)) + wos2(i) = gwt(i) /( sinc(i)*sinc(i) ) + END DO + + ! If NLAT is odd, set values at the equator + + IF( MOD(nlat,2) .NE. 0 ) THEN + i = nzero+1 + cosc(i) = 0.0_HIGH + c = 2.0_HIGH + CALL lgord( d, cosc(i), nlat-1 ) + d = d*d*fi*fi + gwt(i) = c *( fi-0.5_HIGH ) / d + colat(i)= ppi*0.5_HIGH + sinc(i) = 1.0_HIGH + wos2(i) = gwt(i) + END IF + + ! Determine the southern hemisphere values by symmetry + + DO i=nlat-nzero+1,nlat + cosc(i) =-cosc(nlat+1-i) + gwt(i) = gwt(nlat+1-i) + colat(i)= ppi-colat(nlat+1-i) + sinc(i) = sinc(nlat+1-i) + wos2(i) = wos2(nlat+1-i) + END DO + + END SUBROUTINE lggaus + + + SUBROUTINE lgord( f, cosc, n ) + + IMPLICIT NONE + + ! LGORD calculates the value of an ordinary Legendre polynomial at a + ! specific latitude. + + ! On entry: + ! cosc - COS(colatitude) + ! n - the degree of the polynomial + + ! On exit: + ! f - the value of the Legendre polynomial of degree N at + ! latitude ASIN(cosc) + + REAL (KIND=HIGH) :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang + INTEGER :: n, k + + ! Determine the colatitude + + colat = ACOS(cosc) + + c1 = SQRT(2.0_HIGH) + DO k=1,n + c1 = c1 * SQRT( 1.0_HIGH - 1.0_HIGH/REAL((4*k*k),HIGH) ) + END DO + + fn = n + ang= fn * colat + s1 = 0.0_HIGH + c4 = 1.0_HIGH + a =-1.0_HIGH + b = 0.0_HIGH + DO k=0,n,2 + IF (k.eq.n) c4 = 0.5_HIGH * c4 + s1 = s1 + c4 * COS(ang) + a = a + 2.0_HIGH + b = b + 1.0_HIGH + fk = k + ang= colat * (fn-fk-2.0_HIGH) + c4 = ( a * (fn-b+1.0_HIGH) / ( b * (fn+fn-a) ) ) * c4 + END DO + + f = s1 * c1 + + END SUBROUTINE lgord + + + SUBROUTINE llij_gauss (lat, lon, proj, i, j) + + IMPLICIT NONE + + REAL(R8) , INTENT(IN) :: lat, lon + REAL(R8) , INTENT(OUT) :: i, j + TYPE (proj_info), INTENT(IN) :: proj + + INTEGER :: n , n_low + LOGICAL :: found = .FALSE. + REAL(R8) :: diff_1 , diff_nlat + + ! The easy one first, get the x location. The calling routine has already made + ! sure that the necessary assumptions concerning the sign of the deltalon and the + ! relative east/west'ness of the longitude and the starting longitude are consistent + ! to allow this easy computation. + + i = ( lon - proj%lon1 ) / proj%loninc + 1.0_r8 + + ! Since this is a global data set, we need to be concerned about wrapping the + ! fields around the globe. + +! IF ( ( proj%loninc .GT. 0 ) .AND. & +! ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. & +! ( lon + proj%loninc .GE. proj%lon1 + 360 ) ) THEN +!! BUG: We may need to set proj%wrap, but proj is intent(in) +!! WHAT IS THIS USED FOR? +!! proj%wrap = .TRUE. +! i = proj%ixdim +! ELSE IF ( ( proj%loninc .LT. 0 ) .AND. & +! ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. & +! ( lon + proj%loninc .LE. proj%lon1 - 360 ) ) THEN +! ! BUG: We may need to set proj%wrap, but proj is intent(in) +! ! WHAT IS THIS USED FOR? +! ! proj%wrap = .TRUE. +! i = proj%ixdim +! END IF + + ! Yet another quicky test, can we find bounding values? If not, then we may be + ! dealing with putting data to a polar projection, so just give them them maximal + ! value for the location. This is an OK assumption for the interpolation across the + ! top of the pole, given how close the longitude lines are. + + IF ( ABS(lat) .GT. ABS(proj%gauss_lat(1)) ) THEN + + diff_1 = lat - proj%gauss_lat(1) + diff_nlat = lat - proj%gauss_lat(proj%nlat*2) + + IF ( ABS(diff_1) .LT. ABS(diff_nlat) ) THEN + j = 1 + ELSE + j = proj%nlat*2 + END IF + + ! If the latitude is between the two bounding values, we have to search and interpolate. + + ELSE + + DO n = 1 , proj%nlat*2 -1 + IF ( ( proj%gauss_lat(n) - lat ) * ( proj%gauss_lat(n+1) - lat ) .LE. 0.0_r8 ) THEN + found = .TRUE. + n_low = n + EXIT + END IF + END DO + + ! Everything still OK? + + IF ( .NOT. found ) THEN + PRINT '(A)','Troubles in river city. No bounding values of latitude found in the Gaussian routines.' + write(6,*) 'ERROR: Gee_no_bounding_lats_Gaussian' + END IF + + j = ( ( proj%gauss_lat(n_low) - lat ) * REAL( (n_low + 1),r8 ) + & + ( lat - proj%gauss_lat(n_low+1) ) * REAL( n_low,r8 ) ) / & + ( proj%gauss_lat(n_low) - proj%gauss_lat(n_low+1) ) + + END IF + + END SUBROUTINE llij_gauss + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!nc -- This subroutine was in the original module_map_utils.f90 (before PROJ_CASSINI), and +! the model_mod.f90 would still like it to be, hence we are including it. + SUBROUTINE gridwind_to_truewind(lon,proj,ugrid,vgrid,utrue,vtrue) + + ! Subroutine to convert a wind from grid north to true north. + + IMPLICIT NONE + + ! Arguments + REAL(r8), INTENT(IN) :: lon ! Longitude of point in degrees + TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure + REAL(r8), INTENT(IN) :: ugrid ! U-component, grid-relative + REAL(r8), INTENT(IN) :: vgrid ! V-component, grid-relative + REAL(r8), INTENT(OUT) :: utrue ! U-component, earth-relative + REAL(r8), INTENT(OUT) :: vtrue ! V-component, earth-relative + + ! Locals + REAL(r8) :: alpha + REAL(r8) :: diff + + IF ((proj%code .EQ. PROJ_PS).OR.(proj%code .EQ. PROJ_LC))THEN + diff = lon - proj%stdlon + IF (diff .GT. 180.0_r8) diff = diff - 360.0_r8 + IF (diff .LT.-180.0_r8) diff = diff + 360.0_r8 + alpha = diff * proj%cone * rad_per_deg * SIGN(1.0_r8,proj%truelat1) + utrue = vgrid * SIN(alpha) + ugrid * COS(alpha) + vtrue = vgrid * COS(alpha) - ugrid * SIN(alpha) +!nc -- we added in a case structure for CASSINI and CYL + ELSEIF ((proj%code .EQ. PROJ_MERC).OR.(proj%code .EQ. PROJ_LATLON) & + .OR.(proj%code .EQ. PROJ_CASSINI).OR.(proj%code .EQ. PROJ_CYL))THEN + utrue = ugrid + vtrue = vgrid + ELSE + PRINT '(A)', 'Unrecognized projection.' + STOP 'GRIDWIND_TO_TRUEWIND' + ENDIF + + RETURN + END SUBROUTINE gridwind_to_truewind +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!nc -- This subroutine was in the original module_map_utils.f90 (before PROJ_CASSINI), and +! the model_mod.f90 would still like it to be, hence we are including it. + SUBROUTINE truewind_to_gridwind(lon,proj,utrue,vtrue,ugrid,vgrid) + + ! Subroutine to compute grid-relative u/v wind components from the earth- + ! relative values for a given projection. + + IMPLICIT NONE + + ! Arguments + REAL(r8), INTENT(IN) :: lon ! Longitude of point in degrees + TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure + REAL(r8), INTENT(IN) :: utrue ! U-component, earth-relative + REAL(r8), INTENT(IN) :: vtrue ! V-component, earth-relative + REAL(r8), INTENT(OUT) :: ugrid ! U-component, grid-relative + REAL(r8), INTENT(OUT) :: vgrid ! V-component, grid-relative + + ! Locals + REAL(r8) :: alpha + REAL(r8) :: diff + + IF ((proj%code .EQ. PROJ_PS).OR.(proj%code .EQ. PROJ_LC))THEN + + diff = proj%stdlon - lon + IF (diff .GT. 180.0_r8) diff = diff - 360.0_r8 + IF (diff .LT.-180.0_r8) diff = diff + 360.0_r8 + alpha = diff * proj%cone * rad_per_deg * SIGN(1.0_r8,proj%truelat1) + ugrid = vtrue * SIN(alpha) + utrue * COS(alpha) + vgrid = vtrue * COS(alpha) - utrue * SIN(alpha) +!nc -- we added in a case structure for CASSINI and CYL + ELSEIF ((proj%code .EQ. PROJ_MERC).OR.(proj%code .EQ. PROJ_LATLON) & + .OR.(proj%code .EQ. PROJ_CASSINI).OR.(proj%code .EQ. PROJ_CYL))THEN + ugrid = utrue + vgrid = vtrue + ELSE + PRINT '(A)', 'Unrecognized map projection.' + STOP 'TRUEWIND_TO_GRIDWIND' + ENDIF + RETURN + END SUBROUTINE truewind_to_gridwind +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +END MODULE map_utils + diff --git a/models/wrf_unified/readme.rst b/models/wrf_unified/readme.rst new file mode 100644 index 000000000..393b41fbd --- /dev/null +++ b/models/wrf_unified/readme.rst @@ -0,0 +1,140 @@ +.. index:: wrf_chem, WRF-Chem, WRF CHEM + +.. _wrf_unified: + +WRF WRF-Chem Unified Model Interface +==================================== + +DART interface module for the Weather Research and Forecasting +`(WRF) `__ +model including the WRF-Chem extension. + +The model interface code supports WRF configurations with multiple domains. Data +for all domains is read into the DART state vector. During the computation of +the forward operators (getting the estimated observation values from each +ensemble member), the search starts in the domain with the highest number, which +is generally the finest nest or one of multiple finer nests. The search stops as +soon as a domain contains the observation location, working its way from largest +number to smallest number domain, ending with domain 1. For example, in a 4 +domain case the data in the state vector that came from ``wrfinput_d04`` is +searched first, then ``wrfinput_d03``, ``wrfinput_d02``, and finally +``wrfinput_d01``. + + +The forward operator is computed from the first (highest resolution) domain that contains the +lat/lon of the observation. During the assimilation phase, when the state values +are adjusted based on the correlations and assimilation increments, all points +in all domains that are within the localization radius are adjusted, regardless +of domain. + +The fields from WRF that are read into the DART state vector are controlled by +namelist. See below for the documentation on the &model_nml entries. The state +vector should include all fields needed to restart a WRF run. There may be +additional fields needed depending on the microphysics scheme selected. + +Namelist +-------- + +The ``&model_nml`` namelist is read from the ``input.nml`` file. 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:: text + + &model_nml + wrf_state_variables = 'NULL' + wrf_state_bounds = 'NULL' + chemistry_separate_file = .false + chem_state_variables = 'NULL' + chem_state_bounds = 'NULL' + num_domains = 1 + calendar_type = 3 # GREGORIAN + assimilation_period_seconds = 216001600 + sfc_elev_max_diff = - 1.0 + vert_localization_coord = 3 # VERTISHEIGHT + allow_perturbed_ics = .false. # testing purposes only + allow_obs_below_vol = .false. # Allow observations above the surface but below the lowest sigma level. + log_vert_interp = .true. # Do the interpolation of pressure values only after taking the log + log_horz_interpM = .false. + log_horz_interpQ = .false. + / + + + +Description of each namelist entry +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :header-rows: 1 + :widths: 25 15 60 + + * - Item + - Type + - Description + * - wrf_state_variables + - character(:,4) + - A 2D array of strings, 4 per wrf field to be added to the dart state vector. The 4 strings are: + + #. WRF field name - must match netcdf name exactly + #. DART QTY name - must match a valid DART QTY_xxx exactly + #. 'UPDATE' or 'NO_COPY_BACK'. If 'UPDATE', the data is written to netcdf file after the assimilation. If 'NO_COPY_BACK', the data is not written back to the wrf netcdf file after the assimilation. + #. A numeric string listing the domain numbers this array is part of. The special string 999 means all domains. For example, '12' means domains 1 and 2, '13' means 1 and 3. + + * - wrf_state_bounds + - character(:,3) + - A 2D array of strings, 3 per wrf array. During the writing of data to/from the wrf netcdf file, + variables listed here will have minimum and maximum values enforced. The 3 strings are: + + #. WRF field name - must match netcdf name exactly + #. Minimum -- specified as a string but must be a numeric value (e.g. '0.1'). Can be 'NULL' to allow any minimum value. + #. Maximum -- specified as a string but must be a numeric value (e.g. '0.1'). Can be 'NULL' to allow any maximum value. + + * - chemistry_separate_file + - logical + - If .true., chemistry fields are stored in a separate file. If .false., chemistry fields are included in the main state file. + * - chem_state_variables + - character(:,4) + - Chemistry state variables, same format as wrf_state_variables. 'NULL' disables chemistry variables. + * - chem_state_bounds + - character(:,3) + - Chemistry state bounds, same format as wrf_state_bounds. 'NULL' disables chemistry bounds. + * - num_domains + - integer + - Total number of WRF domains, including nested domains. + * - calendar_type + - integer + - Calendar type. Should be 3 (GREGORIAN) for WRF. + * - assimilation_period_seconds + - integer + - The time (in seconds) between assimilations. This is modified if necessary to be an integer multiple of the underlying model timestep. + * - sfc_elev_max_diff + - real(r8) + - TODO If > 0, the maximum difference, in meters, between an observation marked as a 'surface obs' as the vertical type (with the surface elevation, in meters, as the numerical vertical location), and the surface elevation as defined by the model. Observations further away from the surface than this threshold are rejected and not assimilated. If the value is negative, this test is skipped. + * - vert_localization_coord + - integer + - Vertical coordinate for vertical localization. + - 1 = model level + - 2 = pressure (in pascals) + - 3 = height (in meters) + - 4 = scale height (unitless) + * - allow_perturbed_ics + - logical + - Should not be used in most cases. Provided only for testing purposes to create a tiny ensemble for non-advancing tests. + * - allow_obs_below_vol + - logical + - If .false., observations above the surface but below the lowest sigma level are rejected. If .true., code will extrapolate downward from data values at levels 1 and 2. + * - log_vert_interp + - logical + - If .true., interpolation of pressure values is done after taking the log. + * - log_horz_interpM + - logical + - If .true., horizontal interpolation for M grid points is done after taking the log. + * - log_horz_interpQ + - logical + - If .true., horizontal interpolation for quad points is done after taking the log. + +References +---------- + +https://www2.mmm.ucar.edu/wrf/users/docs/user_guide_v4/contents.html diff --git a/models/wrf_unified/work/input.nml b/models/wrf_unified/work/input.nml new file mode 100644 index 000000000..dc5a27e3d --- /dev/null +++ b/models/wrf_unified/work/input.nml @@ -0,0 +1,557 @@ +&probit_transform_nml + / + +&algorithm_info_nml + qceff_table_filename = '' + / + +&perfect_model_obs_nml + read_input_state_from_file = .true. + single_file_in = .false. + input_state_files = 'wrfinput_d01' + init_time_days = -1 + init_time_seconds = -1 + + write_output_state_to_file = .false. + single_file_out = .false. + output_state_files = 'perfect_output_d01.nc' + output_interval = 1 + + obs_seq_in_file_name = "obs_seq.in" + obs_seq_out_file_name = "obs_seq.out" + first_obs_days = -1 + first_obs_seconds = -1 + last_obs_days = -1 + last_obs_seconds = -1 + + async = 0 + adv_ens_command = "../shell_scripts/advance_model.csh" + + trace_execution = .true. + output_timestamps = .false. + print_every_nth_obs = -1 + output_forward_op_errors = .true. + silence = .false. + / + +&filter_nml + async = 0, + adv_ens_command = "../shell_scripts/advance_model.csh", + ens_size = 3, + obs_sequence_in_name = "obs_seq.10", + obs_sequence_out_name = "obs_seq.final", + input_state_file_list = "input_wrf_d01.txt", "input_wrf_d02.txt", "input_wrf_d03.txt", "input_chem_d01.txt", "input_chem_d02.txt", "input_chem_d03.txt" + output_state_file_list = "output_wrf_d01.txt", "output_wrf_d02.txt", "output_wrf_d03.txt", "output_chem_d01.txt", "output_chem_d02.txt", "output_chem_d03.txt" + init_time_days = 154289, + init_time_seconds = 10801, + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + num_output_state_members = 3, + num_output_obs_members = 3, + output_interval = 1, + num_groups = 1, + distributed_state = .true. + compute_posterior = .true. + output_forward_op_errors = .true., + output_timestamps = .true., + trace_execution = .true., + + stages_to_write = 'output' + output_members = .false. + output_mean = .true. + output_sd = .true. + write_all_stages_at_end = .false. + + inf_flavor = 0, 0, + inf_initial_from_restart = .true., .false., + inf_sd_initial_from_restart = .true., .false., + inf_initial = 1.0, 1.00, + inf_lower_bound = 1.0, 1.0, + inf_upper_bound = 1000000.0, 1000000.0, + inf_damping = 0.9, 1.0, + inf_sd_initial = 0.6, 0.0, + inf_sd_lower_bound = 0.6, 0.0, + inf_sd_max_change = 1.05, 1.05, + / + +&quality_control_nml + input_qc_threshold = 3.0, + outlier_threshold = 3.0, + enable_special_outlier_code = .false. + / + +&fill_inflation_restart_nml + write_prior_inf = .true. + prior_inf_mean = 1.00 + prior_inf_sd = 0.6 + + write_post_inf = .false. + post_inf_mean = 1.00 + post_inf_sd = 0.6 + + input_state_files = 'wrfinput_d01', 'wrfinput_d02' + single_file = .false. + verbose = .false. + / + + +# cutoff is in radians; for the earth, 0.05 is about 300 km. +# cutoff is defined to be the half-width of the localization radius, +# so 0.05 radians for cutoff is about an 600 km effective +# localization radius, where the influence of an obs decreases +# to ~half at 300 km, and ~0 at the edges of the area. +&assim_tools_nml + cutoff = 0.05, + sort_obs_inc = .false., + spread_restoration = .false., + sampling_error_correction = .false., + adaptive_localization_threshold = -1, + output_localization_diagnostics = .false., + localization_diagnostics_file = 'localization_diagnostics', + convert_all_state_verticals_first = .true. + convert_all_obs_verticals_first = .true. + print_every_nth_obs = 1, + distribute_mean = .true. + / + +&cov_cutoff_nml + select_localization = 1, + / + +&obs_sequence_nml + write_binary_obs_sequence = .false., + / + +&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' + quantity_files = '../../../assimilation_code/modules/observations/default_quantities_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_radar_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_gps_mod.f90', + '../../../observations/forward_operators/obs_def_vortex_mod.f90', + '../../../observations/forward_operators/obs_def_gts_mod.f90' + / + +&obs_kind_nml + assimilate_these_obs_types = 'RADIOSONDE_TEMPERATURE', + 'RADIOSONDE_U_WIND_COMPONENT', + 'RADIOSONDE_V_WIND_COMPONENT', + 'SAT_U_WIND_COMPONENT', + 'SAT_V_WIND_COMPONENT', + 'AIRCRAFT_U_WIND_COMPONENT', + 'AIRCRAFT_V_WIND_COMPONENT', + 'AIRCRAFT_TEMPERATURE', + 'ACARS_U_WIND_COMPONENT', + 'ACARS_V_WIND_COMPONENT', + 'ACARS_TEMPERATURE', + 'GPSRO_REFRACTIVITY', + 'LAND_SFC_ALTIMETER' + evaluate_these_obs_types = 'RADIOSONDE_SPECIFIC_HUMIDITY', + / + +# Notes for obs_def_radar_mod_nml: +# (1) Reflectivity limit can be applied to observations and/or forward operator. +# (2) The default constants below match the WRF defaults. They will need to +# be changed for other cases, depending on which microphysics scheme is used. +# + +&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 = 2 , + allow_dbztowt_conv = .false., + 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, + / + + +# Notes for model_nml: +# (1) vert_localization_coord must be one of: +# 1 = model level +# 2 = pressure +# 3 = height +# 4 = scale height + +# supported in the DART state vector, and what settings should be used here. +# 'UPDATE' and 'NO_COPY_BACK' are supported in the 4th column; 'NO_UPDATE' is +# not yet supported. + +&model_nml + wrf_state_variables = 'U', 'QTY_U_WIND_COMPONENT', 'UPDATE','999', + 'V', 'QTY_V_WIND_COMPONENT', 'UPDATE','999', + 'W', 'QTY_VERTICAL_VELOCITY', 'UPDATE','999', + 'PH', 'QTY_GEOPOTENTIAL_HEIGHT', 'UPDATE','999', + 'T', 'QTY_POTENTIAL_TEMPERATURE','UPDATE','999', + 'MU', 'QTY_PRESSURE', 'UPDATE','999', + 'QVAPOR','QTY_VAPOR_MIXING_RATIO', 'UPDATE','999', + 'PSFC', 'QTY_SURFACE_PRESSURE', 'UPDATE','999', + wrf_state_bounds = 'QVAPOR','0.0','NULL', + 'QRAIN', '0.0','NULL', + 'QCLOUD','0.0','NULL', + chem_state_variables = 'o3', 'QTY_O3', 'UPDATE', '999', + 'no', 'QTY_NO', 'UPDATE', '999' + chem_state_bounds = 'o3', '0.0', 'NULL' + chemistry_separate_file = .true., + num_domains = 3, + calendar_type = 3, + assimilation_period_seconds = 21600, + vert_localization_coord = 3, + sfc_elev_max_diff = -1.0, + allow_obs_below_vol = .false. + / + +# vert_normalization_X is amount of X equiv to 1 radian in horiz. +# vert localization is 'cutoff' times the pressure/height/levels, +# only if horiz_dist_only is set to .false. in the namelist below. +# the default nlon/nlat should be good for most experiments. it sets +# an internal structure that speeds up searches. don't change it +# based on your grid size. nlon must be an odd number. +&location_nml + horiz_dist_only = .true., + vert_normalization_pressure = 6666666.7, + vert_normalization_height = 5000000.0, + vert_normalization_level = 2666.7, + vert_normalization_scale_height = 10.0, + approximate_distance = .false., + nlon = 71, + nlat = 36, + output_box_info = .false., + / + +&utilities_nml + TERMLEVEL = 1, + module_details = .false., + logfilename = 'dart_log.out', + nmlfilename = 'dart_log.nml', + write_nml = 'file', + / + +&mpi_utilities_nml + / + +®_factor_nml + select_regression = 1, + input_reg_file = "time_mean_reg", + save_reg_diagnostics = .false., + reg_diagnostics_file = "reg_diagnostics", + / + +# layout = 2 spreads the IO tasks across the nodes. +# This can greatly improve the performance in IO if +# tasks_per_node is set to match your hardware +&ensemble_manager_nml + layout = 2, + tasks_per_node = 16 + / + +&obs_def_gps_nml + max_gpsro_obs = 100000, + / + +&obs_def_tpw_nml + / + +# The times in the namelist for the obs_diag program are vectors +# that follow the following sequence: +# year month day hour minute second +# max_num_bins can be used to specify a fixed number of bins, +# in which case last_bin_center should be safely in the future. +# +# Acceptable latitudes range from [-90, 90] +# Acceptable longitudes range from [ 0, Inf] + +&obs_diag_nml + obs_sequence_name = 'obs_seq.final', + obs_sequence_list = '', + first_bin_center = 2003, 1, 1, 0, 0, 0 , + last_bin_center = 2003, 1, 2, 0, 0, 0 , + bin_separation = 0, 0, 0,12, 0, 0 , + bin_width = 0, 0, 0, 6, 0, 0 , + time_to_skip = 0, 0, 0, 0, 0, 0 , + max_num_bins = 1000, + trusted_obs = 'null', + Nregions = 4, + lonlim1 = 0.0, 0.0, 0.0, 235.0, + lonlim2 = 360.0, 360.0, 360.0, 295.0, + latlim1 = 20.0, -80.0, -20.0, 25.0, + latlim2 = 80.0, -20.0, 20.0, 55.0, + reg_names = 'Northern Hemisphere', 'Southern Hemisphere', 'Tropics', 'North America', + print_mismatched_locs = .false., + create_rank_histogram = .true., + outliers_in_histogram = .true., + use_zero_error_obs = .false., + verbose = .false. + / + +&schedule_nml + calendar = 'Gregorian', + first_bin_start = 1601, 1, 1, 0, 0, 0, + first_bin_end = 2999, 1, 1, 0, 0, 0, + last_bin_end = 2999, 1, 1, 0, 0, 0, + bin_interval_days = 1000000, + bin_interval_seconds = 0, + max_num_bins = 1000, + print_table = .true., + / + +&obs_seq_to_netcdf_nml + obs_sequence_name = 'obs_seq.final', + obs_sequence_list = '', + append_to_netcdf = .false., + lonlim1 = 0.0, + lonlim2 = 360.0, + latlim1 = -90.0, + latlim2 = 90.0, + verbose = .false., + / + +# There is one GIGANTIC difference between the obsdef_mask.txt and .nc +# The netCDF file intentionally ignores the effect of nTmin/nTmax. +# The netCDF file has ALL matching stations, regardless of temporal coverage. + +&obs_seq_coverage_nml + obs_sequences = '' + obs_sequence_list = 'obs_coverage_list.txt' + obs_of_interest = 'METAR_U_10_METER_WIND' + textfile_out = 'METAR_U_10_METER_WIND_obsdef_mask.txt' + netcdf_out = 'METAR_U_10_METER_WIND_obsdef_mask.nc' + first_analysis = 2003, 1, 1, 0, 0, 0 + last_analysis = 2003, 1, 2, 0, 0, 0 + forecast_length_days = 1 + forecast_length_seconds = 0 + verification_interval_seconds = 21600 + temporal_coverage_percent = 100.0 + lonlim1 = 0.0 + lonlim2 = 360.0 + latlim1 = -90.0 + latlim2 = 90.0 + verbose = .true. + / + +# selections_file is a list of obs_defs output +# from the obs_seq_coverage utility. + +&obs_selection_nml + filename_seq = 'obs_seq.out', + filename_seq_list = '', + filename_out = 'obs_seq.processed', + selections_file = 'obsdef_mask.txt', + selections_is_obs_seq = .false., + print_only = .false., + calendar = 'gregorian', + / + +&obs_seq_verify_nml + obs_sequences = '' + obs_sequence_list = 'obs_verify_list.txt' + input_template = 'obsdef_mask.nc' + netcdf_out = 'forecast.nc' + obtype_string = 'METAR_U_10_METER_WIND' + print_every = 10000 + verbose = .true. + debug = .false. + / + +&obs_sequence_tool_nml + num_input_files = 1, + filename_seq = 'obs_seq.out', + filename_out = 'obs_seq.processed', + first_obs_days = -1, + first_obs_seconds = -1, + last_obs_days = -1, + last_obs_seconds = -1, + obs_types = '', + keep_types = .false., + print_only = .false., + min_lat = -90.0, + max_lat = 90.0, + min_lon = 0.0, + max_lon = 360.0, + / + +&replace_wrf_fields_nml + debug = .false., + fail_on_missing_field = .false., + fieldnames = "SNOWC", + "ALBBCK", + "TMN", + "TSK", + "SH2O", + "SMOIS", + "SEAICE", + "HGT_d01", + "TSLB", + "SST", + "SNOWH", + "SNOW", + fieldlist_file = '', + / + +&obs_common_subset_nml + num_to_compare_at_once = 2, + filename_seq = 'obs_seq1.final', 'obs_seq2.final', + filename_seq_list = '', + filename_out_suffix = '.common' , + calendar = 'Gregorian', + print_every = 1000, + dart_qc_threshold = 3, + print_only = .false., + / + +&wrf_dart_to_fields_nml + include_slp = .true., + include_wind_components = .true., + include_height_on_pres = .true., + include_temperature = .true., + include_rel_humidity = .true., + include_surface_fields = .false., + include_sat_ir_temp = .false., + pres_levels = 70000., + / + +&ncepobs_nml + year = 2010, + month = 06, + day = 00, + tot_days = 1, + max_num = 1000000, + ObsBase = 'temp_obs.', + select_obs = 0, + ADPUPA = .false., + AIRCAR = .false., + AIRCFT = .false., + SATEMP = .false., + SFCSHP = .false., + ADPSFC = .false., + SATWND = .true., + obs_U = .false., + obs_V = .false., + obs_T = .false., + obs_PS = .false., + obs_QV = .false., + daily_file = .true., + obs_time = .false., + lat1 = 10.00, + lat2 = 60.00, + lon1 = 210.0, + lon2 = 300.0 + / + +&prep_bufr_nml + obs_window_upa = 1.0, + obs_window_air = 1.0, + obs_window_cw = 1.0, + otype_use = 242.0, 243.0, 245.0, 246.0, 251.0, 252.0, 253.0, 257.0, 259.0 + qctype_use = 0, 1, 2, 3, 4, 9, 15 + / + +&convert_cosmic_gps_nml + gpsro_netcdf_file = '', + gpsro_netcdf_filelist = 'flist', + gpsro_out_file = 'obs_seq.gpsro', + local_operator = .true., + obs_levels = 0.22, 0.55, 1.1, 1.8, 2.7, 3.7, 4.9, + 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, + ray_ds = 5000.0, + ray_htop = 13000.1, + / + +&wrf_obs_preproc_nml + + file_name_input = 'obs_seq20110901' + file_name_output = 'obs_seq.europe.prev' + + overwrite_obs_time = .false. + + obs_boundary = 0.0 + increase_bdy_error = .false. + maxobsfac = 2.5 + obsdistbdy = 1.0 + + sfc_elevation_check = .false. + sfc_elevation_tol = 3000.0 + obs_pressure_top = 0.0 + obs_height_top = 2.0e10 + + include_sig_data = .true. + tc_sonde_radii = -1.0 + + superob_aircraft = .true. + aircraft_horiz_int = 800.0 + aircraft_pres_int = 25000.0 + + superob_sat_winds = .true. + sat_wind_horiz_int = 800.0 + sat_wind_pres_int = 25000.0 + + overwrite_ncep_satwnd_qc = .false. + overwrite_ncep_sfc_qc = .false. + / + +! sonde_extra = 'obs_seq.rawin' +! land_sfc_extra = 'obs_seq.land_sfc' +! metar_extra = 'obs_seq.metar' +! marine_sfc_extra = 'obs_seq.marine' +! sat_wind_extra = 'obs_seq.satwnd' +! profiler_extra = 'obs_seq.profiler' +! gpsro_extra = 'obs_seq.gpsro' +! acars_extra = 'obs_seq.acars' +! trop_cyclone_extra = 'obs_seq.tc' + +&state_vector_io_nml + single_precision_output = .true., + / + +&compare_states_nml + / + +&closest_member_tool_nml + input_restart_file_list = 'input_file_list_d01.txt', + output_file_name = 'closest_results.txt' + ens_size = 3, + single_restart_file_in = .false., + difference_method = 4, + use_only_qtys = 'QTY_U_WIND_COMPONENT' + / + +# To test both domains, you must change 'model_nml:num_domains = 2' + +&model_mod_check_nml + input_state_files = 'wrfinput_d01', 'wrfinput_d02', 'wrfinput_d03', 'wrfchem_d01', 'wrfchem_d02','wrfchem_d03' + output_state_files = 'mmc_output1.nc', 'mmc_output2.nc', 'mmc_output3.nc', 'mmc_output4.nc', 'mmc_output5.nc', 'mmc_output6.nc' + test1thru = 4 + run_tests = 1,2,3,4,5 + x_ind = 87370 + loc_of_interest = 231.0, 40.0, 500.0 + quantity_of_interest = 'QTY_POTENTIAL_TEMPERATURE' + interp_test_dlon = 0.1 + interp_test_dlat = 0.1 + interp_test_dvert = 1000.0 + interp_test_lonrange = 250.0, 260.0 + interp_test_latrange = 30.0, 45.0 + interp_test_vertrange = 2000.0, 4000.0 + interp_test_vertcoord = 'VERTISHEIGHT' + verbose = .false. + / + diff --git a/models/wrf_unified/work/quickbuild.sh b/models/wrf_unified/work/quickbuild.sh new file mode 100755 index 000000000..8853c48f7 --- /dev/null +++ b/models/wrf_unified/work/quickbuild.sh @@ -0,0 +1,49 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL=wrf_unified +LOCATION=threed_sphere + +programs=( +closest_member_tool +filter +model_mod_check +perfect_model_obs +) + +serial_programs=( +#radiance_obs_seq_to_netcdf # needs rttov +fill_inflation_restart +obs_diag +obs_sequence_tool +) + + +model_serial_programs=( +) + +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build DART +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@"