Skip to content

Commit

Permalink
Improve specification of grid definition file,variable used for
Browse files Browse the repository at this point in the history 
diagnostic remapping. The format is now:
        DIAG_REMAP_Z_GRID_DEF = "FILE:OM3_zgrid.nc,zw"

Also use library to read netCDF files. mom-ocean#62
  • Loading branch information
@nichannah
nichannah committed Jun 5, 2015
1 parent 70e5686 commit d556620
Showing 1 changed file with 99 additions and 159 deletions.
258 changes: 99 additions & 159 deletions src/framework/MOM_diag_mediator.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,9 @@ module MOM_diag_mediator
use MOM_error_handler, only : MOM_error, FATAL, is_root_pe
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_io, only : slasher, vardesc
use MOM_io, only : file_exists, field_exists, field_size
use MOM_io, only : slasher, vardesc, mom_read_data
use MOM_string_functions, only : extractWord
use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc
use MOM_string_functions, only : lowercase
use MOM_time_manager, only : time_type
Expand Down Expand Up @@ -138,10 +140,9 @@ module MOM_diag_mediator
real :: missing_value = 1.0e+20

! Interface locations for Z remapping
real, dimension(:), pointer :: z_remap_int => null()
! Difference between model interfaces and remap interfaces, this only needs to
! be calculated when thicknesses change.
real, dimension(:,:,:), allocatable :: dz_remap_int
real, dimension(:), allocatable :: zi_remap
! Layer locations for Z remapping
real, dimension(:), allocatable :: zl_remap
! Number of z levels used for remapping
integer :: nz_remap

Expand Down Expand Up @@ -170,10 +171,11 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical)

integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_zzl, id_zzi
integer :: k, nz
integer :: nzi(4)
real :: zlev(G%ks:G%ke), zinter(G%ks:G%ke+1)
logical :: set_vert, Cartesian_grid
character(len=80) :: grid_config, units_temp
character(len=200) :: in_dir, zgrid_file ! strings for directory/file
character(len=200) :: inputdir, string, filename, varname, dimname

! This include declares and sets the variable "version".
#include "version_variable.h"
Expand Down Expand Up @@ -244,24 +246,60 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical)
endif

