From e401fdfcf0b24785fb47d3f85442eb547bc57225 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 7 Jul 2016 11:15:58 +1000 Subject: [PATCH 01/16] Diagnostic regridding for zstar, sigma and rho vertical coordinates. #334 --- src/ALE/MOM_ALE.F90 | 2 +- src/ALE/MOM_regridding.F90 | 9 +- src/ALE/regrid_consts.F90 | 41 +- src/core/MOM.F90 | 21 +- src/core/MOM_dynamics_legacy_split.F90 | 6 +- src/core/MOM_dynamics_split_RK2.F90 | 4 +- src/core/MOM_dynamics_unsplit.F90 | 4 +- src/framework/MOM_diag_mediator.F90 | 935 +++++------------- src/framework/MOM_diag_remap.F90 | 567 +++++++++++ src/framework/MOM_error_handler.F90 | 12 + .../lateral/MOM_mixed_layer_restrat.F90 | 6 +- .../lateral/MOM_thickness_diffuse.F90 | 4 +- .../vertical/MOM_bulk_mixed_layer.F90 | 4 +- .../vertical/MOM_diabatic_aux.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 6 +- 15 files changed, 895 insertions(+), 728 deletions(-) create mode 100644 src/framework/MOM_diag_remap.F90 diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a97fdd7dce..57434dcc1a 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -952,7 +952,7 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz ) call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, & compressibility_fraction=compress_fraction ) - if (coordMode(1:2) == 'Z*') then + if (coordMode(1:5) == 'ZSTAR') then call get_param(param_file, mod, "ZSTAR_RIGID_SURFACE_THRESHOLD", height_of_rigid_surface, & "A threshold height used to detect the presence of a rigid-surface\n"//& "depressing the upper-surface of the model, such as an ice-shelf.\n"//& diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 0e9fca8edf..604244859d 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -863,9 +863,8 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS ) ! Local depth (G%bathyT is positive) nominalDepth = G%bathyT(i,j)*GV%m_to_H - call build_rho_column(CS, remapCS, tv%eqn_of_state, nz, nominalDepth, & - h(i, j, :)*GV%H_to_m, & - tv%T(i, j, :), tv%S(i, j, :), zNew) + call build_rho_column(CS, remapCS, nz, nominalDepth, h(i, j, :)*GV%H_to_m, & + tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew) if (CS%integrate_downward_for_e) then zOld(1) = 0. @@ -926,7 +925,7 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS ) end subroutine build_rho_grid -subroutine build_rho_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, zInterface) +subroutine build_rho_column(CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface) ! The algorithn operates as follows within each ! column: ! 1. Given T & S within each layer, the layer densities are computed. @@ -941,11 +940,11 @@ subroutine build_rho_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, zInte type(regridding_CS), intent(in) :: CS !< Regridding control structure type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options - type(EOS_type), pointer :: eqn_of_state !< Equation of state structure integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom (positive in m) real, dimension(nz), intent(in) :: h !< Layer thicknesses, in m real, dimension(nz), intent(in) :: T, S !< T and S for column + type(EOS_type), pointer :: eqn_of_state !< Equation of state structure real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces ! Local variables diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 0042efc223..415dadaf41 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -8,22 +8,36 @@ module regrid_consts implicit none ; public -! List of regridding types -integer, parameter :: REGRIDDING_LAYER = -1 !< Layer mode -integer, parameter :: REGRIDDING_ZSTAR = 0 !< z* coordinates -integer, parameter :: REGRIDDING_RHO = 1 !< Target interface densities -integer, parameter :: REGRIDDING_SIGMA = 2 !< Sigma coordinates -integer, parameter :: REGRIDDING_ARBITRARY = 3 !< Arbitrary coordinates -integer, parameter :: REGRIDDING_HYCOM1 = 4 !< Simple HyCOM coordinates without BBL -integer, parameter :: REGRIDDING_SLIGHT = 5 !< Stretched coordinates in the +integer, parameter :: REGRIDDING_NUM_TYPES = 7 + +! List of regridding types. These should be consecutive and starting at 1. +! This allows them to be used as array indices. +integer, parameter :: REGRIDDING_LAYER = 1 !< Layer mode +integer, parameter :: REGRIDDING_ZSTAR = 2 !< z* coordinates +integer, parameter :: REGRIDDING_RHO = 3 !< Target interface densities +integer, parameter :: REGRIDDING_SIGMA = 4 !< Sigma coordinates +integer, parameter :: REGRIDDING_ARBITRARY = 5 !< Arbitrary coordinates +integer, parameter :: REGRIDDING_HYCOM1 = 6 !< Simple HyCOM coordinates without BBL +integer, parameter :: REGRIDDING_SLIGHT = REGRIDDING_NUM_TYPES !< Stretched coordinates in the !! lightest water, isopycnal below -character(len=5), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string -character(len=2), parameter :: REGRIDDING_ZSTAR_STRING = "Z*" !< z* string -character(len=3), parameter :: REGRIDDING_RHO_STRING = "RHO" !< Rho string -character(len=5), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string +character(len=6), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string +character(len=6), parameter :: REGRIDDING_ZSTAR_STRING = "ZSTAR" !< z* string +character(len=6), parameter :: REGRIDDING_RHO_STRING = "RHO" !< Rho string +character(len=6), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string +character(len=6), parameter :: REGRIDDING_ARBITRARY_STRING = "ARB" !< Arbitrary coordinates character(len=6), parameter :: REGRIDDING_HYCOM1_STRING = "HYCOM1" !< Hycom string character(len=6), parameter :: REGRIDDING_SLIGHT_STRING = "SLIGHT" !< Hybrid S-rho string -character(len=5), parameter :: DEFAULT_COORDINATE_MODE = REGRIDDING_LAYER_STRING !< Default coordinate mode +character(len=6), parameter :: DEFAULT_COORDINATE_MODE = REGRIDDING_LAYER_STRING !< Default coordinate mode + +integer, dimension(REGRIDDING_NUM_TYPES), parameter :: vertical_coords = & + (/ REGRIDDING_LAYER, REGRIDDING_ZSTAR, REGRIDDING_RHO, & + REGRIDDING_SIGMA, REGRIDDING_ARBITRARY, & + REGRIDDING_HYCOM1, REGRIDDING_SLIGHT /) + +character(len=*), dimension(REGRIDDING_NUM_TYPES), parameter :: vertical_coord_strings = & + (/ REGRIDDING_LAYER_STRING, REGRIDDING_ZSTAR_STRING, REGRIDDING_RHO_STRING, & + REGRIDDING_SIGMA_STRING, REGRIDDING_ARBITRARY_STRING, & + REGRIDDING_HYCOM1_STRING, REGRIDDING_SLIGHT_STRING /) interface coordinateUnits module procedure coordinateUnitsI @@ -44,6 +58,7 @@ function coordinateMode(string) case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1 case (trim(REGRIDDING_SLIGHT_STRING)); coordinateMode = REGRIDDING_SLIGHT + case (trim(REGRIDDING_ARBITRARY_STRING)); coordinateMode = REGRIDDING_ARBITRARY case default ; call MOM_error(FATAL, "coordinateMode: "//& "Unrecognized choice of coordinate ("//trim(string)//").") end select diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4baee857da..cef723db8c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -33,7 +33,7 @@ module MOM use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : diag_mediator_infrastructure_init use MOM_diag_mediator, only : diag_register_area_ids -use MOM_diag_mediator, only : diag_set_thickness_ptr, diag_update_target_grids +use MOM_diag_mediator, only : diag_set_state_ptrs, diag_update_remap_grids 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 @@ -737,7 +737,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) call post_diags_TS_vardec(G, CS, dtdia) @@ -819,7 +819,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) endif endif @@ -1064,7 +1064,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. This needs to ! happen after the H update and before the next post_data. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) call post_diags_TS_vardec(G, CS, CS%dt_trans) @@ -1942,17 +1942,18 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! and before MOM_diagnostics_init call diag_masks_set(G, GV%ke, CS%missing, 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 up pointers within diag mediator control structure, + ! this needs to occur _after_ CS%h etc. have been allocated. + call diag_set_state_ptrs(CS%h, CS%T, CS%S, CS%tv%eqn_of_state, diag) ! This call sets up the diagnostic axes. These are needed, ! e.g. to generate the target grids below. call set_axes_info(G, GV, param_file, diag) - ! Whenever thickness changes let the diag manager know, target grids - ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(diag) + ! Whenever thickness/T/S changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! FIXME: are h, T, S updated at the same time? Review these for T, S updates. + call diag_update_remap_grids(diag) ! Diagnose static fields AND associate areas/volumes with axes call write_static_fields(G, CS%diag) diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index bce371cb19..a1eae25771 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -79,7 +79,7 @@ module MOM_dynamics_legacy_split use MOM_diag_mediator, only : diag_mediator_init, enable_averaging 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 : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -937,7 +937,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & call cpu_clock_end(id_clock_continuity) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) if (G%nonblocking_updates) then call cpu_clock_begin(id_clock_pass) pid_h = pass_var_start(h, G%Domain) @@ -978,7 +978,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & call cpu_clock_end(id_clock_continuity) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) call cpu_clock_end(id_clock_pass) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 4bec4fec97..0b2064634a 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -15,7 +15,7 @@ module MOM_dynamics_split_RK2 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging 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 : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids use MOM_domains, only : MOM_domains_init use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -837,7 +837,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_end(id_clock_pass) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") call cpu_clock_begin(id_clock_pass) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index fc392cbb88..f95084a851 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_unsplit use MOM_diag_mediator, only : diag_mediator_init, enable_averaging 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 : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -448,7 +448,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & call cpu_clock_end(id_clock_pass) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) call enable_averaging(0.5*dt, Time_local, CS%diag) ! Here the second half of the thickness fluxes are offered for averaging. diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 06c616f938..e8e59b1186 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -30,20 +30,22 @@ module MOM_diag_mediator use MOM_coms, only : PE_here use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE -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_error_handler, only : MOM_error, FATAL, is_root_pe, assert +use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, field_exists, field_size use MOM_io, only : slasher, vardesc, query_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 -use MOM_remapping, only : remapping_CS, initialize_remapping, dzFromH1H2 -use MOM_remapping, only : remapping_core_h -use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution -use MOM_regridding, only : build_zstar_column, set_regrid_params use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : EOS_type +use MOM_diag_remap, only : diag_remap_ctrl +use MOM_diag_remap, only : diag_remap_update, diag_remap_set_diag_axes +use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap +use MOM_diag_remap, only : diag_remap_set_vertical_axes, diag_remap_get_nz +use MOM_diag_remap, only : diag_remap_axes_setup_done +use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE, REGRIDDING_NUM_TYPES +use regrid_consts, only : vertical_coords, vertical_coord_strings use diag_axis_mod, only : get_diag_axis_name use diag_manager_mod, only : diag_manager_init, diag_manager_end @@ -53,18 +55,12 @@ module MOM_diag_mediator 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 +#define __DO_SAFETY_CHECKS__ #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) - -#define __DO_SAFETY_CHECKS__ public set_axes_info, post_data, register_diag_field, time_type public post_data_1d_k @@ -76,10 +72,10 @@ module MOM_diag_mediator public diag_axis_init, ocean_register_diag, register_static_field public register_scalar_field public define_axes_group, diag_masks_set -public diag_set_thickness_ptr -public diag_update_target_grids public diag_register_area_ids public diag_register_volume_ids +public diag_set_state_ptrs, diag_update_remap_grids +public post_data_3d_low interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -114,6 +110,9 @@ module MOM_diag_mediator integer :: fms_diag_id !< Underlying FMS diag_manager id. character(16) :: debug_str = '' !< For FATAL errors and debugging. type(axes_grp), pointer :: remap_axes => null() + integer :: vertical_coord + integer, dimension(:), allocatable :: axes + real, pointer, dimension(:,:) :: mask2d => null() real, pointer, dimension(:,:,:) :: mask3d => null() type(diag_type), pointer :: next => null() !< Pointer to the next diag. @@ -140,7 +139,6 @@ module MOM_diag_mediator type(axes_grp) :: axesBi, axesTi, axesCui, axesCvi type(axes_grp) :: axesB1, axesT1, axesCu1, axesCv1 type(axes_grp) :: axesZi, axesZL - type(axes_grp) :: axesTzi, axesTZL, axesBZL, axesCuZL, axesCvZL ! Mask arrays for diagnostics real, dimension(:,:), pointer :: mask2dT => null() @@ -165,28 +163,13 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations real :: missing_value = 1.0e+20 - ! Needed to diagnostic regridding using ALE - type(remapping_CS) :: remap_cs - type(regridding_CS) :: regrid_cs - ! Nominal interface locations for Z remapping - real, dimension(:), allocatable :: zi_remap - ! Nominal layer locations for Z remapping - real, dimension(:), allocatable :: zl_remap - ! Number of z levels used for remapping - integer :: nz_remap + type(diag_remap_ctrl), dimension(REGRIDDING_NUM_TYPES) :: diag_remaps - ! Output grid thicknesses - real, dimension(:,:,:), allocatable :: h_zoutput - - ! Keep track of which remapping is needed for diagnostic output - logical :: do_z_remapping_on_u, do_z_remapping_on_v, do_z_remapping_on_T - logical :: remapping_initialized - - !> String appended to module name for z*-remapped diagnostics - character(len=6) :: z_remap_suffix = '_z_new' - - ! Pointer to H and G for remapping + ! Pointer to H, G and T&S needed for remapping real, dimension(:,:,:), pointer :: h => null() + real, dimension(:,:,:), pointer :: T => null() + real, dimension(:,:,:), pointer :: S => null() + type(EOS_type), pointer :: eqn_of_state => null() type(ocean_grid_type), pointer :: G => null() #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) @@ -198,7 +181,7 @@ module MOM_diag_mediator end type diag_ctrl ! CPU clocks -integer :: id_clock_diag_mediator, id_clock_diag_z_remap, id_clock_diag_grid_updates +integer :: id_clock_diag_mediator, id_clock_diag_remap, id_clock_diag_grid_updates contains @@ -211,12 +194,10 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) logical, optional, intent(in) :: set_vertical !< If true or missing, set up !! vertical axes ! Local variables - integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_zzl, id_zzi - integer :: k, nz - integer :: nzi(4) + integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh + integer :: i, k, nz real :: zlev(GV%ke), zinter(GV%ke+1) logical :: set_vert - character(len=200) :: inputdir, string, filename, varname, dimname ! This include declares and sets the variable "version". #include "version_variable.h" @@ -257,108 +238,6 @@ subroutine set_axes_info(G, GV, 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, "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"//& - " UNIFORM - vertical grid is uniform\n"//& - " between surface and max depth.\n",& - default="") - if (len_trim(string) > 0) then - if (trim(string) == 'UNIFORM') then - ! initialise a uniform coordinate with depth - nzi(1) = GV%ke + 1 - allocate(diag_cs%zi_remap(nzi(1))) - allocate(diag_cs%zl_remap(nzi(1) - 1)) - - diag_cs%zi_remap(1) = 0 - do k = 2,nzi(1) - diag_cs%zi_remap(K) = diag_cs%zi_remap(K - 1) + G%max_depth / real(nzi(1) - 1) - enddo - elseif (index(trim(string), 'FILE:') == 1) then - ! read coordinate information from a file - if (string(6:6)=='.' .or. string(6:6)=='/') then - inputdir = "." - filename = trim(extractWord(trim(string(6:200)), 1)) - else - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - filename = trim(inputdir) // trim(extractWord(trim(string(6:200)), 1)) - endif - varname = trim(extractWord(trim(string(1:200)), 2)) - dimname = trim(extractWord(trim(string(1: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 - ! Log the expanded result as a comment since it cannot be read back in - call log_param(param_file, mod, "! Remapping z diagnostics", & - 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) - else - ! unsupported method - call MOM_error(FATAL,"set_axes_info: "//& - "Incorrect remapping grid specification. Only 'FILE:file,var' and"//& - "'UNIFORM' are currently supported."//& - "Found '"//trim(string)//"'") - endif - - diag_cs%zi_remap(:) = abs(diag_cs%zi_remap(:)) ! Always convert heights into depths - ! 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) - call get_param(param_file, mod, "DIAG_REMAP_Z_MODULE_SUFFIX", diag_cs%z_remap_suffix, & - 'This is the string attached to the end of "ocean_model"\n'// & - 'for use in the model column of the diag_table to indicate\n'// & - 'a diagnostic should be remapped to z*-coordinates.', & - default='_z_new') - if (trim(diag_cs%z_remap_suffix) == '_z') then - ! This will conflict with the older MOM_diag_to_Z module for z-output - call get_param(param_file, mod, "Z_OUTPUT_GRID_FILE", string, default="", do_not_log=.true.) - if (len(trim(string))>0) call MOM_error(FATAL,"MOM_diag_mediator, set_axes_info: "// & - "Z_OUTPUT_GRID_FILE must be blank to use DIAG_REMAP_Z_MODULE_SUFFIX='_z'") - endif - 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 - - ! Axes for z remapping - call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & - x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', is_h_point=.true.) - call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & - x_cell_method='point', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) - call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL, & - x_cell_method='point', y_cell_method='mean', v_cell_method='mean', is_h_point=.false.) - call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL, & - x_cell_method='mean', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) - ! Vertical axes for the interfaces and layers call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, & v_cell_method='point') @@ -395,6 +274,10 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & x_cell_method='mean', y_cell_method='point', is_h_point=.false.) + do i=1, size(diag_cs%diag_remaps) + call diag_remap_set_vertical_axes(diag_cs%diag_remaps(i), G, GV, param_file) + enddo + end subroutine set_axes_info !> Attaches the id of cell areas to axes groupsfor use with cell_measures @@ -432,36 +315,6 @@ subroutine diag_register_volume_ids(diag_cs, id_vol_t) endif end subroutine diag_register_volume_ids -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 - !> Defines a group of "axes" from list of handles subroutine define_axes_group(diag_cs, handles, axes, & x_cell_method, y_cell_method, v_cell_method, is_h_point) @@ -613,6 +466,7 @@ subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) type(diag_type), pointer :: diag => null() if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) + ! Iterate over list of diag 'variants' (e.g. CMOR aliases) and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & 'post_data_2d: Unregistered diagnostic id') @@ -735,9 +589,12 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) ! (in,opt) mask - If present, use this real array as the data mask. type(diag_type), pointer :: diag => null() - real, allocatable :: remapped_field(:,:,:) + integer :: nz, i, j, k + real, dimension(:,:,:), allocatable :: remapped_field + real :: missing_value if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) + ! 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, & @@ -745,208 +602,68 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) diag => diag_cs%diags(diag_field_id) do while (associated(diag)) - if (associated(diag%remap_axes)) then - ! Remap this field to another vertical coordinate. + if (diag%vertical_coord == coordinateMode(DEFAULT_COORDINATE_MODE)) then + ! No remapping + call post_data_3d_low(diag, field, diag_cs, is_static, mask) + else + ! Remap this field to another vertical coordinate. if (present(mask)) then call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif - if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) - 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 (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) - 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) + if (.not. diag_remap_axes_setup_done(diag_cs%diag_remaps(diag%vertical_coord))) then + call MOM_error(FATAL,"post_data_3d: attempt to do diagnostic vertical "// & + "remapping but vertical axes not configured. "// & + "See option DIAG_REMAP__GRID_DEF") endif - if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) - deallocate(remapped_field) - if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) - else - call post_data_3d_low(diag, field, diag_cs, is_static, mask) - endif - diag => diag%next - enddo - if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) -end subroutine post_data_3d + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + + nz = diag_remap_get_nz(diag_cs%diag_remaps(diag%vertical_coord)) + allocate(remapped_field(size(field, 1), size(field, 2), nz)) -!> Remap diagnostic field to z-star vertical grid. -!! \note This uses grids generated by diag_update_target_grids() -subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information - type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure - real, dimension(:,:,:), intent(inout) :: remapped_field !< Field argument remapped to z star coordinate - ! Local variables - real, dimension(diag_cs%nz_remap+1) :: dz - real, dimension(diag_cs%nz_remap) :: h_dest - real, dimension(size(diag_cs%h, 3)) :: h_src - integer :: nz_src, nz_dest - integer :: i, j, k - - call assert(diag_cs%remapping_initialized, & - 'remap_diag_to_z: Remmaping not initialized.') - call assert(size(field, 3) == size(diag_cs%h, 3), & - 'remap_diag_to_z: Remap field and thickness z-axes do not match.') - - remapped_field = diag_cs%missing_value - nz_src = size(field, 3) - nz_dest = diag_cs%nz_remap - - if (is_u_axes(diag%remap_axes, diag_cs)) then - do j=diag_cs%G%jsc, diag_cs%G%jec - do I=diag_cs%G%iscB, diag_cs%G%iecB - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j,1)+diag%mask3d(i+1,j,1) == 0.) cycle - endif #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - ! Check that H is up-to-date. - do k=RANGE_K(diag_cs%h) - if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & - "remap_diag_to_z: H has changed since remapping grids were updated."//& - " diag debug hint: "//diag%debug_str) - enddo -#endif - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) - h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i+1,j,:) ) - call remapping_core_h(nz_src, h_src(:), field(I,j,:), & - nz_dest, h_dest(:), remapped_field(I,j,:), diag_cs%remap_cs) - do k=1, nz_dest - if (h_dest(k) == 0.) then ! This only works for z-like output - remapped_field(i, j, k:nz_dest) = diag_cs%missing_value - exit - endif - enddo - enddo - enddo - elseif (is_v_axes(diag%remap_axes, diag_cs)) then - do J=diag_cs%G%jscB, diag_cs%G%jecB - do i=diag_cs%G%isc, diag_cs%G%iec - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j,1)+diag%mask3d(i,j+1,1) == 0.) cycle + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + do j=RANGE_J(diag_cs%h) + do i=RANGE_I(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) then + call MOM_error(FATAL, & + "post_data_3d: H has changed since remapping grids were updated") endif -#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - ! Check that H is up-to-date. - do k=RANGE_K(diag_cs%h) - if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & - "remap_diag_to_z: H has changed since remapping grids were updated."//& - " diag debug hint: "//diag%debug_str) - enddo -#endif - h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) - h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i,j+1,:) ) - call remapping_core_h(nz_src, h_src(:), field(i,J,:), & - nz_dest, h_dest(:), remapped_field(i,J,:), diag_cs%remap_cs) - do k=1, nz_dest - if (h_dest(k) == 0.) then ! This only works for z-like output - remapped_field(i,J,k:nz_dest) = diag_cs%missing_value - exit - endif - enddo enddo enddo - else - do j=diag_cs%G%jsc, diag_cs%G%jec - do i=diag_cs%G%isc, diag_cs%G%iec - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j, 1) == 0.) cycle - endif -#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - ! Check that H is up-to-date. - do k=RANGE_K(diag_cs%h) - if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & - "remap_diag_to_z: H has changed since remapping grids were updated."//& - " diag debug hint: "//diag%debug_str) - enddo + enddo #endif - h_dest(:) = diag_cs%h_zoutput(i,j,:) - call remapping_core_h(nz_src, diag_cs%h(i,j,:), field(i,j,:), & - nz_dest, h_dest(:), remapped_field(i,j,:), diag_cs%remap_cs) - do k=1, nz_dest - if (h_dest(k) == 0.) then ! This only works for z-like output - remapped_field(i,j,k:nz_dest) = diag_cs%missing_value - exit - endif - enddo - enddo - enddo - endif -end subroutine remap_diag_to_z + call diag_remap_do_remap(diag_cs%diag_remaps(diag%vertical_coord), & + diag_cs%G, diag_cs%h, diag%axes, diag%mask3d, & + diag_cs%missing_value, field, remapped_field) -!> Build/update target vertical grids for diagnostic remapping. -!! \note The target grids need to be updated whenever sea surface -!! height changes. -subroutine diag_update_target_grids(diag_cs) - type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure - - ! Local variables - type(ocean_grid_type), pointer :: G - real :: depth - integer :: nz_src, nz_dest - integer :: i, j, k - logical :: force, h_changed - real, dimension(size(diag_cs%h, 3)) :: h_src_1d - real, dimension(diag_cs%nz_remap) :: h_zout_1d - real, dimension(diag_cs%nz_remap+1) :: z_out_1d - - nz_dest = diag_cs%nz_remap - nz_src = size(diag_cs%h, 3) - G => diag_cs%G - - ! The interface positions for z remapping were not provided, so don't do - ! anything. If z remapping of diagnostics is requested then an error will - ! be triggered on post(). See param DIAG_REMAP_Z_GRID_DEF - if (.not. allocated(diag_cs%zi_remap)) then - return - endif - if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) - - if (.not. diag_cs%remapping_initialized) then - call assert(allocated(diag_cs%zi_remap), & - 'update_diag_target_grids: Remapping axis not initialized') + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) - ! Initialize remapping system, on the first call - call initialize_regridding(nz_dest, 'Z*', 'PPM_IH4', diag_cs%regrid_cs) - call initialize_remapping(diag_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false.) - call set_regrid_params(diag_cs%regrid_cs, min_thickness=0.) - call setCoordinateResolution(diag_cs%zi_remap(2:) - & - diag_cs%zi_remap(:nz_dest), diag_cs%regrid_cs) - - allocate(diag_cs%h_zoutput(G%isd:G%ied,G%jsd:G%jed,nz_dest)) - endif - - ! Store z-star thicknesses on h-points - do j=G%jsc, G%jec+1 - do i=G%isc, G%iec+1 - if (G%mask2dT(i,j)==0.) then - diag_cs%h_zoutput(i,j,:) = 0. - cycle + 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(:,:,:nz)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) endif - h_src_1d(:) = diag_cs%h(i,j,:) - call build_zstar_column(diag_cs%regrid_cs, nz_dest, G%bathyT(i,j), & - sum(h_src_1d(:)), z_out_1d) - h_zout_1d(:) = z_out_1d(1:nz_dest) - z_out_1d(2:nz_dest+1) - diag_cs%h_zoutput(i,j,:) = h_zout_1d(:) - enddo + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + deallocate(remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) + endif + + diag => diag%next enddo - diag_cs%remapping_initialized = .true. -#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - ! Keep a copy of H - used to check whether grids are up-to-date - ! when doing remapping. - diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:) -#endif - if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates) + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) + +end subroutine post_data_3d -end subroutine diag_update_target_grids subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag @@ -1138,180 +855,122 @@ function register_diag_field(module_name, field_name, axes, init_time, & character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction. Use '' have no method. character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction. Use '' have no method. real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file + ! Local variables real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs - type(diag_type), pointer :: diag => null(), cmor_diag => null(), z_remap_diag => null(), cmor_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, cm_string, msg + type(diag_type), pointer :: diag + integer, dimension(3) :: remap_handles + integer :: i, primary_id, fms_id + character(len=256) :: posted_cmor_units, posted_cmor_standard_name + character(len=256) :: posted_cmor_long_name, cm_string, msg + character(len=256) :: new_module_name MOM_missing_value = axes%diag_cs%missing_value if(present(missing_value)) MOM_missing_value = missing_value diag_cs => axes%diag_cs - primary_id = -1 + primary_id = DIAG_FIELD_NOT_FOUND diag => null() - cmor_diag => null() - z_remap_diag => null() - cmor_z_remap_diag => null() - ! Set up the 'primary' diagnostic, first get an underlying FMS id - fms_id = register_diag_field_expand_axes(module_name, field_name, axes, & - 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) - call attach_cell_methods(fms_id, axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) - if (fms_id /= DIAG_FIELD_NOT_FOUND) then - ! If the diagnostic is needed then allocate and id and space - primary_id = get_new_diag_id(diag_cs) - call alloc_diag_with_id(primary_id, diag_cs, diag) - call assert(associated(diag), 'register_diag_field: diag allocation failed') - diag%fms_diag_id = fms_id - diag%debug_str = trim(field_name) - call set_diag_mask(diag, diag_cs, axes) - if (present(conversion)) diag%conversion_factor = conversion - endif - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(diag), module_name, field_name, cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif - - ! For the CMOR variation of the both the native and _z diagnostics - if (present(cmor_field_name)) then - ! Fallback values for strings set to "NULL" - posted_cmor_units = "not provided" ! - posted_cmor_standard_name = "not provided" ! Values might be able to be replaced with a CS%missing field? - posted_cmor_long_name = "not provided" ! - - ! If attributes are present for MOM variable names, use them first for the register_diag_field - ! call for CMOR verison of the variable - if (present(units)) posted_cmor_units = units - if (present(standard_name)) posted_cmor_standard_name = standard_name - if (present(long_name)) posted_cmor_long_name = long_name - - ! If specified in the call to register_diag_field, override attributes with the CMOR versions - if (present(cmor_units)) posted_cmor_units = cmor_units - if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name - if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name - endif - - ! Set up the CMOR variation of the native diagnostic - if (present(cmor_field_name)) then - fms_id = register_diag_field_expand_axes(module_name, cmor_field_name, axes, init_time, & - long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), & - missing_value=MOM_missing_value, range=range, mask_variant=mask_variant, & - standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, & - err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) - call attach_cell_methods(fms_id, axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) - if (fms_id /= DIAG_FIELD_NOT_FOUND) then - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - ! This will add the cmore variation to the 'primary' diagnostic. - ! In the case where there is no primary, it will become the primary. - call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) - cmor_diag%fms_diag_id = fms_id - cmor_diag%debug_str = trim(cmor_field_name) - call set_diag_mask(cmor_diag, diag_cs, axes) - if (present(conversion)) cmor_diag%conversion_factor = conversion - endif - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, cm_string, & - msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) + do i=1,size(vertical_coords) + if (vertical_coord_strings(i) /= 'LAYER') then + new_module_name = trim(module_name//'_'//lowercase(trim(vertical_coord_strings(i)))) + else + new_module_name = module_name endif - endif - ! Remap to z vertical coordinate, note that only diagnostics on layers - ! (not interfaces) are supported, also B axes are not supported yet - if (is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) .and. axes%rank == 3) then - if (get_diag_field_id_fms(trim(module_name)//trim(diag_cs%z_remap_suffix), field_name) /= DIAG_FIELD_NOT_FOUND) 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 DIAG_REMAP_Z_GRID_DEF') + if (get_diag_field_id_fms(new_module_name, field_name) /= DIAG_FIELD_NOT_FOUND) then + if (primary_id == DIAG_FIELD_NOT_FOUND) then + primary_id = get_new_primary_diag_id(diag_cs) endif - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) + call alloc_diag_with_id(primary_id, diag_cs, diag) + call assert(associated(diag), 'register_diag_field: diag allocation failed') + call set_diag_mask_and_axes(diag, diag_cs, axes) + diag%vertical_coord = vertical_coords(i) + + ! Register for vertical remapping + if (axes%rank == 3 .and. & + vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then + call diag_remap_set_diag_axes(diag_cs%diag_remaps(i), diag%axes) endif - call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) - call set_diag_mask(z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(z_remap_diag, diag_cs, axes) - if (present(conversion)) z_remap_diag%conversion_factor = conversion - call assert(associated(z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_expand_axes(module_name//trim(diag_cs%z_remap_suffix), field_name, & - z_remap_diag%remap_axes, & - 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) - call attach_cell_methods(fms_id, z_remap_diag%remap_axes, cm_string, & + + fms_id = register_diag_field_expand_axes(new_module_name, field_name, diag%axes, & + 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) + call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method) - z_remap_diag%fms_diag_id = fms_id - z_remap_diag%debug_str = trim(field_name) + call assert(fms_id /= DIAG_FIELD_NOT_FOUND, & + 'register_diag_field: Could not register diag') + diag%fms_diag_id = fms_id - if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. - elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. - else - diag_cs%do_z_remapping_on_T = .true. + if (is_root_pe() .and. diag_CS%doc_unit > 0) then + msg = '' + if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' + call log_available_diag(associated(diag), module_name, field_name, cm_string, & + msg, diag_CS, long_name, units, standard_name) endif endif - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), field_name, & - cm_string, msg, diag_CS, long_name, units, standard_name) + + if (.not. present(cmor_field_name)) then + cycle endif - ! Remap to z vertical coordinate with CMOR names and attributes - if (present(cmor_field_name)) then - if (get_diag_field_id_fms(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name) /= DIAG_FIELD_NOT_FOUND) 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 DIAG_REMAP_Z_GRID_DEF') - endif - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - call alloc_diag_with_id(primary_id, diag_cs, cmor_z_remap_diag) - call set_diag_mask(cmor_z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(cmor_z_remap_diag, diag_cs, axes) - if (present(conversion)) cmor_z_remap_diag%conversion_factor = conversion - call assert(associated(cmor_z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_expand_axes(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & - cmor_z_remap_diag%remap_axes, & - init_time, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=trim(posted_cmor_standard_name), & - verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & - interp_method=interp_method, tile_count=tile_count) - call attach_cell_methods(fms_id, cmor_z_remap_diag%remap_axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) - cmor_z_remap_diag%fms_diag_id = fms_id - cmor_z_remap_diag%debug_str = trim(cmor_field_name) - - if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. - elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. - else - diag_cs%do_z_remapping_on_T = .true. - endif + ! Set up the CMOR variation of the native diagnostic + if (get_diag_field_id_fms(new_module_name, cmor_field_name) /= DIAG_FIELD_NOT_FOUND) then + ! Fallback values for strings set to "NULL" + ! Values might be able to be replaced with a CS%missing field? + posted_cmor_units = "not provided" + posted_cmor_standard_name = "not provided" + posted_cmor_long_name = "not provided" + + ! If attributes are present for MOM variable names, use them first for the register_diag_field + ! call for CMOR verison of the variable + if (present(units)) posted_cmor_units = units + if (present(standard_name)) posted_cmor_standard_name = standard_name + if (present(long_name)) posted_cmor_long_name = long_name + + ! If specified in the call to register_diag_field, override attributes with the CMOR versions + if (present(cmor_units)) posted_cmor_units = cmor_units + if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name + if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name + + if (primary_id == DIAG_FIELD_NOT_FOUND) then + primary_id = get_new_primary_diag_id(diag_cs) endif + ! This will add the CMOR variation to the 'primary' diagnostic. + ! In the case where there is no primary, it will become the primary. + call alloc_diag_with_id(primary_id, diag_cs, diag) + call set_diag_mask_and_axes(diag, diag_cs, axes) + diag%vertical_coord = vertical_coords(i) + + ! Register for vertical remapping + if (axes%rank == 3 .and. & + vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then + call diag_remap_set_diag_axes(diag_cs%diag_remaps(i), diag%axes) + endif + + fms_id = register_diag_field_expand_axes(new_module_name, cmor_field_name, diag%axes, & + init_time, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), & + missing_value=MOM_missing_value, range=range, mask_variant=mask_variant, & + standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log,& + err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) + call attach_cell_methods(fms_id, axes, cm_string, & + cell_methods, x_cell_method, y_cell_method, v_cell_method) + call assert(fms_id /= DIAG_FIELD_NOT_FOUND, & + 'register_diag_field: Could not register diag') + diag%fms_diag_id = fms_id + if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(cmor_z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & - cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & + call log_available_diag(associated(diag), module_name, cmor_field_name, cm_string, & + msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif endif - endif + enddo register_diag_field = primary_id @@ -1501,7 +1160,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & MOM_missing_value = diag_cs%missing_value if(present(missing_value)) MOM_missing_value = missing_value - primary_id = -1 + primary_id = DIAG_FIELD_NOT_FOUND diag => null() cmor_diag => null() @@ -1509,7 +1168,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, standard_name=standard_name, do_not_log=do_not_log, err_msg=err_msg) if (fms_id /= DIAG_FIELD_NOT_FOUND) then - primary_id = get_new_diag_id(diag_cs) + primary_id = get_new_primary_diag_id(diag_cs) call alloc_diag_with_id(primary_id, diag_cs, diag) call assert(associated(diag), 'register_scalar_field: diag allocation failed') diag%fms_diag_id = fms_id @@ -1538,8 +1197,8 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & missing_value=MOM_missing_value, range=range, & standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, err_msg=err_msg) if (fms_id /= DIAG_FIELD_NOT_FOUND) then - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) + if (primary_id == DIAG_FIELD_NOT_FOUND) then + primary_id = get_new_primary_diag_id(diag_cs) endif call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id @@ -1606,7 +1265,7 @@ function register_static_field(module_name, field_name, axes, & if(present(missing_value)) MOM_missing_value = missing_value diag_cs => axes%diag_cs - primary_id = -1 + primary_id = DIAG_FIELD_NOT_FOUND diag => null() cmor_diag => null() @@ -1616,7 +1275,7 @@ function register_static_field(module_name, field_name, axes, & do_not_log=do_not_log, & interp_method=interp_method, tile_count=tile_count, area=area) if (fms_id /= DIAG_FIELD_NOT_FOUND) then - primary_id = get_new_diag_id(diag_cs) + primary_id = get_new_primary_diag_id(diag_cs) call alloc_diag_with_id(primary_id, diag_cs, diag) call assert(associated(diag), 'register_static_field: diag allocation failed') diag%fms_diag_id = fms_id @@ -1646,8 +1305,8 @@ function register_static_field(module_name, field_name, axes, & standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, & interp_method=interp_method, tile_count=tile_count, area=area) if (fms_id /= DIAG_FIELD_NOT_FOUND) then - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) + if (primary_id == DIAG_FIELD_NOT_FOUND) then + primary_id = get_new_primary_diag_id(diag_cs) endif call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id @@ -1817,7 +1476,7 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) character(len=40) :: mod = "MOM_diag_mediator" ! This module's name. id_clock_diag_mediator = cpu_clock_id('(Ocean diagnostics framework)', grain=CLOCK_MODULE) - id_clock_diag_z_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=CLOCK_ROUTINE) + id_clock_diag_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=CLOCK_ROUTINE) id_clock_diag_grid_updates = cpu_clock_id('(Ocean diagnostics grid updates)', grain=CLOCK_ROUTINE) ! Allocate and initialize list of all diagnostics (and variants) @@ -1827,14 +1486,19 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) call initialize_diag_type(diag_cs%diags(i)) enddo - ! Keep a pointer to the grid, this is needed for regridding + ! Initialise vertical coordinates and test assumptions needed for diagnostic + ! remapping. + do i=1, size(diag_cs%diag_remaps) + call diag_remap_init(diag_cs%diag_remaps(i), vertical_coords(i)) + enddo + + ! Keep pointers grid, h, T, S needed diagnostic remapping diag_cs%G => G diag_cs%h => null() - diag_cs%nz_remap = -1 - diag_cs%do_z_remapping_on_u = .false. - diag_cs%do_z_remapping_on_v = .false. - diag_cs%do_z_remapping_on_T = .false. - diag_cs%remapping_initialized = .false. + diag_cs%T => null() + diag_cs%S => null() + diag_cs%eqn_of_state => null() + #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz)) diag_cs%h_old(:,:,:) = 0.0 @@ -1885,19 +1549,51 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) end subroutine diag_mediator_init -subroutine diag_set_thickness_ptr(h, diag_cs) +subroutine diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs) - real, dimension(:,:,:), target, intent(in) :: h + real, dimension(:,:,:), target, intent(in) :: h, T, S + type(EOS_type), pointer, intent(in) :: eqn_of_state !< Equation of state structure type(diag_ctrl), intent(inout) :: diag_cs ! (inout) diag_cs - diag mediator control structure ! (in) h - a pointer to model thickness + ! (in) T - a pointer to model temperature + ! (in) S - a pointer to model salinity - ! Keep pointer to h, needed for the z axis remapping + ! Keep pointers to h, T, S needed for the diagnostic remapping diag_cs%h => h + diag_cs%T => T + diag_cs%S => S + diag_cs%eqn_of_state => eqn_of_state end subroutine +!> Build/update vertical grids for diagnostic remapping. +!! \note The target grids need to be updated whenever sea surface +!! height changes. +subroutine diag_update_remap_grids(diag_cs) + type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure + + integer :: i + + if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) + + do i=1, size(diag_cs%diag_remaps) + call diag_remap_update(diag_cs%diag_remaps(i), & + diag_cs%G, diag_cs%h, diag_cs%T, diag_cs%S, & + diag_cs%eqn_of_state) + enddo + +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Keep a copy of H - used to check whether grids are up-to-date + ! when doing remapping. + diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:) +#endif + + if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates) + +end subroutine diag_update_remap_grids + !> diag_masks_set sets up the 2d and 3d masks for diagnostics subroutine diag_masks_set(G, nz, missing_value, diag_cs) type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type. @@ -1951,16 +1647,19 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) type(diag_ctrl), intent(inout) :: diag_cs logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() + ! Local variables + integer :: i + if (diag_CS%doc_unit > -1) then close(diag_CS%doc_unit) ; diag_CS%doc_unit = -3 endif deallocate(diag_cs%diags) - if (allocated(diag_cs%zi_remap)) deallocate(diag_cs%zi_remap) - if (allocated(diag_cs%zl_remap)) deallocate(diag_cs%zl_remap) + do i=1, size(diag_cs%diag_remaps) + call diag_remap_end(diag_cs%diag_remaps(i)) + enddo - if (allocated(diag_cs%h_zoutput)) deallocate(diag_cs%h_zoutput) deallocate(diag_cs%mask3dTL) deallocate(diag_cs%mask3dBuL) deallocate(diag_cs%mask3dCuL) @@ -2000,33 +1699,20 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s -subroutine set_diag_remap_axes(diag, diag_cs, axes) - ! Associate a remapping axes with the a diagnostic based on the axes info. - type(diag_ctrl), target, intent(in) :: diag_cs - type(diag_type), pointer, intent(inout) :: diag - type(axes_grp), intent(in) :: axes - - diag%remap_axes => null() - if (axes%rank .eq. 3) then - if ((axes%id .eq. diag_cs%axesTL%id)) then - diag%remap_axes => diag_cs%axesTZL - elseif(axes%id .eq. diag_cs%axesBL%id) then - diag%remap_axes => diag_cs%axesBZL - elseif(axes%id .eq. diag_cs%axesCuL%id ) then - diag%remap_axes => diag_cs%axesCuZL - elseif(axes%id .eq. diag_cs%axesCvL%id) then - diag%remap_axes => diag_cs%axesCvZL - endif - endif - -end subroutine set_diag_remap_axes - -subroutine set_diag_mask(diag, diag_cs, axes) +subroutine set_diag_mask_and_axes(diag, diag_cs, axes) ! Associate a mask with the a diagnostic based on the axes info. + ! Also keep a pointer to the axes. This is needed when doing vertical + ! remapping of the diagnostic. type(diag_ctrl), target, intent(in) :: diag_cs type(diag_type), pointer, intent(inout) :: diag - type(axes_grp), intent(in) :: axes + type(axes_grp), intent(in) :: axes + + ! Local variables + integer :: i + + allocate(diag%axes(size(axes%handles))) + diag%axes(:) = axes%handles(:) diag%mask2d => null() diag%mask3d => null() @@ -2065,95 +1751,23 @@ subroutine set_diag_mask(diag, diag_cs, axes) !call assert(associated(diag%mask2d), "MOM_diag_mediator.F90: Invalid 2d axes id") endif -end subroutine set_diag_mask - -function is_layer_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_layer_axes - - is_layer_axes = .false. - - if (axes%id == diag_cs%axesTZL%id .or. & - axes%id == diag_cs%axesBZL%id .or. & - axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCvZL%id .or. & - axes%id == diag_cs%axesZL%id .or. & - axes%id == diag_cs%axesTL%id .or. & - axes%id == diag_cs%axesBL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCvL%id) then - is_layer_axes = .true. - endif - -end function is_layer_axes - -function is_u_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_u_axes - - is_u_axes = .false. - - if (axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCui%id .or. & - axes%id == diag_cs%axesCu1%id) then - is_u_axes = .true. - endif - -end function is_u_axes - -function is_v_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_v_axes - - is_v_axes = .false. - - if (axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCui%id .or. & - axes%id == diag_cs%axesCu1%id) then - is_v_axes = .true. - endif - -end function is_v_axes - -function is_B_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_B_axes +end subroutine set_diag_mask_and_axes - is_B_axes = .false. +!> Allocate a new diagnostic id, it may be necessary to expand the diagnostics +!> array. +function get_new_primary_diag_id(diag_cs) - if (axes%id == diag_cs%axesBZL%id .or. & - axes%id == diag_cs%axesBL%id .or. & - axes%id == diag_cs%axesBi%id .or. & - axes%id == diag_cs%axesB1%id) then - is_B_axes = .true. - endif - -end function is_B_axes + integer :: get_new_primary_diag_id + type(diag_ctrl), intent(inout) :: diag_cs + ! Arguments: + ! (inout) diag_cs - diagnostics control structure -!> Returns a new diagnostic id, it may be necessary to expand the diagnostics array. -integer function get_new_diag_id(diag_cs) - type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure - ! Local variables type(diag_type), dimension(:), allocatable :: tmp integer :: i if (diag_cs%next_free_diag_id > size(diag_cs%diags)) then call assert(diag_cs%next_free_diag_id - size(diag_cs%diags) == 1, & - 'get_new_diag_id: inconsistent diag id') + 'get_primary_new_diag_id: inconsistent diag id') ! Increase the size of diag_cs%diags and copy data over. ! Do not use move_alloc() because it is not supported by Fortran 90 @@ -2170,10 +1784,10 @@ integer function get_new_diag_id(diag_cs) enddo endif - get_new_diag_id = diag_cs%next_free_diag_id + get_new_primary_diag_id = diag_cs%next_free_diag_id diag_cs%next_free_diag_id = diag_cs%next_free_diag_id + 1 -end function get_new_diag_id +end function get_new_primary_diag_id !> Initializes a diag_type (used after allocating new memory) subroutine initialize_diag_type(diag) @@ -2186,12 +1800,14 @@ subroutine initialize_diag_type(diag) diag%mask3d => null() diag%next => null() diag%conversion_factor = 0. + diag%vertical_coord = coordinateMode(DEFAULT_COORDINATE_MODE) + end subroutine initialize_diag_type ! Make a new diagnostic. Either use memory which is in the array of 'primary' ! diagnostics, or if that is in use, insert it to the list of secondary diags. -subroutine alloc_diag_with_id(diag_id, diag_cs, diag) - integer, intent(in) :: diag_id +subroutine alloc_diag_with_id(primary_diag_id, diag_cs, diag) + integer, intent(in) :: primary_diag_id type(diag_ctrl), target, intent(inout) :: diag_cs type(diag_type), pointer, intent(out) :: diag @@ -2202,12 +1818,12 @@ subroutine alloc_diag_with_id(diag_id, diag_cs, diag) type(diag_type), pointer :: tmp - if (.not. diag_cs%diags(diag_id)%in_use) then - diag => diag_cs%diags(diag_id) + if (.not. diag_cs%diags(primary_diag_id)%in_use) then + diag => diag_cs%diags(primary_diag_id) else allocate(diag) - tmp => diag_cs%diags(diag_id)%next - diag_cs%diags(diag_id)%next => diag + tmp => diag_cs%diags(primary_diag_id)%next + diag_cs%diags(primary_diag_id)%next => diag diag%next => tmp endif diag%in_use = .true. @@ -2248,47 +1864,4 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string end subroutine log_available_diag -subroutine check_field_and_mask_shape_2d(diag, field, mask) - type(diag_type), intent(in) :: diag - real, intent(in) :: field(:,:) - real, intent(in) :: mask(:,:) - - integer :: i - - do i=1, size(shape(field)) - if (size(field, i) /= size(mask, i)) then - call MOM_error(FATAL,"check_field_and_mask_2d: field and mask have "//& - "different shape, diag debug hint: "//diag%debug_str) - endif - enddo - -end subroutine - -subroutine check_field_and_mask_shape_3d(diag, field, mask) - type(diag_type), intent(in) :: diag - real, intent(in) :: field(:, :, :) - real, intent(in) :: mask(:, :, :) - - integer :: i - - do i=1, size(shape(field)) - if (size(field, i) /= size(mask, i)) then - call MOM_error(FATAL,"check_field_and_mask_3d: field and mask have "//& - "different shape, diag debug hint: "//diag%debug_str) - endif - enddo - -end subroutine - -subroutine assert(logical_arg, msg) - - logical, intent(in) :: logical_arg - character(len=*), intent(in) :: msg - - if (.not. logical_arg) then - call MOM_error(FATAL, msg) - endif - -end subroutine assert - end module MOM_diag_mediator diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 new file mode 100644 index 0000000000..e9b3b61016 --- /dev/null +++ b/src/framework/MOM_diag_remap.F90 @@ -0,0 +1,567 @@ +module MOM_diag_remap + +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +!********+*********+*********+*********+*********+*********+*********+** +!* * +!* The subroutines here are used for runtime remapping of * +!* diagnostics. * +!* * +!********+*********+*********+*********+*********+*********+*********+** + +use MOM_error_handler, only : MOM_error, FATAL, assert +use MOM_file_parser, only : get_param, log_param, param_file_type +use MOM_io, only : slasher, mom_read_data +use MOM_io, only : file_exists, field_size +use MOM_string_functions, only : lowercase, extractWord +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : EOS_type +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : remapping_core_h +use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution +use MOM_regridding, only : build_zstar_column, build_rho_column, build_sigma_column +use MOM_regridding, only : set_regrid_params, uniformResolution +use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE, vertical_coord_strings + +use diag_axis_mod, only : get_diag_axis_name +use diag_manager_mod, only : diag_axis_init + +use netcdf + +implicit none ; private + +public diag_remap_ctrl, diag_remap_set_diag_axes +public diag_remap_init, diag_remap_end, diag_remap_update, diag_remap_do_remap +public diag_remap_set_vertical_axes, diag_remap_axes_setup_done, diag_remap_get_nz + +type :: diag_remap_ctrl + ! Whether remappping initialized + logical :: initialized + ! The vertical coordinate that we remap to + integer :: vertical_coord + ! Remap and regrid types needed for remaping using ALE + type(remapping_CS) :: remap_cs + type(regridding_CS) :: regrid_cs + ! Remap grid thicknesses + real, dimension(:,:,:), allocatable :: h + ! Number of vertical levels used for remapping + integer :: nz + ! Nominal layer thicknesses + real, dimension(:), allocatable :: dz + ! Axes id's for the above + integer :: interface_axes_id + integer :: layer_axes_id +end type diag_remap_ctrl + +contains + +subroutine diag_remap_init(remap_cs, vertical_coord) + type(diag_remap_ctrl), intent(inout) :: remap_cs + integer, intent(in) :: vertical_coord + + remap_cs%vertical_coord = vertical_coord + remap_cs%initialized = .false. + +end subroutine diag_remap_init + +subroutine diag_remap_end(remap_cs) + type(diag_remap_ctrl), intent(inout) :: remap_cs + + if (allocated(remap_cs%h)) deallocate(remap_cs%h) + if (allocated(remap_cs%dz)) deallocate(remap_cs%dz) + +end subroutine diag_remap_end + +subroutine diag_remap_set_vertical_axes(remap_cs, G, GV, param_file) + type(diag_remap_ctrl), intent(inout) :: remap_cs + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file structure + + if (remap_cs%vertical_coord /= coordinateMode('LAYER')) then + call setup_axes(remap_cs, G, GV, param_file) + endif + +end subroutine diag_remap_set_vertical_axes + +subroutine setup_axes(remap_cs, G, GV, param_file) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remap control structure + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file structure + + integer :: nzi(4), nzl(4), k + character(len=200) :: inputdir, string, filename, int_varname, layer_varname + character(len=40) :: mod = "MOM_diag_mediator" ! This module's name. + character(len=8) :: units, expected_units + character(len=34) :: longname + + real, allocatable, dimension(:) :: interfaces, layers + + ! Read info needed for z-space remapping + call get_param(param_file, mod, & + 'DIAG_REMAP_'//trim(vertical_coord_strings(remap_cs%vertical_coord))//'_GRID_DEF', & + string, & + "This sets the file and variable names that define the\n"//& + "vertical grid used for diagnostic output remapping. \n"//& + " 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"//& + " UNIFORM - vertical grid is uniform\n"//& + " between surface and max depth.\n",& + default="") + if (len_trim(string) > 0) then + + if (trim(string) == 'UNIFORM') then + nzi(1) = GV%ke + 1 + nzl(1) = GV%ke + allocate(remap_cs%dz(nzl(1))) + allocate(interfaces(nzi(1))) + allocate(layers(nzl(1))) + + remap_cs%dz(:) = uniformResolution(GV%ke, vertical_coord_strings(remap_cs%vertical_coord), & + G%max_depth, GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2)), & + GV%Rlay(GV%ke)+0.5*(GV%Rlay(GV%ke)-GV%Rlay(GV%ke-1))) + + if (remap_cs%vertical_coord == coordinateMode('RHO')) then + interfaces(1) = GV%Rlay(1) + else + interfaces(1) = 0 + endif + + ! Calculate interface and layer positions + do k=2,nzi(1) + interfaces(k) = interfaces(k - 1) + remap_cs%dz(k - 1) + enddo + layers(:) = interfaces(1:nzl(1)) + remap_cs%dz(:) / 2 + + elseif (index(trim(string), 'FILE:') == 1) then + ! read coordinate information from a file + if (string(6:6)=='.' .or. string(6:6)=='/') then + inputdir = "." + filename = trim(extractWord(trim(string(6:200)), 1)) + else + call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + filename = trim(inputdir) // trim(extractWord(trim(string(6:200)), 1)) + endif + int_varname = trim(extractWord(trim(string(1:200)), 2)) + layer_varname = trim(extractWord(trim(string(1: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 + + if (remap_cs%vertical_coord == coordinateMode('SIGMA')) then + expected_units = 'nondim' + elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then + expected_units = 'kg m-3' + else + expected_units = 'meters' + endif + + ! Check that the vars have expected format, units etc. + if (.not. check_grid_def(trim(filename), trim(int_varname), & + trim(expected_units))) then + call MOM_error(FATAL,"set_axes_info: Bad grid definition in "//& + "'"//trim(filename)//"'") + endif + if (.not. check_grid_def(trim(filename), trim(layer_varname), & + trim(expected_units))) then + call MOM_error(FATAL,"set_axes_info: Bad grid definition in "//& + "'"//trim(filename)//"'") + endif + + ! Log the expanded result as a comment since it cannot be read back in + call log_param(param_file, mod, "! Remapping diagnostics", & + trim(inputdir)//"/"//trim(filename)//","//trim(int_varname)//","//trim(layer_varname)) + + ! Get interfaces + call field_size(filename, int_varname, nzi) + call assert(nzi(1) /= 0, 'set_axes_info: bad interface dimension size') + allocate(interfaces(nzi(1))) + call MOM_read_data(filename, int_varname, interfaces) + ! Always convert heights into depths + interfaces(:) = abs(interfaces(:)) + + ! Get layer dimensions + call field_size(filename, layer_varname, nzl) + call assert(nzl(1) /= 0 .and. nzl(1) == nzi(1) - 1, 'set_axes_info: bad layer dimension size') + allocate(layers(nzl(1))) + call MOM_read_data(filename, layer_varname, layers) + ! Always convert heights into depths + layers(:) = abs(layers(:)) + + allocate(remap_cs%dz(nzl(1))) + remap_cs%dz(:) = interfaces(2:) - interfaces(1:nzi(1)-1) + + else + ! unsupported method + call MOM_error(FATAL,"set_axes_info: "//& + "Incorrect remapping grid specification. Only 'FILE:file,var' and"//& + "'UNIFORM' are currently supported."//& + "Found '"//trim(string)//"'") + endif + + remap_cs%nz = nzl(1) + + if (remap_cs%vertical_coord == coordinateMode('SIGMA')) then + units = 'nondim' + longname = 'Fraction' + elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then + units = 'kg m-3' + longname = 'Target Potential Density' + else + units = 'meters' + longname = 'Depth' + endif + + ! Make axes objects + remap_cs%layer_axes_id = diag_axis_init(lowercase(trim(vertical_coord_strings(remap_cs%vertical_coord)))//'_l', & + layers, trim(units), 'z', & + trim(longname)//' at cell center', direction=-1) + remap_cs%interface_axes_id = diag_axis_init(lowercase(trim(vertical_coord_strings(remap_cs%vertical_coord)))//'_i', & + interfaces, trim(units), 'z', & + trim(longname)//' Depth at interface', direction=-1) + deallocate(interfaces) + deallocate(layers) + else + ! In this case the axes associated with these will never be used, however + ! they need to be positive otherwise FMS complains. + remap_cs%layer_axes_id = 1 + remap_cs%interface_axes_id = 1 + endif + +end subroutine setup_axes + +function check_grid_def(filename, varname, expected_units) + ! Do some basic checks on the vertical grid definition file, variable + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + character(len=*), intent(in) :: expected_units + logical :: check_grid_def, check_units + + character (len=200) :: units, long_name + integer :: ncid, status, intid, vid + + check_units = .false. ! FIXME: remove this. + 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 + + if (check_units) then + status = NF90_GET_ATT(ncid, vid, "units", units) + if (status /= NF90_NOERR) then + check_grid_def = .false. + endif + if (trim(units) /= trim(expected_units)) then + if (trim(expected_units) == "meters") then + if (trim(units) /= "m") then + check_grid_def = .false. + endif + else + check_grid_def = .false. + endif + endif + endif + +end function check_grid_def + +subroutine diag_remap_set_diag_axes(remap_cs, axes) + type(diag_remap_ctrl), intent(inout) :: remap_cs + integer, dimension(:), intent(inout) :: axes + + character(len=2) :: axis_name + integer :: vertical_axis + + call assert(size(axes) == 3, & + 'diag_remap_update_diag_axes: Unexpected number of axes') + + call get_diag_axis_name(axes(3), axis_name) + + if (axis_name == 'zl') then + axes(3) = remap_cs%layer_axes_id + elseif (axis_name == 'zi') then + axes(3) = remap_cs%interface_axes_id + else + call assert(.false., & + 'diag_remap_update_diag_axes: Unexpected vertical axis name') + endif + +end subroutine + +function diag_remap_get_nz(remap_cs) + type(diag_remap_ctrl), intent(in) :: remap_cs + integer :: diag_remap_get_nz + + diag_remap_get_nz = remap_cs%nz + +end function + +function diag_remap_axes_setup_done(remap_cs) + type(diag_remap_ctrl), intent(in) :: remap_cs + logical :: diag_remap_axes_setup_done + + if (allocated(remap_cs%dz)) then + diag_remap_axes_setup_done = .true. + else + diag_remap_axes_setup_done = .false. + endif + +end function + +!> Build/update target vertical grids for diagnostic remapping. +!! \note The target grids need to be updated whenever sea surface +!! height changes. +subroutine diag_remap_update(remap_cs, G, h, T, S, eqn_of_state) + type(diag_remap_ctrl), intent(inout) :: remap_cs + type(ocean_grid_type), pointer :: G + real, dimension(:, :, :), intent(in) :: h, T, S + type(EOS_type), pointer, intent(in) :: eqn_of_state + + ! Local variables + integer :: i, j, k, nz + logical :: checked + real, dimension(remap_cs%nz + 1) :: zInterfaces + real, dimension(remap_cs%nz) :: resolution + + if (remap_cs%vertical_coord == coordinateMode('LAYER') .or. & + .not. diag_remap_axes_setup_done(remap_cs)) then + return + endif + + checked = .false. + nz = remap_cs%nz + + if (.not. remap_cs%initialized) then + ! Initialize remapping and regridding on the first call + call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false.) + + call initialize_regridding(remap_cs%nz, & + vertical_coord_strings(remap_cs%vertical_coord), 'PPM_IH4', remap_cs%regrid_cs) + call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.) + call setCoordinateResolution(remap_cs%dz, remap_cs%regrid_cs) + + allocate(remap_cs%h(G%isd:G%ied,G%jsd:G%jed,remap_cs%nz)) + endif + + ! Calculate remapping thicknesses for different target grids based on + ! nominal/target interface locations. This happens for every call on the + ! assumption that h, T, S has changed. + do j=G%jsd, G%jed + do i=G%isd, G%ied + if (G%mask2dT(i,j)==0.) then + remap_cs%h(i,j,:) = 0. + cycle + endif + + if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then + call build_zstar_column(remap_cs%regrid_cs, nz, & + G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) + elseif (remap_cs%vertical_coord == coordinateMode('SIGMA')) then + call build_sigma_column(remap_cs%regrid_cs, nz, & + G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) + elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then + call build_rho_column(remap_cs%regrid_cs, remap_cs%remap_cs, nz, & + G%bathyT(i,j), h(i,j,:), T(i, j, :), S(i, j, :), & + eqn_of_state, zInterfaces) + endif + remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) + enddo + enddo + + remap_cs%initialized = .true. + +end subroutine diag_remap_update + +subroutine diag_remap_do_remap(remap_cs, G, h, axes, mask, missing_value, field, remapped_field) + type(diag_remap_ctrl), intent(in) :: remap_cs + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + integer, dimension(3), intent(in) :: axes !< Axes for the field + real, dimension(:,:,:), intent(in) :: mask !< A mask for the field + real, intent(in) :: missing_value + real, dimension(:,:,:), intent(in) :: field(:,:,:) + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate + + ! Local variables + real, dimension(G%isd:G%ied, G%jsd:G%jed, remap_cs%nz) :: h_dest + real, dimension(size(h, 3)) :: h_src + character(len=8) :: x_axis, y_axis, z_axis + + integer :: nz_src, nz_dest + integer :: i, j, k + + nz_src = size(field, 3) + nz_dest = remap_cs%nz + + call assert(remap_cs%initialized, 'diag_remap_do_remap: remap_cs not initialized.') + call assert(size(field, 3) == size(h, 3), & + 'diag_remap_do_remap: Remap field and thickness z-axes do not match.') + + call get_diag_axis_name(axes(1), x_axis) + call get_diag_axis_name(axes(2), y_axis) + call get_diag_axis_name(axes(3), z_axis) + + call assert(.not. is_layer_axis(remap_cs, z_axis), & + 'diag_remap_do_remap: Remapping on interfaces not supported') + + remapped_field(:,:,:) = missing_value + h_dest(:,:,:) = 0. + + if (is_u_axis(x_axis, y_axis)) then + do j=G%jsc, G%jec + do I=G%iscB, G%iecB + if (mask(i,j,1)+mask(i+1,j,1) == 0.) cycle + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(i, j, :) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) + call remapping_core_h(nz_src, h_src(:), field(I,j,:), & + nz_dest, h_dest(i, j, :), remapped_field(I,j,:), & + remap_cs%remap_cs) + enddo + enddo + elseif (is_v_axis(x_axis, y_axis)) then + do J=G%jscB, G%jecB + do i=G%isc, G%iec + if (mask(i,j,1)+mask(i,j+1,1) == 0.) cycle + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(i, j, :) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:) ) + call remapping_core_h(nz_src, h_src(:), field(i,J,:), & + nz_dest, h_dest(i, j, :), remapped_field(i,J,:), & + remap_cs%remap_cs) + enddo + enddo + elseif (is_B_axis(x_axis, y_axis)) then + do j=G%jscB, G%jecB + do I=G%iscB, G%iecB + if (mask(i,j+1,1)+mask(i+1,j,1) == 0.) cycle + h_src(:) = 0.5 * (h(i,j+1,:) + h(i+1,j,:)) + h_dest(i,j,:) = 0.5 * (remap_cs%h(i,j+1,:) + remap_cs%h(i+1,j,:)) + call remapping_core_h(nz_src, h_src(:), field(I,j,:), & + nz_dest, h_dest(i,j,:), remapped_field(I,j,:), & + remap_cs%remap_cs) + enddo + enddo + elseif (is_T_axis(x_axis, y_axis)) then + do j=G%jsc, G%jec + do i=G%isc, G%iec + if (mask(i,j, 1) == 0.) cycle + h_dest(i,j,:) = remap_cs%h(i,j,:) + call remapping_core_h(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, h_dest(i,j,:), remapped_field(i,j,:), & + remap_cs%remap_cs) + enddo + enddo + else + call assert(.false., 'diag_remap_do_remap: Unsupported axis combination') + endif + + ! Now mask out zero thickness layers, only applicable to zstar output. + ! Notice we keep h_dest just for this purpose. + if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then + do j=G%jsd, G%jed + do i=G%isd, G%ied + do k=1, nz_dest + ! Find first zero-thickness layer and mask from there to the bottom. + if (h_dest(i, j, k) == 0.) then + remapped_field(i,j,k:nz_dest) = missing_value + exit + endif + enddo + enddo + enddo + endif + +end subroutine diag_remap_do_remap + +function is_u_axis(x_axis, y_axis) + character(len=*), intent(in) :: x_axis, y_axis + logical :: is_u_axis + + if (trim(x_axis) == 'xq' .and. trim(y_axis) == 'yh') then + is_u_axis = .true. + else + is_u_axis = .false. + endif + +end function is_u_axis + +function is_v_axis(x_axis, y_axis) + character(len=*), intent(in) :: x_axis, y_axis + logical :: is_v_axis + + if (trim(x_axis) == 'xh' .and. trim(y_axis) == 'yq') then + is_v_axis = .true. + else + is_v_axis = .false. + endif + +end function is_v_axis + +function is_B_axis(x_axis, y_axis) + character(len=*), intent(in) :: x_axis, y_axis + logical :: is_B_axis + + if (trim(x_axis) == 'xq' .and. trim(y_axis) == 'yq') then + is_B_axis = .true. + else + is_B_axis = .false. + endif + +end function is_B_axis + +function is_T_axis(x_axis, y_axis) + character(len=*), intent(in) :: x_axis, y_axis + logical :: is_T_axis + + if (trim(x_axis) == 'xh' .and. trim(y_axis) == 'yh') then + is_T_axis = .true. + else + is_T_axis = .false. + endif + +end function is_T_axis + +function is_layer_axis(remap_cs, z_axis) + type(diag_remap_ctrl), intent(in) :: remap_cs + character(len=*), intent(in) :: z_axis + logical :: is_layer_axis + + if (z_axis == trim(vertical_coord_strings(remap_cs%vertical_coord))//'_zl') then + is_layer_axis = .true. + else + is_layer_axis = .false. + endif + +end function is_layer_axis + +end module MOM_diag_remap diff --git a/src/framework/MOM_error_handler.F90 b/src/framework/MOM_error_handler.F90 index 9d24e35036..cc1b305729 100644 --- a/src/framework/MOM_error_handler.F90 +++ b/src/framework/MOM_error_handler.F90 @@ -18,6 +18,7 @@ module MOM_error_handler public MOM_error, MOM_mesg, NOTE, WARNING, FATAL, is_root_pe, stdlog, stdout public MOM_set_verbosity, MOM_get_verbosity, MOM_verbose_enough public callTree_showQuery, callTree_enter, callTree_leave, callTree_waypoint +public assert ! Verbosity level: ! 0 - FATAL messages only @@ -172,4 +173,15 @@ subroutine callTree_waypoint(mesg,n) endif end subroutine callTree_waypoint +subroutine assert(logical_arg, msg) + + logical, intent(in) :: logical_arg + character(len=*), intent(in) :: msg + + if (.not. logical_arg) then + call MOM_error(FATAL, msg) + endif + +end subroutine assert + end module MOM_error_handler diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index f8bf55b215..0398851a6b 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -11,7 +11,7 @@ module MOM_mixed_layer_restrat use MOM_checksums, only : hchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type -use MOM_diag_mediator, only : diag_update_target_grids +use MOM_diag_mediator, only : diag_update_remap_grids use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type @@ -372,7 +372,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, MLD, G, GV, ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. ! This needs to happen after the H update and before the next post_data. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) ! Offer diagnostic fields for averaging. if (query_averaging_enabled(CS%diag)) then @@ -612,7 +612,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, GV, CS) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) ! Offer diagnostic fields for averaging. if (query_averaging_enabled(CS%diag) .and. & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index bfdc74a14a..1502b4301f 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -48,7 +48,7 @@ module MOM_thickness_diffuse use MOM_checksums, only : hchksum, uchksum, vchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type -use MOM_diag_mediator, only : diag_update_target_grids +use MOM_diag_mediator, only : diag_update_remap_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_file_parser, only : get_param, log_version, param_file_type @@ -359,7 +359,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. ! This needs to happen after the H update and before the next post_data. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index e98e487116..8adafbde89 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -55,7 +55,7 @@ module MOM_bulk_mixed_layer use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc -use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_target_grids +use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_remap_grids use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -811,7 +811,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, CS, & ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. ! This needs to happen after the H update and before the next post_data. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) !$OMP end parallel diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 587e55f6ce..d15d18a511 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -71,7 +71,7 @@ module MOM_diabatic_aux use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr -use MOM_diag_mediator, only : diag_ctrl, time_type! , diag_update_target_grids +use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_EOS, only : calculate_density, calculate_TFreeze use MOM_EOS, only : calculate_specific_vol_derivs use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 93e2bfb4b6..66000afed4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -14,7 +14,7 @@ module MOM_diabatic_driver use MOM_diabatic_aux, only : make_frazil, adjust_salt, insert_brine, differential_diffuse_T_S, triDiagTS use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr -use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_target_grids +use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diffConvection, only : diffConvection_CS, diffConvection_init @@ -405,7 +405,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) ! Whenever thickness changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) ! Set_opacity estimates the optical properties of the water column. ! It will need to be modified later to include information about the @@ -1031,7 +1031,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) ! Whenever thickness changes let the diag manager know, as the ! target grids for vertical remapping may need to be regenerated. - call diag_update_target_grids(CS%diag) + call diag_update_remap_grids(CS%diag) ! diagnostics if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & From df6cd1720f1238cc2cf05319b09e89303c12eb11 Mon Sep 17 00:00:00 2001 From: Nic Hannah Date: Wed, 20 Jul 2016 18:27:51 +1000 Subject: [PATCH 02/16] Correct available_diags for vertically remapped diagnostics. #334 --- src/framework/MOM_diag_mediator.F90 | 34 ++++++++++++++++++----------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index e8e59b1186..50f1e018b7 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -871,15 +871,19 @@ function register_diag_field(module_name, field_name, axes, init_time, & diag_cs => axes%diag_cs primary_id = DIAG_FIELD_NOT_FOUND - diag => null() do i=1,size(vertical_coords) if (vertical_coord_strings(i) /= 'LAYER') then + ! For now we don't support remapping diagnostics onto interfaces. + if (is_interface_axis(diag_cs, axes)) then + cycle + endif new_module_name = trim(module_name//'_'//lowercase(trim(vertical_coord_strings(i)))) else new_module_name = module_name endif + diag => null() if (get_diag_field_id_fms(new_module_name, field_name) /= DIAG_FIELD_NOT_FOUND) then if (primary_id == DIAG_FIELD_NOT_FOUND) then primary_id = get_new_primary_diag_id(diag_cs) @@ -906,12 +910,15 @@ function register_diag_field(module_name, field_name, axes, init_time, & 'register_diag_field: Could not register diag') diag%fms_diag_id = fms_id - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(diag), module_name, field_name, cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif + endif + ! This diag is being registered, so it is available. + print*, 'registering ', new_module_name, field_name + if (is_root_pe() .and. diag_CS%doc_unit > 0) then + msg = '' + print*, 'calling log_available_diag' + if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' + call log_available_diag(associated(diag), new_module_name, field_name, cm_string, & + msg, diag_CS, long_name, units, standard_name) endif if (.not. present(cmor_field_name)) then @@ -919,6 +926,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif ! Set up the CMOR variation of the native diagnostic + diag => null() if (get_diag_field_id_fms(new_module_name, cmor_field_name) /= DIAG_FIELD_NOT_FOUND) then ! Fallback values for strings set to "NULL" ! Values might be able to be replaced with a CS%missing field? @@ -963,12 +971,12 @@ function register_diag_field(module_name, field_name, axes, init_time, & 'register_diag_field: Could not register diag') diag%fms_diag_id = fms_id - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(diag), module_name, cmor_field_name, cm_string, & - msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif + endif + if (is_root_pe() .and. diag_CS%doc_unit > 0) then + msg = 'native name is "'//trim(field_name)//'"' + call log_available_diag(associated(diag), new_module_name, cmor_field_name, cm_string, & + msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & + posted_cmor_standard_name) endif enddo From 52405299cc6b35d8fba5deb8209a369f60689fe4 Mon Sep 17 00:00:00 2001 From: Nic Hannah Date: Wed, 20 Jul 2016 18:53:38 +1000 Subject: [PATCH 03/16] Extend vertical remapping to handle SLIGHT and HYCOM1. #334 --- src/framework/MOM_diag_mediator.F90 | 9 ++++----- src/framework/MOM_diag_remap.F90 | 6 ++++++ 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 50f1e018b7..2cf1f9ac9d 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -614,8 +614,8 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) if (.not. diag_remap_axes_setup_done(diag_cs%diag_remaps(diag%vertical_coord))) then call MOM_error(FATAL,"post_data_3d: attempt to do diagnostic vertical "// & - "remapping but vertical axes not configured. "// & - "See option DIAG_REMAP__GRID_DEF") + "remapping but vertical axes not configured. See option "// & + "DIAG_REMAP_"//vertical_coord_strings(diag%vertical_coord)//"_GRID_DEF") endif if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) @@ -875,7 +875,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & do i=1,size(vertical_coords) if (vertical_coord_strings(i) /= 'LAYER') then ! For now we don't support remapping diagnostics onto interfaces. - if (is_interface_axis(diag_cs, axes)) then + ! Also don't supprot arbitarary remapping. + if (is_interface_axis(diag_cs, axes) .or. vertical_coord_strings(i) == 'ARB') then cycle endif new_module_name = trim(module_name//'_'//lowercase(trim(vertical_coord_strings(i)))) @@ -912,10 +913,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif ! This diag is being registered, so it is available. - print*, 'registering ', new_module_name, field_name if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = '' - print*, 'calling log_available_diag' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' call log_available_diag(associated(diag), new_module_name, field_name, cm_string, & msg, diag_CS, long_name, units, standard_name) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index e9b3b61016..a4f586bc12 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -394,6 +394,12 @@ subroutine diag_remap_update(remap_cs, G, h, T, S, eqn_of_state) call build_rho_column(remap_cs%regrid_cs, remap_cs%remap_cs, nz, & G%bathyT(i,j), h(i,j,:), T(i, j, :), S(i, j, :), & eqn_of_state, zInterfaces) + elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then + call build_slight_column(remap_cs%regrid_cs, nz, & + G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) + elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then + call build_hycom1_column(remap_cs%regrid_cs, nz, & + G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) endif remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) enddo From c3ce7bcb668adf2f1b6910031cdca65b39bd0447 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 5 Sep 2016 15:21:50 -0400 Subject: [PATCH 04/16] Use conversion=GV%H_to_kg_m2 for umo,vmo diagnostics - The umo,vmo CMOR diagnostics were using in-place global multiplications which is known to lead to unitialized variable errors within the diag manager. - Registering with conversion=GV%H_to_kg_m2 instead so that the diag_mediator can create a local array with which to do the conversion. - No answer changes. --- src/core/MOM_dynamics_split_RK2.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 4bec4fec97..b8bcdbfe67 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -891,8 +891,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (CS%id_vav > 0) call post_data(CS%id_vav, v_av, CS%diag) if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag) if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag) - if (CS%id_umo > 0) call post_data(CS%id_umo, uh*GV%H_to_kg_m2, CS%diag) - if (CS%id_vmo > 0) call post_data(CS%id_vmo, vh*GV%H_to_kg_m2, CS%diag) + if (CS%id_umo > 0) call post_data(CS%id_umo, uh, CS%diag) + if (CS%id_vmo > 0) call post_data(CS%id_vmo, vh, CS%diag) ! depth summed zonal mass transport if (CS%id_umo_2d > 0) then @@ -1199,10 +1199,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil 'Meridional Thickness Flux', flux_units) CS%id_umo = register_diag_field('ocean_model', 'umo', & diag%axesCuL, Time,'Zonal Mass Transport (including SGS param)', 'kg/s',& - cmor_standard_name='ocean_mass_x_transport', cmor_long_name='Ocean Mass X Transport') + cmor_standard_name='ocean_mass_x_transport', cmor_long_name='Ocean Mass X Transport', & + conversion=GV%H_to_kg_m2) CS%id_vmo = register_diag_field('ocean_model', 'vmo', & diag%axesCvL, Time,'Meridional Mass Transport (including SGS param)', 'kg/s',& - cmor_standard_name='ocean_mass_y_transport', cmor_long_name='Ocean Mass Y Transport') + cmor_standard_name='ocean_mass_y_transport', cmor_long_name='Ocean Mass Y Transport', & + conversion=GV%H_to_kg_m2) CS%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', & diag%axesCu1, Time,'Zonal Mass Transport (including SGS param) Vertical Sum', 'kg/s',& cmor_standard_name='ocean_mass_x_transport_vertical_sum', & From b0a60b8a8b0bfbaf9b79b72f5764d40837df5cd6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 5 Sep 2016 16:39:21 -0400 Subject: [PATCH 05/16] Re-factor of MOM_diag_mediator + cell_methods fix - register_diag_field now only has a two-way split and (used to be four for native, CMOR, z, z-CMOR). Call tree is now: - register_diag_field() calls register_diag_field_expand_cmor() for both the native and z-remapped field. - register_diag_field_expand_cmor() calls register_diag_field_expand_axes() for both the native names and CMOR names. This is where the diag_type is allocated and filled. - register_diag_field_expand_axes() calls register_diag_field_fms() one of four ways depending on presence of area interp_method optional arguments. - Fix: The z-diagnostic cell_methods was incorrectly using the native vertical dimension name for some fields when the diagnostic was not requested. This led to the wrong names being provided in available_diags. The re-factor corrected this unexpectedly. - No answer changes. --- src/framework/MOM_diag_mediator.F90 | 375 +++++++++++++++------------- 1 file changed, 199 insertions(+), 176 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 06c616f938..5dd28b9c25 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -98,6 +98,8 @@ module MOM_diag_mediator character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group includes vertical direction. ! For detecting position on the grid logical :: is_h_point = .false. !< If true, indicates that this axes group is for an h-point located field. + logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid. False for any other grid. Used for rank>2. + logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located fields that must be remapped to these axes. Used for rank>2. ! ID's for cell_measures integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp. integer :: id_volume = -1 !< The diag_manager id for volume to be used for cell_measure of variables with this axes_grp. @@ -113,7 +115,7 @@ module MOM_diag_mediator logical :: in_use !< True if this entry is being used. integer :: fms_diag_id !< Underlying FMS diag_manager id. character(16) :: debug_str = '' !< For FATAL errors and debugging. - type(axes_grp), pointer :: remap_axes => null() + type(axes_grp), pointer :: axes => null() real, pointer, dimension(:,:) :: mask2d => null() real, pointer, dimension(:,:,:) :: mask3d => null() type(diag_type), pointer :: next => null() !< Pointer to the next diag. @@ -345,19 +347,25 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) 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 + id_zzl = diag_axis_init('zl_remap', (/0.,1./), "meters", "z", & + "Depth of cell center", direction=-1) + id_zzi = diag_axis_init('zi_remap', (/0.,1./), "meters", "z", & + 'Depth of interfaces', direction=-1) endif ! Axes for z remapping call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & - x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', is_h_point=.true.) + x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & + is_h_point=.true., is_native=.false., needs_remapping=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & - x_cell_method='point', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + x_cell_method='point', y_cell_method='point', v_cell_method='mean', & + is_native=.false.) call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL, & - x_cell_method='point', y_cell_method='mean', v_cell_method='mean', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & + is_native=.false., needs_remapping=.true.) call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL, & - x_cell_method='mean', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + x_cell_method='mean', y_cell_method='point', v_cell_method='mean', & + is_native=.false., needs_remapping=.true.) ! Vertical axes for the interfaces and layers call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, & @@ -367,33 +375,35 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) ! Axis groupings for the model layers call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%axesTL, & - x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', is_h_point=.true.) + x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & + is_h_point=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%axesBL, & - x_cell_method='point', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + x_cell_method='point', y_cell_method='point', v_cell_method='mean') call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%axesCuL, & - x_cell_method='point', y_cell_method='mean', v_cell_method='mean', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean', v_cell_method='mean') call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%axesCvL, & - x_cell_method='mean', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + x_cell_method='mean', y_cell_method='point', v_cell_method='mean') ! Axis groupings for the model interfaces call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%axesTi, & - x_cell_method='mean', y_cell_method='mean', v_cell_method='point', is_h_point=.true.) + x_cell_method='mean', y_cell_method='mean', v_cell_method='point', & + is_h_point=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%axesCui, & - x_cell_method='point', y_cell_method='mean', v_cell_method='point', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean', v_cell_method='point') call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%axesCvi, & - x_cell_method='mean', y_cell_method='point', v_cell_method='point', is_h_point=.false.) + x_cell_method='mean', y_cell_method='point', v_cell_method='point') call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, & - x_cell_method='point', y_cell_method='point', v_cell_method='point', is_h_point=.false.) + x_cell_method='point', y_cell_method='point', v_cell_method='point') ! Axis groupings for 2-D arrays call define_axes_group(diag_cs, (/ id_xh, id_yh /), diag_cs%axesT1, & x_cell_method='mean', y_cell_method='mean', is_h_point=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1, & - x_cell_method='point', y_cell_method='point', is_h_point=.false.) + x_cell_method='point', y_cell_method='point') call define_axes_group(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1, & - x_cell_method='point', y_cell_method='mean', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean') call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & - x_cell_method='mean', y_cell_method='point', is_h_point=.false.) + x_cell_method='mean', y_cell_method='point') end subroutine set_axes_info @@ -464,7 +474,7 @@ end function check_grid_def !> Defines a group of "axes" from list of handles subroutine define_axes_group(diag_cs, handles, axes, & - x_cell_method, y_cell_method, v_cell_method, is_h_point) + x_cell_method, y_cell_method, v_cell_method, is_h_point, is_native, needs_remapping) type(diag_ctrl), target, intent(in) :: diag_cs !< Diagnostics control structure integer, dimension(:), intent(in) :: handles !< A list of 1D axis handles type(axes_grp), intent(out) :: axes !< The group of 1D axes @@ -472,6 +482,9 @@ subroutine define_axes_group(diag_cs, handles, axes, & character(len=*), optional, intent(in) :: y_cell_method !< A y-direction cell method used to construct the "cell_methods" attribute in CF convention character(len=*), optional, intent(in) :: v_cell_method !< A vertical direction cell method used to construct the "cell_methods" attribute in CF convention logical, optional, intent(in) :: is_h_point !< If true, indicates this axes group for h-point located fields + logical, optional, intent(in) :: is_native !< If true, indicates that this axes group is for a native model grid. False for any other grid. + logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is for a intensive layer-located fields + !! that must be remapped to these axes. Used for rank>2. ! Local variables integer :: n n = size(handles) @@ -503,6 +516,8 @@ subroutine define_axes_group(diag_cs, handles, axes, & axes%v_cell_method = '' endif if (present(is_h_point)) axes%is_h_point = is_h_point + if (present(is_native)) axes%is_native = is_native + if (present(needs_remapping)) axes%needs_remapping = needs_remapping end subroutine define_axes_group subroutine set_diag_mediator_grid(G, diag_cs) @@ -744,8 +759,9 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) 'post_data_3d: Unregistered diagnostic id') diag => diag_cs%diags(diag_field_id) do while (associated(diag)) + call assert(associated(diag%axes), 'post_data_3d: axes is not associated') - if (associated(diag%remap_axes)) then + if (diag%axes%needs_remapping) then ! Remap this field to another vertical coordinate. if (present(mask)) then @@ -799,7 +815,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) nz_src = size(field, 3) nz_dest = diag_cs%nz_remap - if (is_u_axes(diag%remap_axes, diag_cs)) then + if (is_u_axes(diag%axes, diag_cs)) then do j=diag_cs%G%jsc, diag_cs%G%jec do I=diag_cs%G%iscB, diag_cs%G%iecB if (associated(diag%mask3d)) then @@ -825,7 +841,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo enddo - elseif (is_v_axes(diag%remap_axes, diag_cs)) then + elseif (is_v_axes(diag%axes, diag_cs)) then do J=diag_cs%G%jscB, diag_cs%G%jecB do i=diag_cs%G%isc, diag_cs%G%iec if (associated(diag%mask3d)) then @@ -1106,12 +1122,123 @@ function get_diag_time_end(diag_cs) get_diag_time_end = diag_cs%time_end end function get_diag_time_end -function register_diag_field(module_name, field_name, axes, init_time, & +!> Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived from one field. +integer function register_diag_field(module_name, field_name, axes, init_time, & + long_name, units, missing_value, range, mask_variant, standard_name, & + verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & + cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & + x_cell_method, y_cell_method, v_cell_method, conversion) + character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes for this field + type(time_type), intent(in) :: init_time !< Time at which a field is first available? + character(len=*), optional, intent(in) :: long_name !< Long name of a field. + character(len=*), optional, intent(in) :: units !< Units of a field. + character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field + real, optional, intent(in) :: missing_value !< A value that indicates missing values. + real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?) + logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided with post_data calls (not used in MOM?) + logical, optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?) + logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?) + character(len=*), optional, intent(out):: err_msg !< String into which an error message might be placed (not used in MOM?) + character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not be interpolated as a scalar + integer, optional, intent(in) :: tile_count !< no clue (not used in MOM?) + character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field + character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field + character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field + character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field + character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute. Use '' to have no attribute. + !! If present, this overrides the default constructed from the default for + !! each individual axis direction. + character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction. Use '' have no method. + character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction. Use '' have no method. + character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction. Use '' have no method. + real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file + ! Local variables + real :: MOM_missing_value + type(diag_ctrl), pointer :: diag_cs + type(axes_grp), pointer :: remap_axes => null() + integer :: dm_id + character(len=256) :: cm_string, msg + logical :: active + + MOM_missing_value = axes%diag_cs%missing_value + if(present(missing_value)) MOM_missing_value = missing_value + + diag_cs => axes%diag_cs + dm_id = -1 + + ! Register the native diagnostic + active = register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, & + 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, & + cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & + cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & + cell_methods=cell_methods, x_cell_method=x_cell_method, & + y_cell_method=y_cell_method, v_cell_method=v_cell_method, & + conversion=conversion) + + ! Register diagnostics remapped to z vertical coordinate, note that only diagnostics on layers + ! (not interfaces) are supported, also B axes are not supported yet + if (axes%rank == 3 .and. (.not. is_B_axes(axes, diag_cs))) then + if (is_layer_axes(axes, diag_cs)) then + remap_axes => null() + if (axes%rank .eq. 3) then + if ((axes%id .eq. diag_cs%axesTL%id)) then + remap_axes => diag_cs%axesTZL + elseif(axes%id .eq. diag_cs%axesBL%id) then + remap_axes => diag_cs%axesBZL + elseif(axes%id .eq. diag_cs%axesCuL%id ) then + remap_axes => diag_cs%axesCuZL + elseif(axes%id .eq. diag_cs%axesCvL%id) then + remap_axes => diag_cs%axesCvZL + endif + endif + call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + active = register_diag_field_expand_cmor(dm_id, & + module_name//trim(diag_cs%z_remap_suffix), field_name, remap_axes, & + 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, & + cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & + cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & + cell_methods=cell_methods, x_cell_method=x_cell_method, & + y_cell_method=y_cell_method, v_cell_method=v_cell_method, & + conversion=conversion) + + if (active .and. .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 DIAG_REMAP_Z_GRID_DEF') + endif + + ! Record whether any remapping is needed at all + if (active) then + if (is_u_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_u = .true. + elseif (is_v_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_v = .true. + else + diag_cs%do_z_remapping_on_T = .true. + endif + endif + endif ! is_layer_axes + endif + + register_diag_field = dm_id + +end function register_diag_field + +!> Returns True if either the native of CMOr version of the diagnostic were registered. Updates 'primary_id' +!! after calling register_diag_field_expand_axes() for both native and CMOR variants of the field. +logical function register_diag_field_expand_cmor(primary_id, module_name, field_name, axes, init_time, & long_name, units, missing_value, range, mask_variant, standard_name, & verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & x_cell_method, y_cell_method, v_cell_method, conversion) - integer :: register_diag_field !< An integer handle for a diagnostic array. + integer, intent(inout) :: primary_id !< The diag_mediator ID for this diagnostic group character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model" character(len=*), intent(in) :: field_name !< Name of the diagnostic field type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes for this field @@ -1141,42 +1268,43 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Local variables real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs - type(diag_type), pointer :: diag => null(), cmor_diag => null(), z_remap_diag => null(), cmor_z_remap_diag => null() - integer :: primary_id, fms_id, remap_id + type(diag_type), pointer :: this_diag => null() + integer :: fms_id character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name, cm_string, msg MOM_missing_value = axes%diag_cs%missing_value if(present(missing_value)) MOM_missing_value = missing_value + register_diag_field_expand_cmor = .false. diag_cs => axes%diag_cs - primary_id = -1 - diag => null() - cmor_diag => null() - z_remap_diag => null() - cmor_z_remap_diag => null() ! Set up the 'primary' diagnostic, first get an underlying FMS id - fms_id = register_diag_field_expand_axes(module_name, field_name, axes, & - 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) + fms_id = register_diag_field_expand_axes(module_name, field_name, axes, 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) call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method) + this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND) then - ! If the diagnostic is needed then allocate and id and space - primary_id = get_new_diag_id(diag_cs) - call alloc_diag_with_id(primary_id, diag_cs, diag) - call assert(associated(diag), 'register_diag_field: diag allocation failed') - diag%fms_diag_id = fms_id - diag%debug_str = trim(field_name) - call set_diag_mask(diag, diag_cs, axes) - if (present(conversion)) diag%conversion_factor = conversion + ! If the diagnostic is needed obtain a diag_mediator ID (if needed) + if (primary_id == -1) primary_id = get_new_diag_id(diag_cs) + ! Create a new diag_type to store links in + call alloc_diag_with_id(primary_id, diag_cs, this_diag) + call assert(associated(this_diag), 'register_diag_field_expand_cmor: diag allocation failed') + ! Record FMS id, masks and conversion factor, in diag_type + this_diag%fms_diag_id = fms_id + this_diag%debug_str = trim(field_name) + call set_diag_mask(this_diag, diag_cs, axes) + this_diag%axes => axes + if (present(conversion)) this_diag%conversion_factor = conversion + register_diag_field_expand_cmor = .true. endif if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(diag), module_name, field_name, cm_string, & + call log_available_diag(associated(this_diag), module_name, field_name, cm_string, & msg, diag_CS, long_name, units, standard_name) endif @@ -1202,123 +1330,40 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Set up the CMOR variation of the native diagnostic if (present(cmor_field_name)) then fms_id = register_diag_field_expand_axes(module_name, cmor_field_name, axes, init_time, & - long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), & - missing_value=MOM_missing_value, range=range, mask_variant=mask_variant, & - standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, & - err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) + long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), & + missing_value=MOM_missing_value, range=range, mask_variant=mask_variant, & + standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, & + err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method) + this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND) then - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - ! This will add the cmore variation to the 'primary' diagnostic. - ! In the case where there is no primary, it will become the primary. - call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) - cmor_diag%fms_diag_id = fms_id - cmor_diag%debug_str = trim(cmor_field_name) - call set_diag_mask(cmor_diag, diag_cs, axes) - if (present(conversion)) cmor_diag%conversion_factor = conversion + ! If the diagnostic is needed obtain a diag_mediator ID (if needed) + if (primary_id == -1) primary_id = get_new_diag_id(diag_cs) + ! Create a new diag_type to store links in + call alloc_diag_with_id(primary_id, diag_cs, this_diag) + call assert(associated(this_diag), 'register_diag_field_expand_cmor: cmor_diag allocation failed') + ! Record FMS id, masks and conversion factor, in diag_type + this_diag%fms_diag_id = fms_id + this_diag%debug_str = trim(cmor_field_name) + call set_diag_mask(this_diag, diag_cs, axes) + this_diag%axes => axes + if (present(conversion)) this_diag%conversion_factor = conversion + register_diag_field_expand_cmor = .true. endif if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, cm_string, & + call log_available_diag(associated(this_diag), module_name, cmor_field_name, cm_string, & msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif endif - ! Remap to z vertical coordinate, note that only diagnostics on layers - ! (not interfaces) are supported, also B axes are not supported yet - if (is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) .and. axes%rank == 3) then - if (get_diag_field_id_fms(trim(module_name)//trim(diag_cs%z_remap_suffix), field_name) /= DIAG_FIELD_NOT_FOUND) 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 DIAG_REMAP_Z_GRID_DEF') - 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) - call set_diag_mask(z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(z_remap_diag, diag_cs, axes) - if (present(conversion)) z_remap_diag%conversion_factor = conversion - call assert(associated(z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_expand_axes(module_name//trim(diag_cs%z_remap_suffix), field_name, & - z_remap_diag%remap_axes, & - 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) - call attach_cell_methods(fms_id, z_remap_diag%remap_axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) - z_remap_diag%fms_diag_id = fms_id - z_remap_diag%debug_str = trim(field_name) - - if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. - elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. - else - diag_cs%do_z_remapping_on_T = .true. - endif - endif - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), field_name, & - cm_string, msg, diag_CS, long_name, units, standard_name) - endif - ! Remap to z vertical coordinate with CMOR names and attributes - if (present(cmor_field_name)) then - if (get_diag_field_id_fms(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name) /= DIAG_FIELD_NOT_FOUND) 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 DIAG_REMAP_Z_GRID_DEF') - endif - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - call alloc_diag_with_id(primary_id, diag_cs, cmor_z_remap_diag) - call set_diag_mask(cmor_z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(cmor_z_remap_diag, diag_cs, axes) - if (present(conversion)) cmor_z_remap_diag%conversion_factor = conversion - call assert(associated(cmor_z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_expand_axes(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & - cmor_z_remap_diag%remap_axes, & - init_time, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=trim(posted_cmor_standard_name), & - verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & - interp_method=interp_method, tile_count=tile_count) - call attach_cell_methods(fms_id, cmor_z_remap_diag%remap_axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) - cmor_z_remap_diag%fms_diag_id = fms_id - cmor_z_remap_diag%debug_str = trim(cmor_field_name) - - if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. - elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. - else - diag_cs%do_z_remapping_on_T = .true. - endif - endif - if (is_root_pe() .and. diag_CS%doc_unit > 0) then - msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(cmor_z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & - cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif - endif - endif - - register_diag_field = primary_id - -end function register_diag_field +end function register_diag_field_expand_cmor -!> Returns ID from register_diag_field_fms (the diag_manager routine) after expanding axes (axes-group) into handles -!! and conditionally adding an FMS area_id for cell_measures. +!> Returns an FMS id from register_diag_field_fms (the diag_manager routine) after expanding axes (axes-group) +!! into handles and conditionally adding an FMS area_id for cell_measures. integer function register_diag_field_expand_axes(module_name, field_name, axes, init_time, & long_name, units, missing_value, range, mask_variant, standard_name, & verbose, do_not_log, err_msg, interp_method, tile_count) @@ -2000,33 +2045,11 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s -subroutine set_diag_remap_axes(diag, diag_cs, axes) - ! Associate a remapping axes with the a diagnostic based on the axes info. - type(diag_ctrl), target, intent(in) :: diag_cs - type(diag_type), pointer, intent(inout) :: diag - type(axes_grp), intent(in) :: axes - - diag%remap_axes => null() - if (axes%rank .eq. 3) then - if ((axes%id .eq. diag_cs%axesTL%id)) then - diag%remap_axes => diag_cs%axesTZL - elseif(axes%id .eq. diag_cs%axesBL%id) then - diag%remap_axes => diag_cs%axesBZL - elseif(axes%id .eq. diag_cs%axesCuL%id ) then - diag%remap_axes => diag_cs%axesCuZL - elseif(axes%id .eq. diag_cs%axesCvL%id) then - diag%remap_axes => diag_cs%axesCvZL - endif - endif - -end subroutine set_diag_remap_axes - +!> Associates the mask pointers within diag with the appropriate mask based on the axes group. subroutine set_diag_mask(diag, diag_cs, axes) - ! Associate a mask with the a diagnostic based on the axes info. - - type(diag_ctrl), target, intent(in) :: diag_cs - type(diag_type), pointer, intent(inout) :: diag - type(axes_grp), intent(in) :: axes + type(diag_ctrl), target, intent(in) :: diag_cs !< Diag_mediator control structure + type(diag_type), pointer, intent(inout) :: diag !< This diag type + type(axes_grp), intent(in) :: axes !< Axes group diag%mask2d => null() diag%mask3d => null() @@ -2181,7 +2204,7 @@ subroutine initialize_diag_type(diag) diag%in_use = .false. diag%fms_diag_id = -1 - diag%remap_axes => null() + diag%axes => null() diag%mask2d => null() diag%mask3d => null() diag%next => null() From 54b5724cc5030a62a3c8fd9a1568ef1c97cd3c3d Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 6 Sep 2016 10:24:54 -0400 Subject: [PATCH 06/16] Re-factor of diag_mediator: removed functions is_u_axes, is_v_axes, ... - Removed functions is_u_axes(), is_v_axes(), is_B_axes(), is_layer() by storing logicals in the axes_group type. - Also added needs_remapping logical for non-native 3D axes groups to indicate to do remapping. - No answer changes. --- src/framework/MOM_diag_mediator.F90 | 206 +++++++++++----------------- 1 file changed, 82 insertions(+), 124 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 5dd28b9c25..595554ef34 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -98,8 +98,17 @@ module MOM_diag_mediator character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group includes vertical direction. ! For detecting position on the grid logical :: is_h_point = .false. !< If true, indicates that this axes group is for an h-point located field. - logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid. False for any other grid. Used for rank>2. - logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located fields that must be remapped to these axes. Used for rank>2. + logical :: is_q_point = .false. !< If true, indicates that this axes group is for a q-point located field. + logical :: is_u_point = .false. !< If true, indicates that this axes group is for a u-point located field. + logical :: is_v_point = .false. !< If true, indicates that this axes group is for a v-point located field. + logical :: is_layer = .false. !< If true, indicates that this axes group is for a layer vertically-located field. + logical :: is_interface = .false. !< If true, indicates that this axes group is for an interface vertically-located field. + logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid. False for any other + !! grid. Used for rank>2. + logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located field + !! that must be remapped to these axes. Used for rank>2. + logical :: needs_interpolation = .false. !< If true, indicates that this axes group is for a sampled interface-located field + !! that must be interpolated to these axes. Used for rank>2. ! ID's for cell_measures integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp. integer :: id_volume = -1 !< The diag_manager id for volume to be used for cell_measure of variables with this axes_grp. @@ -353,57 +362,64 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) 'Depth of interfaces', direction=-1) endif - ! Axes for z remapping + ! Axes for z layers call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & - is_h_point=.true., is_native=.false., needs_remapping=.true.) + is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & x_cell_method='point', y_cell_method='point', v_cell_method='mean', & - is_native=.false.) + is_q_point=.true., is_layer=.true., is_native=.false.) + !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for axesBZL call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL, & x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & - is_native=.false., needs_remapping=.true.) + is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL, & x_cell_method='mean', y_cell_method='point', v_cell_method='mean', & - is_native=.false., needs_remapping=.true.) + is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) ! Vertical axes for the interfaces and layers call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, & - v_cell_method='point') + v_cell_method='point', is_interface=.true.) call define_axes_group(diag_cs, (/ id_zL /), diag_cs%axesZL, & - v_cell_method='mean') + v_cell_method='mean', is_layer=.true.) ! Axis groupings for the model layers call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%axesTL, & x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & - is_h_point=.true.) + is_h_point=.true., is_layer=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%axesBL, & - x_cell_method='point', y_cell_method='point', v_cell_method='mean') + x_cell_method='point', y_cell_method='point', v_cell_method='mean', & + is_q_point=.true., is_layer=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%axesCuL, & - x_cell_method='point', y_cell_method='mean', v_cell_method='mean') + x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & + is_u_point=.true., is_layer=.true.) call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%axesCvL, & - x_cell_method='mean', y_cell_method='point', v_cell_method='mean') + x_cell_method='mean', y_cell_method='point', v_cell_method='mean', & + is_v_point=.true., is_layer=.true.) ! Axis groupings for the model interfaces call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%axesTi, & x_cell_method='mean', y_cell_method='mean', v_cell_method='point', & - is_h_point=.true.) + is_h_point=.true., is_interface=.true.) + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, & + x_cell_method='point', y_cell_method='point', v_cell_method='point', & + is_q_point=.true., is_interface=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%axesCui, & - x_cell_method='point', y_cell_method='mean', v_cell_method='point') + x_cell_method='point', y_cell_method='mean', v_cell_method='point', & + is_u_point=.true., is_interface=.true.) call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%axesCvi, & - x_cell_method='mean', y_cell_method='point', v_cell_method='point') - call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, & - x_cell_method='point', y_cell_method='point', v_cell_method='point') + x_cell_method='mean', y_cell_method='point', v_cell_method='point', & + is_v_point=.true., is_interface=.true.) ! Axis groupings for 2-D arrays call define_axes_group(diag_cs, (/ id_xh, id_yh /), diag_cs%axesT1, & x_cell_method='mean', y_cell_method='mean', is_h_point=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1, & - x_cell_method='point', y_cell_method='point') + x_cell_method='point', y_cell_method='point', is_q_point=.true.) call define_axes_group(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1, & - x_cell_method='point', y_cell_method='mean') + x_cell_method='point', y_cell_method='mean', is_u_point=.true.) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & - x_cell_method='mean', y_cell_method='point') + x_cell_method='mean', y_cell_method='point', is_v_point=.true.) end subroutine set_axes_info @@ -474,7 +490,10 @@ end function check_grid_def !> Defines a group of "axes" from list of handles subroutine define_axes_group(diag_cs, handles, axes, & - x_cell_method, y_cell_method, v_cell_method, is_h_point, is_native, needs_remapping) + x_cell_method, y_cell_method, v_cell_method, & + is_h_point, is_q_point, is_u_point, is_v_point, & + is_layer, is_interface, & + is_native, needs_remapping, needs_interpolation) type(diag_ctrl), target, intent(in) :: diag_cs !< Diagnostics control structure integer, dimension(:), intent(in) :: handles !< A list of 1D axis handles type(axes_grp), intent(out) :: axes !< The group of 1D axes @@ -482,9 +501,16 @@ subroutine define_axes_group(diag_cs, handles, axes, & character(len=*), optional, intent(in) :: y_cell_method !< A y-direction cell method used to construct the "cell_methods" attribute in CF convention character(len=*), optional, intent(in) :: v_cell_method !< A vertical direction cell method used to construct the "cell_methods" attribute in CF convention logical, optional, intent(in) :: is_h_point !< If true, indicates this axes group for h-point located fields + logical, optional, intent(in) :: is_q_point !< If true, indicates this axes group for q-point located fields + logical, optional, intent(in) :: is_u_point !< If true, indicates this axes group for u-point located fields + logical, optional, intent(in) :: is_v_point !< If true, indicates this axes group for v-point located fields + logical, optional, intent(in) :: is_layer !< If true, indicates that this axes group is for a layer vertically-located field. + logical, optional, intent(in) :: is_interface !< If true, indicates that this axes group is for an interface vertically-located field. logical, optional, intent(in) :: is_native !< If true, indicates that this axes group is for a native model grid. False for any other grid. - logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is for a intensive layer-located fields + logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is for a intensive layer-located field !! that must be remapped to these axes. Used for rank>2. + logical, optional, intent(in) :: needs_interpolation !< If true, indicates that this axes group is for a sampled interface-located field + !! that must be interpolated to these axes. Used for rank>2. ! Local variables integer :: n n = size(handles) @@ -516,8 +542,14 @@ subroutine define_axes_group(diag_cs, handles, axes, & axes%v_cell_method = '' endif if (present(is_h_point)) axes%is_h_point = is_h_point + if (present(is_q_point)) axes%is_q_point = is_q_point + if (present(is_u_point)) axes%is_u_point = is_u_point + if (present(is_v_point)) axes%is_v_point = is_v_point + if (present(is_layer)) axes%is_layer = is_layer + if (present(is_interface)) axes%is_interface = is_interface if (present(is_native)) axes%is_native = is_native if (present(needs_remapping)) axes%needs_remapping = needs_remapping + if (present(needs_interpolation)) axes%needs_interpolation = needs_interpolation end subroutine define_axes_group subroutine set_diag_mediator_grid(G, diag_cs) @@ -815,7 +847,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) nz_src = size(field, 3) nz_dest = diag_cs%nz_remap - if (is_u_axes(diag%axes, diag_cs)) then + if (diag%axes%is_u_point) then do j=diag_cs%G%jsc, diag_cs%G%jec do I=diag_cs%G%iscB, diag_cs%G%iecB if (associated(diag%mask3d)) then @@ -841,7 +873,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo enddo - elseif (is_v_axes(diag%axes, diag_cs)) then + elseif (diag%axes%is_v_point) then do J=diag_cs%G%jscB, diag_cs%G%jecB do i=diag_cs%G%isc, diag_cs%G%iec if (associated(diag%mask3d)) then @@ -867,7 +899,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo enddo - else + elseif (diag%axes%is_h_point) then do j=diag_cs%G%jsc, diag_cs%G%jec do i=diag_cs%G%isc, diag_cs%G%iec if (associated(diag%mask3d)) then @@ -892,6 +924,10 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo enddo enddo + else + call MOM_error(FATAL, & + "remap_diag_to_z: Q point remapping is not coded yet. "//& + " debug hint: "//diag%debug_str) endif end subroutine remap_diag_to_z @@ -1180,23 +1216,22 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion) - ! Register diagnostics remapped to z vertical coordinate, note that only diagnostics on layers - ! (not interfaces) are supported, also B axes are not supported yet - if (axes%rank == 3 .and. (.not. is_B_axes(axes, diag_cs))) then - if (is_layer_axes(axes, diag_cs)) then - remap_axes => null() - if (axes%rank .eq. 3) then - if ((axes%id .eq. diag_cs%axesTL%id)) then - remap_axes => diag_cs%axesTZL - elseif(axes%id .eq. diag_cs%axesBL%id) then - remap_axes => diag_cs%axesBZL - elseif(axes%id .eq. diag_cs%axesCuL%id ) then - remap_axes => diag_cs%axesCuZL - elseif(axes%id .eq. diag_cs%axesCvL%id) then - remap_axes => diag_cs%axesCvZL - endif + ! Register diagnostics remapped to z vertical coordinate + if (axes%rank == 3 .and. axes%is_layer) then + remap_axes => null() + if (axes%rank .eq. 3) then + if ((axes%id .eq. diag_cs%axesTL%id)) then + remap_axes => diag_cs%axesTZL + elseif(axes%id .eq. diag_cs%axesBL%id) then + remap_axes => diag_cs%axesBZL + elseif(axes%id .eq. diag_cs%axesCuL%id ) then + remap_axes => diag_cs%axesCuZL + elseif(axes%id .eq. diag_cs%axesCvL%id) then + remap_axes => diag_cs%axesCvZL endif - call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + endif + call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + if (remap_axes%needs_remapping) then active = register_diag_field_expand_cmor(dm_id, & module_name//trim(diag_cs%z_remap_suffix), field_name, remap_axes, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & @@ -1216,16 +1251,16 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & ! Record whether any remapping is needed at all if (active) then - if (is_u_axes(axes, diag_cs)) then + if (remap_axes%is_u_point) then diag_cs%do_z_remapping_on_u = .true. - elseif (is_v_axes(axes, diag_cs)) then + elseif (remap_axes%is_v_point) then diag_cs%do_z_remapping_on_v = .true. - else + elseif (remap_axes%is_h_point) then diag_cs%do_z_remapping_on_T = .true. endif endif - endif ! is_layer_axes - endif + endif ! remap_axes%needs_remapping + endif ! axes%rank == 3 register_diag_field = dm_id @@ -2090,83 +2125,6 @@ subroutine set_diag_mask(diag, diag_cs, axes) end subroutine set_diag_mask -function is_layer_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_layer_axes - - is_layer_axes = .false. - - if (axes%id == diag_cs%axesTZL%id .or. & - axes%id == diag_cs%axesBZL%id .or. & - axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCvZL%id .or. & - axes%id == diag_cs%axesZL%id .or. & - axes%id == diag_cs%axesTL%id .or. & - axes%id == diag_cs%axesBL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCvL%id) then - is_layer_axes = .true. - endif - -end function is_layer_axes - -function is_u_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_u_axes - - is_u_axes = .false. - - if (axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCui%id .or. & - axes%id == diag_cs%axesCu1%id) then - is_u_axes = .true. - endif - -end function is_u_axes - -function is_v_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_v_axes - - is_v_axes = .false. - - if (axes%id == diag_cs%axesCuZL%id .or. & - axes%id == diag_cs%axesCuL%id .or. & - axes%id == diag_cs%axesCui%id .or. & - axes%id == diag_cs%axesCu1%id) then - is_v_axes = .true. - endif - -end function is_v_axes - -function is_B_axes(axes, diag_cs) - - type(axes_grp), intent(in) :: axes - type(diag_ctrl), intent(in) :: diag_cs - - logical :: is_B_axes - - is_B_axes = .false. - - if (axes%id == diag_cs%axesBZL%id .or. & - axes%id == diag_cs%axesBL%id .or. & - axes%id == diag_cs%axesBi%id .or. & - axes%id == diag_cs%axesB1%id) then - is_B_axes = .true. - endif - -end function is_B_axes - !> Returns a new diagnostic id, it may be necessary to expand the diagnostics array. integer function get_new_diag_id(diag_cs) type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure From 11141a2efaf392252875a03709b0e93f654566eb Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 6 Sep 2016 16:12:20 -0400 Subject: [PATCH 07/16] Implemented z-coordinate diagnostics of interface variables - Interface variables cannot be remapped so our z-coordinate diagnostics have been limited to layer-centered fields. This implements interpolation of interface-centered fields to the interfaces of the target grid (same as used by the diagnostic remapping). - This added a new group of axes for the z-coordinate interfaces. - Note that the (horizontal) B-grid located diagnostics are not implemented (as they are not implemented in remapping code either). - Had to provide a work-around for the current existence of MOM_diag_to_Z which calls register_diag_field for 3D fields that should not be remapped. Todo: defaults for is_native should be flipped. - No answer changes. --- src/framework/MOM_diag_mediator.F90 | 243 ++++++++++++++++++++++------ 1 file changed, 190 insertions(+), 53 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 595554ef34..cd029c970c 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -31,6 +31,7 @@ module MOM_diag_mediator use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE use MOM_error_handler, only : MOM_error, FATAL, is_root_pe +use MOM_diag_vkernels, only : interpolate_column 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 : file_exists, field_exists, field_size @@ -39,7 +40,7 @@ module MOM_diag_mediator 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, initialize_remapping, dzFromH1H2 +use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution use MOM_regridding, only : build_zstar_column, set_regrid_params @@ -107,7 +108,7 @@ module MOM_diag_mediator !! grid. Used for rank>2. logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located field !! that must be remapped to these axes. Used for rank>2. - logical :: needs_interpolation = .false. !< If true, indicates that this axes group is for a sampled interface-located field + logical :: needs_interpolating = .false. !< If true, indicates that this axes group is for a sampled interface-located field !! that must be interpolated to these axes. Used for rank>2. ! ID's for cell_measures integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp. @@ -151,7 +152,8 @@ module MOM_diag_mediator type(axes_grp) :: axesBi, axesTi, axesCui, axesCvi type(axes_grp) :: axesB1, axesT1, axesCu1, axesCv1 type(axes_grp) :: axesZi, axesZL - type(axes_grp) :: axesTzi, axesTZL, axesBZL, axesCuZL, axesCvZL + type(axes_grp) :: axesTZL, axesBZL, axesCuZL, axesCvZL + type(axes_grp) :: axesTzi, axesBZi, axesCuZi, axesCvZi ! Mask arrays for diagnostics real, dimension(:,:), pointer :: mask2dT => null() @@ -362,21 +364,6 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) 'Depth of interfaces', direction=-1) endif - ! Axes for z layers - call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & - x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & - is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) - call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & - x_cell_method='point', y_cell_method='point', v_cell_method='mean', & - is_q_point=.true., is_layer=.true., is_native=.false.) - !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for axesBZL - call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL, & - x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & - is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) - call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL, & - x_cell_method='mean', y_cell_method='point', v_cell_method='mean', & - is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) - ! Vertical axes for the interfaces and layers call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, & v_cell_method='point', is_interface=.true.) @@ -421,6 +408,36 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & x_cell_method='mean', y_cell_method='point', is_v_point=.true.) + ! Axes for z layers + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & + x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & + is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & + x_cell_method='point', y_cell_method='point', v_cell_method='mean', & + is_q_point=.true., is_layer=.true., is_native=.false.) + !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for axesBZL + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL, & + x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & + is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL, & + x_cell_method='mean', y_cell_method='point', v_cell_method='mean', & + is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) + + ! Axes for z interfaces + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzi /), diag_cs%axesTZi, & + x_cell_method='mean', y_cell_method='mean', v_cell_method='point', & + is_h_point=.true., is_layer=.true., is_native=.false., needs_interpolating=.true.) + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzi /), diag_cs%axesBZi, & + x_cell_method='point', y_cell_method='point', v_cell_method='point', & + is_q_point=.true., is_layer=.true., is_native=.false.) + !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for axesBZL + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zzi /), diag_cs%axesCuZi, & + x_cell_method='point', y_cell_method='mean', v_cell_method='point', & + is_u_point=.true., is_layer=.true., is_native=.false., needs_interpolating=.true.) + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zzi /), diag_cs%axesCvZi, & + x_cell_method='mean', y_cell_method='point', v_cell_method='point', & + is_v_point=.true., is_layer=.true., is_native=.false., needs_interpolating=.true.) + end subroutine set_axes_info !> Attaches the id of cell areas to axes groupsfor use with cell_measures @@ -493,7 +510,7 @@ subroutine define_axes_group(diag_cs, handles, axes, & x_cell_method, y_cell_method, v_cell_method, & is_h_point, is_q_point, is_u_point, is_v_point, & is_layer, is_interface, & - is_native, needs_remapping, needs_interpolation) + is_native, needs_remapping, needs_interpolating) type(diag_ctrl), target, intent(in) :: diag_cs !< Diagnostics control structure integer, dimension(:), intent(in) :: handles !< A list of 1D axis handles type(axes_grp), intent(out) :: axes !< The group of 1D axes @@ -509,7 +526,7 @@ subroutine define_axes_group(diag_cs, handles, axes, & logical, optional, intent(in) :: is_native !< If true, indicates that this axes group is for a native model grid. False for any other grid. logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is for a intensive layer-located field !! that must be remapped to these axes. Used for rank>2. - logical, optional, intent(in) :: needs_interpolation !< If true, indicates that this axes group is for a sampled interface-located field + logical, optional, intent(in) :: needs_interpolating !< If true, indicates that this axes group is for a sampled interface-located field !! that must be interpolated to these axes. Used for rank>2. ! Local variables integer :: n @@ -549,7 +566,7 @@ subroutine define_axes_group(diag_cs, handles, axes, & if (present(is_interface)) axes%is_interface = is_interface if (present(is_native)) axes%is_native = is_native if (present(needs_remapping)) axes%needs_remapping = needs_remapping - if (present(needs_interpolation)) axes%needs_interpolation = needs_interpolation + if (present(needs_interpolating)) axes%needs_interpolating = needs_interpolating end subroutine define_axes_group subroutine set_diag_mediator_grid(G, diag_cs) @@ -795,7 +812,6 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) if (diag%axes%needs_remapping) 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 @@ -815,6 +831,26 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) deallocate(remapped_field) if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + elseif (diag%axes%needs_interpolating) then + ! Interpolate this field to another vertical coordinate. + if (present(mask)) then + call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") + endif + + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap+1)) + call vertically_interpolate_diag_field(diag_cs, diag, field, remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + 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) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + deallocate(remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) else call post_data_3d_low(diag, field, diag_cs, is_static, mask) endif @@ -824,15 +860,14 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) end subroutine post_data_3d -!> Remap diagnostic field to z-star vertical grid. +!> Remap diagnostic field to alternative vertical grid. !! \note This uses grids generated by diag_update_target_grids() subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure - real, dimension(:,:,:), intent(inout) :: remapped_field !< Field argument remapped to z star coordinate + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field argument remapped to alternative coordinate ! Local variables - real, dimension(diag_cs%nz_remap+1) :: dz real, dimension(diag_cs%nz_remap) :: h_dest real, dimension(size(diag_cs%h, 3)) :: h_src integer :: nz_src, nz_dest @@ -932,6 +967,96 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z +!> Vertically interpolate diagnostic field to alternative vertical grid. +!! \note This uses grids generated by diag_update_target_grids() +subroutine vertically_interpolate_diag_field(diag_cs, diag, field, interpolated_field) + type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure + type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped + real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate + ! Local variables + real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(diag_cs%h, 3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(diag_cs%remapping_initialized, & + 'vertically_interpolate_diag_field: Remmaping not initialized.') + call assert(size(field, 3) == size(diag_cs%h, 3)+1, & + 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.') + + interpolated_field(:,:,:) = diag_cs%missing_value + nz_src = size(diag_cs%h, 3) + nz_dest = diag_cs%nz_remap + + if (diag%axes%is_u_point) then + do j=diag_cs%G%jsc, diag_cs%G%jec + do I=diag_cs%G%iscB, diag_cs%G%iecB + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j,1)+diag%mask3d(i+1,j,1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_interpolate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i+1,j,:) ) + call interpolate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, diag_cs%missing_value, interpolated_field(I,j,:)) + enddo + enddo + elseif (diag%axes%is_v_point) then + do J=diag_cs%G%jscB, diag_cs%G%jecB + do i=diag_cs%G%isc, diag_cs%G%iec + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j,1)+diag%mask3d(i,j+1,1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_interpolate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i,j+1,:) ) + call interpolate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, diag_cs%missing_value, interpolated_field(i,J,:)) + enddo + enddo + elseif (diag%axes%is_h_point) then + do j=diag_cs%G%jsc, diag_cs%G%jec + do i=diag_cs%G%isc, diag_cs%G%iec + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j, 1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_interpolate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = diag_cs%h(i,j,:) + h_dest(:) = diag_cs%h_zoutput(i,j,:) + call interpolate_column(nz_src, h_src, field(i,j,:), & + nz_dest, h_dest, diag_cs%missing_value, interpolated_field(i,j,:)) + enddo + enddo + else + call MOM_error(FATAL, & + "vertically_interpolate_diag_field: Q point remapping is not coded yet. "//& + " debug hint: "//diag%debug_str) + endif + +end subroutine vertically_interpolate_diag_field + !> Build/update target vertical grids for diagnostic remapping. !! \note The target grids need to be updated whenever sea surface !! height changes. @@ -1217,7 +1342,7 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & conversion=conversion) ! Register diagnostics remapped to z vertical coordinate - if (axes%rank == 3 .and. axes%is_layer) then + if (axes%rank == 3) then ! .and. axes%is_layer) then remap_axes => null() if (axes%rank .eq. 3) then if ((axes%id .eq. diag_cs%axesTL%id)) then @@ -1228,38 +1353,50 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & remap_axes => diag_cs%axesCuZL elseif(axes%id .eq. diag_cs%axesCvL%id) then remap_axes => diag_cs%axesCvZL + elseif(axes%id .eq. diag_cs%axesTi%id) then + remap_axes => diag_cs%axesTZi + elseif(axes%id .eq. diag_cs%axesBi%id) then + remap_axes => diag_cs%axesBZi + elseif(axes%id .eq. diag_cs%axesCui%id ) then + remap_axes => diag_cs%axesCuZi + elseif(axes%id .eq. diag_cs%axesCvi%id) then + remap_axes => diag_cs%axesCvZi endif endif - call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') - if (remap_axes%needs_remapping) then - active = register_diag_field_expand_cmor(dm_id, & - module_name//trim(diag_cs%z_remap_suffix), field_name, remap_axes, & - 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, & - cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & - cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & - cell_methods=cell_methods, x_cell_method=x_cell_method, & - y_cell_method=y_cell_method, v_cell_method=v_cell_method, & - conversion=conversion) + ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will + ! always exist but in the mean-time we have to do this check: + !call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + if (associated(remap_axes)) then + if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then + active = register_diag_field_expand_cmor(dm_id, & + module_name//trim(diag_cs%z_remap_suffix), field_name, remap_axes, & + 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, & + cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & + cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & + cell_methods=cell_methods, x_cell_method=x_cell_method, & + y_cell_method=y_cell_method, v_cell_method=v_cell_method, & + conversion=conversion) - if (active .and. .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 DIAG_REMAP_Z_GRID_DEF') - endif + if (active .and. .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 DIAG_REMAP_Z_GRID_DEF') + endif - ! Record whether any remapping is needed at all - if (active) then - if (remap_axes%is_u_point) then - diag_cs%do_z_remapping_on_u = .true. - elseif (remap_axes%is_v_point) then - diag_cs%do_z_remapping_on_v = .true. - elseif (remap_axes%is_h_point) then - diag_cs%do_z_remapping_on_T = .true. + ! Record whether any remapping is needed at all + if (active) then + if (remap_axes%is_u_point) then + diag_cs%do_z_remapping_on_u = .true. + elseif (remap_axes%is_v_point) then + diag_cs%do_z_remapping_on_v = .true. + elseif (remap_axes%is_h_point) then + diag_cs%do_z_remapping_on_T = .true. + endif endif - endif - endif ! remap_axes%needs_remapping + endif ! remap_axes%needs_remapping + endif ! associated(remap_axes) endif ! axes%rank == 3 register_diag_field = dm_id From bfd21a4e80ded4c61c94e5e5108694dce4c604ee Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 7 Sep 2016 11:51:37 -0400 Subject: [PATCH 08/16] Implemented remapped diagnostics of vertically extrinsic fields - Some layer data is already vertically integrated (e.g. hT, uh, vh, ...) and cannot be remapped as if it were a concentration (intrinsic variable). Setting optional argument v_extrinsic=.true. to register_diag_field() now indicates that the field should be re-gridded as an extrinsic field. - No answer changes. --- src/framework/MOM_diag_mediator.F90 | 162 +++++++++++++++++++++++++--- 1 file changed, 146 insertions(+), 16 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index cd029c970c..accabf5a28 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -31,7 +31,7 @@ module MOM_diag_mediator use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE use MOM_error_handler, only : MOM_error, FATAL, is_root_pe -use MOM_diag_vkernels, only : interpolate_column +use MOM_diag_vkernels, only : interpolate_column, reintegrate_column 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 : file_exists, field_exists, field_size @@ -130,6 +130,7 @@ module MOM_diag_mediator real, pointer, dimension(:,:,:) :: mask3d => null() type(diag_type), pointer :: next => null() !< Pointer to the next diag. real :: conversion_factor = 0. !< A factor to multiply data by before posting to FMS, if non-zero. + logical :: v_extrinsic = .false. !< True for vertically extrinsic fields (vertically integrated). False for intrinsic (concentrations). end type diag_type !> The following data type a list of diagnostic fields an their variants, @@ -810,7 +811,27 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) do while (associated(diag)) call assert(associated(diag%axes), 'post_data_3d: axes is not associated') - if (diag%axes%needs_remapping) then + if (diag%v_extrinsic .and. .not.diag%axes%is_native) then + ! The field is vertically integrated and needs to be re-gridded + if (present(mask)) then + call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") + endif + + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call vertically_reintegrate_diag_field(diag_cs, diag, field, remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + 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) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + deallocate(remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + elseif (diag%axes%needs_remapping) then ! Remap this field to another vertical coordinate. if (present(mask)) then call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") @@ -967,6 +988,96 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) end subroutine remap_diag_to_z +!> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. +!! \note This uses grids generated by diag_update_target_grids() +subroutine vertically_reintegrate_diag_field(diag_cs, diag, field, reintegrated_field) + type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure + type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped + real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate + ! Local variables + real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(diag_cs%h, 3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(diag_cs%remapping_initialized, & + 'vertically_reintegrate_diag_field: Remmaping not initialized.') + call assert(size(field, 3) == size(diag_cs%h, 3), & + 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.') + + reintegrated_field(:,:,:) = diag_cs%missing_value + nz_src = size(diag_cs%h, 3) + nz_dest = diag_cs%nz_remap + + if (diag%axes%is_u_point) then + do j=diag_cs%G%jsc, diag_cs%G%jec + do I=diag_cs%G%iscB, diag_cs%G%iecB + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j,1)+diag%mask3d(i+1,j,1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_reintegrate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i+1,j,:) ) + call reintegrate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, diag_cs%missing_value, reintegrated_field(I,j,:)) + enddo + enddo + elseif (diag%axes%is_v_point) then + do J=diag_cs%G%jscB, diag_cs%G%jecB + do i=diag_cs%G%isc, diag_cs%G%iec + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j,1)+diag%mask3d(i,j+1,1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_reintegrate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i,j+1,:) ) + call reintegrate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, diag_cs%missing_value, reintegrated_field(i,J,:)) + enddo + enddo + elseif (diag%axes%is_h_point) then + do j=diag_cs%G%jsc, diag_cs%G%jec + do i=diag_cs%G%isc, diag_cs%G%iec + if (associated(diag%mask3d)) then + if (diag%mask3d(i,j, 1) == 0.) cycle + endif +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "vertically_reintegrate_diag_field: H has changed since remapping grids were updated."//& + " diag debug hint: "//diag%debug_str) + enddo +#endif + h_src(:) = diag_cs%h(i,j,:) + h_dest(:) = diag_cs%h_zoutput(i,j,:) + call reintegrate_column(nz_src, h_src, field(i,j,:), & + nz_dest, h_dest, diag_cs%missing_value, reintegrated_field(i,j,:)) + enddo + enddo + else + call MOM_error(FATAL, & + "vertically_reintegrate_diag_field: Q point remapping is not coded yet. "//& + " debug hint: "//diag%debug_str) + endif + +end subroutine vertically_reintegrate_diag_field + !> Vertically interpolate diagnostic field to alternative vertical grid. !! \note This uses grids generated by diag_update_target_grids() subroutine vertically_interpolate_diag_field(diag_cs, diag, field, interpolated_field) @@ -1288,7 +1399,7 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & long_name, units, missing_value, range, mask_variant, standard_name, & verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & - x_cell_method, y_cell_method, v_cell_method, conversion) + x_cell_method, y_cell_method, v_cell_method, conversion, v_extrinsic) character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model" character(len=*), intent(in) :: field_name !< Name of the diagnostic field type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes for this field @@ -1315,6 +1426,7 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction. Use '' have no method. character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction. Use '' have no method. real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file + logical, optional, intent(in) :: v_extrinsic !< True for vertically extrinsic fields (vertically integrated). Default/absent for intrinsic. ! Local variables real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs @@ -1339,7 +1451,7 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & - conversion=conversion) + conversion=conversion, v_extrinsic=v_extrinsic) ! Register diagnostics remapped to z vertical coordinate if (axes%rank == 3) then ! .and. axes%is_layer) then @@ -1378,7 +1490,7 @@ integer function register_diag_field(module_name, field_name, axes, init_time, & cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & - conversion=conversion) + conversion=conversion, v_extrinsic=v_extrinsic) if (active .and. .not. allocated(diag_cs%zi_remap)) then call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & @@ -1409,7 +1521,7 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ long_name, units, missing_value, range, mask_variant, standard_name, & verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & - x_cell_method, y_cell_method, v_cell_method, conversion) + x_cell_method, y_cell_method, v_cell_method, conversion, v_extrinsic) integer, intent(inout) :: primary_id !< The diag_mediator ID for this diagnostic group character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model" character(len=*), intent(in) :: field_name !< Name of the diagnostic field @@ -1437,6 +1549,7 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction. Use '' have no method. character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction. Use '' have no method. real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file + logical, optional, intent(in) :: v_extrinsic !< True for vertically extrinsic fields (vertically integrated). Default/absent for intrinsic. ! Local variables real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs @@ -1457,7 +1570,8 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_id, axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) + cell_methods, x_cell_method, y_cell_method, v_cell_method, & + v_extrinsic=v_extrinsic) this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND) then ! If the diagnostic is needed obtain a diag_mediator ID (if needed) @@ -1470,6 +1584,7 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ this_diag%debug_str = trim(field_name) call set_diag_mask(this_diag, diag_cs, axes) this_diag%axes => axes + if (present(v_extrinsic)) this_diag%v_extrinsic = v_extrinsic if (present(conversion)) this_diag%conversion_factor = conversion register_diag_field_expand_cmor = .true. endif @@ -1507,7 +1622,8 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, & err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_id, axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method) + cell_methods, x_cell_method, y_cell_method, v_cell_method, & + v_extrinsic=v_extrinsic) this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND) then ! If the diagnostic is needed obtain a diag_mediator ID (if needed) @@ -1520,6 +1636,7 @@ logical function register_diag_field_expand_cmor(primary_id, module_name, field_ this_diag%debug_str = trim(cmor_field_name) call set_diag_mask(this_diag, diag_cs, axes) this_diag%axes => axes + if (present(v_extrinsic)) this_diag%v_extrinsic = v_extrinsic if (present(conversion)) this_diag%conversion_factor = conversion register_diag_field_expand_cmor = .true. endif @@ -1598,7 +1715,8 @@ integer function register_diag_field_expand_axes(module_name, field_name, axes, end function register_diag_field_expand_axes !> Attaches "cell_methods" attribute to a variable based on defaults for axes_grp or optional arguments. -subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y_cell_method, v_cell_method) +subroutine attach_cell_methods(id, axes, ostring, cell_methods, & + x_cell_method, y_cell_method, v_cell_method, v_extrinsic) integer, intent(in) :: id !< Handle to diagnostic type(axes_grp), intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes for this field character(len=*), intent(out) :: ostring !< The cell_methods strings that would appear in the file @@ -1608,12 +1726,14 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction. Use '' have no method. character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction. Use '' have no method. character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction. Use '' have no method. + logical, optional, intent(in) :: v_extrinsic !< True for vertically extrinsic fields (vertically integrated). Default/absent for intrinsic. ! Local variables character(len=9) :: axis_name ostring = '' if (present(cell_methods)) then - if (present(x_cell_method) .or. present(y_cell_method) .or. present(v_cell_method)) then + if (present(x_cell_method) .or. present(y_cell_method) .or. present(v_cell_method) & + .or. present(v_extrinsic)) then call MOM_error(FATAL, "attach_cell_methods: " // & 'Individual direction cell method was specified along with a "cell_methods" string.') endif @@ -1625,39 +1745,49 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y if (present(x_cell_method)) then if (len(trim(x_cell_method))>0) then call get_diag_axis_name(axes%handles(1), axis_name) - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(x_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(x_cell_method) endif else if (len(trim(axes%x_cell_method))>0) then call get_diag_axis_name(axes%handles(1), axis_name) - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(axes%x_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%x_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%x_cell_method) endif endif if (present(y_cell_method)) then if (len(trim(y_cell_method))>0) then call get_diag_axis_name(axes%handles(2), axis_name) - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(y_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(y_cell_method) endif else if (len(trim(axes%y_cell_method))>0) then call get_diag_axis_name(axes%handles(2), axis_name) - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(axes%y_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%y_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%y_cell_method) endif endif if (present(v_cell_method)) then + if (present(v_extrinsic)) call MOM_error(FATAL, "attach_cell_methods: " // & + 'Vertical cell method was specified along with the vertically extrinsic flag.') if (len(trim(v_cell_method))>0) then if (axes%rank==1) then call get_diag_axis_name(axes%handles(1), axis_name) elseif (axes%rank==3) then call get_diag_axis_name(axes%handles(3), axis_name) endif - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(v_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(v_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(v_cell_method) endif + elseif (present(v_extrinsic)) then + if (axes%rank==1) then + call get_diag_axis_name(axes%handles(1), axis_name) + elseif (axes%rank==3) then + call get_diag_axis_name(axes%handles(3), axis_name) + endif + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':sum') + ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':sum' else if (len(trim(axes%v_cell_method))>0) then if (axes%rank==1) then @@ -1665,7 +1795,7 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y elseif (axes%rank==3) then call get_diag_axis_name(axes%handles(3), axis_name) endif - call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//': '//trim(axes%v_cell_method)) + call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%v_cell_method)) ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%v_cell_method) endif endif From 3a7efe50a04d298f88967c895cf1565497169f6a Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 7 Sep 2016 11:54:13 -0400 Subject: [PATCH 09/16] Labeled some diagnostics as vertically integrated by layer - 3D diagnostics h, uh, vh have been labeled as vertically extrinsic (meaning integrated over the layer) so that remapping does not treat the fields as concentrations. - Changes availabl_diags. - No answer changes. --- src/core/MOM.F90 | 10 +++++----- src/core/MOM_dynamics_split_RK2.F90 | 12 ++++++------ src/core/MOM_dynamics_unsplit.F90 | 4 ++-- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 ++-- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4baee857da..aa25839912 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2216,7 +2216,7 @@ subroutine register_diags(Time, G, GV, CS, ADp) 'Meridional velocity', 'meter second-1', cmor_field_name='vo', cmor_units='m s-1', & cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity') CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, & - 'Layer Thickness', thickness_units, v_cell_method='sum') + 'Layer Thickness', thickness_units, v_extrinsic=.true.) CS%id_volo = register_scalar_field('ocean_model', 'volo', Time, diag,& long_name='Total volume of liquid ocean', units='m3', & @@ -2366,7 +2366,7 @@ subroutine register_diags(Time, G, GV, CS, ADp) CS%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, Time, & 'Meridional velocity before diabatic forcing', 'meter second-1') CS%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, Time, & - 'Layer Thickness before diabatic forcing', thickness_units, v_cell_method='sum') + 'Layer Thickness before diabatic forcing', thickness_units, v_extrinsic=.true.) CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, & 'Interface Heights before diabatic forcing', 'meter') if (CS%diabatic_first .and. (.not. CS%adiabatic)) then @@ -2375,7 +2375,7 @@ subroutine register_diags(Time, G, GV, CS, ADp) CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, & 'Meridional velocity before remapping', 'meter second-1') CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, & - 'Layer Thickness before remapping', thickness_units, v_cell_method='sum') + 'Layer Thickness before remapping', thickness_units, v_extrinsic=.true.) CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, & 'Temperature before remapping', 'degC') CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, & @@ -2428,7 +2428,7 @@ subroutine register_diags_TS_tendency(Time, G, CS) cmor_field_name="opottemptend", cmor_units="W m-2", & cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content", & cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content", & - v_cell_method='sum') + v_extrinsic=.true.) CS%id_Th_tendency_2d = register_diag_field('ocean_model', 'Th_tendency_2d', diag%axesT1, Time, & 'Vertical sum of net time tendency for heat', 'W/m2', & cmor_field_name="opottemptend_2d", cmor_units="W m-2", & @@ -2468,7 +2468,7 @@ subroutine register_diags_TS_tendency(Time, G, CS) cmor_field_name="osalttend", cmor_units="kg m-2 s-1", & cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content", & cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content", & - v_cell_method='sum') + v_extrinsic=.true.) CS%id_Sh_tendency_2d = register_diag_field('ocean_model', 'Sh_tendency_2d', diag%axesT1, Time, & 'Vertical sum of net time tendency for salt', 'kg/(m2 * s)', & cmor_field_name="osalttend_2d", cmor_units="kg m-2 s-1", & diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b8bcdbfe67..219a7c361e 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1194,25 +1194,25 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & - 'Zonal Thickness Flux', flux_units) + 'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extrinsic=.true.) CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, & - 'Meridional Thickness Flux', flux_units) + 'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extrinsic=.true.) CS%id_umo = register_diag_field('ocean_model', 'umo', & diag%axesCuL, Time,'Zonal Mass Transport (including SGS param)', 'kg/s',& cmor_standard_name='ocean_mass_x_transport', cmor_long_name='Ocean Mass X Transport', & - conversion=GV%H_to_kg_m2) + conversion=GV%H_to_kg_m2, y_cell_method='sum', v_extrinsic=.true.) CS%id_vmo = register_diag_field('ocean_model', 'vmo', & diag%axesCvL, Time,'Meridional Mass Transport (including SGS param)', 'kg/s',& cmor_standard_name='ocean_mass_y_transport', cmor_long_name='Ocean Mass Y Transport', & - conversion=GV%H_to_kg_m2) + conversion=GV%H_to_kg_m2, x_cell_method='sum', v_extrinsic=.true.) CS%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', & diag%axesCu1, Time,'Zonal Mass Transport (including SGS param) Vertical Sum', 'kg/s',& cmor_standard_name='ocean_mass_x_transport_vertical_sum', & - cmor_long_name='Ocean Mass X Transport Vertical Sum') + cmor_long_name='Ocean Mass X Transport Vertical Sum', y_cell_method='sum') CS%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', & diag%axesCv1, Time,'Meridional Mass Transport (including SGS param) Vertical Sum', 'kg/s',& cmor_standard_name='ocean_mass_y_transport_vertical_sum', & - cmor_long_name='Ocean Mass Y Transport Vertical Sum') + cmor_long_name='Ocean Mass Y Transport Vertical Sum', x_cell_method='sum') CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, & 'Zonal Coriolis and Advective Acceleration', 'meter second-2') diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index fc392cbb88..ab6290c6b1 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -682,9 +682,9 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, & flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & - 'Zonal Thickness Flux', flux_units) + 'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extrinsic=.true.) CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, & - 'Meridional Thickness Flux', flux_units) + 'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extrinsic=.true.) CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, & 'Zonal Coriolis and Advective Acceleration', 'meter second-2') CS%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, Time, & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 8e228a9189..e750a68643 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -646,9 +646,9 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, param_file, diag, CS flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & - 'Zonal Thickness Flux', flux_units) + 'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extrinsic=.true.) CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, & - 'Meridional Thickness Flux', flux_units) + 'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extrinsic=.true.) CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, & 'Zonal Coriolis and Advective Acceleration', 'meter second-2') CS%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, Time, & From ccdc9b617e9692aee6fc0d3c45b1f95bb56cfc1f Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 8 Sep 2016 12:23:32 -0400 Subject: [PATCH 10/16] Doxygenized assert() - New subroutines and functions should be doxygenized at inception. --- src/framework/MOM_error_handler.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/framework/MOM_error_handler.F90 b/src/framework/MOM_error_handler.F90 index cc1b305729..3194ced107 100644 --- a/src/framework/MOM_error_handler.F90 +++ b/src/framework/MOM_error_handler.F90 @@ -173,10 +173,10 @@ subroutine callTree_waypoint(mesg,n) endif end subroutine callTree_waypoint +!> Issues a FATAL error if the assertion fails, i.e. the first argument is false. subroutine assert(logical_arg, msg) - - logical, intent(in) :: logical_arg - character(len=*), intent(in) :: msg + logical, intent(in) :: logical_arg !< If false causes a FATAL error + character(len=*), intent(in) :: msg !< Message to issue in case of failed assertion if (.not. logical_arg) then call MOM_error(FATAL, msg) From b77344b0eab3cb01f9f3be781a56d3a84cd17930 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 8 Sep 2016 12:43:02 -0400 Subject: [PATCH 11/16] Clarification: renamed vector of diag_remap to diag_remap_cs - The vector diag_remap was a vector of control structures. To help me understand the new code I'm renaming by adding _cs. - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master. --- src/framework/MOM_diag_mediator.F90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 2cf1f9ac9d..8993330c3f 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -163,7 +163,7 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations real :: missing_value = 1.0e+20 - type(diag_remap_ctrl), dimension(REGRIDDING_NUM_TYPES) :: diag_remaps + type(diag_remap_ctrl), dimension(REGRIDDING_NUM_TYPES) :: diag_remap_cs ! Pointer to H, G and T&S needed for remapping real, dimension(:,:,:), pointer :: h => null() @@ -274,8 +274,8 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & x_cell_method='mean', y_cell_method='point', is_h_point=.false.) - do i=1, size(diag_cs%diag_remaps) - call diag_remap_set_vertical_axes(diag_cs%diag_remaps(i), G, GV, param_file) + do i=1, size(diag_cs%diag_remap_cs) + call diag_remap_set_vertical_axes(diag_cs%diag_remap_cs(i), G, GV, param_file) enddo end subroutine set_axes_info @@ -612,7 +612,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif - if (.not. diag_remap_axes_setup_done(diag_cs%diag_remaps(diag%vertical_coord))) then + if (.not. diag_remap_axes_setup_done(diag_cs%diag_remap_cs(diag%vertical_coord))) then call MOM_error(FATAL,"post_data_3d: attempt to do diagnostic vertical "// & "remapping but vertical axes not configured. See option "// & "DIAG_REMAP_"//vertical_coord_strings(diag%vertical_coord)//"_GRID_DEF") @@ -620,7 +620,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) - nz = diag_remap_get_nz(diag_cs%diag_remaps(diag%vertical_coord)) + nz = diag_remap_get_nz(diag_cs%diag_remap_cs(diag%vertical_coord)) allocate(remapped_field(size(field, 1), size(field, 2), nz)) #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) @@ -637,7 +637,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) enddo #endif - call diag_remap_do_remap(diag_cs%diag_remaps(diag%vertical_coord), & + call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%vertical_coord), & diag_cs%G, diag_cs%h, diag%axes, diag%mask3d, & diag_cs%missing_value, field, remapped_field) @@ -897,7 +897,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Register for vertical remapping if (axes%rank == 3 .and. & vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then - call diag_remap_set_diag_axes(diag_cs%diag_remaps(i), diag%axes) + call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes) endif fms_id = register_diag_field_expand_axes(new_module_name, field_name, diag%axes, & @@ -956,7 +956,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Register for vertical remapping if (axes%rank == 3 .and. & vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then - call diag_remap_set_diag_axes(diag_cs%diag_remaps(i), diag%axes) + call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes) endif fms_id = register_diag_field_expand_axes(new_module_name, cmor_field_name, diag%axes, & @@ -1495,8 +1495,8 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) ! Initialise vertical coordinates and test assumptions needed for diagnostic ! remapping. - do i=1, size(diag_cs%diag_remaps) - call diag_remap_init(diag_cs%diag_remaps(i), vertical_coords(i)) + do i=1, size(diag_cs%diag_remap_cs) + call diag_remap_init(diag_cs%diag_remap_cs(i), vertical_coords(i)) enddo ! Keep pointers grid, h, T, S needed diagnostic remapping @@ -1585,8 +1585,8 @@ subroutine diag_update_remap_grids(diag_cs) if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) - do i=1, size(diag_cs%diag_remaps) - call diag_remap_update(diag_cs%diag_remaps(i), & + do i=1, size(diag_cs%diag_remap_cs) + call diag_remap_update(diag_cs%diag_remap_cs(i), & diag_cs%G, diag_cs%h, diag_cs%T, diag_cs%S, & diag_cs%eqn_of_state) enddo @@ -1663,8 +1663,8 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) deallocate(diag_cs%diags) - do i=1, size(diag_cs%diag_remaps) - call diag_remap_end(diag_cs%diag_remaps(i)) + do i=1, size(diag_cs%diag_remap_cs) + call diag_remap_end(diag_cs%diag_remap_cs(i)) enddo deallocate(diag_cs%mask3dTL) From 4804926aad757009429bacdb63f7ae6360e9eae9 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 8 Sep 2016 13:21:50 -0400 Subject: [PATCH 12/16] Added axes_groups for each diagnostic coordinate - For each diag_remap_ctrl element there are now also eight axes_groups just as there are eight native axes_groups. - Not being used yet! - Still doesn't compile after a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master. --- src/framework/MOM_diag_mediator.F90 | 43 +++++++++++++++++++++++++---- src/framework/MOM_diag_remap.F90 | 12 ++++++++ 2 files changed, 49 insertions(+), 6 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 8993330c3f..d3db2fb178 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -43,7 +43,7 @@ module MOM_diag_mediator use MOM_diag_remap, only : diag_remap_update, diag_remap_set_diag_axes use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap use MOM_diag_remap, only : diag_remap_set_vertical_axes, diag_remap_get_nz -use MOM_diag_remap, only : diag_remap_axes_setup_done +use MOM_diag_remap, only : diag_remap_axes_setup_done, diag_remap_get_vertical_ids use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE, REGRIDDING_NUM_TYPES use regrid_consts, only : vertical_coords, vertical_coord_strings @@ -163,8 +163,13 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations real :: missing_value = 1.0e+20 + !> Control structure for each possible coordinate type(diag_remap_ctrl), dimension(REGRIDDING_NUM_TYPES) :: diag_remap_cs + !> Axes groups for each possible coordinate (these will all be 3D groups) + type(axes_grp), dimension(REGRIDDING_NUM_TYPES) :: remap_axesTL, remap_axesBL, remap_axesCuL, remap_axesCvL + type(axes_grp), dimension(REGRIDDING_NUM_TYPES) :: remap_axesTi, remap_axesBi, remap_axesCui, remap_axesCvi + ! Pointer to H, G and T&S needed for remapping real, dimension(:,:,:), pointer :: h => null() real, dimension(:,:,:), pointer :: T => null() @@ -275,35 +280,61 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) x_cell_method='mean', y_cell_method='point', is_h_point=.false.) do i=1, size(diag_cs%diag_remap_cs) + ! For each possible diagnostic coordinate call diag_remap_set_vertical_axes(diag_cs%diag_remap_cs(i), G, GV, param_file) + ! This fetches the 1D-axis id for layers and interfaces and overwrite id_zl and id_zi from above + call diag_remap_get_vertical_ids(diag_cs%diag_remap_cs(i), id_zl, id_zi) + ! Axis groupings for the model layers + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%remap_axesTL(i), & + x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', is_h_point=.true.) + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%remap_axesBL(i), & + x_cell_method='point', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%remap_axesCuL(i), & + x_cell_method='point', y_cell_method='mean', v_cell_method='mean', is_h_point=.false.) + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%remap_axesCvL(i), & + x_cell_method='mean', y_cell_method='point', v_cell_method='mean', is_h_point=.false.) + ! Axis groupings for the model interfaces + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), & + x_cell_method='mean', y_cell_method='mean', v_cell_method='point', is_h_point=.true.) + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%remap_axesCui(i), & + x_cell_method='point', y_cell_method='mean', v_cell_method='point', is_h_point=.false.) + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%remap_axesCvi(i), & + x_cell_method='mean', y_cell_method='point', v_cell_method='point', is_h_point=.false.) + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), & + x_cell_method='point', y_cell_method='point', v_cell_method='point', is_h_point=.false.) enddo end subroutine set_axes_info -!> Attaches the id of cell areas to axes groupsfor use with cell_measures +!> Attaches the id of cell areas to axes groups for use with cell_measures subroutine diag_register_area_ids(diag_cs, id_area_t, id_area_q) type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure integer, optional, intent(in) :: id_area_t !< Diag_mediator id for area of h-cells integer, optional, intent(in) :: id_area_q !< Diag_mediator id for area of q-cells ! Local variables - integer :: fms_id + integer :: fms_id, i if (present(id_area_t)) then fms_id = diag_cs%diags(id_area_t)%fms_diag_id diag_cs%axesT1%id_area = fms_id diag_cs%axesTi%id_area = fms_id diag_cs%axesTL%id_area = fms_id - diag_cs%axesTZL%id_area = fms_id + do i=1, size(diag_cs%diag_remap_cs) + diag_cs%remap_axesTL(i)%id_area = fms_id + ! Note to AJA: why am I not doing TZi too? + enddo endif if (present(id_area_q)) then fms_id = diag_cs%diags(id_area_q)%fms_diag_id diag_cs%axesB1%id_area = fms_id diag_cs%axesBi%id_area = fms_id diag_cs%axesBL%id_area = fms_id - diag_cs%axesBZL%id_area = fms_id + do i=1, size(diag_cs%diag_remap_cs) + diag_cs%remap_axesBL(i)%id_area = fms_id + enddo endif end subroutine diag_register_area_ids -!> Attaches the id of cell volumes to axes groupsfor use with cell_measures +!> Attaches the id of cell volumes to axes groups for use with cell_measures subroutine diag_register_volume_ids(diag_cs, id_vol_t) type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure integer, optional, intent(in) :: id_vol_t !< Diag_manager id for volume of h-cells diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index a4f586bc12..da8db145b8 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -52,6 +52,7 @@ module MOM_diag_remap public diag_remap_ctrl, diag_remap_set_diag_axes public diag_remap_init, diag_remap_end, diag_remap_update, diag_remap_do_remap public diag_remap_set_vertical_axes, diag_remap_axes_setup_done, diag_remap_get_nz +public diag_remap_get_vertical_ids type :: diag_remap_ctrl ! Whether remappping initialized @@ -319,6 +320,17 @@ subroutine diag_remap_set_diag_axes(remap_cs, axes) end subroutine +!> Get layer and interface axes ids for this coordinate +subroutine diag_remap_get_vertical_ids(remap_cs, id_layer, id_interface) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + integer, intent(out) :: id_layer !< 1D-axes id for layer points + integer, intent(out) :: id_interface !< 1D-axes id for interface points + + id_layer = remap_cs%layer_axes_id + id_interface = remap_cs%interface_axes_id + +end subroutine diag_remap_get_vertical_ids + function diag_remap_get_nz(remap_cs) type(diag_remap_ctrl), intent(in) :: remap_cs integer :: diag_remap_get_nz From 9338ee552362053b541dd8a134d366865d7f5c5e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 8 Sep 2016 13:28:17 -0400 Subject: [PATCH 13/16] Clarification: renamed diag_type%axes to diag_type%axes_ids - The coexistence of axes_groups and axes ids is confusing. I've renamed %axes to $axes_ids to help me understand the code. - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master. --- src/framework/MOM_diag_mediator.F90 | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index d3db2fb178..4c2db6fa88 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -111,8 +111,7 @@ module MOM_diag_mediator character(16) :: debug_str = '' !< For FATAL errors and debugging. type(axes_grp), pointer :: remap_axes => null() integer :: vertical_coord - integer, dimension(:), allocatable :: axes - + integer, dimension(:), allocatable :: axes_ids !< 1D-axis ids real, pointer, dimension(:,:) :: mask2d => null() real, pointer, dimension(:,:,:) :: mask3d => null() type(diag_type), pointer :: next => null() !< Pointer to the next diag. @@ -669,7 +668,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) #endif call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%vertical_coord), & - diag_cs%G, diag_cs%h, diag%axes, diag%mask3d, & + diag_cs%G, diag_cs%h, diag%axes_ids, diag%mask3d, & diag_cs%missing_value, field, remapped_field) if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) @@ -928,10 +927,10 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Register for vertical remapping if (axes%rank == 3 .and. & vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then - call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes) + call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes_ids) endif - fms_id = register_diag_field_expand_axes(new_module_name, field_name, diag%axes, & + fms_id = register_diag_field_expand_axes(new_module_name, field_name, diag%axes_ids, & 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, & @@ -987,10 +986,10 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Register for vertical remapping if (axes%rank == 3 .and. & vertical_coords(i) /= coordinateMode(DEFAULT_COORDINATE_MODE)) then - call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes) + call diag_remap_set_diag_axes(diag_cs%diag_remap_cs(i), diag%axes_ids) endif - fms_id = register_diag_field_expand_axes(new_module_name, cmor_field_name, diag%axes, & + fms_id = register_diag_field_expand_axes(new_module_name, cmor_field_name, diag%axes_ids,& init_time, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), & missing_value=MOM_missing_value, range=range, mask_variant=mask_variant, & standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log,& @@ -1749,8 +1748,8 @@ subroutine set_diag_mask_and_axes(diag, diag_cs, axes) ! Local variables integer :: i - allocate(diag%axes(size(axes%handles))) - diag%axes(:) = axes%handles(:) + allocate(diag%axes_ids(size(axes%handles))) + diag%axes_ids(:) = axes%handles(:) diag%mask2d => null() diag%mask3d => null() From d393a082af57454f95a9dc0577c7eeec4cdc9a26 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 8 Sep 2016 13:38:17 -0400 Subject: [PATCH 14/16] Re-introduced name "Z*" for ZSTAR coordinate mode - This keeps the new name "ZSTAR" but also allows use of the old coordinate name "Z*". - Note this code does not compile due to a broken merge of nicjhan-334-rho-sigma-diagnostic-regridding on dev/master. --- src/ALE/MOM_ALE.F90 | 2 +- src/ALE/regrid_consts.F90 | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 57434dcc1a..5379e6f927 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -952,7 +952,7 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz ) call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, & compressibility_fraction=compress_fraction ) - if (coordMode(1:5) == 'ZSTAR') then + if (coordinateMode(coordMode) == REGRIDDING_ZSTAR) then call get_param(param_file, mod, "ZSTAR_RIGID_SURFACE_THRESHOLD", height_of_rigid_surface, & "A threshold height used to detect the presence of a rigid-surface\n"//& "depressing the upper-surface of the model, such as an ice-shelf.\n"//& diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 415dadaf41..b5d9dd139a 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -21,6 +21,7 @@ module regrid_consts integer, parameter :: REGRIDDING_SLIGHT = REGRIDDING_NUM_TYPES !< Stretched coordinates in the !! lightest water, isopycnal below character(len=6), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string +character(len=6), parameter :: REGRIDDING_ZSTAR_STRING_OLD = "Z*" !< z* string (legacy name) character(len=6), parameter :: REGRIDDING_ZSTAR_STRING = "ZSTAR" !< z* string character(len=6), parameter :: REGRIDDING_RHO_STRING = "RHO" !< Rho string character(len=6), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string @@ -54,6 +55,7 @@ function coordinateMode(string) select case ( uppercase(trim(string)) ) case (trim(REGRIDDING_LAYER_STRING)); coordinateMode = REGRIDDING_LAYER case (trim(REGRIDDING_ZSTAR_STRING)); coordinateMode = REGRIDDING_ZSTAR + case (trim(REGRIDDING_ZSTAR_STRING_OLD)); coordinateMode = REGRIDDING_ZSTAR case (trim(REGRIDDING_RHO_STRING)); coordinateMode = REGRIDDING_RHO case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1 From 47c82b8dccc48681a2e09cbdb683f128de2e12b3 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 9 Sep 2016 11:54:05 -0400 Subject: [PATCH 15/16] Fixed reintegrate_column() for missing data cases - When the source or destination columns were completely vanished the reintegrate_column() in MOM_diag_vkernels was getting stuck in a never ending loop. - Added exit tests for these special cases. - Also added unit tests to cover these scenarios. --- src/framework/MOM_diag_vkernels.F90 | 41 ++++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/src/framework/MOM_diag_vkernels.F90 b/src/framework/MOM_diag_vkernels.F90 index 692a33faf2..18d283bffc 100644 --- a/src/framework/MOM_diag_vkernels.F90 +++ b/src/framework/MOM_diag_vkernels.F90 @@ -105,11 +105,16 @@ subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, real :: uh_src_rem, uh_dest_rem, duh ! Incremental amounts of stuff integer :: k_src, k_dest ! Index of cell in src and dest columns integer :: iter + logical :: src_ran_out, src_exists + + uh_dest(:) = missing_value k_src = 0 k_dest = 0 h_dest_rem = 0. h_src_rem = 0. + src_ran_out = .false. + src_exists = .false. do while(.true.) if (h_src_rem==0. .and. k_src0.) duh = uh_src_rem h_src_rem = 0. uh_src_rem = 0. h_dest_rem = max(0., h_dest_rem - dh) @@ -147,9 +160,11 @@ subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, h_dest_rem = 0. endif uh_dest(k_dest) = uh_dest(k_dest) + duh - if (k_dest+k_src==nsrc+ndest) exit + if (k_dest==ndest .and. (k_src==nsrc .or. h_dest_rem==0.)) exit enddo + if (.not. src_exists) uh_dest(1:ndest) = missing_value + end subroutine reintegrate_column !> Returns true if any unit tests for module MOM_diag_vkernels fail @@ -159,7 +174,7 @@ logical function diag_vkernels_unit_tests() logical :: fail write(0,*) '============ MOM_diag_kernels: diag_vkernels_unit_tests ==================' - write(0,*) '- - - - - - - - - - reintegration tests - - - - - - - - - - - - - - - - -' + write(0,*) '- - - - - - - - - - interpolation tests - - - - - - - - - - - - - - - - -' fail = test_interp('Identity: 3 layer', & 3, (/1.,2.,3./), (/1.,2.,3.,4./), & @@ -238,6 +253,26 @@ logical function diag_vkernels_unit_tests() 5, (/0.,3.,0.,3.,0./), (/0.,-4.,0.,2.,0./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail + fail = test_reintegrate('D: 3 layer to 3 (vanished)', & + 3, (/2.,2.,2./), (/-5.,2.,1./), & + 3, (/0.,0.,0./), (/0.,0.,0./) ) + diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail + + fail = test_reintegrate('D: 3 layer (vanished) to 3', & + 3, (/0.,0.,0./), (/-5.,2.,1./), & + 3, (/2.,2.,2./), (/missing_value, missing_value, missing_value/) ) + diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail + + fail = test_reintegrate('D: 3 layer (vanished) to 3 (vanished)', & + 3, (/0.,0.,0./), (/-5.,2.,1./), & + 3, (/0.,0.,0./), (/missing_value, missing_value, missing_value/) ) + diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail + + fail = test_reintegrate('D: 3 layer (vanished) to 3 (vanished)', & + 3, (/0.,0.,0./), (/0.,0.,0./), & + 3, (/0.,0.,0./), (/missing_value, missing_value, missing_value/) ) + diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail + write(0,*) '==========================================================================' contains From 1f0a1d755e8d6afa82efb1f308271c2831cb9451 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 9 Sep 2016 11:55:52 -0400 Subject: [PATCH 16/16] Renamed old diag-to-z diagnostics module - To avoid a conflict with the general coordinate diagnostics, I've renamed the diagnostics module registered within MOM_diag_to_z to ocean_model_zold. --- src/diagnostics/MOM_diag_to_Z.F90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 4c019bd9d8..4a3adefbfc 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -926,14 +926,14 @@ 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_zold', name, CS%axesTz, Time, & long_name, units, missing_value=CS%missing_tr(m), & standard_name=standard_name) - CS%id_tr_xyave(m) = register_diag_field('ocean_model_z', name//'_xyave', CS%axesZ, Time, & + CS%id_tr_xyave(m) = register_diag_field('ocean_model_zold', name//'_xyave', CS%axesZ, 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_zold', 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, & @@ -1038,24 +1038,24 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) call define_axes_group(diag, (/ diag%axesCv1%handles(1), diag%axesCv1%handles(2), zint_axis /), CS%axesCvzi) call define_axes_group(diag, (/ z_axis /), CS%axesZ) - CS%id_u_z = register_diag_field('ocean_model_z', 'u', CS%axesCuz, Time, & + CS%id_u_z = register_diag_field('ocean_model_zold', 'u', CS%axesCuz, Time, & 'Zonal Velocity in Depth Space', 'meter second-1', & missing_value=CS%missing_vel, cmor_field_name='uo', cmor_units='m s-1',& cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity') if (CS%id_u_z>0) call safe_alloc_ptr(CS%u_z,IsdB,IedB,jsd,jed,CS%nk_zspace) - CS%id_v_z = register_diag_field('ocean_model_z', 'v', CS%axesCvz, Time, & + CS%id_v_z = register_diag_field('ocean_model_zold', 'v', CS%axesCvz, Time, & 'Meridional Velocity in Depth Space', 'meter second-1', & missing_value=CS%missing_vel, cmor_field_name='vo', cmor_units='m s-1',& cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity') if (CS%id_v_z>0) call safe_alloc_ptr(CS%v_z,isd,ied,JsdB,JedB,CS%nk_zspace) - CS%id_uh_z = register_diag_field('ocean_model_z', 'uh', CS%axesCuz, Time, & + CS%id_uh_z = register_diag_field('ocean_model_zold', 'uh', CS%axesCuz, Time, & 'Zonal Mass Transport (including SGS param) in Depth Space', flux_units, & missing_value=CS%missing_trans) if (CS%id_uh_z>0) call safe_alloc_ptr(CS%uh_z,IsdB,IedB,jsd,jed,CS%nk_zspace) - CS%id_vh_z = register_diag_field('ocean_model_z', 'vh', CS%axesCvz, Time, & + CS%id_vh_z = register_diag_field('ocean_model_zold', 'vh', CS%axesCvz, Time, & 'Meridional Mass Transport (including SGS param) in Depth Space', flux_units,& missing_value=CS%missing_trans) if (CS%id_vh_z>0) call safe_alloc_ptr(CS%vh_z,isd,ied,JsdB,JedB,CS%nk_zspace) @@ -1064,25 +1064,25 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) ! Check whether diag-table is requesting any z-space files; issue a warning if it is. - id_test = register_diag_field('ocean_model_z', 'u', diag%axesCu1, Time, & + id_test = register_diag_field('ocean_model_zold', 'u', diag%axesCu1, Time, & 'Zonal Velocity in Depth Space', 'meter second-1', cmor_field_name='uo') if (id_test>0) call MOM_error(WARNING, & "MOM_diag_to_Z_init: u cannot be output without "//& "an appropriate depth-space target file.") - id_test = register_diag_field('ocean_model_z', 'v', diag%axesCv1, Time, & + id_test = register_diag_field('ocean_model_zold', 'v', diag%axesCv1, Time, & 'Meridional Velocity in Depth Space', 'meter second-1', cmor_field_name='vo') if (id_test>0) call MOM_error(WARNING, & "MOM_diag_to_Z_init: v cannot be output without "//& "an appropriate depth-space target file.") - id_test = register_diag_field('ocean_model_z', 'uh', diag%axesCu1, Time, & + id_test = register_diag_field('ocean_model_zold', 'uh', diag%axesCu1, Time, & 'Meridional Volume Transport in Depth Space', flux_units) if (id_test>0) call MOM_error(WARNING, & "MOM_diag_to_Z_init: uh cannot be output without "//& "an appropriate depth-space target file.") - id_test = register_diag_field('ocean_model_z', 'vh', diag%axesCv1, Time, & + id_test = register_diag_field('ocean_model_zold', 'vh', diag%axesCv1, Time, & 'Meridional Volume Transport in Depth Space', flux_units) if (id_test>0) call MOM_error(WARNING, & "MOM_diag_to_Z_init: vh cannot be output without "//& @@ -1337,7 +1337,7 @@ function register_Z_diag(var_desc, CS, day, missing) "register_Z_diag: unknown z_grid component "//trim(z_grid)) end select - register_Z_diag = register_diag_field("ocean_model_z", trim(var_name), axes, & + register_Z_diag = register_diag_field("ocean_model_zold", trim(var_name), axes, & day, trim(longname), trim(units), missing_value=missing) end function register_Z_diag @@ -1377,7 +1377,7 @@ function register_Zint_diag(var_desc, CS, day) "register_Z_diag: unknown hor_grid component "//trim(hor_grid)) end select - register_Zint_diag = register_diag_field("ocean_model_z", trim(var_name),& + register_Zint_diag = register_diag_field("ocean_model_zold", trim(var_name),& axes, day, trim(longname), trim(units), missing_value=CS%missing_value) end function register_Zint_diag