Skip to content

Commit

Permalink
Merge pull request #25 from matt-frey/add-statistics-netcdf
Browse files Browse the repository at this point in the history
Add field diagnostics
  • Loading branch information
matt-frey authored Nov 16, 2022
2 parents 4403444 + 95864e7 commit 3206ca1
Show file tree
Hide file tree
Showing 5 changed files with 271 additions and 43 deletions.
89 changes: 47 additions & 42 deletions src/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -6,76 +6,81 @@ AM_FCFLAGS = \
-I $(top_builddir)/src/netcdf

bin_PROGRAMS = ps3d genspec genspec2d beltrami production
ps3d_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
utils/merge_sort.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
advance.f90 \
ps3d_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
utils/merge_sort.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
fields/field_diagnostics_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
advance.f90 \
ps3d.f90

ps3d_LDADD = \
$(top_builddir)/src/utils/libps_utils.la \
$(top_builddir)/src/fft/libps_fft.la \
$(top_builddir)/src/netcdf/libps_netcdf.la

genspec_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
genspec_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
fields/field_diagnostics_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
genspec.f90

genspec_LDADD = \
$(top_builddir)/src/utils/libps_utils.la \
$(top_builddir)/src/fft/libps_fft.la \
$(top_builddir)/src/netcdf/libps_netcdf.la

genspec2d_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
genspec2d_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
fields/field_diagnostics_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
genspec2d.f90

genspec2d_LDADD = \
$(top_builddir)/src/utils/libps_utils.la \
$(top_builddir)/src/fft/libps_fft.la \
$(top_builddir)/src/netcdf/libps_netcdf.la

beltrami_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
beltrami_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
fields/field_diagnostics_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
beltrami.f90

beltrami_LDADD = \
$(top_builddir)/src/utils/libps_utils.la \
$(top_builddir)/src/fft/libps_fft.la \
$(top_builddir)/src/netcdf/libps_netcdf.la

production_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
production_SOURCES = \
utils/options.f90 \
utils/parameters.f90 \
inversion/inversion_utils.f90 \
fields/fields.f90 \
fields/field_netcdf.f90 \
fields/field_diagnostics_netcdf.f90 \
inversion/inversion.f90 \
utils/utils.f90 \
production.f90

production_LDADD = \
Expand Down
205 changes: 205 additions & 0 deletions src/fields/field_diagnostics_netcdf.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
! =============================================================================
! Write field diagnostics to NetCDF.
! =============================================================================
module field_diagnostics_netcdf
use fields
use inversion_utils
use netcdf_utils
use netcdf_writer
use netcdf_reader
use constants, only : one
use parameters, only : lower, extent, nx, ny, nz
use config, only : package_version, cf_version
use timer, only : start_timer, stop_timer
use options, only : write_netcdf_options
use physics, only : write_physical_quantities
implicit none

private

character(len=512) :: ncfname
integer :: ncid
integer :: t_axis_id, t_dim_id, n_writes, ke_id, en_id
#ifdef ENABLE_BUOYANCY
integer :: pe_id, bmax_id, bmin_id
#endif

double precision :: restart_time

integer :: field_stats_io_timer

public :: create_netcdf_field_stats_file, &
write_netcdf_field_stats, &
field_stats_io_timer


contains

! Create the NetCDF field diagnostic file.
! @param[in] basename of the file
! @param[in] overwrite the file
subroutine create_netcdf_field_stats_file(basename, overwrite)
character(*), intent(in) :: basename
logical, intent(in) :: overwrite
logical :: l_exist

ncfname = basename // '_field_stats.nc'

restart_time = -one
n_writes = 1

call exist_netcdf_file(ncfname, l_exist)

if (l_exist) then
call open_netcdf_file(ncfname, NF90_NOWRITE, ncid)
call get_num_steps(ncid, n_writes)
call get_time(ncid, restart_time)
call read_netcdf_field_stats_content
call close_netcdf_file(ncid)
n_writes = n_writes + 1
return
endif

call create_netcdf_file(ncfname, overwrite, ncid)

call write_netcdf_info(ncid=ncid, &
ps3d_version=package_version, &
file_type='field_stats', &
cf_version=cf_version)

call write_netcdf_box(ncid, lower, extent, (/nx, ny, nz/))

call write_physical_quantities(ncid)

call write_netcdf_options(ncid)

call define_netcdf_temporal_dimension(ncid, t_dim_id, t_axis_id)

call define_netcdf_dataset( &
ncid=ncid, &
name='ke', &
long_name='kinetic energy', &
std_name='', &
unit='m^5/s^2', &
dtype=NF90_DOUBLE, &
dimids=(/t_dim_id/), &
varid=ke_id)

call define_netcdf_dataset( &
ncid=ncid, &
name='en', &
long_name='enstrophy', &
std_name='', &
unit='m^3/s^2', &
dtype=NF90_DOUBLE, &
dimids=(/t_dim_id/), &
varid=en_id)

