From 0dda96cf4f5d8ac9a2844c06fbbcd4fcf5b64f4d Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:37:32 -0700 Subject: [PATCH 01/12] Split Zstar grid creation up so can be done column-by-column. #62 --- src/ALE/MOM_regridding.F90 | 68 ++++++++++++++++++++++++-------------- 1 file changed, 44 insertions(+), 24 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index c9f14c6c9e..391036447d 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -98,6 +98,7 @@ module MOM_regridding public getCoordinateUnits public getCoordinateShortName public getStaticThickness +public buildGridZstarColumn public DEFAULT_COORDINATE_MODE character(len=158), parameter, public :: regriddingCoordinateModeDoc = & @@ -143,7 +144,7 @@ module MOM_regridding ! Deviation tolerance between succesive grids in regridding iterations real, parameter :: DEVIATION_TOLERANCE = 1e-10 ! Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are -! used to build the new grid by finding the coordinates associated with +! used to build the new grid by finding the coordinates associated with ! target densities and interpolations of degree larger than 1. integer, parameter :: NR_ITERATIONS = 8 ! Tolerance for Newton-Raphson iterations (stop when increment falls below this) @@ -355,9 +356,8 @@ subroutine checkGridsMatch( G, h, dzInterface ) 'Non-zero dzInterface at bottom!') enddo enddo - -end subroutine checkGridsMatch +end subroutine checkGridsMatch !------------------------------------------------------------------------------ ! Build uniform z*-ccordinate grid with partial steps @@ -405,31 +405,11 @@ subroutine buildGridZstar( CS, G, h, dzInterface ) do k = 1,nz totalThickness = totalThickness + h(i,j,k) end do - minThickness = min( CS%min_thickness, totalThickness/float(nz) ) - ! Position of free-surface - eta = totalThickness - nominalDepth - - ! z* = (z-eta) / stretching where stretching = (H+eta)/H - ! z = eta + stretching * z* - stretching = totalThickness / nominalDepth - - ! Integrate down from the top for a notional new grid, ignoring topography - zNew(1) = eta - do k = 1,nz - dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing - zNew(k+1) = zNew(k) - dh - enddo + call buildGridZStarColumn(CS, nz, nominalDepth, totalThickness, zNew) - ! The rest of the model defines grids integrating up from the bottom zOld(nz+1) = - nominalDepth - zNew(nz+1) = - nominalDepth do k = nz,1,-1 - ! Adjust interface position to accomodate inflating layers - ! without disturbing the interface above - if ( zNew(k) < (zNew(k+1) + minThickness) ) then - zNew(k) = zNew(k+1) + minThickness - endif zOld(k) = zOld(k+1) + h(i,j,k) enddo @@ -463,6 +443,46 @@ subroutine buildGridZstar( CS, G, h, dzInterface ) end subroutine buildGridZstar +subroutine buildGridZstarColumn( CS, nz, depth, totalThickness, zInterface) + + ! Arguments + type(regridding_CS), intent(in) :: CS + integer, intent(in) :: nz + real, intent(in) :: depth + real, intent(in) :: totalThickness + real, dimension(nz+1), intent(inout) :: zInterface + + real :: eta, stretching, dh + real :: minThickness + integer :: k + + minThickness = min( CS%min_thickness, totalThickness/float(nz) ) + + ! Position of free-surface + eta = totalThickness - depth + + ! z* = (z-eta) / stretching where stretching = (H+eta)/H + ! z = eta + stretching * z* + stretching = totalThickness / depth + + ! Integrate down from the top for a notional new grid, ignoring topography + zInterface(1) = eta + do k = 1,nz + dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing + zInterface(k+1) = zInterface(k) - dh + enddo + + ! Integrating up from the bottom adjusting interface position to accomodate + ! inflating layers without disturbing the interface above + zInterface(nz+1) = -depth + do k = nz,1,-1 + if ( zInterface(k) < (zInterface(k+1) + minThickness) ) then + zInterface(k) = zInterface(k+1) + minThickness + endif + enddo + +end subroutine buildGridZstarColumn + !------------------------------------------------------------------------------ ! Build sigma grid From 691109751d794774b7d88aac3151af401bc932ee Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:38:47 -0700 Subject: [PATCH 02/12] Fix output from remapping conservation check. #62 --- src/ALE/MOM_remapping.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 4385cae22c..487d12ce8c 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -400,7 +400,7 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) if (k <= n0) then hTmp = h0(k) + ( dx(k+1) - dx(k) ) if (hTmp < 0.) then - write(0,*) 'k,h0(k),hNew,dx(+1),dx(0)=',k,h0(k),dx(k+1),dx(k) + write(0,*) 'k,h0(k),hTmp,dx(k+1),dx(k)=',k,h0(k),hTmp,dx(k+1),dx(k) call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& 'negative h implied by fluxes' ) endif @@ -534,9 +534,9 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) if (k<=n0) then; hTmp = h0(k); else; hTmp = 0.; endif z0 = z0 + hTmp ; z1 = z1 + ( hTmp + ( dx(k+1) - dx(k) ) ) enddo - if (abs(totalHU2-totalHU0) > (err0+err2)*real(n1) .and. (err0+err2)/=0.) then + if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then write(0,*) 'h0=',h0 - write(0,*) 'hf=',h0+dx(2:n1+1)-dx(1:n1) + write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1) write(0,*) 'u0=',u0 write(0,*) 'u1=',u1 write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 From b10da826eb39c75ed03db0a8673111710e548b46 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:40:11 -0700 Subject: [PATCH 03/12] Only initialise diag_mediator after h and G have been set up. #62 --- src/core/MOM.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c78182c60d..899d83c58b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1638,8 +1638,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - call diag_mediator_init(G, param_file, diag) - ! Read relevant parameters and write them to the model log. call log_version(param_file, "MOM", version, "") call get_param(param_file, "MOM", "VERBOSITY", verbosity, & @@ -1920,6 +1918,10 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif call callTree_waypoint("state variables allocated (initialize_MOM)") + ! Initialise the diag mediator, this needs to occur _after_ CS%h has been + ! allocated. + call diag_mediator_init(G, CS%h, param_file, diag) + ! Set the fields that are needed for bitwise identical restarting ! the time stepping scheme. call restart_init(param_file, CS%restart_CSp) @@ -2379,7 +2381,7 @@ subroutine register_diags(Time, G, CS, ADp) endif if (CS%use_temperature) then CS%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, Time, & - 'Potential Temperature', 'Celsius') + long_name='Potential Temperature', units='Celsius') CS%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, Time, & 'Salinity', 'PPT') endif From 51c9ec6816b3341987ce831a7f019dc07cd4b093 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:42:03 -0700 Subject: [PATCH 04/12] Disable old Z space diagnostic remapping for tracers. Will eventually be removed. #62 --- src/diagnostics/MOM_diag_to_Z.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 511a21f61a..c255cb22ae 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -862,11 +862,11 @@ subroutine register_Z_tracer_low(tr_ptr, name, long_name, units, standard_name, CS%missing_tr(m) = CS%missing_value ! This could be changed later, if desired. if (CS%nk_zspace > 0) then - CS%id_tr(m) = register_diag_field('ocean_model_z', name, CS%axesTz, Time, & + CS%id_tr(m) = register_diag_field('ocean_model_old_z', name, CS%axesTz, Time, & long_name, units, missing_value=CS%missing_tr(m), & standard_name=standard_name) else - id_test = register_diag_field('ocean_model_z', name, CS%diag%axesT1, Time, & + id_test = register_diag_field('ocean_model_old_z', name, CS%diag%axesT1, Time, & long_name, units, missing_value=CS%missing_tr(m), & standard_name=standard_name) if (id_test>0) call MOM_error(WARNING, & @@ -1092,7 +1092,7 @@ subroutine get_Z_depths(depth_file, int_depth_name, int_depth, cell_depth_name, 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" - edge_index = diag_axis_init(int_depth_name, int_depth, units, 'z', & + edge_index = diag_axis_init(int_depth_name//'_old', 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. @@ -1109,7 +1109,7 @@ subroutine get_Z_depths(depth_file, int_depth_name, int_depth, cell_depth_name, 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, cell_depth, units, 'z',& + z_axis_index = diag_axis_init(cell_depth_name//'_old', cell_depth, units, 'z',& long_name, edges = edge_index, direction=-1) deallocate(cell_depth) From f320f0a8c21c3fa4dd20bf4bb781dff433578300 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:43:31 -0700 Subject: [PATCH 05/12] Create diag variation that remaps T+layer diagnostics automatically. #62 --- src/framework/MOM_diag_mediator.F90 | 347 ++++++++++++++++++++++++++-- 1 file changed, 322 insertions(+), 25 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a6d3ac850b..e0e8e89cbf 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -31,19 +31,33 @@ 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 : vardesc +use MOM_io, only : slasher, vardesc use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type +use MOM_remapping, only : remapping_CS, remapping_core, initialize_remapping, dzFromH1H2 +use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution, buildGridZStarColumn use diag_manager_mod, only : diag_manager_init, diag_manager_end use diag_manager_mod, only : send_data, diag_axis_init use diag_manager_mod, only : register_diag_field_fms=>register_diag_field use diag_manager_mod, only : register_static_field_fms=>register_static_field +use diag_manager_mod, only : get_diag_field_id_fms=>get_diag_field_id use diag_manager_mod, only : DIAG_FIELD_NOT_FOUND +use netcdf + implicit none ; private +#include + +#define RANGE_I(a) lbound(a, 1),ubound(a, 1) +#define RANGE_J(a) lbound(a, 2),ubound(a, 2) +#define RANGE_K(a) lbound(a, 3),ubound(a, 3) +#define DIM_I(a) lbound(a, 1):ubound(a, 1) +#define DIM_J(a) lbound(a, 2):ubound(a, 2) +#define DIM_K(a) lbound(a, 3):ubound(a, 3) + public set_axes_info, post_data, register_diag_field, time_type public post_data_1d_k public safe_alloc_ptr, safe_alloc_alloc @@ -75,6 +89,7 @@ module MOM_diag_mediator type, private :: diag_type logical :: in_use integer :: fms_diag_id ! underlying FMS diag id + type(axesType), pointer :: remap_axes => null() real, pointer, dimension(:,:) :: mask2d => null() real, pointer, dimension(:,:,:) :: mask3d => null() type(diag_type), pointer :: next => null() ! pointer to the next diag @@ -98,6 +113,7 @@ module MOM_diag_mediator type(axesType) :: axesBi, axesTi, axesCui, axesCvi type(axesType) :: axesB1, axesT1, axesCu1, axesCv1 type(axesType) :: axesZi, axesZL + type(axesType) :: axesTzi, axesTZL ! Mask arrays for diagnostics real, dimension(:,:), pointer :: mask2dT => null() @@ -122,6 +138,18 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations 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 + ! Number of z levels used for remapping + integer :: nz_remap + + ! Pointer to H and G for Z remapping + real, dimension(:,:,:), pointer :: h => null() + type(ocean_grid_type), pointer :: G => null() + end type diag_ctrl integer :: doc_unit = -1 @@ -141,17 +169,17 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) ! (inout) diag_cs - structure used to regulate diagnostic output. ! (in,opt) set_vertical - If true (or missing), set up the vertical axes. - integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, k, nz + integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_zzl, id_zzi + integer :: k, nz 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 ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_diag_mediator" ! This module's name. - nz = G%ke - set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical ! Read all relevant parameters and write them to the model log. @@ -185,12 +213,6 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) else Cartesian_grid = .false. endif - - zInter(1:nz+1) = G%GV%sInterface(1:nz+1) - zLev(1:nz) = G%GV%sLayer(1:nz) - -! do i=1,nz ; zlev(i) = real(i) ; enddo -! do i=1,nz+1 ; zinter(i) = real(i) - 0.5 ; enddo if(G%symmetric) then id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', & @@ -202,13 +224,16 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) 'q point nominal longitude', Domain2=G%Domain%mpp_domain) id_yq = diag_axis_init('yq', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', & 'q point nominal latitude', Domain2=G%Domain%mpp_domain) - endif + endif id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', & 'h point nominal longitude', Domain2=G%Domain%mpp_domain) id_yh = diag_axis_init('yh', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', & 'h point nominal latitude', Domain2=G%Domain%mpp_domain) if (set_vert) then + nz = G%ke + zinter(1:nz+1) = G%GV%sInterface(1:nz+1) + zlev(1:nz) = G%GV%sLayer(1:nz) id_zl = diag_axis_init('zl', zlev, trim(G%GV%zAxisUnits), 'z', & 'Layer '//trim(G%GV%zAxisLongName), & direction=G%GV%direction) @@ -219,6 +244,25 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) id_zl = -1 ; id_zi = -1 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") + endif + ! 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) @@ -240,7 +284,10 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) call defineAxes(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1) 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 subroutine defineAxes(diag_cs, handles, axes) @@ -261,6 +308,131 @@ 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, 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, 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 @@ -362,7 +534,6 @@ subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) ! (in,opt) is_static - If true, this is a static field that is always offered. ! (in,opt) mask - If present, use this real array as the data mask. - integer :: altId type(diag_type), pointer :: diag => null() ! Iterate over list of diag 'variants' (e.g. CMOR aliases) and post each. @@ -370,7 +541,7 @@ subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) 'post_data_2d: Unregistered diagnostic id') diag => diag_cs%diags(diag_field_id) do while (associated(diag)) - call post_data_2d_low(diag, field, diag_cs, is_static, mask) + call post_data_2d_low(diag, field, diag_cs, is_static, mask) diag => diag%next enddo @@ -463,6 +634,7 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) end subroutine post_data_2d_low subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) + integer, intent(in) :: diag_field_id real, intent(in) :: field(:,:,:) type(diag_ctrl), target, intent(in) :: diag_cs @@ -477,20 +649,105 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) ! (in) static - If true, this is a static field that is always offered. ! (in,opt) mask - If present, use this real array as the data mask. - integer :: altId type(diag_type), pointer :: diag => null() + real, allocatable :: remapped_field(:,:,:) - ! Iterate over list of diag 'variants', e.g. CMOR aliases. + ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical + ! grids, and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & 'post_data_3d: Unregistered diagnostic id') diag => diag_cs%diags(diag_field_id) do while (associated(diag)) - call post_data_3d_low(diag, field, diag_cs, is_static, mask) + + if (associated(diag%remap_axes)) then + ! Remap this field to another vertical coordinate. + + if (present(mask)) then + call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") + endif + + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call remap_diag_to_z(field, diag, diag_cs, remapped_field) + if (associated(diag%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%mask3d(:,:,:diag_cs%nz_remap)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + deallocate(remapped_field) + else + call post_data_3d_low(diag, field, diag_cs, is_static, mask) + endif diag => diag%next enddo end subroutine post_data_3d +subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) + + real, dimension(:,:,:), intent(in) :: field + type(diag_type), intent(in) :: diag + type(diag_ctrl), intent(in) :: diag_cs + real, dimension(:,:,:), intent(inout) :: remapped_field + +! Arguments: +! (in) field - The diagnostic field to be remapped +! (in) diag - structure used to regulate diagnostic output +! (out) remapped_field - the field argument remapped to z star coordinate + + type(remapping_CS) :: remap_cs + 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 + + 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') + + 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) + call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', regrid_cs) + call setCoordinateResolution(h_remap, regrid_cs) + call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) + do j=RANGE_J(field) + do i=RANGE_I(field) + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j, 1) == 0.) cycle + endif + + ! 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) + + ! 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) + 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, :)) + + ! Lower levels of the remapped data get squashed to follow bathymetry, + ! 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 + remapped_field(i, j, k:diag_cs%nz_remap) = diag_cs%missing_value + exit + endif + enddo + enddo + enddo + +end subroutine remap_diag_to_z + subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, intent(in) :: field(:,:,:) @@ -647,7 +904,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & cmor_long_name, cmor_units, cmor_standard_name) integer :: register_diag_field character(len=*), intent(in) :: module_name, field_name - type(axesType), intent(in) :: axes + type(axesType), target, intent(in) :: axes type(time_type), intent(in) :: init_time character(len=*), optional, intent(in) :: long_name, units, standard_name real, optional, intent(in) :: missing_value, range(2) @@ -685,8 +942,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs - type(diag_type), pointer :: diag => null(), cmor_diag => null() - integer :: primary_id, fms_id + type(diag_type), pointer :: diag => null(), cmor_diag => null(), z_remap_diag => null() + integer :: primary_id, fms_id, remap_id character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name MOM_missing_value = axes%diag_cs%missing_value @@ -696,6 +953,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & primary_id = -1 diag => null() cmor_diag => null() + z_remap_diag => null() ! Set up the 'primary' diagnostic, first get an underlying FMS id fms_id = register_diag_field_fms(module_name, field_name, axes%handles, & @@ -748,8 +1006,31 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif endif - ! Set up any other variations here, e.g. remapped vertical dimension, - ! or decimation. + ! Remap to z vertical coordinate, note that only + ! diagnostics on T points and layers (not interfaces) are supported, + ! for other diagnostics the old remapping approach can still be used + if (axes%id == diag_cs%axesTL%id) then + fms_id = register_diag_field_fms(module_name//'_z', field_name, diag_cs%axesTZL%handles, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + 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 + call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & + 'destination grid spec provided, see param Z_OUTPUT_GRID_FILE') + endif + + if (primary_id == -1) then + primary_id = get_new_diag_id(diag_cs) + endif + call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) + z_remap_diag%fms_diag_id = fms_id + z_remap_diag%remap_axes => diag_cs%axesTZL + call set_diag_mask(z_remap_diag, diag_cs, diag_cs%axesTZL) + endif + endif ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. doc_unit > 0) then @@ -760,6 +1041,10 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif + if (axes%id == diag_cs%axesTL%id) then + call log_available_diag(associated(z_remap_diag), module_name//'_z', field_name, & + long_name, units, standard_name) + endif endif register_diag_field = primary_id @@ -1072,8 +1357,9 @@ function ocean_register_diag(var_desc, G, diag_cs, day) end function ocean_register_diag -subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) - type(ocean_grid_type), intent(inout) :: G +subroutine diag_mediator_init(G, h, param_file, diag_cs, err_msg) + type(ocean_grid_type), target, intent(inout) :: G + real, dimension(NIMEM_,NJMEM_,NKMEM_), target, intent(in) :: h type(param_file_type), intent(in) :: param_file type(diag_ctrl), intent(inout) :: diag_cs character(len=*), optional, intent(out) :: err_msg @@ -1095,8 +1381,16 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) do i=1, DIAG_ALLOC_CHUNK_SIZE diag_cs%diags(i)%in_use = .false. diag_cs%diags(i)%next => null() + diag_cs%diags(i)%remap_axes => null() + diag_cs%diags(i)%mask2d => null() + diag_cs%diags(i)%mask3d => null() enddo + ! Keep pointers to these, needed for the z axis remapping + diag_cs%h => h + diag_cs%G => G + diag_cs%nz_remap = -1 + diag_cs%is = G%isc - (G%isd-1) ; diag_cs%ie = G%iec - (G%isd-1) diag_cs%js = G%jsc - (G%jsd-1) ; diag_cs%je = G%jec - (G%jsd-1) diag_cs%isd = G%isd ; diag_cs%ied = G%ied @@ -1225,7 +1519,8 @@ subroutine set_diag_mask(diag, diag_cs, axes) if (axes%rank .eq. 3) then diag%mask2d => null() diag%mask3d => null() - if (axes%id .eq. diag_cs%axesTL%id) then + if ((axes%id .eq. diag_cs%axesTL%id) .or. & + (axes%id .eq. diag_cs%axesTZL%id)) then diag%mask3d => diag_cs%mask3dTL elseif(axes%id .eq. diag_cs%axesBL%id) then diag%mask3d => diag_cs%mask3dBuL @@ -1289,6 +1584,8 @@ function get_new_diag_id(diag_cs) do i=diag_cs%next_free_diag_id, size(diag_cs%diags) diag_cs%diags(i)%in_use = .false. diag_cs%diags(i)%next => null() + diag_cs%diags(i)%mask2d => null() + diag_cs%diags(i)%mask3d => null() enddo endif From d38f43404fb13ef15f8a7e675cfcb0a72c321edc Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Tue, 2 Jun 2015 13:27:56 -0700 Subject: [PATCH 06/12] Reduce the minimum regridding thickness. Necessary for conservation when there are many thin layers. #62 --- src/ALE/MOM_regridding.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 391036447d..518a032522 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -119,7 +119,7 @@ module MOM_regridding " PQM_IH6IH5 (5th-order accurate)" character(len=6), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2" logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false. -real, parameter, public :: regriddingDefaultMinThickness = 1.e-3 +real, parameter, public :: regriddingDefaultMinThickness = 1.e-10 ! ----------------------------------------------------------------------------- ! The following are private constants From 9ecfe58e0817e4e8f983cf53f78068e2f1e1daaf Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 4 Jun 2015 11:01:05 -0700 Subject: [PATCH 07/12] Set the minimum layer thickness for remapping only. #62 --- src/ALE/MOM_regridding.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 518a032522..391036447d 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -119,7 +119,7 @@ module MOM_regridding " PQM_IH6IH5 (5th-order accurate)" character(len=6), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2" logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false. -real, parameter, public :: regriddingDefaultMinThickness = 1.e-10 +real, parameter, public :: regriddingDefaultMinThickness = 1.e-3 ! ----------------------------------------------------------------------------- ! The following are private constants diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index e0e8e89cbf..3bd1592a5e 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -36,7 +36,7 @@ module MOM_diag_mediator use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_remapping, only : remapping_CS, remapping_core, initialize_remapping, dzFromH1H2 -use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution, buildGridZStarColumn +use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution, buildGridZStarColumn, setRegriddingMinimumThickness use diag_manager_mod, only : diag_manager_init, diag_manager_end use diag_manager_mod, only : send_data, diag_axis_init @@ -713,8 +713,10 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! Nominal thicknesses to remap to h_remap(:) = diag_cs%z_remap_int(2:) - diag_cs%z_remap_int(: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) call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) + do j=RANGE_J(field) do i=RANGE_I(field) if (associated(diag%mask3d)) then @@ -725,7 +727,6 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! 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) - ! Calculate how much thicknesses change between source and dest grids, do ! remapping z_int_tmp = -z_int_tmp From 9ab9b684a0c45cd336a04ed707fb395275d0c7df Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 5 Jun 2015 07:17:41 -0700 Subject: [PATCH 08/12] Don't move call to diag_mediator_init, setup pointer to H with setter method. #62 --- src/core/MOM.F90 | 10 ++++++---- src/framework/MOM_diag_mediator.F90 | 22 +++++++++++++++++----- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 899d83c58b..3fc1616c88 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -337,7 +337,7 @@ module MOM use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_coms, only : reproducing_sum -use MOM_diag_mediator, only : diag_mediator_init, enable_averaging +use MOM_diag_mediator, only : diag_mediator_init, enable_averaging, diag_set_thickness_ptr use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : register_scalar_field @@ -1638,6 +1638,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + call diag_mediator_init(G, param_file, diag) + ! Read relevant parameters and write them to the model log. call log_version(param_file, "MOM", version, "") call get_param(param_file, "MOM", "VERBOSITY", verbosity, & @@ -1918,9 +1920,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif call callTree_waypoint("state variables allocated (initialize_MOM)") - ! Initialise the diag mediator, this needs to occur _after_ CS%h has been - ! allocated. - call diag_mediator_init(G, CS%h, param_file, diag) + ! Set up a pointers h within diag mediator control structure, + ! this needs to occur _after_ CS%h has been allocated. + call diag_set_thickness_ptr(CS%h, diag) ! Set the fields that are needed for bitwise identical restarting ! the time stepping scheme. diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 3bd1592a5e..157312162a 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -67,6 +67,7 @@ module MOM_diag_mediator public diag_axis_init, ocean_register_diag, register_static_field public register_scalar_field public defineAxes, diag_masks_set +public diag_set_thickness_ptr interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -715,7 +716,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) 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) - call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) + call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) do j=RANGE_J(field) do i=RANGE_I(field) @@ -1358,9 +1359,8 @@ function ocean_register_diag(var_desc, G, diag_cs, day) end function ocean_register_diag -subroutine diag_mediator_init(G, h, param_file, diag_cs, err_msg) +subroutine diag_mediator_init(G, param_file, diag_cs, err_msg) type(ocean_grid_type), target, intent(inout) :: G - real, dimension(NIMEM_,NJMEM_,NKMEM_), target, intent(in) :: h type(param_file_type), intent(in) :: param_file type(diag_ctrl), intent(inout) :: diag_cs character(len=*), optional, intent(out) :: err_msg @@ -1387,8 +1387,7 @@ subroutine diag_mediator_init(G, h, param_file, diag_cs, err_msg) diag_cs%diags(i)%mask3d => null() enddo - ! Keep pointers to these, needed for the z axis remapping - diag_cs%h => h + ! Keep a pointer to the grid, this is needed for regridding diag_cs%G => G diag_cs%nz_remap = -1 @@ -1431,6 +1430,19 @@ subroutine diag_mediator_init(G, h, param_file, diag_cs, err_msg) end subroutine diag_mediator_init +subroutine diag_set_thickness_ptr(h, diag_cs) + + real, dimension(NIMEM_,NJMEM_,NKMEM_), target, intent(in) :: h + type(diag_ctrl), intent(inout) :: diag_cs + + ! (inout) diag_cs - diag mediator control structure + ! (in) h - a pointer to model thickness + + ! Keep pointer to h, needed for the z axis remapping + diag_cs%h => h + +end subroutine + subroutine diag_masks_set(G, missing_value, diag_cs) ! Setup the 2d masks for diagnostics type(ocean_grid_type), target, intent(in) :: G From 423b0ee2b507a054b786383ca127743ac2f7866b Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 5 Jun 2015 07:20:35 -0700 Subject: [PATCH 09/12] Remove unnecessary addition of optional arguments. #62 --- src/core/MOM.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3fc1616c88..4f9d0a974a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2383,7 +2383,7 @@ subroutine register_diags(Time, G, CS, ADp) endif if (CS%use_temperature) then CS%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, Time, & - long_name='Potential Temperature', units='Celsius') + 'Potential Temperature', 'Celsius') CS%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, Time, & 'Salinity', 'PPT') endif From d4e25ec770903f83eee46f0cd5b3fba6dbee881b Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 5 Jun 2015 07:51:31 -0700 Subject: [PATCH 10/12] Change name of ALE remapped Z level diagnostics to have _z_new suffix, avoids conflict with existing Z level diagnostics. #62 --- src/diagnostics/MOM_diag_to_Z.F90 | 8 ++++---- src/framework/MOM_diag_mediator.F90 | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index c255cb22ae..511a21f61a 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -862,11 +862,11 @@ subroutine register_Z_tracer_low(tr_ptr, name, long_name, units, standard_name, CS%missing_tr(m) = CS%missing_value ! This could be changed later, if desired. if (CS%nk_zspace > 0) then - CS%id_tr(m) = register_diag_field('ocean_model_old_z', name, CS%axesTz, Time, & + CS%id_tr(m) = register_diag_field('ocean_model_z', name, CS%axesTz, Time, & long_name, units, missing_value=CS%missing_tr(m), & standard_name=standard_name) else - id_test = register_diag_field('ocean_model_old_z', name, CS%diag%axesT1, Time, & + id_test = register_diag_field('ocean_model_z', name, CS%diag%axesT1, Time, & long_name, units, missing_value=CS%missing_tr(m), & standard_name=standard_name) if (id_test>0) call MOM_error(WARNING, & @@ -1092,7 +1092,7 @@ subroutine get_Z_depths(depth_file, int_depth_name, int_depth, cell_depth_name, 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" - edge_index = diag_axis_init(int_depth_name//'_old', int_depth, units, 'z', & + edge_index = diag_axis_init(int_depth_name, 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. @@ -1109,7 +1109,7 @@ subroutine get_Z_depths(depth_file, int_depth_name, int_depth, cell_depth_name, 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//'_old', cell_depth, units, 'z',& + z_axis_index = diag_axis_init(cell_depth_name, cell_depth, units, 'z',& long_name, edges = edge_index, direction=-1) deallocate(cell_depth) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 157312162a..94e06b80cc 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -394,7 +394,7 @@ subroutine get_Z_output_grid(depth_file, int_depth_name, int_depth, cell_depth_n 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, int_depth, units, 'z', & + 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 @@ -416,7 +416,7 @@ subroutine get_Z_output_grid(depth_file, int_depth_name, int_depth, cell_depth_n 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, cell_depth, units, 'z',& + 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) @@ -1012,7 +1012,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! diagnostics on T points and layers (not interfaces) are supported, ! for other diagnostics the old remapping approach can still be used if (axes%id == diag_cs%axesTL%id) then - fms_id = register_diag_field_fms(module_name//'_z', field_name, diag_cs%axesTZL%handles, & + fms_id = register_diag_field_fms(module_name//'_z_new', field_name, diag_cs%axesTZL%handles, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & @@ -1044,7 +1044,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_standard_name) endif if (axes%id == diag_cs%axesTL%id) then - call log_available_diag(associated(z_remap_diag), module_name//'_z', field_name, & + call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & long_name, units, standard_name) endif endif From 70e568647b150bb1840313948de79293c98b6d14 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 5 Jun 2015 07:53:34 -0700 Subject: [PATCH 11/12] Remove unnecessary include. #62 --- src/framework/MOM_diag_mediator.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 94e06b80cc..5f6fe8c66b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -49,8 +49,6 @@ module MOM_diag_mediator implicit none ; private -#include - #define RANGE_I(a) lbound(a, 1),ubound(a, 1) #define RANGE_J(a) lbound(a, 2),ubound(a, 2) #define RANGE_K(a) lbound(a, 3),ubound(a, 3) @@ -1432,7 +1430,7 @@ end subroutine diag_mediator_init subroutine diag_set_thickness_ptr(h, diag_cs) - real, dimension(NIMEM_,NJMEM_,NKMEM_), target, intent(in) :: h + real, dimension(:,:,:), target, intent(in) :: h type(diag_ctrl), intent(inout) :: diag_cs ! (inout) diag_cs - diag mediator control structure From d556620c0916b159cff0a374be33a1c4c110b81a Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Fri, 5 Jun 2015 13:26:33 -0700 Subject: [PATCH 12/12] Improve specification of grid definition file,variable used for diagnostic remapping. The format is now: DIAG_REMAP_Z_GRID_DEF = "FILE:OM3_zgrid.nc,zw" Also use library to read netCDF files. #62 --- src/framework/MOM_diag_mediator.F90 | 258 +++++++++++----------------- 1 file changed, 99 insertions(+), 159 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 5f6fe8c66b..a6022ab3de 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -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 @@ -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 @@ -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" @@ -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:, - where is a file within\n"//& + " the INPUTDIR, 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) @@ -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 @@ -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 @@ -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) @@ -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, :)) @@ -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 @@ -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