Skip to content

Commit

Permalink
Obs Converter Clean up
Browse files Browse the repository at this point in the history 
Did a cleanup of the obs converter:
- Removed comments
- Removed unused variables
- Improved style
- Better use of DART's modules
- Improved initialization and functions
- Ran and tested the converter after the changes were made
  • Loading branch information
@mgharamti
mgharamti committed Sep 4, 2024
1 parent 48d6d59 commit e3e1b93
Show file tree
Hide file tree
Showing 5 changed files with 79 additions and 341 deletions.
15 changes: 0 additions & 15 deletions assimilation_code/modules/observations/default_quantities_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,22 +387,7 @@
! QTY_WIND_TURBINE_POWER
! QTY_W_CURRENT_COMPONENT
! QTY_X_LAMBDA
! QTY_NITRATE_CONCENTRATION
! QTY_PHOSPHATE_CONCENTRATION
! QTY_DISSOLVED_OXYGEN
! QTY_PHYTOPLANKTON_BIOMASS
! QTY_ALKALINITY
! QTY_DISSOLVED_INORGANIC_CARBON
! QTY_DISSOLVED_ORGANIC_CARBON
! QTY_DISSOLVED_ORGANIC_NITROGEN
! QTY_DISSOLVED_ORGANIC_P
! QTY_DISSOLVED_INORGANIC_IRON
! QTY_SURFACE_CHLOROPHYLL
! QTY_LAYER_THICKNESS
! QTY_COLUMN_DEPTH
! QTY_DISSOLVED_INORGANIC_SIO3
! QTY_MESOZOOPLANKTON_CARBON
! QTY_MICROZOOPLANKTON_CARBON
! QTY_BGC_PARAM
!
! END DART PREPROCESS QUANTITY DEFINITIONS
Expand Down
80 changes: 33 additions & 47 deletions observations/obs_converters/BATS/bats_climatology.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import os
#!/usr/bin/env python3
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
Expand All @@ -8,23 +8,23 @@
#################################################################

# input and output datafiles
input_file_path = 'data/bats_bottle.txt'
output_txt_path = 'data/bats_climatology.txt'
output_nc_path = 'data/bats_climatology.nc'
input_file_path = 'bats_bottle.txt'
output_txt_path = 'bats_climatology.txt'
output_nc_path = 'bats_climatology.nc'

first_data_line = 61 # file is read starting at this line (so that header information is ignored)
month_columns = [18, 19] # first and last columns in the data file containing month values
depth_columns = [58, 63] # first and last columns in the data file containing depth values
depth_columns = [64, 69] # first and last columns in the data file containing depth values
num_obs_types = 6 # number of variables being read from the datafile
max_num_data = 16660 # upper bound on the number of data points to be processed

# first and last columns in the data file containing values for each variable
obs_val_columns = [[102, 107], # O2
[126, 131], # CO2
[134, 139], # Alk
[141, 147], # NO31
[158, 164], # PO41
[166, 172]] # Si1
obs_val_columns = [[113, 119], # O2
[137, 143], # CO2
[145, 151], # Alk
[153, 159], # NO31
[170, 176], # PO41
[178, 184]] # Si1