#ifdef ENABLE_BUOYANCY
call define_netcdf_dataset( &
ncid=ncid, &
name='pe', &
long_name='potential energy', &
std_name='', &
unit='m^5/s^2', &
dtype=NF90_DOUBLE, &
dimids=(/t_dim_id/), &
varid=pe_id)

call define_netcdf_dataset( &
ncid=ncid, &
name='min_buoyancy', &
long_name='minimum gridded buoyancy', &
std_name='', &
unit='m/s^2', &
dtype=NF90_DOUBLE, &
dimids=(/t_dim_id/), &
varid=bmin_id)

call define_netcdf_dataset( &
ncid=ncid, &
name='max_buoyancy', &
long_name='maximum gridded buoyancy', &
std_name='', &
unit='m/s^2', &
dtype=NF90_DOUBLE, &
dimids=(/t_dim_id/), &
varid=bmax_id)
#endif
call close_definition(ncid)

end subroutine create_netcdf_field_stats_file

! Pre-condition: Assumes an open file
subroutine read_netcdf_field_stats_content

call get_dim_id(ncid, 't', t_dim_id)

call get_var_id(ncid, 't', t_axis_id)

call get_var_id(ncid, 'ke', ke_id)

call get_var_id(ncid, 'en', en_id)

#ifdef ENABLE_BUOYANCY
call get_var_id(ncid, 'pe', pe_id)

call get_var_id(ncid, 'min_buoyancy', bmin_id)

call get_var_id(ncid, 'max_buoyancy', bmax_id)
#endif
end subroutine read_netcdf_field_stats_content

! Write a step in the field diagnostic file.
! @param[in] t is the time
! @param[in] dt is the time step
subroutine write_netcdf_field_stats(t)
double precision, intent(in) :: t
double precision :: ke, en
#ifdef ENABLE_BUOYANCY
double precision :: bmin, bmax, pe
#endif

call start_timer(field_stats_io_timer)

if (t <= restart_time) then
call stop_timer(field_stats_io_timer)
return
endif

ke = get_kinetic_energy()
en = get_enstrophy()

#ifdef ENABLE_BUOYANCY
call field_combine_physical(sbuoy, buoy)
pe = get_potential_energy()
bmin = minval(buoy)
bmax = maxval(buoy)
#endif

call open_netcdf_file(ncfname, NF90_WRITE, ncid)

! write time
call write_netcdf_scalar(ncid, t_axis_id, t, n_writes)

!
! write diagnostics
!
call write_netcdf_scalar(ncid, ke_id, ke, n_writes)
call write_netcdf_scalar(ncid, en_id, en, n_writes)
#ifdef ENABLE_BUOYANCY
call write_netcdf_scalar(ncid, pe_id, pe, n_writes)
call write_netcdf_scalar(ncid, bmin_id, bmin, n_writes)
call write_netcdf_scalar(ncid, bmax_id, bmax, n_writes)
#endif

! increment counter
n_writes = n_writes + 1

call close_netcdf_file(ncid)

call stop_timer(field_stats_io_timer)

end subroutine write_netcdf_field_stats

end module field_diagnostics_netcdf
2 changes: 2 additions & 0 deletions src/ps3d.f90
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ program ps3d
use timer
use fields
use field_netcdf, only : field_io_timer, read_netcdf_fields
use field_diagnostics_netcdf, only : field_stats_io_timer
use inversion_mod, only : vor2vel_timer &
, vtend_timer &
, vor2vel &
Expand Down Expand Up @@ -48,6 +49,7 @@ subroutine pre_run

call register_timer('ps', ps_timer)
call register_timer('field I/O', field_io_timer)
call register_timer('field diagnostics I/O', field_stats_io_timer)
call register_timer('vor2vel', vor2vel_timer)
call register_timer('vorticity tendency', vtend_timer)
call register_timer('advance', advance_timer)
Expand Down
4 changes: 4 additions & 0 deletions src/utils/options.f90
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ module options
type info
double precision :: field_freq = one
logical :: write_fields = .true.
double precision :: field_stats_freq = one
logical :: write_field_stats = .true.
logical :: overwrite = .false.
character(len=512) :: basename = ''
end type info
Expand Down Expand Up @@ -129,6 +131,8 @@ subroutine write_netcdf_options(ncid)

call write_netcdf_attribute(ncid, "field_freq", output%field_freq)
call write_netcdf_attribute(ncid, "write_fields", output%write_fields)
call write_netcdf_attribute(ncid, "field_stats_freq", output%field_stats_freq)
call write_netcdf_attribute(ncid, "write_field_stats", output%write_field_stats)
call write_netcdf_attribute(ncid, "overwrite", output%overwrite)
call write_netcdf_attribute(ncid, "basename", trim(output%basename))

Expand Down
Loading

0 comments on commit 3206ca1

Please sign in to comment.