! Read info needed for z-space remapping
call get_param(param_file, mod, "Z_OUTPUT_GRID_FILE", zgrid_file, &
"The file that specifies the vertical grid for \n"//&
"depth-space diagnostics, or blank to disable \n"//&
"depth-space output.", default="")
if (len_trim(zgrid_file) > 0) then
call get_param(param_file, mod, "INPUTDIR", in_dir, &
"The directory in which input files are found.", default=".")
in_dir = slasher(in_dir)
call get_Z_output_grid(trim(in_dir)//trim(zgrid_file), "zw", diag_cs%z_remap_int, &
"zt", id_zzl, id_zzi, diag_cs%nz_remap)
call log_param(param_file, mod, "!INPUTDIR/Z_OUTPUT_GRID_FILE", &
trim(in_dir)//trim(zgrid_file))
call log_param(param_file, mod, "!NK_ZSPACE (from file)", diag_cs%nz_remap, &
"The number of depth-space levels. This is determined \n"//&
"from the size of the variable zw in the output grid file.", &
units="nondim")
call get_param(param_file, mod, "DIAG_REMAP_Z_GRID_DEF", string, &
"This sets the file and variable names that define the\n"//&
"vertical grid used for diagnostic output remapping to\n"//&
"Z space. It should look like:\n"//&
" FILE:<file>,<variable> - where <file> is a file within\n"//&
" the INPUTDIR, <variable> is\n"//&
" the name of the variable that\n"//&
" contains interface positions.\n",&
default="")
if (len_trim(string) > 0) then
call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
filename = trim(inputdir) // trim(extractWord(trim(string(6:200)), 1))
varname = trim(extractWord(trim(string(6:200)), 2))
dimname = trim(extractWord(trim(string(6:200)), 3))

if (.not. file_exists(trim(filename))) then
call MOM_error(FATAL,"set_axes_info: Specified file not found: "//&
"Looking for '"//trim(filename)//"'")
endif
! Check that the grid has expected format, units etc.
if (.not. check_grid_def(trim(filename), trim(varname))) then
call MOM_error(FATAL,"set_axes_info: Bad grid definition in "//&
"'"//trim(filename)//"'")
endif
call log_param(param_file, mod, "DIAG_REMAP_Z_GRID_DEF", &
trim(inputdir)//"/"//trim(filename)//","//trim(varname))

! Get interface dimensions
call field_size(filename, varname, nzi)
call assert(nzi(1) /= 0, 'set_axes_info: bad z-axis dimension size')
allocate(diag_cs%zi_remap(nzi(1)))
allocate(diag_cs%zl_remap(nzi(1) - 1))
call MOM_read_data(filename, varname, diag_cs%zi_remap)
! Calculate layer positions
diag_cs%zl_remap(:) = diag_cs%zi_remap(1:nzi(1)-1) + &
(diag_cs%zi_remap(2:) - diag_cs%zi_remap(:nzi(1)-1)) / 2
diag_cs%nz_remap = nzi(1) - 1

! Make axes objects
id_zzl = diag_axis_init('zl_remap', diag_cs%zl_remap, "meters", "z", &
"Depth of cell center", direction=-1)
id_zzi = diag_axis_init('zi_remap', diag_cs%zi_remap, "meters", "z", &
'Depth of interfaces', direction=-1)
else
! In this case the axes associated with these will never be used, however
! they need to be positive otherwise FMS complains.
id_zzl = 1
id_zzi = 1
endif

! Axis for z remapping
call defineAxes(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL)

! Vertical axes for the interfaces and layers
call defineAxes(diag_cs, (/ id_zi /), diag_cs%axesZi)
call defineAxes(diag_cs, (/ id_zL /), diag_cs%axesZL)
Expand All @@ -284,11 +322,38 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical)
call defineAxes(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1)
call defineAxes(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1)

! Axis for z remapping
call defineAxes(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL)

end subroutine set_axes_info

function check_grid_def(filename, varname)
! Do some basic checks on the vertical grid definition file, variable
character(len=*), intent(in) :: filename
character(len=*), intent(in) :: varname
logical :: check_grid_def

character (len=200) :: units, long_name
integer :: ncid, status, intid, vid

check_grid_def = .true.
status = NF90_OPEN(filename, NF90_NOWRITE, ncid);
if (status /= NF90_NOERR) then
check_grid_def = .false.
endif

status = NF90_INQ_VARID(ncid, varname, vid)
if (status /= NF90_NOERR) then
check_grid_def = .false.
endif

status = NF90_GET_ATT(ncid, vid, "units", units)
if (status /= NF90_NOERR) then
check_grid_def = .false.
endif
if (trim(units) /= "meters" .and. trim(units) /= "m") then
check_grid_def = .false.
endif

end function check_grid_def

subroutine defineAxes(diag_cs, handles, axes)
! Defines "axes" from list of handle and associates mask
type(diag_ctrl), target, intent(in) :: diag_cs
Expand All @@ -307,131 +372,6 @@ subroutine defineAxes(diag_cs, handles, axes)

end subroutine defineAxes

subroutine get_Z_output_grid(depth_file, int_depth_name, int_depth, cell_depth_name, &
z_axis_index, z_int_axis_index, nz)
character(len=*), intent(in) :: depth_file
character(len=*), intent(in) :: int_depth_name
real, dimension(:), pointer :: int_depth
character(len=*), intent(in) :: cell_depth_name
integer, intent(out) :: z_axis_index
integer, intent(out) :: z_int_axis_index
integer, intent(out) :: nz

! This subroutine reads the depths of the interfaces bounding the intended
! layers from a NetCDF file. If no appropriate file is found, -1 is returned
! as the number of layers in the output file. Also, a diag_manager axis is set
! up with the same information as this axis.

real, allocatable :: cell_depth(:)
character (len=200) :: units, long_name
integer :: ncid, status, intid, intvid, layid, layvid, k, ni

nz = -1

status = NF90_OPEN(depth_file, NF90_NOWRITE, ncid);
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
" Difficulties opening "//trim(depth_file)//" - "//&
trim(NF90_STRERROR(status)))
nz = -1 ; return
endif

status = NF90_INQ_DIMID(ncid, int_depth_name, intid)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting ID of dimension "//&
trim(int_depth_name)//" in "//trim(depth_file))
endif

status = nf90_Inquire_Dimension(ncid, intid, len=ni)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting number of interfaces of "//&
trim(int_depth_name)//" in "//trim(depth_file))
endif

if (ni < 2) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
"At least two interface depths must be specified in "//trim(depth_file))
endif

status = NF90_INQ_DIMID(ncid, cell_depth_name, layid)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting ID of dimension "//&
trim(cell_depth_name)//" in "//trim(depth_file))
endif
status = nf90_Inquire_Dimension(ncid, layid, len=nz)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting number of interfaces of "//&
trim(cell_depth_name)//" in "//trim(depth_file))
endif
if (ni /= nz+1) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
"The interface depths must have one more point than cell centers in "//&
trim(depth_file))
endif

allocate(int_depth(nz+1))
allocate(cell_depth(nz))

status = NF90_INQ_VARID(ncid, int_depth_name, intvid)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting ID of variable "//&
trim(int_depth_name)//" in "//trim(depth_file))
endif
status = NF90_GET_VAR(ncid, intvid, int_depth)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Reading variable "//&
trim(int_depth_name)//" in "//trim(depth_file))
endif
status = NF90_GET_ATT(ncid, intvid, "units", units)
if (status /= NF90_NOERR) units = "meters"
status = NF90_GET_ATT(ncid, intvid, "long_name", long_name)
if (status /= NF90_NOERR) long_name = "Depth of edges"
z_int_axis_index = diag_axis_init(int_depth_name//'_new', int_depth, units, 'z', &
long_name, direction=-1)

! Create an fms axis with the same data as the cell_depth array in the
! input file.
status = NF90_INQ_VARID(ncid, cell_depth_name, layvid)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Getting ID of variable "//&
trim(cell_depth_name)//" in "//trim(depth_file))
endif
status = NF90_GET_VAR(ncid, layvid, cell_depth)
if (status /= NF90_NOERR) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
trim(NF90_STRERROR(status))//" Reading variable "//&
trim(cell_depth_name)//" in "//trim(depth_file))
endif
status = NF90_GET_ATT(ncid, layvid, "units", units)
if (status /= NF90_NOERR) units = "meters"
status = NF90_GET_ATT(ncid, layvid, "long_name", long_name)
if (status /= NF90_NOERR) long_name = "Depth of cell center"

z_axis_index = diag_axis_init(cell_depth_name//'_new', cell_depth, units, 'z',&
long_name, edges = z_int_axis_index, direction=-1)
deallocate(cell_depth)
status = nf90_close(ncid)

! Check the sign convention and change to the MOM "height" convention.
!if (int_depth(1) < int_depth(2)) then
! do k=1,nz+1 ; int_depth(k) = -1*int_depth(k) ; enddo
!endif

! Check for inversions in grid.
do k=1,nz ; if (int_depth(k) > int_depth(k+1)) then
call MOM_error(FATAL,"MOM_diag_to_Z get_Z_depths: "//&
"Inverted interface depths in output grid in "//depth_file)
endif ; enddo

end subroutine get_Z_output_grid

subroutine set_diag_mediator_grid(G, diag_cs)
type(ocean_grid_type), intent(inout) :: G
type(diag_ctrl), intent(inout) :: diag_cs
Expand Down Expand Up @@ -700,17 +640,17 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field)
type(regridding_CS) :: regrid_cs
integer :: nz_src, i, j, k
real, dimension(diag_cs%nz_remap) :: h_new, h_remap
real, dimension(diag_cs%nz_remap+1) :: dz, z_int_tmp
real, dimension(diag_cs%nz_remap+1) :: dz, zi_tmp

call assert(size(field, 3) == size(diag_cs%h, 3), &
'Remap field and thickness z-axes do not match.')
call assert(associated(diag_cs%z_remap_int), 'Remapping axis not initialized')
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

remapped_field = diag_cs%missing_value
nz_src = size(field, 3)

! Nominal thicknesses to remap to
h_remap(:) = diag_cs%z_remap_int(2:) - diag_cs%z_remap_int(:diag_cs%nz_remap)
h_remap(:) = diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap)
call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', regrid_cs)
call setRegriddingMinimumThickness(diag_cs%G%Angstrom, regrid_cs)
call setCoordinateResolution(h_remap, regrid_cs)
Expand All @@ -725,11 +665,11 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field)
! Calculate thicknesses for new Z* using nominal thicknesses from
! h_remap, current bathymetry and total thickness.
call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, diag_cs%G%bathyT(i, j), &
sum(diag_cs%h(i, j, :)), z_int_tmp)
sum(diag_cs%h(i, j, :)), zi_tmp)
! Calculate how much thicknesses change between source and dest grids, do
! remapping
z_int_tmp = -z_int_tmp
h_new(:) = z_int_tmp(2:) - z_int_tmp(:size(z_int_tmp)-1)
zi_tmp = -zi_tmp
h_new(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1)
call dzFromH1H2(nz_src, diag_cs%h(i, j, :), diag_cs%nz_remap, h_new, dz)
call remapping_core(remap_cs, nz_src, diag_cs%h(i, j, :), field(i,j,:), &
diag_cs%nz_remap, dz, remapped_field(i, j, :))
Expand All @@ -738,7 +678,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field)
! their depth does not corrospond to the nominal depth at that level
! (and the nominal depth is below the bathymetry), so remove
do k=1, diag_cs%nz_remap
if (diag_cs%z_remap_int(k) >= diag_cs%G%bathyT(i, j)) then
if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then
remapped_field(i, j, k:diag_cs%nz_remap) = diag_cs%missing_value
exit
endif
Expand Down Expand Up @@ -1017,9 +957,9 @@ function register_diag_field(module_name, field_name, axes, init_time, &
interp_method=interp_method, tile_count=tile_count)
if (fms_id /= DIAG_FIELD_NOT_FOUND) then
! Check that we have read in necessary remapping information from file
if (diag_cs%nz_remap == -1) then
if (.not. allocated(diag_cs%zi_remap)) then
call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// &
'destination grid spec provided, see param Z_OUTPUT_GRID_FILE')
'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF')
endif

if (primary_id == -1) then
Expand Down

0 comments on commit d556620

Please sign in to comment.