# layer interface depths for data in prog_z.nc output file from MARBL, run "ncdump -v "z_i" prog_z.nc" to generate.
marbl_depths = [0, 5, 15, 25, 40, 62.5, 87.5, 112.5, 137.5, 175, 225, 275, 350, 450,
Expand All @@ -45,57 +45,43 @@
sq_means = np.zeros((num_obs_types, 12, len(marbl_depths)))
sample_counts = np.zeros((num_obs_types, 12, len(marbl_depths)), dtype=int)

bats_data = open(input_file_path, "r")
line_number = 0

for line in bats_data.readlines():
line_number += 1

if(line_number < first_data_line):
continue
with open(input_file_path, "r") as bats_data:
for line in bats_data.readlines()[first_data_line-1:]:

# extracting the depth and month where this observation was taken

depth_val = float(line[depth_columns[0]-1:depth_columns[1]])
month_index = int(line[month_columns[0]-1:month_columns[1]]) - 1
# extracting the depth and month where this observation was taken
depth_val = float(line[depth_columns[0]-1:depth_columns[1]])
month_index = int(line[month_columns[0]-1:month_columns[1]]) - 1

# assigning this observation to one of the MARBL layers

depth_index = len(marbl_depths) - 1
# assigning this observation to one of the MARBL layers
depth_index = len(marbl_depths) - 1

while((depth_val < marbl_depths[depth_index]) and (depth_index > 0)):
depth_index -= 1
while((depth_val < marbl_depths[depth_index]) and (depth_index > 0)):
depth_index -= 1

# adding this observation into the climatology

for obs_type_index in range(num_obs_types):
obs_val = float(line[obs_val_columns[obs_type_index][0]-1:obs_val_columns[obs_type_index][1]])
# adding this observation into the climatology
for obs_type_index in range(num_obs_types):
obs_val = float(line[obs_val_columns[obs_type_index][0]-1:obs_val_columns[obs_type_index][1]])

if(abs(obs_val + 999.0) > 1e-8): # this filters out 'missing' values
prev_mean = clim_means[obs_type_index, month_index, depth_index]
prev_sqmean = sq_means[obs_type_index, month_index, depth_index]
prev_samples = sample_counts[obs_type_index, month_index, depth_index]
if(abs(obs_val + 999.0) > 1e-8): # this filters out 'missing' values
prev_mean = clim_means[obs_type_index, month_index, depth_index]
prev_sqmean = sq_means[obs_type_index, month_index, depth_index]
prev_samples = sample_counts[obs_type_index, month_index, depth_index]

new_mean = (prev_samples*prev_mean + obs_val)/(prev_samples + 1)
new_sqmean = (prev_samples*prev_sqmean + obs_val**2)/(prev_samples + 1)
new_mean = (prev_samples*prev_mean + obs_val)/(prev_samples + 1)
new_sqmean = (prev_samples*prev_sqmean + obs_val**2)/(prev_samples + 1)

clim_means[obs_type_index, month_index, depth_index] = new_mean
sq_means[obs_type_index, month_index, depth_index] = new_sqmean
clim_means[obs_type_index, month_index, depth_index] = new_mean
sq_means[obs_type_index, month_index, depth_index] = new_sqmean

sample_counts[obs_type_index, month_index, depth_index] += 1
sample_counts[obs_type_index, month_index, depth_index] += 1

print("month =",month_index,", type =",obs_type_index,", depth =",marbl_depths[depth_index],", value =",obs_val)
print("month =",month_index,", type =",obs_type_index,", depth =",marbl_depths[depth_index],", value =",obs_val)

bats_data.close()

# setting up a CSV file to write the climatology into; this will be read by the DART data converter.

os.system("rm -f "+output_txt_path)
output_txt = open(output_txt_path, "w")

# setting up a netCDF file to write the climatology into; this will be useful for creating MARBL-DART diagnostics.

os.system("rm -f "+output_nc_path)
output_nc = nc.Dataset(output_nc_path, "w")

output_nc.createDimension("Month")
Expand Down
151 changes: 15 additions & 136 deletions observations/obs_converters/BATS/bats_to_clim_obs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,30 +2,29 @@
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

! This file is meant to read a text file containing bottle data from the
! Bermuda Atlantic Time-Series Study (https://bats.bios.asu.edu/), which
! is converted to climatological estimates with monthly resolution.

program bats_to_clim_obs

use types_mod, only : r8, PI, DEG2RAD
use types_mod, only : r8
use utilities_mod, only : initialize_utilities, finalize_utilities, &
open_file, close_file, &
find_namelist_in_file, check_namelist_read, &
error_handler, E_ERR, E_MSG, nmlfileunit, &
do_nml_file, do_nml_term
use time_manager_mod, only : time_type, set_calendar_type, set_date, &
operator(>=), increment_time, get_time, &
operator(-), NOLEAP, GREGORIAN, operator(+), &
operator(-), GREGORIAN, operator(+), &
print_date
use location_mod, only : VERTISHEIGHT, VERTISPRESSURE
use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
static_init_obs_sequence, init_obs, write_obs_seq, &
init_obs_sequence, get_num_obs, set_copy_meta_data, &
set_qc_meta_data, destroy_obs_sequence

use obs_utilities_mod, only : create_3d_obs, add_obs_to_seq
use obs_kind_mod, only : BATS_OXYGEN, BATS_INORGANIC_CARBON, BATS_ALKALINITY, &
BATS_NITRATE, BATS_PHOSPHATE, BATS_SILICATE

Expand All @@ -40,31 +39,29 @@ program bats_to_clim_obs
BATS_PHOSPHATE, BATS_SILICATE/)

real(r8), parameter :: MIN_OBS_ERROR = 0.1_r8
real(r8), parameter :: bats_lon = 360.0_r8 - 64.0_r8
real(r8), parameter :: bats_lat = 31.0_r8
real(r8), parameter :: BATS_LON = 360.0_r8 - 64.0_r8
real(r8), parameter :: BATS_LAT = 31.0_r8

! namelist variables, changeable at runtime
character(len=256) :: text_input_file, obs_out_dir
integer :: max_lines
real(r8) :: obs_err_var_inflation
logical :: debug
character(len=256) :: text_input_file = 'bats_climatology.txt'
character(len=256) :: obs_out_dir = 'obs_seq_files'
integer :: max_lines = 3000
real(r8) :: obs_err_var_inflation = 10.0
logical :: debug = .true.

namelist /bats_to_clim_obs_nml/ text_input_file, max_lines, obs_err_var_inflation, obs_out_dir, debug

! local variables
character (len=294) :: input_line, obs_out_file
character (len=6) :: month_str

integer :: oday, osec, rcio, iunit, otype, char_index, comma_index, line_number, otype_index, &
year, month, month_old, day, hour, minute, second, num_copies, num_qc, max_obs
integer :: comma_locs(4)

integer :: oday, osec, rcio, iunit, char_index, comma_index, line_number, otype_index, &
month, month_old, num_copies, num_qc, max_obs
integer :: comma_locs(4)
logical :: first_obs, new_obs_seq

real(r8) :: qc, obs_err
real(r8) :: lat, lon, vert, ovalue

! the uncertainties corresponding to the observations above
real(r8) :: vert, ovalue

type(obs_sequence_type) :: obs_seq
type(obs_type) :: obs, prev_obs
Expand Down Expand Up @@ -229,7 +226,7 @@ program bats_to_clim_obs
end if

! adding the observation
call create_3d_obs(bats_lat, bats_lon, vert, VERTISHEIGHT, &
call create_3d_obs(BATS_LAT, BATS_LON, vert, VERTISHEIGHT, &
ovalue, OTYPE_ORDERING(otype_index), obs_err, &
oday, osec, qc, obs)

Expand All @@ -247,122 +244,4 @@ program bats_to_clim_obs
! end of main program
call finalize_utilities()

contains


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! create_3d_obs - subroutine that is used to create an observation
! type from observation data.
!
! NOTE: assumes the code is using the threed_sphere locations module,
! that the observation has a single data value and a single
! qc value, and that this obs type has no additional required
! data (e.g. gps and radar obs need additional data per obs)
!
! inputs:
! lat - latitude of observation
! lon - longitude of observation
! vval - vertical coordinate
! vkind - kind of vertical coordinate (pressure, level, etc)
! obsv - observation value
! otype - observation type
! oerr - observation error
! day - gregorian day
! sec - gregorian second
! qc - quality control value
! outputs:
! obs - observation type
!
! created Oct. 2007 Ryan Torn, NCAR/MMM
! adapted for more generic use 11 Mar 2010, nancy collins, ncar/image
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_3d_obs(lat, lon, vval, vkind, obsv, otype, oerr, day, sec, qc, obs)
use types_mod, only : r8
use obs_def_mod, only : obs_def_type, set_obs_def_time, set_obs_def_type_of_obs, &
set_obs_def_error_variance, set_obs_def_location
use obs_sequence_mod, only : obs_type, set_obs_values, set_qc, set_obs_def
use time_manager_mod, only : time_type, set_time
use location_mod, only : set_location

integer, intent(in) :: otype, vkind, day, sec
real(r8), intent(in) :: lat, lon, vval, obsv, oerr, qc
type(obs_type), intent(inout) :: obs

real(r8) :: obs_val(1), qc_val(1)
type(obs_def_type) :: obs_def

call set_obs_def_location(obs_def, set_location(lon, lat, vval, vkind))
call set_obs_def_type_of_obs(obs_def, otype)
call set_obs_def_time(obs_def, set_time(sec, day))
call set_obs_def_error_variance(obs_def, oerr * oerr)
call set_obs_def(obs, obs_def)

obs_val(1) = obsv
call set_obs_values(obs, obs_val)
qc_val(1) = qc
call set_qc(obs, qc_val)

end subroutine create_3d_obs


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! add_obs_to_seq -- adds an observation to a sequence. inserts if first
! obs, inserts with a prev obs to save searching if that's possible.
!
! seq - observation sequence to add obs to
! obs - observation, already filled in, ready to add
! obs_time - time of this observation, in dart time_type format
! prev_obs - the previous observation that was added to this sequence
! (will be updated by this routine)
! prev_time - the time of the previously added observation
! (will also be updated by this routine)
! first_obs - should be initialized to be .true., and then will be
! updated by this routine to be .false. after the first obs
! has been added to this sequence.
!
! created Mar 8, 2010 nancy collins, ncar/image
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine add_obs_to_seq(seq, obs, obs_time, prev_obs, prev_time, first_obs)
use types_mod, only : r8
use obs_sequence_mod, only : obs_sequence_type, obs_type, insert_obs_in_seq
use time_manager_mod, only : time_type, operator(>=)

type(obs_sequence_type), intent(inout) :: seq
type(obs_type), intent(inout) :: obs, prev_obs
type(time_type), intent(in) :: obs_time
type(time_type), intent(inout) :: prev_time
logical, intent(inout) :: first_obs

! insert(seq,obs) always works (i.e. it inserts the obs in
! proper time format) but it can be slow with a long file.
! supplying a previous observation that is older (or the same
! time) as the new one speeds up the searching a lot.

if(first_obs) then ! for the first observation, no prev_obs
call insert_obs_in_seq(seq, obs)
first_obs = .false.
else
if(obs_time >= prev_time) then ! same time or later than previous obs
call insert_obs_in_seq(seq, obs, prev_obs)
else ! earlier, search from start of seq
call insert_obs_in_seq(seq, obs)
endif
endif

! update for next time
prev_obs = obs
prev_time = obs_time

end subroutine add_obs_to_seq

end program bats_to_clim_obs

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$
Loading

0 comments on commit e3e1b93

Please sign in to comment.