diff --git a/.doxygen b/.doxygen index 8b1e2aafab..6e58dd5075 100644 --- a/.doxygen +++ b/.doxygen @@ -68,7 +68,7 @@ OUTPUT_DIRECTORY = # performance problems for the file system. # The default value is: NO. -CREATE_SUBDIRS = YES +CREATE_SUBDIRS = NO # If the ALLOW_UNICODE_NAMES tag is set to YES, doxygen will allow non-ASCII # characters to appear in the names of generated files. If set to NO, non-ASCII diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 04369646e4..2442f6399b 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -14,9 +14,8 @@ module MOM_ALE use MOM_error_handler, only : callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, param_file_type, log_param -use MOM_io, only : file_exists, field_exists, MOM_read_data use MOM_io, only : vardesc, var_desc, fieldtype, SINGLE_FILE -use MOM_io, only : create_file, write_field, close_file, slasher +use MOM_io, only : create_file, write_field, close_file use MOM_regridding, only : initialize_regridding, regridding_main, end_regridding use MOM_regridding, only : uniformResolution use MOM_regridding, only : inflate_vanished_layers_old, setCoordinateResolution @@ -35,14 +34,14 @@ module MOM_ALE use MOM_remapping, only : remapping_core_h, remapping_core_w use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme use MOM_remapping, only : remapping_CS, dzFromH1H2 -use MOM_string_functions, only : uppercase, extractWord +use MOM_string_functions, only : uppercase, extractWord, extract_integer use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use regrid_defs, only : PRESSURE_RECONSTRUCTION_PLM !use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE -use regrid_consts, only : coordinateUnits, coordinateMode +use regrid_consts, only : coordinateUnits, coordinateMode, state_dependent use regrid_consts, only : REGRIDDING_ZSTAR, REGRIDDING_RHO use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_SLIGHT use regrid_edge_values, only : edge_values_implicit_h4 @@ -192,9 +191,7 @@ subroutine ALE_init( param_file, GV, max_depth, CS) default=.true.) ! Initialize and configure regridding - allocate( dz(GV%ke) ) - call ALE_initRegridding( GV, max_depth, param_file, mod, CS%regridCS, dz ) - deallocate( dz ) + call ALE_initRegridding( GV, max_depth, param_file, mod, CS%regridCS) ! Initialize and configure remapping call get_param(param_file, mod, "REMAPPING_SCHEME", string, & @@ -1028,386 +1025,26 @@ integer function pressureReconstructionScheme(CS) end function pressureReconstructionScheme - -!> Initialize regridding module -subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz ) +!> Initializes regridding for the main ALE algorithm +subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS) type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure real, intent(in) :: max_depth !< The maximum depth of the ocean, in m. type(param_file_type), intent(in) :: param_file !< parameter file character(len=*), intent(in) :: mod !< Name of calling module type(regridding_CS), intent(out) :: regridCS !< Regridding parameters and work arrays - real, dimension(:), intent(out) :: dz !< Resolution (thickness) in units of coordinate - ! Local variables - character(len=80) :: string, varName ! Temporary strings - character(len=40) :: coordMode, interpScheme, coordUnits ! Temporary strings - character(len=200) :: inputdir, fileName - character(len=320) :: message ! Temporary strings - integer :: K, ke - logical :: tmpLogical, fix_haloclines, set_max, do_sum - real :: filt_len, strat_tol, index_scale - real :: tmpReal, compress_fraction - real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int - integer :: nz_fixed_sfc - real :: rho_target(GV%ke+1) ! Target density used in HYBRID mode - real, dimension(size(dz)) :: h_max ! Maximum layer thicknesses, in m. - real, dimension(size(dz)) :: dz_max ! Thicknesses used to find maximum interface depths, in m. - real, dimension(size(dz)+1) :: z_max ! Maximum tinterface depths, in m. - - ke = size(dz) ! Number of levels in resolution vector - - call get_param(param_file, mod, "REGRIDDING_COORDINATE_MODE", coordMode, & + character(len=30) :: coord_mode + + call get_param(param_file, mod, "REGRIDDING_COORDINATE_MODE", coord_mode, & "Coordinate mode for vertical regridding.\n"//& "Choose among the following possibilities:\n"//& - trim(regriddingCoordinateModeDoc),& + trim(regriddingCoordinateModeDoc), & default=DEFAULT_COORDINATE_MODE, fail_if_missing=.true.) - call get_param(param_file, mod, "REGRIDDING_COORDINATE_UNITS", coordUnits, & - "Units of the regridding coordinuate.",& - default=coordinateUnits(coordMode)) - - call get_param(param_file, mod, "INTERPOLATION_SCHEME", interpScheme, & - "This sets the interpolation scheme to use to\n"//& - "determine the new grid. These parameters are\n"//& - "only relevant when REGRIDDING_COORDINATE_MODE is\n"//& - "set to a function of state. Otherwise, it is not\n"//& - "used. It can be one of the following schemes:\n"//& - trim(regriddingInterpSchemeDoc),& - default=regriddingDefaultInterpScheme) - - call get_param(param_file, mod, "REGRID_COMPRESSIBILITY_FRACTION", compress_fraction, & - "When interpolating potential density profiles we can add\n"//& - "some artificial compressibility solely to make homogenous\n"//& - "regions appear stratified.", default=0.) - - call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, & - compressibility_fraction=compress_fraction ) - - call get_param(param_file, mod, "ALE_COORDINATE_CONFIG", string, & - "Determines how to specify the coordinate\n"//& - "resolution. Valid options are:\n"//& - " PARAM - use the vector-parameter ALE_RESOLUTION\n"//& - " UNIFORM - uniformly distributed\n"//& - " FILE:string - read from a file. The string specifies\n"//& - " the filename and variable name, separated\n"//& - " by a comma or space, e.g. FILE:lev.nc,Z\n"//& - " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//& - " HYBRID:string - read from a file. The string specifies\n"//& - " the filename and two variable names, separated\n"//& - " by a comma or space, for sigma-2 and dz. e.g.\n"//& - " HYBRID:vgrid.nc,sigma2,dz",& - default='UNIFORM') - message = "The distribution of vertical resolution for the target\n"//& - "grid used for Eulerian-like coordinates. For example,\n"//& - "in z-coordinate mode, the parameter is a list of level\n"//& - "thicknesses (in m). In sigma-coordinate mode, the list\n"//& - "is of non-dimensional fractions of the water column." - select case ( trim(string) ) - case ("UNIFORM") - dz(:) = uniformResolution(GV%ke, coordMode, 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)) ) - call log_param(param_file, mod, "!ALE_RESOLUTION", dz, & - trim(message), units=trim(coordUnits)) - case ("PARAM") - call get_param(param_file, mod, "ALE_RESOLUTION", dz, & - trim(message), units=trim(coordUnits), fail_if_missing=.true.) - case default - if (index(trim(string),'FILE:')==1) then - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - - if (string(6:6)=='.' .or. string(6:6)=='/') then - ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path - fileName = trim( extractWord(trim(string(6:80)), 1) ) - else - ! Otherwise assume we should look for the file in INPUTDIR - fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) - endif - if (.not. file_exists(fileName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") - - varName = trim( extractWord(trim(string(6:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") - if (len_trim(varName)==0) then - if (field_exists(fileName,'dz')) then; varName = 'dz' - elseif (field_exists(fileName,'dsigma')) then; varName = 'dsigma' - elseif (field_exists(fileName,'ztest')) then; varName = 'ztest' - else ; call MOM_error(FATAL,"ALE_initRegridding: "// & - "Coordinate variable not specified and none could be guessed.") - endif - endif - call MOM_read_data(trim(fileName), trim(varName), dz) - call log_param(param_file, mod, "!ALE_RESOLUTION", dz, & - trim(message), units=coordinateUnits(coordMode)) - elseif (index(trim(string),'FNC1:')==1) then - call dz_function1( trim(string(6:)), dz ) - call log_param(param_file, mod, "!ALE_RESOLUTION", dz, & - trim(message), units=coordinateUnits(coordMode)) - elseif (index(trim(string),'HYBRID:')==1) then - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - - ! The following assumes the FILE: syntax of above but without "FILE:" in the string - fileName = trim( extractWord(trim(string(8:)), 1) ) - if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) - if (.not. file_exists(fileName)) call MOM_error(FATAL,"ALE_initRegridding: HYBRID "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") - varName = trim( extractWord(trim(string(8:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,"ALE_initRegridding: HYBRID "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") - call MOM_read_data(trim(fileName), trim(varName), rho_target) - call set_target_densities( regridCS, rho_target ) - varName = trim( extractWord(trim(string(8:)), 3) ) - if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz - call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) - else ! Read dz from file - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,"ALE_initRegridding: HYBRID "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") - call MOM_read_data(trim(fileName), trim(varName), dz) - endif - call log_param(param_file, mod, "!ALE_RESOLUTION", dz, & - trim(message), units=coordinateUnits(coordMode)) - call log_param(param_file, mod, "!TARGET_DENSITIES", rho_target, & - 'HYBRID target densities for itnerfaces', units=coordinateUnits(coordMode)) - - else - call MOM_error(FATAL,"ALE_initRegridding: "// & - "Unrecognized coordinate configuraiton"//trim(string)) - endif - end select - if (coordinateMode(coordMode) == REGRIDDING_ZSTAR .or. & - coordinateMode(coordMode) == REGRIDDING_HYCOM1 .or. & - coordinateMode(coordMode) == REGRIDDING_SLIGHT) then - ! Adjust target grid to be consistent with max_depth - ! This is a work around to the from_Z initialization... ??? - tmpReal = sum( dz(:) ) - if (tmpReal < max_depth) then - dz(ke) = dz(ke) + ( max_depth - tmpReal ) - elseif (tmpReal > max_depth) then - if ( dz(ke) + ( max_depth - tmpReal ) > 0. ) then - dz(ke) = dz(ke) + ( max_depth - tmpReal ) - else - call MOM_error(FATAL,"ALE_initRegridding: "// & - "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) - endif - endif - endif - call setCoordinateResolution( dz, regridCS ) - if (coordinateMode(coordMode) == REGRIDDING_RHO) call set_target_densities_from_GV( GV, regridCS ) - - call get_param(param_file, mod, "MIN_THICKNESS", tmpReal, & - "When regridding, this is the minimum layer\n"//& - "thickness allowed.", units="m",& - default=regriddingDefaultMinThickness ) - - call get_param(param_file, mod, "BOUNDARY_EXTRAPOLATION", tmpLogical, & - "When defined, a proper high-order reconstruction\n"//& - "scheme is used within boundary cells rather\n"//& - "than PCM. E.g., if PPM is used for remapping, a\n"//& - "PPM reconstruction will also be used within\n"//& - "boundary cells.", default=regriddingDefaultBoundaryExtrapolation) - call set_regrid_params( regridCS, min_thickness=tmpReal, boundary_extrapolation=tmpLogical ) - - if (coordinateMode(coordMode) == REGRIDDING_SLIGHT) then - ! Set SLight-specific regridding parameters. - call get_param(param_file, mod, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, & - "The nominal thickness of fixed thickness near-surface\n"//& - "layers with the SLight coordinate.", units="m", default=1.0) - call get_param(param_file, mod, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, & - "The number of fixed-depth surface layers with the SLight\n"//& - "coordinate.", units="nondimensional", default=2) - call get_param(param_file, mod, "SLIGHT_SURFACE_AVG_DEPTH", Rho_avg_depth, & - "The thickness of the surface region over which to average\n"//& - "when calculating the density to use to define the interior\n"//& - "with the SLight coordinate.", units="m", default=1.0) - call get_param(param_file, mod, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, & - "The number of layers to offset the surface density when\n"//& - "defining where the interior ocean starts with SLight.", & - units="nondimensional", default=2.0) - call get_param(param_file, mod, "SLIGHT_FIX_HALOCLINES", fix_haloclines, & - "If true, identify regions above the reference pressure\n"//& - "where the reference pressure systematically underestimates\n"//& - "the stratification and use this in the definition of the\n"//& - "interior with the SLight coordinate.", default=.false.) - - call set_regrid_params( regridCS, dz_min_surface=dz_fixed_sfc, & - nz_fixed_surface=nz_fixed_sfc, Rho_ML_avg_depth=Rho_avg_depth, & - nlay_ML_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines) - if (fix_haloclines) then - ! Set additional parameters related to SLIGHT_FIX_HALOCLINES. - call get_param(param_file, mod, "HALOCLINE_FILTER_LENGTH", filt_len, & - "A length scale over which to smooth the temperature and\n"//& - "salinity before identifying erroneously unstable haloclines.", & - units="m", default=2.0) - call get_param(param_file, mod, "HALOCLINE_STRAT_TOL", strat_tol, & - "A tolerance for the ratio of the stratification of the\n"//& - "apparent coordinate stratification to the actual value\n"//& - "that is used to identify erroneously unstable haloclines.\n"//& - "This ratio is 1 when they are equal, and sensible values \n"//& - "are between 0 and 0.5.", units="nondimensional", default=0.2) - call set_regrid_params(regridCS, halocline_filt_len=filt_len, & - halocline_strat_tol=strat_tol) - endif - endif - - call get_param(param_file, mod, "MAXIMUM_INT_DEPTH_CONFIG", string, & - "Determines how to specify the maximum interface depths.\n"//& - "Valid options are:\n"//& - " NONE - there are no maximum interface depths\n"//& - " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//& - " FILE:string - read from a file. The string specifies\n"//& - " the filename and variable name, separated\n"//& - " by a comma or space, e.g. FILE:lev.nc,Z\n"//& - " FNC1:string - FNC1:dz_min,H_total,power,precision",& - default='NONE') - message = "The list of maximum depths for each interface." - if ( trim(string) == "NONE") then - ! Do nothing. - elseif ( trim(string) == "PARAM") then - call get_param(param_file, mod, "MAXIMUM_INTERFACE_DEPTHS", z_max, & - trim(message), units="m", fail_if_missing=.true.) - call set_regrid_max_depths( regridCS, z_max, GV%m_to_H ) - elseif (index(trim(string),'FILE:')==1) then - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - - if (string(6:6)=='.' .or. string(6:6)=='/') then - ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path - fileName = trim( extractWord(trim(string(6:80)), 1) ) - else - ! Otherwise assume we should look for the file in INPUTDIR - fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) - endif - if (.not. file_exists(fileName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") - - do_sum = .false. - varName = trim( extractWord(trim(string(6:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") - if (len_trim(varName)==0) then - if (field_exists(fileName,'z_max')) then; varName = 'z_max' - elseif (field_exists(fileName,'dz')) then; varName = 'dz' ; do_sum = .true. - elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' ; do_sum = .true. - else ; call MOM_error(FATAL,"ALE_initRegridding: "// & - "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") - endif - endif - if (do_sum) then - call MOM_read_data(trim(fileName), trim(varName), dz_max) - z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(k) ; enddo - else - call MOM_read_data(trim(fileName), trim(varName), z_max) - endif - call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, & - trim(message), units=coordinateUnits(coordMode)) - call set_regrid_max_depths( regridCS, z_max, GV%m_to_H ) - elseif (index(trim(string),'FNC1:')==1) then - call dz_function1( trim(string(6:)), dz_max ) - if ((coordinateMode(coordMode) == REGRIDDING_SLIGHT) .and. & - (dz_fixed_sfc > 0.0)) then - do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo - endif - z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(K) ; enddo - call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, & - trim(message), units=coordinateUnits(coordMode)) - call set_regrid_max_depths( regridCS, z_max, GV%m_to_H ) - else - call MOM_error(FATAL,"ALE_initRegridding: "// & - "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string)) - endif - - ! Optionally specify maximum thicknesses for each layer, enforced by moving - ! the interface below a layer downward. - call get_param(param_file, mod, "MAX_LAYER_THICKNESS_CONFIG", string, & - "Determines how to specify the maximum layer thicknesses.\n"//& - "Valid options are:\n"//& - " NONE - there are no maximum layer thicknesses\n"//& - " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//& - " FILE:string - read from a file. The string specifies\n"//& - " the filename and variable name, separated\n"//& - " by a comma or space, e.g. FILE:lev.nc,Z\n"//& - " FNC1:string - FNC1:dz_min,H_total,power,precision",& - default='NONE') - message = "The list of maximum thickness for each layer." - if ( trim(string) == "NONE") then - ! Do nothing. - elseif ( trim(string) == "PARAM") then - call get_param(param_file, mod, "MAX_LAYER_THICKNESS", h_max, & - trim(message), units="m", fail_if_missing=.true.) - call set_regrid_max_thickness( regridCS, h_max, GV%m_to_H ) - elseif (index(trim(string),'FILE:')==1) then - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - - if (string(6:6)=='.' .or. string(6:6)=='/') then - ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path - fileName = trim( extractWord(trim(string(6:80)), 1) ) - else - ! Otherwise assume we should look for the file in INPUTDIR - fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) - endif - if (.not. file_exists(fileName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") - - varName = trim( extractWord(trim(string(6:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,"ALE_initRegridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") - if (len_trim(varName)==0) then - if (field_exists(fileName,'h_max')) then; varName = 'h_max' - elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' - else ; call MOM_error(FATAL,"ALE_initRegridding: "// & - "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") - endif - endif - call MOM_read_data(trim(fileName), trim(varName), h_max) - call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, & - trim(message), units=coordinateUnits(coordMode)) - call set_regrid_max_thickness( regridCS, h_max, GV%m_to_H ) - elseif (index(trim(string),'FNC1:')==1) then - call dz_function1( trim(string(6:)), h_max ) - call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, & - trim(message), units=coordinateUnits(coordMode)) - call set_regrid_max_thickness( regridCS, h_max, GV%m_to_H ) - else - call MOM_error(FATAL,"ALE_initRegridding: "// & - "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string)) - endif + call initialize_regridding(regridCS, GV, max_depth, param_file, mod, coord_mode, '', '') end subroutine ALE_initRegridding - -!> Parses a string and generates a dz(:) profile -subroutine dz_function1( string, dz ) - character(len=*), intent(in) :: string !< String with list of parameters - real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses - - ! Local variables - integer :: nk, k - real :: dz_min, power, prec, H_total - - nk = size(dz) ! Number of cells - prec = -1024. - read( string, *) dz_min, H_total, power, prec - if (prec == -1024.) call MOM_error(FATAL,"dz_function1: "// & - "Problem reading FNC1: string ="//trim(string)) - ! Create profile of ( dz - dz_min ) - do k = 1, nk - dz(k) = (real(k-1)/real(nk-1))**power - enddo - dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total - dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec - dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total - dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec - dz(nk) = dz(nk) + ( H_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer - dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec - dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min - -end subroutine dz_function1 - - !> Query the target coordinate interfaces positions function ALE_getCoordinate( CS ) type(ALE_CS), pointer :: CS !< module control structure @@ -1449,7 +1086,7 @@ subroutine ALE_update_regrid_weights( dt, CS ) if (CS%regrid_time_scale > 0.0) then w = CS%regrid_time_scale / (CS%regrid_time_scale + dt) endif - call set_regrid_params( CS%regridCS, old_grid_weight=w ) + call set_regrid_params(CS%regridCS, old_grid_weight=w) endif end subroutine ALE_update_regrid_weights diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9180cbd70e..c2eedf49f8 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -4,11 +4,15 @@ module MOM_regridding ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : param_file_type, get_param, log_param +use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data +use MOM_io, only : vardesc, var_desc, fieldtype, SINGLE_FILE +use MOM_io, only : create_file, write_field, close_file, slasher use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : EOS_type, calculate_density, calculate_density_derivs use MOM_EOS, only : calculate_compress -use MOM_string_functions,only : uppercase +use MOM_string_functions,only : uppercase, extractWord, extract_integer, extract_real use regrid_edge_values, only : edge_values_explicit_h2, edge_values_explicit_h4 use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6 @@ -22,12 +26,15 @@ module MOM_regridding use P3M_functions, only : P3M_interpolation, P3M_boundary_extrapolation use MOM_remapping, only : remapping_core_h use MOM_remapping, only : remapping_CS +use regrid_consts, only : state_dependent, coordinateUnits use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA use regrid_consts, only : REGRIDDING_ARBITRARY use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_SLIGHT +use netcdf ! Used by check_grid_def() + implicit none ; private #include @@ -139,7 +146,7 @@ module MOM_regridding public initialize_regridding, end_regridding, regridding_main public inflate_vanished_layers_old, check_remapping_grid, check_grid_column public adjust_interface_motion -public set_regrid_params +public set_regrid_params, get_regrid_size public uniformResolution, setCoordinateResolution public build_zstar_column, build_rho_column, build_sigma_column public set_target_densities_from_GV, set_target_densities @@ -210,52 +217,566 @@ module MOM_regridding contains -!> Initialization of regridding -subroutine initialize_regridding( nk, coordMode, interpScheme, CS, compressibility_fraction ) - integer, intent(in) :: nk !< Number of levels - character(len=*), intent(in) :: coordMode !< Coordinate mode to use - character(len=*), intent(in) :: interpScheme !< Interpolation mode to use - type(regridding_CS), intent(inout) :: CS !< Regridding control structure - real, optional, intent(in) :: compressibility_fraction !< Fraction of compressibility - !! to add to potential density profiles (nondim) - - CS%nk = nk - - CS%regridding_scheme = coordinateMode(coordMode) - - select case ( uppercase(trim(interpScheme)) ) - case ("P1M_H2"); CS%interpolation_scheme = INTERPOLATION_P1M_H2 - case ("P1M_H4"); CS%interpolation_scheme = INTERPOLATION_P1M_H4 - case ("P1M_IH2"); CS%interpolation_scheme = INTERPOLATION_P1M_IH4 - case ("PLM"); CS%interpolation_scheme = INTERPOLATION_PLM - case ("PPM_H4"); CS%interpolation_scheme = INTERPOLATION_PPM_H4 - case ("PPM_IH4"); CS%interpolation_scheme = INTERPOLATION_PPM_IH4 - case ("P3M_IH4IH3"); CS%interpolation_scheme = INTERPOLATION_P3M_IH4IH3 - case ("P3M_IH6IH5"); CS%interpolation_scheme = INTERPOLATION_P3M_IH6IH5 - case ("PQM_IH4IH3"); CS%interpolation_scheme = INTERPOLATION_PQM_IH4IH3 - case ("PQM_IH6IH5"); CS%interpolation_scheme = INTERPOLATION_PQM_IH6IH5 - case default ; call MOM_error(FATAL, "read_regridding_options: "//& - "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interpScheme)//").") - end select +!> Initialization and configures a regridding control structure based on customizable run-time parameters +subroutine initialize_regridding(CS, GV, max_depth, param_file, mod, coord_mode, param_prefix, param_suffix) + type(regridding_CS), intent(inout) :: CS !< Regridding control structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, intent(in) :: max_depth !< The maximum depth of the ocean, in m. + type(param_file_type), intent(in) :: param_file !< Parameter file + character(len=*), intent(in) :: mod !< Name of calling module. + character(len=*), intent(in) :: coord_mode !< Coordinate mode + character(len=*), intent(in) :: param_prefix !< String to prefix to parameter names. + !! If empty, causes main model parameters to be used. + character(len=*), intent(in) :: param_suffix !< String to append to parameter names. + ! Local variables + integer :: ke ! Number of levels + character(len=80) :: string, string2, varName ! Temporary strings + character(len=40) :: coord_units, param_name, coord_res_param ! Temporary strings + character(len=200) :: inputdir, fileName + character(len=320) :: message ! Temporary strings + character(len=12) :: expected_units ! Temporary strings + logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters + logical :: coord_is_state_dependent, ierr + real :: filt_len, strat_tol, index_scale, tmpReal + real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int + integer :: nz_fixed_sfc, k, nzf(4) + real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate + real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses, in m. + real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths, in m. + real, dimension(:), allocatable :: z_max ! Maximum interface depths, in m. + real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode + ! Thicknesses that give level centers corresponding to table 2 of WOA09 + real, dimension(40) :: woa09_dz = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., & + 37.5, 50., 50., 75., 100., 100., 100., 100., & + 100., 100., 100., 100., 100., 100., 100., 175., & + 250., 375., 500., 500., 500., 500., 500., 500., & + 500., 500., 500., 500., 500., 500., 500., 500. /) + + call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + + main_parameters=.false. + if (len_trim(param_prefix)==0) main_parameters=.true. + if (main_parameters .and. len_trim(param_suffix)>0) call MOM_error(FATAL,trim(mod)//', initialize_regridding: '// & + 'Suffix provided without prefix for parameter names!') + + CS%nk = 0 + CS%regridding_scheme = coordinateMode(coord_mode) + coord_is_state_dependent = state_dependent(coord_mode) + + if (main_parameters) then + ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS) + call get_param(param_file, mod, "REGRIDDING_COORDINATE_UNITS", coord_units, & + "Units of the regridding coordinuate.",& + default=coordinateUnits(coord_mode)) + else + coord_units=coordinateUnits(coord_mode) + endif + + if (coord_is_state_dependent) then + if (main_parameters) then + param_name = "INTERPOLATION_SCHEME" + string2 = regriddingDefaultInterpScheme + else + param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix) + string2 = 'PPM_H4' ! Default for diagnostics + endif + call get_param(param_file, mod, "INTERPOLATION_SCHEME", string, & + "This sets the interpolation scheme to use to\n"//& + "determine the new grid. These parameters are\n"//& + "only relevant when REGRIDDING_COORDINATE_MODE is\n"//& + "set to a function of state. Otherwise, it is not\n"//& + "used. It can be one of the following schemes:\n"//& + trim(regriddingInterpSchemeDoc), default=trim(string2)) + call set_regrid_params(CS, interp_scheme=string) + else + CS%interpolation_scheme = -1 ! Cause error if ever used + endif + + if (main_parameters .and. coord_is_state_dependent) then + call get_param(param_file, mod, "REGRID_COMPRESSIBILITY_FRACTION", tmpReal, & + "When interpolating potential density profiles we can add\n"//& + "some artificial compressibility solely to make homogenous\n"//& + "regions appear stratified.", default=0.) + call set_regrid_params(CS, compress_fraction=tmpReal) + endif + + ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG) + if (main_parameters) then + param_name = "ALE_COORDINATE_CONFIG" + coord_res_param = "ALE_RESOLUTION" + string2 = 'UNIFORM' + else + param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix) + coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix) + string2 = 'UNIFORM' + if (max_depth>3000.) string2='WOA09' ! For convenience + endif + call get_param(param_file, mod, param_name, string, & + "Determines how to specify the coordinate\n"//& + "resolution. Valid options are:\n"//& + " PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//& + " UNIFORM[:N] - uniformly distributed\n"//& + " FILE:string - read from a file. The string specifies\n"//& + " the filename and variable name, separated\n"//& + " by a comma or space, e.g. FILE:lev.nc,dz\n"//& + " or FILE:lev.nc,interfaces=zw\n"//& + " WOA09[:N] - the WOA09 vertical grid (approximately)\n"//& + " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//& + " HYBRID:string - read from a file. The string specifies\n"//& + " the filename and two variable names, separated\n"//& + " by a comma or space, for sigma-2 and dz. e.g.\n"//& + " HYBRID:vgrid.nc,sigma2,dz",& + default=trim(string2)) + message = "The distribution of vertical resolution for the target\n"//& + "grid used for Eulerian-like coordinates. For example,\n"//& + "in z-coordinate mode, the parameter is a list of level\n"//& + "thicknesses (in m). In sigma-coordinate mode, the list\n"//& + "is of non-dimensional fractions of the water column." + if (index(trim(string),'UNIFORM')==1) then + if (len_trim(string)==7) then + ke = GV%ke ! Use model nk by default + tmpReal = max_depth + elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then + ! Format is "UNIFORM:N" or "UNIFORM:N,dz" + ke = extract_integer(string(9:len_trim(string)),'',1) + tmpReal = extract_real(string(9:len_trim(string)),',',2,missing_value=max_depth) + else + call MOM_error(FATAL,trim(mod)//', initialize_regridding: '// & + 'Unable to interpret "'//trim(string)//'".') + endif + allocate(dz(ke)) + dz(:) = uniformResolution(ke, coord_mode, tmpReal, & + GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2)), & + GV%Rlay(ke)+0.5*(GV%Rlay(ke)-GV%Rlay(ke-1)) ) + if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, & + trim(message), units=trim(coord_units)) + elseif (trim(string)=='PARAM') then + ! Read coordinate resolution (main model = ALE_RESOLUTION) + ke = GV%ke ! Use model nk by default + allocate(dz(ke)) + call get_param(param_file, mod, coord_res_param, dz, & + trim(message), units=trim(coord_units), fail_if_missing=.true.) + elseif (index(trim(string),'FILE:')==1) then + ! FILE:filename,var_name is assumed to be reading level thickness variables + ! FILE:filename,interfaces=var_name reads positions + if (string(6:6)=='.' .or. string(6:6)=='/') then + ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path + fileName = trim( extractWord(trim(string(6:80)), 1) ) + else + ! Otherwise assume we should look for the file in INPUTDIR + fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) + endif + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + + varName = trim( extractWord(trim(string(6:)), 2) ) + if (len_trim(varName)==0) then + if (field_exists(fileName,'dz')) then; varName = 'dz' + elseif (field_exists(fileName,'dsigma')) then; varName = 'dsigma' + elseif (field_exists(fileName,'ztest')) then; varName = 'ztest' + else ; call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Coordinate variable not specified and none could be guessed.") + endif + endif + ! This check fails when the variable is a dimension variable! -AJA + !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + ! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + if (CS%regridding_scheme == REGRIDDING_SIGMA) then + expected_units = 'nondim' + elseif (CS%regridding_scheme == REGRIDDING_RHO) then + expected_units = 'kg m-3' + else + expected_units = 'meters' + endif + if (index(trim(varName),'interfaces=')==1) then + varName=trim(varName(12:)) + call check_grid_def(filename, varName, expected_units, message, ierr) + if (ierr) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "//& + "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message)) + call field_size(trim(fileName), trim(varName), nzf) + ke = nzf(1)-1 + allocate(dz(ke)) + allocate(z_max(ke+1)) + call MOM_read_data(trim(fileName), trim(varName), z_max) + dz(:) = abs(z_max(1:ke) - z_max(2:ke+1)) + deallocate(z_max) + else + ! Assume reading resolution + call field_size(trim(fileName), trim(varName), nzf) + ke = nzf(1) + allocate(dz(ke)) + call MOM_read_data(trim(fileName), trim(varName), dz) + endif + if (main_parameters .and. ke/=GV%ke) then + call MOM_error(FATAL,trim(mod)//', initialize_regridding: '// & + 'Mismatch in number of model levels and "'//trim(string)//'".') + endif + if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + elseif (index(trim(string),'FNC1:')==1) then + ke = GV%ke; allocate(dz(ke)) + call dz_function1( trim(string(6:)), dz ) + if (main_parameters) call log_param(param_file, mod, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + elseif (index(trim(string),'HYBRID:')==1) then + ke = GV%ke; allocate(dz(ke)) + ! The following assumes the FILE: syntax of above but without "FILE:" in the string + allocate(rho_target(ke+1)) + fileName = trim( extractWord(trim(string(8:)), 1) ) + if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: HYBRID "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + varName = trim( extractWord(trim(string(8:)), 2) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: HYBRID "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + call MOM_read_data(trim(fileName), trim(varName), rho_target) + varName = trim( extractWord(trim(string(8:)), 3) ) + if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz + call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) + else ! Read dz from file + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: HYBRID "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + call MOM_read_data(trim(fileName), trim(varName), dz) + endif + if (main_parameters) then + call log_param(param_file, mod, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + call log_param(param_file, mod, "!TARGET_DENSITIES", rho_target, & + 'HYBRID target densities for itnerfaces', units=coordinateUnits(coord_mode)) + endif + elseif (index(trim(string),'WOA09')==1) then + if (len_trim(string)==5) then + tmpReal = 0. ; ke = 0 + do while (tmpReal40 .or. ke<1) call MOM_error(FATAL,trim(mod)//', initialize_regridding: '// & + 'For "WOA05:N" N must 0 max_depth) then + if ( dz(ke) + ( max_depth - tmpReal ) > 0. ) then + dz(ke) = dz(ke) + ( max_depth - tmpReal ) + else + call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) + endif + endif + endif + endif + + CS%nk=ke + ! Target resolution (for fixed coordinates) + allocate( CS%coordinateResolution(CS%nk) ); CS%coordinateResolution(:) = -1.E30 + if (state_dependent(CS%regridding_scheme)) then + ! Target values + allocate( CS%target_density(CS%nk+1) ); CS%target_density(:) = -1.E30 + endif + + call setCoordinateResolution(dz, CS) + if (allocated(rho_target)) then + call set_target_densities(CS, rho_target) + deallocate(rho_target) + endif + ! \todo This line looks like it would overwrite the target densities set just above? + if (coordinateMode(coord_mode) == REGRIDDING_RHO) call set_target_densities_from_GV(GV, CS) + + if (main_parameters) then + call get_param(param_file, mod, "MIN_THICKNESS", tmpReal, & + "When regridding, this is the minimum layer\n"//& + "thickness allowed.", units="m",& + default=regriddingDefaultMinThickness ) + call set_regrid_params(CS, min_thickness=tmpReal) + else + call set_regrid_params(CS, min_thickness=0.) + endif - CS%boundary_extrapolation = regriddingDefaultBoundaryExtrapolation + if (main_parameters .and. coord_is_state_dependent) then + call get_param(param_file, mod, "BOUNDARY_EXTRAPOLATION", tmpLogical, & + "When defined, a proper high-order reconstruction\n"//& + "scheme is used within boundary cells rather\n"//& + "than PCM. E.g., if PPM is used for remapping, a\n"//& + "PPM reconstruction will also be used within\n"//& + "boundary cells.", default=regriddingDefaultBoundaryExtrapolation) + call set_regrid_params(CS, boundary_extrapolation=tmpLogical) + else + call set_regrid_params(CS, boundary_extrapolation=.false.) + endif - CS%min_thickness = regriddingDefaultMinThickness + if (coordinateMode(coord_mode) == REGRIDDING_SLIGHT) then + ! Set SLight-specific regridding parameters. + call get_param(param_file, mod, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, & + "The nominal thickness of fixed thickness near-surface\n"//& + "layers with the SLight coordinate.", units="m", default=1.0) + call get_param(param_file, mod, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, & + "The number of fixed-depth surface layers with the SLight\n"//& + "coordinate.", units="nondimensional", default=2) + call get_param(param_file, mod, "SLIGHT_SURFACE_AVG_DEPTH", Rho_avg_depth, & + "The thickness of the surface region over which to average\n"//& + "when calculating the density to use to define the interior\n"//& + "with the SLight coordinate.", units="m", default=1.0) + call get_param(param_file, mod, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, & + "The number of layers to offset the surface density when\n"//& + "defining where the interior ocean starts with SLight.", & + units="nondimensional", default=2.0) + call get_param(param_file, mod, "SLIGHT_FIX_HALOCLINES", fix_haloclines, & + "If true, identify regions above the reference pressure\n"//& + "where the reference pressure systematically underestimates\n"//& + "the stratification and use this in the definition of the\n"//& + "interior with the SLight coordinate.", default=.false.) + + call set_regrid_params(CS, dz_min_surface=dz_fixed_sfc, & + nz_fixed_surface=nz_fixed_sfc, Rho_ML_avg_depth=Rho_avg_depth, & + nlay_ML_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines) + if (fix_haloclines) then + ! Set additional parameters related to SLIGHT_FIX_HALOCLINES. + call get_param(param_file, mod, "HALOCLINE_FILTER_LENGTH", filt_len, & + "A length scale over which to smooth the temperature and\n"//& + "salinity before identifying erroneously unstable haloclines.", & + units="m", default=2.0) + call get_param(param_file, mod, "HALOCLINE_STRAT_TOL", strat_tol, & + "A tolerance for the ratio of the stratification of the\n"//& + "apparent coordinate stratification to the actual value\n"//& + "that is used to identify erroneously unstable haloclines.\n"//& + "This ratio is 1 when they are equal, and sensible values \n"//& + "are between 0 and 0.5.", units="nondimensional", default=0.2) + call set_regrid_params(CS, halocline_filt_len=filt_len, & + halocline_strat_tol=strat_tol) + endif - if (present(compressibility_fraction)) CS%compressibility_fraction = compressibility_fraction + endif - call allocate_regridding( CS ) + if (main_parameters .and. coord_is_state_dependent) then + call get_param(param_file, mod, "MAXIMUM_INT_DEPTH_CONFIG", string, & + "Determines how to specify the maximum interface depths.\n"//& + "Valid options are:\n"//& + " NONE - there are no maximum interface depths\n"//& + " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//& + " FILE:string - read from a file. The string specifies\n"//& + " the filename and variable name, separated\n"//& + " by a comma or space, e.g. FILE:lev.nc,Z\n"//& + " FNC1:string - FNC1:dz_min,H_total,power,precision",& + default='NONE') + message = "The list of maximum depths for each interface." + allocate(z_max(ke+1)) + allocate(dz_max(ke)) + if ( trim(string) == "NONE") then + ! Do nothing. + elseif ( trim(string) == "PARAM") then + call get_param(param_file, mod, "MAXIMUM_INTERFACE_DEPTHS", z_max, & + trim(message), units="m", fail_if_missing=.true.) + call set_regrid_max_depths(CS, z_max, GV%m_to_H) + elseif (index(trim(string),'FILE:')==1) then + if (string(6:6)=='.' .or. string(6:6)=='/') then + ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path + fileName = trim( extractWord(trim(string(6:80)), 1) ) + else + ! Otherwise assume we should look for the file in INPUTDIR + fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) + endif + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + + do_sum = .false. + varName = trim( extractWord(trim(string(6:)), 2) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + if (len_trim(varName)==0) then + if (field_exists(fileName,'z_max')) then; varName = 'z_max' + elseif (field_exists(fileName,'dz')) then; varName = 'dz' ; do_sum = .true. + elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' ; do_sum = .true. + else ; call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") + endif + endif + if (do_sum) then + call MOM_read_data(trim(fileName), trim(varName), dz_max) + z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(k) ; enddo + else + call MOM_read_data(trim(fileName), trim(varName), z_max) + endif + call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, & + trim(message), units=coordinateUnits(coord_mode)) + call set_regrid_max_depths(CS, z_max, GV%m_to_H) + elseif (index(trim(string),'FNC1:')==1) then + call dz_function1( trim(string(6:)), dz_max ) + if ((coordinateMode(coord_mode) == REGRIDDING_SLIGHT) .and. & + (dz_fixed_sfc > 0.0)) then + do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo + endif + z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(K) ; enddo + call log_param(param_file, mod, "!MAXIMUM_INT_DEPTHS", z_max, & + trim(message), units=coordinateUnits(coord_mode)) + call set_regrid_max_depths(CS, z_max, GV%m_to_H) + else + call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string)) + endif + deallocate(z_max) + deallocate(dz_max) + + ! Optionally specify maximum thicknesses for each layer, enforced by moving + ! the interface below a layer downward. + call get_param(param_file, mod, "MAX_LAYER_THICKNESS_CONFIG", string, & + "Determines how to specify the maximum layer thicknesses.\n"//& + "Valid options are:\n"//& + " NONE - there are no maximum layer thicknesses\n"//& + " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//& + " FILE:string - read from a file. The string specifies\n"//& + " the filename and variable name, separated\n"//& + " by a comma or space, e.g. FILE:lev.nc,Z\n"//& + " FNC1:string - FNC1:dz_min,H_total,power,precision",& + default='NONE') + message = "The list of maximum thickness for each layer." + allocate(h_max(ke)) + if ( trim(string) == "NONE") then + ! Do nothing. + elseif ( trim(string) == "PARAM") then + call get_param(param_file, mod, "MAX_LAYER_THICKNESS", h_max, & + trim(message), units="m", fail_if_missing=.true.) + call set_regrid_max_thickness(CS, h_max, GV%m_to_H) + elseif (index(trim(string),'FILE:')==1) then + if (string(6:6)=='.' .or. string(6:6)=='/') then + ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path + fileName = trim( extractWord(trim(string(6:80)), 1) ) + else + ! Otherwise assume we should look for the file in INPUTDIR + fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) + endif + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + + varName = trim( extractWord(trim(string(6:)), 2) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + if (len_trim(varName)==0) then + if (field_exists(fileName,'h_max')) then; varName = 'h_max' + elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' + else ; call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") + endif + endif + call MOM_read_data(trim(fileName), trim(varName), h_max) + call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, & + trim(message), units=coordinateUnits(coord_mode)) + call set_regrid_max_thickness(CS, h_max, GV%m_to_H) + elseif (index(trim(string),'FNC1:')==1) then + call dz_function1( trim(string(6:)), h_max ) + call log_param(param_file, mod, "!MAX_LAYER_THICKNESS", h_max, & + trim(message), units=coordinateUnits(coord_mode)) + call set_regrid_max_thickness(CS, h_max, GV%m_to_H) + else + call MOM_error(FATAL,trim(mod)//", initialize_regridding: "// & + "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string)) + endif + deallocate(h_max) + endif + deallocate(dz) end subroutine initialize_regridding +!> Do some basic checks on the vertical grid definition file, variable +subroutine check_grid_def(filename, varname, expected_units, msg, ierr) + character(len=*), intent(in) :: filename !< File name + character(len=*), intent(in) :: varname !< Variable name + character(len=*), intent(in) :: expected_units !< Expected units of variable + character(len=*), intent(inout) :: msg !< Message to use for errors + logical, intent(out) :: ierr !< True if an error occurs + ! Local variables + character (len=200) :: units, long_name + integer :: ncid, status, intid, vid + + ierr = .false. + status = NF90_OPEN(trim(filename), NF90_NOWRITE, ncid); + if (status /= NF90_NOERR) then + ierr = .true. + msg = 'File not found: '//trim(filename) + return + endif + + status = NF90_INQ_VARID(ncid, trim(varname), vid) + if (status /= NF90_NOERR) then + ierr = .true. + msg = 'Var not found: '//trim(varname) + return + endif + + status = NF90_GET_ATT(ncid, vid, "units", units) + if (status /= NF90_NOERR) then + ierr = .true. + msg = 'Attribute not found: units' + return + endif + + if (trim(units) /= trim(expected_units)) then + if (trim(expected_units) == "meters") then + if (trim(units) /= "m") then + ierr = .true. + endif + else + ierr = .true. + endif + endif + + if (ierr) then + msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units) + endif + +end subroutine check_grid_def + !> Deallocation of regridding memory subroutine end_regridding(CS) type(regridding_CS), intent(inout) :: CS !< Regridding control structure - call regridding_memory_deallocation( CS ) + deallocate( CS%coordinateResolution ) + if (allocated(CS%target_density)) deallocate( CS%target_density ) + if (allocated(CS%max_interface_depths) ) deallocate( CS%max_interface_depths ) + if (allocated(CS%max_layer_thickness) ) deallocate( CS%max_layer_thickness ) end subroutine end_regridding +!> Numeric value of interpolation_scheme corresponding to scheme name +integer function interpolation_scheme(interp_scheme) + character(len=*), intent(in) :: interp_scheme !< Name of interpolation scheme + + select case ( uppercase(trim(interp_scheme)) ) + case ("P1M_H2"); interpolation_scheme = INTERPOLATION_P1M_H2 + case ("P1M_H4"); interpolation_scheme = INTERPOLATION_P1M_H4 + case ("P1M_IH2"); interpolation_scheme = INTERPOLATION_P1M_IH4 + case ("PLM"); interpolation_scheme = INTERPOLATION_PLM + case ("PPM_H4"); interpolation_scheme = INTERPOLATION_PPM_H4 + case ("PPM_IH4"); interpolation_scheme = INTERPOLATION_PPM_IH4 + case ("P3M_IH4IH3"); interpolation_scheme = INTERPOLATION_P3M_IH4IH3 + case ("P3M_IH6IH5"); interpolation_scheme = INTERPOLATION_P3M_IH6IH5 + case ("PQM_IH4IH3"); interpolation_scheme = INTERPOLATION_PQM_IH4IH3 + case ("PQM_IH6IH5"); interpolation_scheme = INTERPOLATION_PQM_IH6IH5 + case default ; call MOM_error(FATAL, "MOM_regridding: "//& + "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").") + end select + +end function interpolation_scheme + !------------------------------------------------------------------------------ ! Dispatching regridding routine: regridding & remapping !------------------------------------------------------------------------------ @@ -876,9 +1397,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. @@ -939,7 +1459,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. @@ -954,11 +1474,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 @@ -2746,9 +3266,9 @@ function getCoordinateShortName( CS ) end function getCoordinateShortName -!> This subroutine can be used to set many of the parameters for MOM_regridding. +!> Can be used to set any of the parameters for MOM_regridding. subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, & - depth_of_time_filter_shallow, depth_of_time_filter_deep, & + interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, & nlay_ML_to_interior, fix_haloclines, halocline_filt_len, & halocline_strat_tol, integrate_downward_for_e) @@ -2756,6 +3276,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the new grid (m) real, optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid + character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates real, optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic (H units) real, optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic (H units) real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density @@ -2775,6 +3296,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri call MOM_error(FATAL,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!') CS%old_grid_weight = old_grid_weight endif + if (present(interp_scheme)) CS%interpolation_scheme = interpolation_scheme(interp_scheme) if (present(depth_of_time_filter_shallow)) CS%depth_of_time_filter_shallow = depth_of_time_filter_shallow if (present(depth_of_time_filter_deep)) CS%depth_of_time_filter_deep = depth_of_time_filter_deep if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then @@ -2797,6 +3319,14 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri end subroutine set_regrid_params +!> Returns the number of levels/layers in the regridding control structure +integer function get_regrid_size(CS) + type(regridding_CS), intent(inout) :: CS !< Regridding control structure + + get_regrid_size = CS%nk + +end function get_regrid_size + !------------------------------------------------------------------------------ ! Return coordinate-derived thicknesses for fixed coordinate systems !------------------------------------------------------------------------------ @@ -2836,44 +3366,33 @@ function getStaticThickness( CS, SSH, depth ) end function getStaticThickness -!------------------------------------------------------------------------------ -! Allocate memory for regridding -!------------------------------------------------------------------------------ -subroutine allocate_regridding( CS ) -!------------------------------------------------------------------------------ -! In this routine, we allocate the memory needed to carry out regridding -! steps in the course of the simulation. - -! For example, to compute implicit edge-value estimates, a tridiagonal system -! must be solved. We allocate the needed memory at the beginning of the -! simulation because the number of layers never changes. -!------------------------------------------------------------------------------ - - ! Arguments - type(regridding_CS), intent(inout) :: CS - +!> Parses a string and generates a dz(:) profile that goes like k**power. +subroutine dz_function1( string, dz ) + character(len=*), intent(in) :: string !< String with list of parameters in form + !! dz_min, H_total, power, precision + real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses ! Local variables - - ! Target values - allocate( CS%target_density(CS%nk+1) ) - - ! Target resolution (for fixed coordinates) - allocate( CS%coordinateResolution(CS%nk) ); CS%coordinateResolution(:) = -1.E30 - -end subroutine allocate_regridding - - -!> Deallocate memory for regridding -subroutine regridding_memory_deallocation( CS ) - type(regridding_CS), intent(inout) :: CS !< Regridding control structure - - ! Target values - deallocate( CS%target_density ) - deallocate( CS%coordinateResolution ) - if (allocated(CS%max_interface_depths) ) deallocate( CS%max_interface_depths ) - if (allocated(CS%max_layer_thickness) ) deallocate( CS%max_layer_thickness ) - -end subroutine regridding_memory_deallocation + integer :: nk, k + real :: dz_min, power, prec, H_total + + nk = size(dz) ! Number of cells + prec = -1024. + read( string, *) dz_min, H_total, power, prec + if (prec == -1024.) call MOM_error(FATAL,"dz_function1: "// & + "Problem reading FNC1: string ="//trim(string)) + ! Create profile of ( dz - dz_min ) + do k = 1, nk + dz(k) = (real(k-1)/real(nk-1))**power + enddo + dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total + dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec + dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total + dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec + dz(nk) = dz(nk) + ( H_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer + dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec + dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min + +end subroutine dz_function1 !> \namespace mom_regridding !! diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 0042efc223..89e3058e80 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -8,28 +8,50 @@ 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 = 2 + +! 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 = 7 !< 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_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 +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_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_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 module procedure coordinateUnitsS end interface +interface state_dependent + module procedure state_dependent_char + module procedure state_dependent_int +end interface + contains !> Parse a string parameter specifying the coordinate mode and @@ -40,10 +62,12 @@ 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 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 @@ -76,4 +100,27 @@ function coordinateUnitsS(string) coordinateUnitsS = coordinateUnitsI(coordMode) end function coordinateUnitsS +!> Returns true if the coordinate is dependent on the state density, returns false otherwise. +logical function state_dependent_char(string) + character(len=*), intent(in) :: string !< String to indicate coordinate mode + + state_dependent_char = state_dependent_int( coordinateMode(string) ) + +end function state_dependent_char + +!> Returns true if the coordinate is dependent on the state density, returns false otherwise. +logical function state_dependent_int(mode) + integer, intent(in) :: mode !< Coordinate mode + select case ( mode ) + case (REGRIDDING_LAYER); state_dependent_int = .true. + case (REGRIDDING_ZSTAR); state_dependent_int = .false. + case (REGRIDDING_RHO); state_dependent_int = .true. + case (REGRIDDING_SIGMA); state_dependent_int = .false. + case (REGRIDDING_HYCOM1); state_dependent_int = .true. + case (REGRIDDING_SLIGHT); state_dependent_int = .true. + case default ; call MOM_error(FATAL, "state_dependent: "//& + "Unrecognized choice of coordinate.") + end select +end function state_dependent_int + end module regrid_consts diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 77a2d2e371..fa5311b9e5 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 @@ -356,8 +356,8 @@ module MOM integer :: id_e_preale = -1 ! Diagnostics for tracer horizontal transport - integer :: id_uhtr = -1 - integer :: id_vhtr = -1 + integer :: id_uhtr = -1, id_umo = -1, id_umo_2d = 1 + integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = 1 ! The remainder provides pointers to child module control structures. type(MOM_dyn_unsplit_CS), pointer :: dyn_unsplit_CSp => NULL() @@ -485,7 +485,8 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) h ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss)) real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_predia, eta_preale - + real, dimension(SZIB_(CS%G), SZJ_(CS%G)) :: umo2d ! Diagnostics + real, dimension(SZI_(CS%G), SZJB_(CS%G)) :: vmo2d ! Diagnostics real :: tot_wt_ssh, Itot_wt_ssh, I_time_int real :: zos_area_mean, volo, ssh_ga @@ -766,7 +767,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) @@ -848,7 +849,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 @@ -969,7 +970,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) if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, & CS%visc, dt, G, GV, CS%MEKE_CSp, CS%uhtr, CS%vhtr) @@ -1116,7 +1117,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) @@ -1208,6 +1209,32 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%id_uhtr > 0) call post_data(CS%id_uhtr, CS%uhtr, CS%diag) if (CS%id_vhtr > 0) call post_data(CS%id_vhtr, CS%vhtr, CS%diag) + if (CS%id_umo_2d > 0) then + umo2d(:,:) = CS%uhtr(:,:,1) + do k = 2, nz + umo2d(:,:) = umo2d(:,:) + CS%uhtr(:,:,k) + enddo + umo2d(:,:) = umo2d(:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) + call post_data(CS%id_umo, umo2d, CS%diag) + endif + if (CS%id_umo > 0) then + ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below + CS%uhtr(:,:,:) = CS%uhtr(:,:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) + call post_data(CS%id_umo, CS%uhtr, CS%diag) + endif + if (CS%id_vmo_2d > 0) then + vmo2d(:,:) = CS%vhtr(:,:,1) + do k = 2, nz + vmo2d(:,:) = vmo2d(:,:) + CS%vhtr(:,:,k) + enddo + vmo2d(:,:) = vmo2d(:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) + call post_data(CS%id_vmo, vmo2d, CS%diag) + endif + if (CS%id_vmo > 0) then + ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below + CS%uhtr(:,:,:) = CS%uhtr(:,:,:) * ( GV%H_to_kg_m2 / CS%dt_trans ) + call post_data(CS%id_vmo, CS%vhtr, CS%diag) + endif call post_diags_TS_tendency(G,GV,CS,dtdia) @@ -2538,17 +2565,18 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo ! 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) @@ -2820,7 +2848,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_extensive=.true.) CS%id_volo = register_scalar_field('ocean_model', 'volo', Time, diag,& long_name='Total volume of liquid ocean', units='m3', & @@ -2970,7 +2998,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_extensive=.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 @@ -2979,7 +3007,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_extensive=.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, & @@ -2997,10 +3025,24 @@ subroutine register_diags(Time, G, GV, CS, ADp) ! Diagnostics related to tracer transport CS%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, Time, & - 'Accumulated zonal thickness fluxes to advect tracers', 'kg') + 'Accumulated zonal thickness fluxes to advect tracers', 'kg', & + y_cell_method='sum', v_extensive=.true.) CS%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, Time, & - 'Accumulated meridional thickness fluxes to advect tracers', 'kg') - + 'Accumulated meridional thickness fluxes to advect tracers', 'kg', & + x_cell_method='sum', v_extensive=.true.) + CS%id_umo = register_diag_field('ocean_model', 'umo', & + diag%axesCuL, Time, 'Ocean Mass X Transport', 'kg/s', & + standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.) + CS%id_vmo = register_diag_field('ocean_model', 'vmo', & + diag%axesCvL, Time, 'Ocean Mass Y Transport', 'kg/s', & + standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.) + CS%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', & + diag%axesCu1, Time, 'Ocean Mass X Transport Vertical Sum', 'kg/s', & + standard_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, 'Ocean Mass Y Transport Vertical Sum', 'kg/s', & + standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum') + end subroutine register_diags @@ -3038,7 +3080,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_extensive=.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", & @@ -3078,7 +3120,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_extensive=.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_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 9ee1d4ccde..31b7b1e889 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -78,7 +78,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 @@ -946,7 +946,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) @@ -987,7 +987,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 ee72632b08..a46eb074d8 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 @@ -283,10 +283,6 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & logical :: do_pass_Ray_uv, do_pass_kv_bbl_thick logical :: showCallTree - ! for diagnostics - real, dimension(SZIB_(G),SZJ_(G)) :: uwork2d - real, dimension(SZI_(G),SZJB_(G)) :: vwork2d - integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -846,7 +842,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) @@ -900,29 +896,6 @@ 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) - - ! depth summed zonal mass transport - if (CS%id_umo_2d > 0) then - do j=js,je ; do I=Isq,Ieq - uwork2d(i,j) = 0.0 - do k=1,nz - uwork2d(i,j) = uwork2d(i,j) + uh(i,j,k)*GV%H_to_kg_m2 - enddo - enddo ; enddo - call post_data(CS%id_umo_2d, uwork2d, CS%diag) - endif - ! depth summed merid mass transport - if (CS%id_vmo_2d > 0) then - do J=Jsq,Jeq ; do i=is,ie - vwork2d(i,j) = 0.0 - do k=1,nz - vwork2d(i,j) = vwork2d(i,j) + vh(i,j,k)*GV%H_to_kg_m2 - enddo - enddo ; enddo - call post_data(CS%id_vmo_2d, vwork2d, CS%diag) - endif if (CS%debug) then call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV) @@ -1203,23 +1176,9 @@ 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_extensive=.true.) CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, & - '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') - 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') - 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') - 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') + 'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true.) 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 1cc93403ff..4dd09df8a8 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 @@ -457,7 +457,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. @@ -695,9 +695,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_extensive=.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_extensive=.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 71a68c286c..2ff81ba218 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -651,9 +651,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_extensive=.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_extensive=.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/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 40c3487134..46694e192b 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -982,14 +982,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, & @@ -1031,7 +1031,6 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) character(len=48) :: flux_units, string integer :: z_axis, zint_axis integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nk, id_test - logical :: diag_mediator_is_using_z isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nk = G%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -1055,18 +1054,7 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) "depth-space diagnostics, or blank to disable \n"//& "depth-space output.", default="") - ! Check that the diag_mediator z-sapce remapping is not using the same module name - string = '' - call get_param(param_file, mod, "DIAG_REMAP_Z_MODULE_SUFFIX", string, & - default='_z_new', do_not_log=.true.) - diag_mediator_is_using_z = .false. - if (trim(string) == '_z') diag_mediator_is_using_z = .true. - if (len_trim(zgrid_file) > 0) then - if (diag_mediator_is_using_z) call MOM_error(FATAL, "MOM_diag_to_Z_init:"// & - "Z_OUTPUT_GRID_FILE can not be used when DIAG_REMAP_Z_MODULE_SUFFIX='_z'." // & - " Z_OUTPUT_GRID_FILE='"//trim(zgrid_file)//"'") - call get_param(param_file, mod, "INPUTDIR", in_dir, & "The directory in which input files are found.", default=".") in_dir = slasher(in_dir) @@ -1094,55 +1082,28 @@ 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) - elseif (.not. diag_mediator_is_using_z) then - - ! 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, & - '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, & - '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, & - '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, & - '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 "//& - "an appropriate depth-space target file.") endif end subroutine MOM_diag_to_Z_init @@ -1356,7 +1317,7 @@ function register_Z_diag(var_desc, CS, day, missing) character(len=48) :: units ! A variable's units. character(len=240) :: longname ! A variable's longname. character(len=8) :: hor_grid, z_grid ! Variable grid info. - type(axes_grp) :: axes + type(axes_grp), pointer :: axes call query_vardesc(var_desc, name=var_name, units=units, longname=longname, & hor_grid=hor_grid, z_grid=z_grid, caller="register_Zint_diag") @@ -1368,21 +1329,21 @@ function register_Z_diag(var_desc, CS, day, missing) case ("z") select case (hor_grid) case ("q") - axes = CS%axesBz + axes => CS%axesBz case ("h") - axes = CS%axesTz + axes => CS%axesTz case ("u") - axes = CS%axesCuz + axes => CS%axesCuz case ("v") - axes = CS%axesCvz + axes => CS%axesCvz case ("Bu") - axes = CS%axesBz + axes => CS%axesBz case ("T") - axes = CS%axesTz + axes => CS%axesTz case ("Cu") - axes = CS%axesCuz + axes => CS%axesCuz case ("Cv") - axes = CS%axesCvz + axes => CS%axesCvz case default call MOM_error(FATAL,& "register_Z_diag: unknown hor_grid component "//trim(hor_grid)) @@ -1393,7 +1354,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 @@ -1408,7 +1369,7 @@ function register_Zint_diag(var_desc, CS, day) character(len=48) :: units ! A variable's units. character(len=240) :: longname ! A variable's longname. character(len=8) :: hor_grid ! Variable grid info. - type(axes_grp) :: axes + type(axes_grp), pointer :: axes call query_vardesc(var_desc, name=var_name, units=units, longname=longname, & hor_grid=hor_grid, caller="register_Zint_diag") @@ -1421,19 +1382,19 @@ function register_Zint_diag(var_desc, CS, day) ! desired axes to register the diagnostic field for. select case (hor_grid) case ("h") - axes = CS%axesTzi + axes => CS%axesTzi case ("q") - axes = CS%axesBzi + axes => CS%axesBzi case ("u") - axes = CS%axesCuzi + axes => CS%axesCuzi case ("v") - axes = CS%axesCvzi + axes => CS%axesCvzi case default call MOM_error(FATAL,& "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 diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index ed17cf5136..6aac92b50a 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -124,6 +124,7 @@ subroutine find_obsolete_params(param_file) call obsolete_char(param_file, "REF_COMPRESS_FILE_TEMP") call obsolete_char(param_file, "REF_COMPRESS_FILE_SALT") call obsolete_char(param_file, "REF_COMPRESS_FILE_DEPTH") + call obsolete_char(param_file, "DIAG_REMAP_Z_GRID_DEF", "Use NUM_DIAG_COORDS, DIAG_COORDS and DIAG_COORD_DEF_Z") call obsolete_logical(param_file, "OLD_RESTRAT_PARAM", .false.) call obsolete_real(param_file, "ML_RESTRAT_COEF", 0.0) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 7bec64561b..6296c8397c 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -30,20 +30,23 @@ 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 +use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap +use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field +use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured +use MOM_diag_remap, only : diag_remap_get_axes_info, diag_remap_set_active +use MOM_diag_remap, only : diag_remap_diag_registration_closed +use MOM_diag_remap, only : horizontally_average_diag_field use diag_axis_mod, only : get_diag_axis_name use diag_manager_mod, only : diag_manager_init, diag_manager_end @@ -55,18 +58,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 @@ -78,10 +75,9 @@ 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 interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -98,8 +94,24 @@ module MOM_diag_mediator character(len=9) :: x_cell_method = '' !< Default nature of data representation, if axes group includes x-direction. character(len=9) :: y_cell_method = '' !< Default nature of data representation, if axes group includes y-direction. character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group includes vertical direction. + ! For remapping + integer :: nz = 0 !< Vertical dimension of diagnostic + integer :: vertical_coordinate_number = 0 !< Index of the corresponding diag_remap_ctrl for this axis group ! 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_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_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. + ! For horizontally averaged diagnositcs (applies to 2d and 3d fields only) + type(axes_grp), pointer :: xyave_axes => null() !< The associated 1d axes for horizontall area-averaged diagnostics ! 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. @@ -114,12 +126,14 @@ module MOM_diag_mediator type, private :: diag_type 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() + integer :: fms_xyave_diag_id = -1 !< For a horizontally area-averaged diagnostic. + character(32) :: debug_str = '' !< For FATAL errors and debugging. + 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. real :: conversion_factor = 0. !< A factor to multiply data by before posting to FMS, if non-zero. + logical :: v_extensive = .false. !< True for vertically extensive fields (vertically integrated). False for intensive (concentrations). end type diag_type !> The following data type a list of diagnostic fields an their variants, @@ -142,7 +156,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() @@ -150,11 +163,11 @@ module MOM_diag_mediator real, dimension(:,:), pointer :: mask2dCu => null() real, dimension(:,:), pointer :: mask2dCv => null() real, dimension(:,:,:), pointer :: mask3dTL => null() - real, dimension(:,:,:), pointer :: mask3dBuL => null() + real, dimension(:,:,:), pointer :: mask3dBL => null() real, dimension(:,:,:), pointer :: mask3dCuL => null() real, dimension(:,:,:), pointer :: mask3dCvL => null() real, dimension(:,:,:), pointer :: mask3dTi => null() - real, dimension(:,:,:), pointer :: mask3dBui => null() + real, dimension(:,:,:), pointer :: mask3dBi => null() real, dimension(:,:,:), pointer :: mask3dCui => null() real, dimension(:,:,:), pointer :: mask3dCvi => null() @@ -165,30 +178,23 @@ module MOM_diag_mediator integer :: next_free_diag_id !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 + real :: missing_value = -1.0e+34 - ! Output grid thicknesses - real, dimension(:,:,:), allocatable :: h_zoutput + !> Number of diagnostic vertical coordinates (remapped) + integer :: num_diag_coords + !> Control structure for each possible coordinate + type(diag_remap_ctrl), dimension(:), allocatable :: diag_remap_cs - ! 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 + !> Axes groups for each possible coordinate (these will all be 3D groups) + type(axes_grp), dimension(:), allocatable :: remap_axesZL, remap_axesZi + type(axes_grp), dimension(:), allocatable :: remap_axesTL, remap_axesBL, remap_axesCuL, remap_axesCvL + type(axes_grp), dimension(:), allocatable :: remap_axesTi, remap_axesBi, remap_axesCui, remap_axesCvi - !> 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__) @@ -200,7 +206,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 @@ -213,22 +219,13 @@ 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" - character(len=40) :: mod = "MOM_diag_mediator" ! This module's name. set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical - ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version, "") - if(G%symmetric) then id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', & 'q point nominal longitude', Domain2=G%Domain%mpp_domain) @@ -259,170 +256,165 @@ 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') + 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.) + x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', & + is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL) 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', & + 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', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean', v_cell_method='mean', & + is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL) 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', & + is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL) ! 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., is_interface=.true., xyave_axes=diag_cs%axesZi) + 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', is_h_point=.false.) + x_cell_method='point', y_cell_method='mean', v_cell_method='point', & + is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi) 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.) - 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='mean', y_cell_method='point', v_cell_method='point', & + is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi) ! 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', 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', is_h_point=.false.) + 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', is_h_point=.false.) + x_cell_method='mean', y_cell_method='point', is_v_point=.true.) + + if (diag_cs%num_diag_coords>0) then + allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesTL(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesBL(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesCuL(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesCvL(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesZi(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesTi(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesBi(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesCui(diag_cs%num_diag_coords)) + allocate(diag_cs%remap_axesCvi(diag_cs%num_diag_coords)) + endif + + do i=1, diag_cs%num_diag_coords + ! For each possible diagnostic coordinate + call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, param_file) + + ! This vertical coordinate has been configured so can be used. + if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then + + ! This fetches the 1D-axis id for layers and interfaces and overwrite + ! id_zl and id_zi from above. It also returns the number of layers. + call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zL, id_zi) + + ! Axes for z layers + call define_axes_group(diag_cs, (/ id_zL /), diag_cs%remap_axesZL(i), & + nz=nz, vertical_coordinate_number=i, & + v_cell_method='mean', & + is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.) + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%remap_axesTL(i), & + nz=nz, vertical_coordinate_number=i, & + 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., & + xyave_axes=diag_cs%remap_axesZL(i)) + + !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBL + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%remap_axesBL(i), & + nz=nz, vertical_coordinate_number=i, & + x_cell_method='point', y_cell_method='point', v_cell_method='mean', & + is_q_point=.true., is_layer=.true., is_native=.false.) + + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%remap_axesCuL(i), & + nz=nz, vertical_coordinate_number=i, & + 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., & + xyave_axes=diag_cs%remap_axesZL(i)) + + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%remap_axesCvL(i), & + nz=nz, vertical_coordinate_number=i, & + 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., & + xyave_axes=diag_cs%remap_axesZL(i)) + + ! Axes for z interfaces + call define_axes_group(diag_cs, (/ id_zi /), diag_cs%remap_axesZi(i), & + nz=nz, vertical_coordinate_number=i, & + v_cell_method='point', & + is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true.) + call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), & + nz=nz, vertical_coordinate_number=i, & + x_cell_method='mean', y_cell_method='mean', v_cell_method='point', & + is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., & + xyave_axes=diag_cs%remap_axesZi(i)) + + !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi + call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), & + nz=nz, vertical_coordinate_number=i, & + x_cell_method='point', y_cell_method='point', v_cell_method='point', & + is_q_point=.true., is_interface=.true., is_native=.false.) + + call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%remap_axesCui(i), & + nz=nz, vertical_coordinate_number=i, & + x_cell_method='point', y_cell_method='mean', v_cell_method='point', & + is_u_point=.true., is_interface=.true., is_native=.false., & + needs_interpolating=.true., & xyave_axes=diag_cs%remap_axesZi(i)) + + call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%remap_axesCvi(i), & + nz=nz, vertical_coordinate_number=i, & + x_cell_method='mean', y_cell_method='point', v_cell_method='point', & + is_v_point=.true., is_interface=.true., is_native=.false., & + needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i)) + endif + 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, diag_cs%num_diag_coords + 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, diag_cs%num_diag_coords + 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 @@ -434,50 +426,38 @@ 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) +subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_number, & + 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_interpolating, & + xyave_axes) 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 + integer, optional, intent(in) :: nz !< Number of layers in this diagnostic grid + integer, optional, intent(in) :: vertical_coordinate_number !< Index number for vertical coordinate character(len=*), optional, intent(in) :: x_cell_method !< A x-direction cell method used to construct the "cell_methods" attribute in CF convention 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 field + !! that must be remapped to these axes. Used for rank>2. + 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. + type(axes_grp), optional, target :: xyave_axes !< The corresponding axes group for horizontally area-average diagnostics ! Local variables integer :: n + n = size(handles) - if (n<1 .or. n>3) call MOM_error(FATAL,"define_axes_group: wrong size for list of handles!") + if (n<1 .or. n>3) call MOM_error(FATAL, "define_axes_group: wrong size for list of handles!") allocate( axes%handles(n) ) axes%id = i2s(handles, n) ! Identifying string axes%rank = n @@ -504,7 +484,19 @@ subroutine define_axes_group(diag_cs, handles, axes, & else axes%v_cell_method = '' endif + if (present(nz)) axes%nz = nz + if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number 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_interpolating)) axes%needs_interpolating = needs_interpolating + if (present(xyave_axes)) axes%xyave_axes => xyave_axes + end subroutine define_axes_group subroutine set_diag_mediator_grid(G, diag_cs) @@ -522,12 +514,11 @@ subroutine set_diag_mediator_grid(G, diag_cs) end subroutine set_diag_mediator_grid -subroutine post_data_0d(diag_field_id, field, diag_cs, is_static, mask) +subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) integer, intent(in) :: diag_field_id real, intent(in) :: field type(diag_ctrl), target, intent(in) :: diag_cs logical, optional, intent(in) :: is_static - real, optional, intent(in) :: mask(:,:) ! Arguments: ! (in) diag_field_id - the id for an output variable returned by a @@ -615,6 +606,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') @@ -677,12 +669,6 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) call MOM_error(FATAL,"post_data_2d_low: peculiar size in j-direction") endif - if (present(mask)) then - call check_field_and_mask_shape_2d(diag, field, mask) - elseif ((diag_cs%ave_enabled) .and. associated(diag%mask2d)) then - call check_field_and_mask_shape_2d(diag, field, diag%mask2d) - endif - if (diag%conversion_factor/=0.) then allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2) ) ) locfield(isv:iev,jsv:jev) = field(isv:iev,jsv:jev) * diag%conversion_factor @@ -692,6 +678,8 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) if (is_stat) then if (present(mask)) then + call assert(size(locfield) == size(mask), & + 'post_data_2d_low is_stat: mask size mismatch: '//diag%debug_str) used = send_data(diag%fms_diag_id, locfield, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=mask) !elseif(associated(diag%mask2d)) then @@ -703,6 +691,8 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) endif elseif (diag_cs%ave_enabled) then if (present(mask)) then + call assert(size(locfield) == size(mask), & + 'post_data_2d_low: mask size mismatch: '//diag%debug_str) used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & weight=diag_cs%time_int, rmask=mask) @@ -737,38 +727,95 @@ 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 + logical :: staggered_in_x, staggered_in_y 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, & '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') + + staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point + staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point - if (associated(diag%remap_axes)) then + if (diag%v_extensive .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_remap>0) call cpu_clock_begin(id_clock_diag_remap) + allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, diag_cs%h, staggered_in_x, staggered_in_y, & + diag%mask3d, diag_cs%missing_value, field, remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_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%axes%nz)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + 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) + 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.") + endif + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) + call diag_remap_do_remap(diag_cs%diag_remap_cs( & + diag%axes%vertical_coordinate_number), & + diag_cs%G, diag_cs%h, staggered_in_x, staggered_in_y, & + diag%mask3d, diag_cs%missing_value, field, remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_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%axes%nz)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + 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) + 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)) - 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 (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz+1)) + call vertically_interpolate_diag_field(diag_cs%diag_remap_cs( & + diag%axes%vertical_coordinate_number), & + diag_cs%G, diag_cs%h, staggered_in_x, staggered_in_y, & + diag%mask3d, diag_cs%missing_value, field, remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_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)) + mask=diag%mask3d(:,:,:diag%axes%nz+1)) 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) + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) deallocate(remapped_field) - if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) else call post_data_3d_low(diag, field, diag_cs, is_static, mask) endif @@ -778,178 +825,6 @@ 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. -!! \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 - 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 -#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 - -!> 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') - - ! 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 - 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 - 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) - -end subroutine diag_update_target_grids - subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, target, intent(in) :: field(:,:,:) @@ -1001,12 +876,6 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) call MOM_error(FATAL,"post_data_3d_low: peculiar size in j-direction") endif - if (present(mask)) then - call check_field_and_mask_shape_3d(diag, field, mask) - elseif ((diag_cs%ave_enabled) .and. associated(diag%mask3d)) then - call check_field_and_mask_shape_3d(diag, field, diag%mask3d) - endif - if (diag%conversion_factor/=0.) then allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2), & lbound(field,3):ubound(field,3) ) ) @@ -1015,36 +884,73 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) locfield => field endif - if (is_stat) then - if (present(mask)) then - used = send_data(diag%fms_diag_id, locfield, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=mask) - !elseif(associated(diag%mask3d)) then - ! used = send_data(diag_field_id, locfield, & - ! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask3d) - else - used = send_data(diag%fms_diag_id, locfield, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) - endif - elseif (diag_cs%ave_enabled) then - if (present(mask)) then - used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int, rmask=mask) - elseif(associated(diag%mask3d)) then - used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int, rmask=diag%mask3d) - else - used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int) + if (diag%fms_diag_id>0) then + if (is_stat) then + if (present(mask)) then + call assert(size(locfield) == size(mask), & + 'post_data_3d_low is_stat: mask size mismatch: '//diag%debug_str) + used = send_data(diag%fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=mask) + !elseif(associated(diag%mask3d)) then + ! used = send_data(diag_field_id, locfield, & + ! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask3d) + else + used = send_data(diag%fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) + endif + elseif (diag_cs%ave_enabled) then + if (present(mask)) then + call assert(size(locfield) == size(mask), & + 'post_data_3d_low: mask size mismatch: '//diag%debug_str) + used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int, rmask=mask) + elseif(associated(diag%mask3d)) then + call assert(size(locfield) == size(diag%mask3d), & + 'post_data_3d_low: mask3d size mismatch: '//diag%debug_str) + used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int, rmask=diag%mask3d) + else + used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int) + endif endif endif + if (diag%fms_xyave_diag_id>0) then + call post_xy_average(diag_cs, diag, locfield) + endif if (diag%conversion_factor/=0.) deallocate( locfield ) end subroutine post_data_3d_low +!> Post the horizontally area-averaged diagnostic +subroutine post_xy_average(diag_cs, diag, field) + type(diag_type), intent(in) :: diag !< This diagnostic + real, target, intent(in) :: field(:,:,:) !< Diagnostic field + type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics mediator control structure + ! Local variable + integer :: nk, i, j, k + real, dimension(size(field,3)) :: averaged_field + logical :: staggered_in_x, staggered_in_y, used + + if (.not.diag%axes%is_h_point) then + call MOM_error(FATAL, 'post_xy_average: Horizontally averaged diagnostic not implemented yet.') + endif + + staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point + staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point + if (diag_cs%ave_enabled) then + call horizontally_average_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, staggered_in_x, staggered_in_y, & + diag%axes%is_layer, diag%v_extensive, & + diag_cs%missing_value, field, averaged_field) + used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, weight=diag_cs%time_int) + endif + +end subroutine post_xy_average + subroutine enable_averaging(time_int_in, time_end_in, diag_cs) real, intent(in) :: time_int_in type(time_type), intent(in) :: time_end_in @@ -1108,12 +1014,12 @@ 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) - integer :: register_diag_field !< An integer handle for a diagnostic array. + x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive) 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 @@ -1140,49 +1046,176 @@ 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_extensive !< True for vertically extensive fields (vertically integrated). Default/absent for intensive. ! 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(axes_grp), pointer :: remap_axes => null() + integer :: dm_id, i + character(len=256) :: new_module_name + 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, v_extensive=v_extensive) + + ! For each diagnostic coordinate register the diagnostic again under a different module name + do i=1,diag_cs%num_diag_coords + new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix) + + ! Register diagnostics remapped to z vertical coordinate + if (axes%rank == 3) then + remap_axes => null() + if ((axes%id .eq. diag_cs%axesTL%id)) then + remap_axes => diag_cs%remap_axesTL(i) + elseif(axes%id .eq. diag_cs%axesBL%id) then + remap_axes => diag_cs%remap_axesBL(i) + elseif(axes%id .eq. diag_cs%axesCuL%id ) then + remap_axes => diag_cs%remap_axesCuL(i) + elseif(axes%id .eq. diag_cs%axesCvL%id) then + remap_axes => diag_cs%remap_axesCvL(i) + elseif(axes%id .eq. diag_cs%axesTi%id) then + remap_axes => diag_cs%remap_axesTi(i) + elseif(axes%id .eq. diag_cs%axesBi%id) then + remap_axes => diag_cs%remap_axesBi(i) + elseif(axes%id .eq. diag_cs%axesCui%id ) then + remap_axes => diag_cs%remap_axesCui(i) + elseif(axes%id .eq. diag_cs%axesCvi%id) then + remap_axes => diag_cs%remap_axesCvi(i) + endif + ! 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, new_module_name, 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, v_extensive=v_extensive) + if (active) then + call diag_remap_set_active(diag_cs%diag_remap_cs(i)) + endif + endif ! remap_axes%needs_remapping + endif ! associated(remap_axes) + endif ! axes%rank == 3 + enddo ! i + + 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 'dm_id' +!! after calling register_diag_field_expand_axes() for both native and CMOR variants of the field. +logical function register_diag_field_expand_cmor(dm_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, v_extensive) + integer, intent(inout) :: dm_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 + 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 + logical, optional, intent(in) :: v_extensive !< True for vertically extensive fields (vertically integrated). Default/absent for intensive. + ! Local variables + real :: MOM_missing_value + type(diag_ctrl), pointer :: diag_cs + type(diag_type), pointer :: this_diag => null() + integer :: fms_id, fms_xyave_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) - 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 + cell_methods, x_cell_method, y_cell_method, v_cell_method, & + v_extensive=v_extensive) 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(fms_id>0, module_name, field_name, cm_string, & msg, diag_CS, long_name, units, standard_name) endif + ! Associated horizontally area-averaged diagnostic + fms_xyave_id = DIAG_FIELD_NOT_FOUND + if (associated(axes%xyave_axes)) then + fms_xyave_id = register_diag_field_expand_axes(module_name, trim(field_name)//'_xyave', & + axes%xyave_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_xyave_id, axes%xyave_axes, cm_string, & + cell_methods, v_cell_method, v_extensive=v_extensive) + 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)//'_xyave"' + call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, & + msg, diag_CS, long_name, units, standard_name) + endif + endif + this_diag => null() + if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then + call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg) + this_diag%fms_xyave_diag_id = fms_xyave_id + + if (present(v_extensive)) this_diag%v_extensive = v_extensive + if (present(conversion)) this_diag%conversion_factor = conversion + register_diag_field_expand_cmor = .true. + endif - ! For the CMOR variation of the both the native and _z diagnostics + ! For the CMOR variation of the above diagnostic if (present(cmor_field_name)) then ! Fallback values for strings set to "NULL" posted_cmor_units = "not provided" ! @@ -1199,128 +1232,54 @@ function register_diag_field(module_name, field_name, axes, init_time, & 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) + 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 + cell_methods, x_cell_method, y_cell_method, v_cell_method, & + v_extensive=v_extensive) 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(fms_id>0, 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 + ! Associated horizontally area-averaged diagnostic + fms_xyave_id = DIAG_FIELD_NOT_FOUND + if (associated(axes%xyave_axes)) then + fms_xyave_id = register_diag_field_expand_axes(module_name, trim(cmor_field_name)//'_xyave', & + axes%xyave_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_xyave_id, axes%xyave_axes, cm_string, & + cell_methods, v_cell_method, v_extensive=v_extensive) 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, & + msg = 'native name is "'//trim(field_name)//'_xyave"' + call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', cm_string, & + msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif endif + this_diag => null() + if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then + call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg) + this_diag%fms_xyave_diag_id = fms_xyave_id + + if (present(v_extensive)) this_diag%v_extensive = v_extensive + if (present(conversion)) this_diag%conversion_factor = conversion + register_diag_field_expand_cmor = .true. + 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) @@ -1382,8 +1341,33 @@ integer function register_diag_field_expand_axes(module_name, field_name, axes, end function register_diag_field_expand_axes +!> Create a diagnostic type and attached to list +subroutine add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg) + type(diag_ctrl), pointer :: diag_cs !< Diagnostics mediator control structure + integer, intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group + integer, intent(in) :: fms_id !< The FMS diag_manager ID for this diagnostic + type(diag_type), pointer :: this_diag !< This diagnostic + type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes for this field + 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 diagnostic + character(len=*), intent(in) :: msg !< Message for errors + + ! If the diagnostic is needed obtain a diag_mediator ID (if needed) + if (dm_id == -1) dm_id = get_new_diag_id(diag_cs) + ! Create a new diag_type to store links in + call alloc_diag_with_id(dm_id, diag_cs, this_diag) + call assert(associated(this_diag), trim(msg)//': diag_type allocation failed') + ! Record FMS id, masks and conversion factor, in diag_type + this_diag%fms_diag_id = fms_id + this_diag%debug_str = trim(module_name)//"-"//trim(field_name) + call set_diag_mask(this_diag, diag_cs, axes) + this_diag%axes => axes + +end subroutine add_diag_to_list + !> 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_extensive) 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 @@ -1393,12 +1377,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_extensive !< True for vertically extensive fields (vertically integrated). Default/absent for intensive. ! 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_extensive)) then call MOM_error(FATAL, "attach_cell_methods: " // & 'Individual direction cell method was specified along with a "cell_methods" string.') endif @@ -1410,39 +1396,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_extensive)) call MOM_error(FATAL, "attach_cell_methods: " // & + 'Vertical cell method was specified along with the vertically extensive 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_extensive)) 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 @@ -1450,7 +1446,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 @@ -1459,8 +1455,8 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y end subroutine attach_cell_methods function register_scalar_field(module_name, field_name, init_time, diag_cs, & - long_name, units, missing_value, range, mask_variant, standard_name, & - verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & + long_name, units, missing_value, range, standard_name, & + do_not_log, err_msg, interp_method, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name) integer :: register_scalar_field character(len=*), intent(in) :: module_name, field_name @@ -1468,10 +1464,9 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & type(diag_ctrl), intent(inout) :: diag_cs character(len=*), optional, intent(in) :: long_name, units, standard_name real, optional, intent(in) :: missing_value, range(2) - logical, optional, intent(in) :: mask_variant, verbose, do_not_log + logical, optional, intent(in) :: do_not_log character(len=*), optional, intent(out):: err_msg character(len=*), optional, intent(in) :: interp_method - integer, optional, intent(in) :: tile_count character(len=*), optional, intent(in) :: cmor_field_name, cmor_long_name character(len=*), optional, intent(in) :: cmor_units, cmor_standard_name @@ -1488,7 +1483,6 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & ! Following params have yet to be used in MOM. ! (in,opt) range - valid range of a variable - ! (in,opt) mask_variant - If true a logical mask must be provided with post_data calls ! (in,opt) verbose - If true, FMS is verbosed ! (in,opt) do_not_log - If true, do not log something ! (out,opt) err_msg - character string into which an error message might be placed @@ -1496,14 +1490,14 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & ! (in,opt) tile_count - no clue real :: MOM_missing_value - integer :: primary_id, fms_id + integer :: dm_id, fms_id type(diag_type), pointer :: diag => null(), cmor_diag => null() character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name MOM_missing_value = diag_cs%missing_value if(present(missing_value)) MOM_missing_value = missing_value - primary_id = -1 + dm_id = -1 diag => null() cmor_diag => null() @@ -1511,11 +1505,11 @@ 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) - call alloc_diag_with_id(primary_id, diag_cs, diag) + dm_id = get_new_diag_id(diag_cs) + call alloc_diag_with_id(dm_id, diag_cs, diag) call assert(associated(diag), 'register_scalar_field: diag allocation failed') diag%fms_diag_id = fms_id - diag%debug_str = trim(field_name) + diag%debug_str = trim(module_name)//"-"//trim(field_name) endif if (present(cmor_field_name)) then @@ -1540,12 +1534,12 @@ 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 (dm_id == -1) then + dm_id = get_new_diag_id(diag_cs) endif - call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) + call alloc_diag_with_id(dm_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id - cmor_diag%debug_str = trim(cmor_field_name) + cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name) endif endif @@ -1560,7 +1554,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & endif endif - register_scalar_field = primary_id + register_scalar_field = dm_id end function register_scalar_field @@ -1601,14 +1595,14 @@ function register_static_field(module_name, field_name, axes, & real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs type(diag_type), pointer :: diag => null(), cmor_diag => null() - integer :: primary_id, fms_id, cmor_id + integer :: dm_id, fms_id, cmor_id character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_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 + dm_id = -1 diag => null() cmor_diag => null() @@ -1618,11 +1612,11 @@ 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) - call alloc_diag_with_id(primary_id, diag_cs, diag) + dm_id = get_new_diag_id(diag_cs) + call alloc_diag_with_id(dm_id, diag_cs, diag) call assert(associated(diag), 'register_static_field: diag allocation failed') diag%fms_diag_id = fms_id - diag%debug_str = trim(field_name) + diag%debug_str = trim(module_name)//"-"//trim(field_name) endif if (present(cmor_field_name)) then @@ -1648,12 +1642,12 @@ 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 (dm_id == -1) then + dm_id = get_new_diag_id(diag_cs) endif - call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) + call alloc_diag_with_id(dm_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id - cmor_diag%debug_str = trim(cmor_field_name) + cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name) endif endif @@ -1668,7 +1662,7 @@ function register_static_field(module_name, field_name, axes, & endif endif - register_static_field = primary_id + register_static_field = dm_id end function register_static_field @@ -1692,14 +1686,14 @@ function ocean_register_diag(var_desc, G, diag_CS, day) integer :: ocean_register_diag !< An integer handle to this diagnostic. type(vardesc), intent(in) :: var_desc !< The vardesc type describing the diagnostic type(ocean_grid_type), intent(in) :: G !< The ocean's grid type - type(diag_ctrl), intent(in) :: diag_CS !< The diagnotic control structure + type(diag_ctrl), intent(in), target :: diag_CS !< The diagnotic control structure type(time_type), intent(in) :: day !< The current model time character(len=64) :: var_name ! A variable's name. character(len=48) :: units ! A variable's units. character(len=240) :: longname ! A variable's longname. character(len=8) :: hor_grid, z_grid ! Variable grid info. - type(axes_grp) :: axes + type(axes_grp), pointer :: axes call query_vardesc(var_desc, units=units, longname=longname, hor_grid=hor_grid, & z_grid=z_grid, caller="ocean_register_diag") @@ -1711,23 +1705,23 @@ function ocean_register_diag(var_desc, G, diag_CS, day) case ("L") select case (hor_grid) case ("q") - axes = diag_cs%axesBL + axes => diag_cs%axesBL case ("h") - axes = diag_cs%axesTL + axes => diag_cs%axesTL case ("u") - axes = diag_cs%axesCuL + axes => diag_cs%axesCuL case ("v") - axes = diag_cs%axesCvL + axes => diag_cs%axesCvL case ("Bu") - axes = diag_cs%axesBL + axes => diag_cs%axesBL case ("T") - axes = diag_cs%axesTL + axes => diag_cs%axesTL case ("Cu") - axes = diag_cs%axesCuL + axes => diag_cs%axesCuL case ("Cv") - axes = diag_cs%axesCvL + axes => diag_cs%axesCvL case ("z") - axes = diag_cs%axeszL + axes => diag_cs%axeszL case default call MOM_error(FATAL, "ocean_register_diag: " // & "unknown hor_grid component "//trim(hor_grid)) @@ -1736,23 +1730,23 @@ function ocean_register_diag(var_desc, G, diag_CS, day) case ("i") select case (hor_grid) case ("q") - axes = diag_cs%axesBi + axes => diag_cs%axesBi case ("h") - axes = diag_cs%axesTi + axes => diag_cs%axesTi case ("u") - axes = diag_cs%axesCui + axes => diag_cs%axesCui case ("v") - axes = diag_cs%axesCvi + axes => diag_cs%axesCvi case ("Bu") - axes = diag_cs%axesBi + axes => diag_cs%axesBi case ("T") - axes = diag_cs%axesTi + axes => diag_cs%axesTi case ("Cu") - axes = diag_cs%axesCui + axes => diag_cs%axesCui case ("Cv") - axes = diag_cs%axesCvi + axes => diag_cs%axesCvi case ("z") - axes = diag_cs%axeszi + axes => diag_cs%axeszi case default call MOM_error(FATAL, "ocean_register_diag: " // & "unknown hor_grid component "//trim(hor_grid)) @@ -1761,21 +1755,21 @@ function ocean_register_diag(var_desc, G, diag_CS, day) case ("1") select case (hor_grid) case ("q") - axes = diag_cs%axesB1 + axes => diag_cs%axesB1 case ("h") - axes = diag_cs%axesT1 + axes => diag_cs%axesT1 case ("u") - axes = diag_cs%axesCu1 + axes => diag_cs%axesCu1 case ("v") - axes = diag_cs%axesCv1 + axes => diag_cs%axesCv1 case ("Bu") - axes = diag_cs%axesB1 + axes => diag_cs%axesB1 case ("T") - axes = diag_cs%axesT1 + axes => diag_cs%axesT1 case ("Cu") - axes = diag_cs%axesCu1 + axes => diag_cs%axesCu1 case ("Cv") - axes = diag_cs%axesCv1 + axes => diag_cs%axesCv1 case default call MOM_error(FATAL, "ocean_register_diag: " // & "unknown hor_grid component "//trim(hor_grid)) @@ -1807,7 +1801,7 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) type(diag_ctrl), intent(inout) :: diag_cs !< A pointer to a type with many variables !! used for diagnostics character(len=*), optional, intent(in) :: doc_file_dir !< A directory in which to create the - !! file + !! file ! This subroutine initializes the diag_mediator and the diag_manager. ! The grid type should have its dimensions set by this point, but it @@ -1816,10 +1810,13 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) logical :: opened, new_file character(len=8) :: this_pe character(len=240) :: doc_file, doc_file_dflt, doc_path + character(len=240), allocatable :: diag_coords(:) +! This include declares and sets the variable "version". +#include "version_variable.h" 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) @@ -1829,14 +1826,43 @@ 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 + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mod, version, "") + + call get_param(param_file, mod, 'NUM_DIAG_COORDS', diag_cs%num_diag_coords, & + 'The number of diagnostic vertical coordinates to use.\n'//& + 'For each coordinate, an entry in DIAG_COORDS must be provided.', & + default=1) + if (diag_cs%num_diag_coords>0) then + allocate(diag_coords(diag_cs%num_diag_coords)) + if (diag_cs%num_diag_coords==1) then ! The default is to provide just one instance of Z* + call get_param(param_file, mod, 'DIAG_COORDS', diag_coords, & + 'A list of string tuples associating diag_table modules to\n'//& + 'a coordinate definition used for diagnostics. Each string\n'//& + 'is of the form "MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME".', & + default='z Z ZSTAR') + else ! If using more than 1 diagnostic coordinate, all must be explicitly defined + call get_param(param_file, mod, 'DIAG_COORDS', diag_coords, & + 'A list of string tuples associating diag_table modules to\n'//& + 'a coordinate definition used for diagnostics. Each string\n'//& + 'is of the form "MODULE_SUFFIX,PARAMETER_SUFFIX,COORDINATE_NAME".', & + fail_if_missing=.true.) + endif + allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords)) + ! Initialize each diagnostic vertical coordinate + do i=1, diag_cs%num_diag_coords + call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i)) + enddo + deallocate(diag_coords) + endif + + ! 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 @@ -1887,19 +1913,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 + ! Local variables + integer :: i + + if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) + + do i=1, diag_cs%num_diag_coords + 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 + +#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. @@ -1915,22 +1973,22 @@ subroutine diag_masks_set(G, nz, missing_value, diag_cs) diag_cs%mask2dCu=> G%mask2dCu diag_cs%mask2dCv=> G%mask2dCv allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:nz)) - allocate(diag_cs%mask3dBuL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz)) + allocate(diag_cs%mask3dBL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz)) allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:nz)) allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:nz)) do k=1,nz diag_cs%mask3dTL(:,:,k) = diag_cs%mask2dT (:,:) - diag_cs%mask3dBuL(:,:,k) = diag_cs%mask2dBu(:,:) + diag_cs%mask3dBL(:,:,k) = diag_cs%mask2dBu(:,:) diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:) diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:) enddo allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:nz+1)) - allocate(diag_cs%mask3dBui(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz+1)) + allocate(diag_cs%mask3dBi(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz+1)) allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:nz+1)) allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:nz+1)) do k=1,nz+1 diag_cs%mask3dTi(:,:,k) = diag_cs%mask2dT (:,:) - diag_cs%mask3dBui(:,:,k) = diag_cs%mask2dBu(:,:) + diag_cs%mask3dBi(:,:,k) = diag_cs%mask2dBu(:,:) diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:) diag_cs%mask3dCvi(:,:,k) = diag_cs%mask2dCv(:,:) enddo @@ -1942,10 +2000,16 @@ end subroutine diag_masks_set subroutine diag_mediator_close_registration(diag_CS) type(diag_ctrl), intent(inout) :: diag_CS + integer :: i + if (diag_CS%doc_unit > -1) then close(diag_CS%doc_unit) ; diag_CS%doc_unit = -2 endif + do i=1, diag_cs%num_diag_coords + call diag_remap_diag_registration_closed(diag_cs%diag_remap_cs(i)) + enddo + end subroutine diag_mediator_close_registration subroutine diag_mediator_end(time, diag_CS, end_diag_manager) @@ -1953,22 +2017,25 @@ 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, diag_cs%num_diag_coords + call diag_remap_end(diag_cs%diag_remap_cs(i)) + enddo - if (allocated(diag_cs%h_zoutput)) deallocate(diag_cs%h_zoutput) deallocate(diag_cs%mask3dTL) - deallocate(diag_cs%mask3dBuL) + deallocate(diag_cs%mask3dBL) deallocate(diag_cs%mask3dCuL) deallocate(diag_cs%mask3dCvL) deallocate(diag_cs%mask3dTi) - deallocate(diag_cs%mask3dBui) + deallocate(diag_cs%mask3dBi) deallocate(diag_cs%mask3dCui) deallocate(diag_cs%mask3dCvi) @@ -2002,150 +2069,58 @@ 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() if (axes%rank .eq. 3) then - if ((axes%id .eq. diag_cs%axesTL%id)) then - diag%mask3d => diag_cs%mask3dTL - elseif(axes%id .eq. diag_cs%axesBL%id) then - diag%mask3d => diag_cs%mask3dBuL - elseif(axes%id .eq. diag_cs%axesCuL%id ) then - diag%mask3d => diag_cs%mask3dCuL - elseif(axes%id .eq. diag_cs%axesCvL%id) then - diag%mask3d => diag_cs%mask3dCvL - elseif(axes%id .eq. diag_cs%axesTi%id) then - diag%mask3d => diag_cs%mask3dTi - elseif(axes%id .eq. diag_cs%axesBi%id) then - diag%mask3d => diag_cs%mask3dBui - elseif(axes%id .eq. diag_cs%axesCui%id ) then - diag%mask3d => diag_cs%mask3dCui - elseif(axes%id .eq. diag_cs%axesCvi%id) then - diag%mask3d => diag_cs%mask3dCvi + if (axes%is_layer) then + if (axes%is_h_point) then + diag%mask3d => diag_cs%mask3dTL + elseif (axes%is_q_point) then + diag%mask3d => diag_cs%mask3dBL + elseif (axes%is_u_point) then + diag%mask3d => diag_cs%mask3dCuL + elseif (axes%is_v_point) then + diag%mask3d => diag_cs%mask3dCvL + endif + elseif (axes%is_interface) then + if (axes%is_h_point) then + diag%mask3d => diag_cs%mask3dTi + elseif (axes%is_q_point) then + diag%mask3d => diag_cs%mask3dBi + elseif (axes%is_u_point) then + diag%mask3d => diag_cs%mask3dCui + elseif (axes%is_v_point) then + diag%mask3d => diag_cs%mask3dCvi + endif endif - !call assert(associated(diag%mask3d), "MOM_diag_mediator.F90: Invalid 3d axes id") + !call assert(associated(diag%mask3d), "set_diag_mask: Invalid 3d axes id."// & + ! " diag:"//diag%debug_str) elseif(axes%rank .eq. 2) then - if (axes%id .eq. diag_cs%axesT1%id) then - diag%mask2d => diag_cs%mask2dT - elseif(axes%id .eq. diag_cs%axesB1%id) then - diag%mask2d => diag_cs%mask2dBu - elseif(axes%id .eq. diag_cs%axesCu1%id) then + + if (axes%is_h_point) then + diag%mask2d => diag_cs%mask2dT + elseif (axes%is_q_point) then + diag%mask2d => diag_cs%mask2dBu + elseif (axes%is_u_point) then diag%mask2d => diag_cs%mask2dCu - elseif(axes%id .eq. diag_cs%axesCv1%id) then + elseif (axes%is_v_point) then diag%mask2d => diag_cs%mask2dCv endif - !call assert(associated(diag%mask2d), "MOM_diag_mediator.F90: Invalid 2d axes id") + !call assert(associated(diag%mask2d), "set_diag_mask.F90: Invalid 2d axes id."// & + ! " diag:"//diag%debug_str) 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 - - 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 @@ -2183,11 +2158,12 @@ 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() diag%conversion_factor = 0. + end subroutine initialize_diag_type ! Make a new diagnostic. Either use memory which is in the array of 'primary' @@ -2250,47 +2226,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..120812e19a --- /dev/null +++ b/src/framework/MOM_diag_remap.F90 @@ -0,0 +1,667 @@ +!> This module is used for runtime remapping of diagnostics to z star, sigma and +!! rho vertical coordinates. It defines the diag_remap_ctrl type which +!! represents a remapping of diagnostics to a particular vertical coordinate. +!! The module is It is used by the diag mediator module in the following way: +!! 1) _init() is called to initialise a diag_remap_ctrl instance. +!! 2) _configure_axes() is called to read the configuration file and set up the +!! vertical coordinate / axes definitions. +!! 3) _get_axes_info() returns information needed for the diag mediator to +!! define new axes for the remapped diagnostics. +!! 4) _update() is called periodically (whenever h, T or S change) to either +!! create or update the target remapping grids. +!! 5) _do_remap() is called from within a diag post() to do the remapping before +!! the diagnostic is written out. + +module MOM_diag_remap + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_coms, only : sum_across_PEs +use MOM_error_handler, only : MOM_error, FATAL, assert +use MOM_diag_vkernels, only : interpolate_column, reintegrate_column +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 +use MOM_regridding, only : build_zstar_column, build_rho_column, build_sigma_column +use MOM_regridding, only : set_regrid_params, get_regrid_size, uniformResolution +use MOM_regridding, only : getCoordinateInterfaces +use regrid_consts, only : coordinateMode + +use diag_axis_mod, only : get_diag_axis_name +use diag_manager_mod, only : diag_axis_init + +implicit none ; private + +public diag_remap_ctrl +public diag_remap_init, diag_remap_end, diag_remap_update, diag_remap_do_remap +public diag_remap_configure_axes, diag_remap_axes_configured +public diag_remap_get_axes_info, diag_remap_set_active +public diag_remap_diag_registration_closed +public vertically_reintegrate_diag_field +public vertically_interpolate_diag_field +public horizontally_average_diag_field + +!> This type represents remapping of diagnostics to a particular vertical +!! coordinate. +!! There is one of these types for each vertical coordinate. The vertical axes +!! of a diagnostic will reference an instance of this type indicating how (or +!! if) the diagnostic should be vertically remapped when being posted. +type :: diag_remap_ctrl + logical :: configured = .false. !< Whether vertical coordinate has been configured + logical :: initialized = .false. !< Whether remappping initialized + logical :: used = .false. !< Whether this coordinate actually gets used. + integer :: vertical_coord = 0 !< The vertical coordinate that we remap to + character(len=10) :: vertical_coord_name ='' !< The coordinate name as understood by ALE + character(len=16) :: diag_coord_name = '' !< A name for the purpose of run-time parameters + character(len=8) :: diag_module_suffix = '' !< The suffix for the module to appear in diag_table + type(remapping_CS) :: remap_cs !< Remapping control structure use for this axes + type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordiantes for this axes + integer :: nz = 0 !< Number of vertical levels used for remapping + real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses + real, dimension(:), allocatable :: dz !< Nominal layer thicknesses + integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces + integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers +end type diag_remap_ctrl + +contains + +!> Initialize a diagnostic remapping type with the given vertical coordinate. +subroutine diag_remap_init(remap_cs, coord_tuple) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure + character(len=*), intent(in) :: coord_tuple !< A string in form of + !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME + + remap_cs%diag_module_suffix = trim(extractWord(coord_tuple, 1)) + remap_cs%diag_coord_name = trim(extractWord(coord_tuple, 2)) + remap_cs%vertical_coord_name = trim(extractWord(coord_tuple, 3)) + remap_cs%vertical_coord = coordinateMode(remap_cs%vertical_coord_name) + remap_cs%configured = .false. + remap_cs%initialized = .false. + remap_cs%used = .false. + remap_cs%nz = 0 + +end subroutine diag_remap_init + +!> De-init a diagnostic remapping type. +!! Free allocated memory. +subroutine diag_remap_end(remap_cs) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure + + if (allocated(remap_cs%h)) deallocate(remap_cs%h) + if (allocated(remap_cs%dz)) deallocate(remap_cs%dz) + remap_cs%configured = .false. + remap_cs%initialized = .false. + remap_cs%used = .false. + remap_cs%nz = 0 + +end subroutine diag_remap_end + +!> Inform that all diagnostics have been registered. +!! If _set_active() has not been called on the remapping control structure +!! will be disabled. This saves time in the case that a vertical coordinate was +!! configured but no diagnostics which use the coordinate appeared in the +!! diag_table. +subroutine diag_remap_diag_registration_closed(remap_cs) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure + + if (.not. remap_cs%used) then + call diag_remap_end(remap_cs) + endif + +end subroutine diag_remap_diag_registration_closed + +!> Indicate that this remapping type is actually used by the diag manager. +!! If this is never called then the type will be disabled to save time. +!! See further explanation with diag_remap_registration_closed. +subroutine diag_remap_set_active(remap_cs) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure + + remap_cs%used = .true. + +end subroutine diag_remap_set_active + +!> Configure the vertical axes for a diagnostic remapping control structure. +!! Reads a configuration parameters to determine coordinate generation. +subroutine diag_remap_configure_axes(remap_cs, GV, param_file) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remap control structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file structure + ! Local variables + integer :: nzi(4), nzl(4), k + character(len=200) :: inputdir, string, filename, int_varname, layer_varname + character(len=40) :: mod = "MOM_diag_remap" ! This module's name. + character(len=8) :: units, expected_units + character(len=34) :: longname, string2 + + character(len=256) :: err_msg + logical :: ierr + + real, allocatable, dimension(:) :: interfaces, layers + + call initialize_regridding(remap_cs%regrid_cs, GV, GV%max_depth, param_file, mod, & + trim(remap_cs%vertical_coord_name), "DIAG_COORD", trim(remap_cs%diag_coord_name)) + call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.) + + remap_cs%nz = get_regrid_size(remap_cs%regrid_cs) + + 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 + allocate(interfaces(remap_cs%nz+1)) + allocate(layers(remap_cs%nz)) + + interfaces(:) = getCoordinateInterfaces(remap_cs%regrid_cs) + layers(:) = 0.5 * ( interfaces(1:remap_cs%nz) + interfaces(2:remap_cs%nz+1) ) + + remap_cs%interface_axes_id = diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_i', & + interfaces, trim(units), 'z', & + trim(longname)//' at interface', direction=-1) + remap_cs%layer_axes_id = diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_l', & + layers, trim(units), 'z', & + trim(longname)//' at cell center', direction=-1, & + edges=remap_cs%interface_axes_id) + + ! Axes have now been configured. + remap_cs%configured = .true. + + deallocate(interfaces) + deallocate(layers) + +end subroutine diag_remap_configure_axes + +!> Get layer and interface axes ids for this coordinate +!! Needed when defining axes groups. +subroutine diag_remap_get_axes_info(remap_cs, nz, id_layer, id_interface) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + integer, intent(out) :: nz !< Number of vertical levels for the coordinate + integer, intent(out) :: id_layer !< 1D-axes id for layer points + integer, intent(out) :: id_interface !< 1D-axes id for interface points + + nz = remap_cs%nz + id_layer = remap_cs%layer_axes_id + id_interface = remap_cs%interface_axes_id + +end subroutine diag_remap_get_axes_info + + +!> Whether or not the axes for this vertical coordinated has been configured. +!! Configuration is complete when diag_remap_configure_axes() has been +!! successfully called. +function diag_remap_axes_configured(remap_cs) + type(diag_remap_ctrl), intent(in) :: remap_cs + logical :: diag_remap_axes_configured + + diag_remap_axes_configured = remap_cs%configured + +end function + +!> Build/update target vertical grids for diagnostic remapping. +!! \note The target grids need to be updated whenever sea surface +!! height or layer thicknesses changes. In the case of density-based +!! coordinates then technically we should also regenerate the +!! target grid whenever T/S change. +subroutine diag_remap_update(remap_cs, G, h, T, S, eqn_of_state) + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), pointer :: G !< The ocean's grid type + real, dimension(:, :, :), intent(in) :: h, T, S !< New thickness, T and S + type(EOS_type), pointer, intent(in) :: eqn_of_state !< A pointer to the equation of state + + ! Local variables + integer :: i, j, k, nz + real, dimension(remap_cs%nz + 1) :: zInterfaces + real, dimension(remap_cs%nz) :: resolution + + ! Note that coordinateMode('LAYER') is never 'configured' so will + ! always return here. + if (.not. remap_cs%configured) then + return + endif + + 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.) + allocate(remap_cs%h(G%isd:G%ied,G%jsd:G%jed, nz)) + remap_cs%initialized = .true. + 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) + elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then +! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & +! G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) + call MOM_error(FATAL,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!") + 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) + call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") + endif + remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) + enddo + enddo + +end subroutine diag_remap_update + +!> Remap diagnostic field to alternative vertical grid. +subroutine diag_remap_do_remap(remap_cs, G, h, staggered_in_x, staggered_in_y, & + mask, missing_value, field, remapped_field) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field + real, intent(in) :: missing_value !< A missing_value to assign land/vanished points + real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate + ! Local variables + real, dimension(remap_cs%nz) :: h_dest + real, dimension(size(h,3)) :: h_src + logical :: mask_vanished_layers + integer :: nz_src, nz_dest + integer :: i, j, k + + 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.') + + nz_src = size(field,3) + nz_dest = remap_cs%nz + mask_vanished_layers = (remap_cs%vertical_coord == coordinateMode('ZSTAR')) + remapped_field(:,:,:) = missing_value + + if (staggered_in_x .and. .not. staggered_in_y) then + ! U-points + do j=G%jsc, G%jec + do I=G%iscB, G%iecB + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 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(:), remapped_field(I,j,:), & + remap_cs%remap_cs) + if (mask_vanished_layers) then ! This only works for z-like output + do k=1, nz_dest + if (h_dest(k) == 0.) remapped_field(i, j, k:nz_dest) = missing_value + enddo + endif + enddo + enddo + elseif (staggered_in_y .and. .not. staggered_in_x) then + ! V-points + do J=G%jscB, G%jecB + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 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(:), remapped_field(i,J,:), & + remap_cs%remap_cs) + if (mask_vanished_layers) then ! This only works for z-like output + do k=1, nz_dest + if (h_dest(k) == 0.) remapped_field(i,j,k) = missing_value + enddo + endif + enddo + enddo + elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then + ! H-points + do j=G%jsc, G%jec + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j, 1) == 0.) cycle + endif + h_dest(:) = remap_cs%h(i,j,:) + call remapping_core_h(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, h_dest(:), remapped_field(i,j,:), & + remap_cs%remap_cs) + if (mask_vanished_layers) then ! This only works for z-like output + do k=1, nz_dest + if (h_dest(k)==0.) remapped_field(i,j,k) = missing_value + enddo + endif + enddo + enddo + else + call assert(.false., 'diag_remap_do_remap: Unsupported axis combination') + endif + +end subroutine diag_remap_do_remap + +!> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. +subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & + mask, missing_value, field, reintegrated_field) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field + real, intent(in) :: missing_value !< A missing_value to assign land/vanished points + 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(remap_cs%nz) :: h_dest + real, dimension(size(h,3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(remap_cs%initialized, 'vertically_reintegrate_diag_field: remap_cs not initialized.') + call assert(size(field, 3) == size(h, 3), & + 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.') + + nz_src = size(field,3) + nz_dest = remap_cs%nz + reintegrated_field(:,:,:) = missing_value + + if (staggered_in_x .and. .not. staggered_in_y) then + ! U-points + do j=G%jsc, G%jec + do I=G%iscB, G%iecB + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * ( remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:) ) + call reintegrate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, missing_value, reintegrated_field(I,j,:)) + enddo + enddo + elseif (staggered_in_y .and. .not. staggered_in_x) then + ! V-points + do J=G%jscB, G%jecB + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * ( remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:) ) + call reintegrate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, missing_value, reintegrated_field(i,J,:)) + enddo + enddo + elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then + ! H-points + do j=G%jsc, G%jec + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j, 1) == 0.) cycle + endif + h_src(:) = h(i,j,:) + h_dest(:) = remap_cs%h(i,j,:) + call reintegrate_column(nz_src, h_src, field(i,j,:), & + nz_dest, h_dest, missing_value, reintegrated_field(i,j,:)) + enddo + enddo + else + call assert(.false., 'vertically_reintegrate_diag_field: Q point remapping is not coded yet.') + endif + +end subroutine vertically_reintegrate_diag_field + +!> Vertically interpolate diagnostic field to alternative vertical grid. +subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & + mask, missing_value, field, interpolated_field) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field + real, intent(in) :: missing_value !< A missing_value to assign land/vanished points + 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(remap_cs%nz) :: h_dest + real, dimension(size(h,3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(remap_cs%initialized, 'vertically_interpolate_diag_field: remap_cs not initialized.') + call assert(size(field, 3) == size(h, 3)+1, & + 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.') + + interpolated_field(:,:,:) = missing_value + + nz_src = size(h,3) + nz_dest = remap_cs%nz + + if (staggered_in_x .and. .not. staggered_in_y) then + ! U-points + do j=G%jsc, G%jec + do I=G%iscB, G%iecB + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * ( remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:) ) + call interpolate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, missing_value, interpolated_field(I,j,:)) + enddo + enddo + elseif (staggered_in_y .and. .not. staggered_in_x) then + ! V-points + do J=G%jscB, G%jecB + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j,1) == 0.) cycle + endif + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * ( remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:) ) + call interpolate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, missing_value, interpolated_field(i,J,:)) + enddo + enddo + elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then + ! H-points + do j=G%jsc, G%jec + do i=G%isc, G%iec + if (associated(mask)) then + if (mask(i,j, 1) == 0.) cycle + endif + h_src(:) = h(i,j,:) + h_dest(:) = remap_cs%h(i,j,:) + call interpolate_column(nz_src, h_src, field(i,j,:), & + nz_dest, h_dest, missing_value, interpolated_field(i,j,:)) + enddo + enddo + else + call assert(.false., 'vertically_interpolate_diag_field: Q point remapping is not coded yet.') + endif + +end subroutine vertically_interpolate_diag_field + +!> Horizontally average field +subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggered_in_y, & + is_layer, is_extensive, & + missing_value, field, averaged_field) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points + logical, intent(in) :: is_layer !< True if the z-axis location is at h points + logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers) + real, intent(in) :: missing_value !< A missing_value to assign land/vanished points + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped + real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged + ! Local variables + real, dimension(remap_cs%nz+1) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages + real :: v1, v2, total_volume, total_stuff, val + integer :: i, j, k + + call assert(remap_cs%initialized, 'horizontally_average_diag_field: remap_cs not initialized.') + if (is_layer) then + call assert(size(field, 3) == remap_cs%nz, & + 'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.') + call assert(size(averaged_field, 1) == remap_cs%nz, & + 'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.') + else + call assert(size(field, 3) == remap_cs%nz+1, & + 'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.') + call assert(size(averaged_field, 1) == remap_cs%nz+1, & + 'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.') + endif + + if (staggered_in_x .and. .not. staggered_in_y) then + if (is_layer) then + ! U-points + do k=1,remap_cs%nz + vol_sum(k) = 0. + stuff_sum(k) = 0. + if (is_extensive) then + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCu(I,j) + v2 = G%areaCu(I-1,j) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j) + enddo ; enddo + else ! Intensive + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCu(I,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i+1,j,k) ) + v2 = G%areaCu(I-1,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i-1,j,k) ) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j) + enddo ; enddo + endif + enddo + else ! Interface + do k=1,remap_cs%nz+1 + vol_sum(k) = 0. + stuff_sum(k) = 0. + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCu(I,j) + v2 = G%areaCu(I-1,j) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j) + enddo ; enddo + enddo + endif + elseif (staggered_in_y .and. .not. staggered_in_x) then + if (is_layer) then + ! V-points + do k=1,remap_cs%nz + vol_sum(k) = 0. + stuff_sum(k) = 0. + if (is_extensive) then + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCv(i,J) + v2 = G%areaCv(i,J-1) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j) + enddo ; enddo + else ! Intensive + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCv(i,J) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j+1,k) ) + v2 = G%areaCv(i,J-1) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j-1,k) ) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j) + enddo ; enddo + endif + enddo + else ! Interface + do k=1,remap_cs%nz+1 + vol_sum(k) = 0. + stuff_sum(k) = 0. + do j=G%jsc, G%jec ; do i=G%isc, G%iec + v1 = G%areaCv(i,J) + v2 = G%areaCv(i,J-1) + vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j) + stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j) + enddo ; enddo + enddo + endif + elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then + if (is_layer) then + ! H-points + do k=1,remap_cs%nz + vol_sum(k) = 0. + stuff_sum(k) = 0. + if (is_extensive) then + do j=G%jsc, G%jec ; do i=G%isc, G%iec + if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then + v1 = G%areaT(i,j) + vol_sum(k) = vol_sum(k) + v1 + stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k) + endif + enddo ; enddo + else ! Intensive + do j=G%jsc, G%jec ; do i=G%isc, G%iec + if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then + v1 = G%areaT(i,j) * remap_cs%h(i,j,k) + vol_sum(k) = vol_sum(k) + v1 + stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k) + endif + enddo ; enddo + endif + enddo + else ! Interface + do k=1,remap_cs%nz+1 + vol_sum(k) = 0. + stuff_sum(k) = 0. + do j=G%jsc, G%jec ; do i=G%isc, G%iec + val = field(i,j,k) + if (G%mask2dT(i,j)>0. .and. val/=missing_value) then + v1 = G%areaT(i,j) + vol_sum(k) = vol_sum(k) + v1 + stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k) + endif + enddo ; enddo + enddo + endif + else + call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.') + endif + + call sum_across_PEs(vol_sum, remap_cs%nz) + call sum_across_PEs(stuff_sum, remap_cs%nz) + + do k=1,remap_cs%nz + if (vol_sum(k)>0.) then + averaged_field(k) = stuff_sum(k) / vol_sum(k) + else + averaged_field(k) = missing_value + endif + enddo + +end subroutine horizontally_average_diag_field + +end module MOM_diag_remap 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 diff --git a/src/framework/MOM_error_handler.F90 b/src/framework/MOM_error_handler.F90 index 9d24e35036..3194ced107 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 +!> 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 !< 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) + endif + +end subroutine assert + end module MOM_error_handler diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index 64c14ad213..63b5063d5a 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -37,6 +37,8 @@ module MOM_string_functions public string_functions_unit_tests public extractWord public extract_word +public extract_integer +public extract_real public remove_spaces public slasher @@ -258,6 +260,54 @@ end function extractWord if (b<=ns .and. nw==n-1) extract_word = trim(string(b:ns)) end function extract_word +!> Returns the integer corresponding to the nth word in the argument. +integer function extract_integer(string, separators, n, missing_value) + character(len=*), intent(in) :: string !< String to scan + character(len=*), intent(in) :: separators !< Characters to use for delineation + integer, intent(in) :: n !< Number of word to extract + integer, optional, intent(in) :: missing_value !< Value to assign if word is missing + ! Local variables + integer :: ns, i, b, e, nw + character(len=20) :: word + + word = extract_word(string, separators, n) + + if (len_trim(word)>0) then + read(word(1:len_trim(word)),*) extract_integer + else + if (present(missing_value)) then + extract_integer = missing_value + else + extract_integer = 0 + endif + endif + +end function extract_integer + +!> Returns the real corresponding to the nth word in the argument. +real function extract_real(string, separators, n, missing_value) + character(len=*), intent(in) :: string !< String to scan + character(len=*), intent(in) :: separators !< Characters to use for delineation + integer, intent(in) :: n !< Number of word to extract + real, optional, intent(in) :: missing_value !< Value to assign if word is missing + ! Local variables + integer :: ns, i, b, e, nw + character(len=20) :: word + + word = extract_word(string, separators, n) + + if (len_trim(word)>0) then + read(word(1:len_trim(word)),*) extract_real + else + if (present(missing_value)) then + extract_real = missing_value + else + extract_real = 0 + endif + endif + +end function extract_real + !> Returns string with all spaces removed. character(len=120) function remove_spaces(string) character(len=*), intent(in) :: string !< String to scan @@ -280,34 +330,74 @@ end function extract_word remove_spaces = trim(remove_spaces) end function remove_spaces +!> Returns true if a unit test of string_functions fails. logical function string_functions_unit_tests() - ! Should only be called from a single/root thread - ! Returns True is a test fails, otherwise False integer :: i(5) = (/ -1, 1, 3, 3, 0 /) real :: r(8) = (/ 0., 1., -2., 1.3, 3.E-11, 3.E-11, 3.E-11, -5.1E12 /) - string_functions_unit_tests = .false. + logical :: fail, verbose + verbose = .false. + fail = .false. write(*,*) '==== MOM_string_functions: string_functions_unit_tests ===' - call localTest(left_int(i(1)),'-1') - call localTest(left_ints(i(:)),'-1, 1, 3, 3, 0') - call localTest(left_real(r(1)),'0.0') - call localTest(left_reals(r(:)),'0.0, 1.0, -2.0, 1.3, 3*3.0E-11, -5.1E+12') - call localTest(left_reals(r(:),sep=' '),'0.0 1.0 -2.0 1.3 3*3.0E-11 -5.1E+12') - call localTest(left_reals(r(:),sep=','),'0.0,1.0,-2.0,1.3,3*3.0E-11,-5.1E+12') - call localTest(extractWord("One Two,Three",1),"One") - call localTest(extractWord("One Two,Three",2),"Two") - call localTest(extractWord("One Two,Three",3),"Three") - call localTest(extractWord("One Two, Three",3),"Three") - call localTest(extractWord(" One Two,Three",1),"One") - call localTest(extract_word("One,Two,Three",",",3),"Three") - call localTest(extract_word("One,Two,Three",",",4),"") + fail = fail .or. localTestS(left_int(-1),'-1') + fail = fail .or. localTestS(left_ints(i(:)),'-1, 1, 3, 3, 0') + fail = fail .or. localTestS(left_real(0.),'0.0') + fail = fail .or. localTestS(left_reals(r(:)),'0.0, 1.0, -2.0, 1.3, 3*3.0E-11, -5.1E+12') + fail = fail .or. localTestS(left_reals(r(:),sep=' '),'0.0 1.0 -2.0 1.3 3*3.0E-11 -5.1E+12') + fail = fail .or. localTestS(left_reals(r(:),sep=','),'0.0,1.0,-2.0,1.3,3*3.0E-11,-5.1E+12') + fail = fail .or. localTestS(extractWord("One Two,Three",1),"One") + fail = fail .or. localTestS(extractWord("One Two,Three",2),"Two") + fail = fail .or. localTestS(extractWord("One Two,Three",3),"Three") + fail = fail .or. localTestS(extractWord("One Two, Three",3),"Three") + fail = fail .or. localTestS(extractWord(" One Two,Three",1),"One") + fail = fail .or. localTestS(extract_word("One,Two,Three",",",3),"Three") + fail = fail .or. localTestS(extract_word("One,Two,Three",",",4),"") + fail = fail .or. localTestS(remove_spaces("1 2 3"),"123") + fail = fail .or. localTestS(remove_spaces(" 1 2 3"),"123") + fail = fail .or. localTestS(remove_spaces("1 2 3 "),"123") + fail = fail .or. localTestS(remove_spaces("123"),"123") + fail = fail .or. localTestS(remove_spaces(" "),"") + fail = fail .or. localTestS(remove_spaces(""),"") + fail = fail .or. localTestI(extract_integer("1","",1),1) + fail = fail .or. localTestI(extract_integer("1,2,3",",",1),1) + fail = fail .or. localTestI(extract_integer("1,2",",",2),2) + fail = fail .or. localTestI(extract_integer("1,2",",",3),0) + fail = fail .or. localTestI(extract_integer("1,2",",",4,4),4) + fail = fail .or. localTestR(extract_real("1.","",1),1.) + fail = fail .or. localTestR(extract_real("1.,2.,3.",",",1),1.) + fail = fail .or. localTestR(extract_real("1.,2.",",",2),2.) + fail = fail .or. localTestR(extract_real("1.,2.",",",3),0.) + fail = fail .or. localTestR(extract_real("1.,2.",",",4,4.),4.) + if (.not. fail) write(*,*) 'Pass' write(*,*) '==========================================================' + string_functions_unit_tests = fail contains - subroutine localTest(str1,str2) + logical function localTestS(str1,str2) character(len=*) :: str1, str2 - write(*,*) '>'//trim(str1)//'<' - if (trim(str1)/=trim(str2)) write(*,*) 'FAIL:',trim(str1),':',trim(str2) - if (trim(str1)/=trim(str2)) string_functions_unit_tests=.true. - end subroutine localTest + localTestS=.false. + if (trim(str1)/=trim(str2)) localTestS=.true. + if (localTestS .or. verbose) then + write(*,*) '>'//trim(str1)//'<' + if (localTestS) write(*,*) trim(str1),':',trim(str2), '<-- FAIL' + endif + end function localTestS + logical function localTestI(i1,i2) + integer :: i1,i2 + localTestI=.false. + if (i1/=i2) localTestI=.true. + if (localTestI .or. verbose) then + write(*,*) i1,i2 + if (localTestI) write(*,*) i1,'!=',i2, '<-- FAIL' + endif + end function localTestI + logical function localTestR(r1,r2) + real :: r1,r2 + localTestR=.false. + if (r1/=r2) localTestR=.true. + if (localTestR .or. verbose) then + write(*,*) r1,r2 + if (localTestR) write(*,*) r1,'!=',r2, '<-- FAIL' + endif + end function localTestR end function string_functions_unit_tests !> Returns a directory name that is terminated with a "/" or "./" if the diff --git a/src/framework/_Diagnostics.dox b/src/framework/_Diagnostics.dox new file mode 100644 index 0000000000..07321cb516 --- /dev/null +++ b/src/framework/_Diagnostics.dox @@ -0,0 +1,226 @@ +/*! \page Diagnostics Diagnostics + +MOM6 diagnostics are orchestrated via the FMS diag_manager, as for previous versions of MOM. +However, because MOM6 is a general coordinate model, the model native-coordinae output can be less familiar to users of earlier generations of MOM. +To alleviate this problem, MOM6 provides both "native" and "remapped" diagnostics; +the former being diagnostics in the actual model coordinate space, and the latter in user-defined coordinates. + +\section diag_table The "diag_table" + +At run-time, diagnostics are controlled by the input file `diag_table` which is interpreted but the FMS package diag_manager. +The diag_table syntax is documented at http://data1.gfdl.noaa.gov/~nnz/MOM/mom5_pubrel_August2012/src/shared/diag_manager/diag_table.html. + +The diag_table file has three kinds of section: Title, File and Field. The title section is mandatory and always the first. +There can be multiple file and field sections, typically either in pairs or grouped in to all files and all fields, +but always with the file section preceding the corresponding field section. + +\subsection diag_table_title Title section + +The first two lines are mandatory and comprise a line with a title and a line with six integers defining a base date against which time will be referenced. + +``` +"My ocean-only test case" +1900 1 1 0 0 0 +``` + + +\subsection diag_table_files File section + +This section defines an arbitrary number of files that will be created. +Each file is limited to a single rate of either sampling or time-averaging. + +``` +"file_name", output_freq, "output_freq_units", file_format, "time_axis_units", "time_axis_name" +``` + +- `file_name` : The name of the file that contains diagnostics at the given frequency (excluding the ".nc" extension). + +- `output_freq` : The period between records in `file_name`, if positive. + Special values of 0 mean write every time step and -1 write only at the end of the run. + +- `output_freq_units` : The units in which `output_freq` is given. + Valid values are "years", "months", "days", "hours", "minutes" or "seconds". + +- `file_format` : Always set to 1, meaning netcdf. + +- `time_axis_units` : The units to use for the time-axis in the file. + Valid values are "years", "months", "days", "hours", "minutes" or "seconds". + +- `time_axis_name` : The name of the time-axis (usually "Time"). + +Optional entries in the file line allow the generation of multiple files are intervals: + +``` +"file_name", output_freq, "output_freq_units", file_format, "time_axis_units", "time_axis_name"[, new_file_freq, "new_file_freq_units"[, "start_time"[, file_duration, "file_duration_units"]]] +``` + +- file_name : The base name of the file that contains diagnostics at the given frequency (excluding the ".nc" extension). + The strings %%yr, %%mo, %%dy, %%hr %%mi, %%sc are expanded to the current year, month, day, hour, minute and second respectively, with new files created every new_file_freq. + +- `new_file_freq` : The period between generation of new files. + +- `new_file_freq_units` : The units in which `new_file_freq` is given. + +- `start_time`, `file_duration`, `file_duration_units` : Even finer grain control of output files. + + +\subsection diag_table_fields Field section + +An arbitrary number of lines, one per diagnostic field: + +``` +"module_name", "field_name", "output_name", "file_name", "time_sampling", "reduction_method", "regional_section", packing +``` + +- `module_name` : Name of the component model. + For native ocean variables this should be "ocean_model". + See \ref remapped_diagnostics for non-native vertical-grid diagnostics in the ocean model. + +- `field_name` : The name of the variable as registered in the model. + +- `output_name` " The name of the variable as it will appear in the file. + This is usually the same as the `field_name` but can be used to rename a diagnostic. + +- `file_name` : One of the files defined above in the section \ref diag_table_files. + +- `time_sampling` : Always set to "all". + +- `reduction_method` : "none" means sample or snapshot. + "average" or "mean" performs a time-average. + "min" or "max" diagnose the minium or maxium over each time period. + +- `regional_section` : "none" means global output. A string of six space separated numbers, "lat_min, lat_max, lon_min, lon_max, vert_min, vert_max", limits the diagnostic to a region. + +- `packing` : Data representation in the file. 1 means "real*8", 2 means "real*4", 4 mean 16-bit integers, 8 means 1-byte. + + +\subsection diag_table_example Example + + +``` +"OM4 1/4 degree" +1900 1 1 0 0 0 + +# Static file +"ocean_static", -1, "months", 1, "days", "time" # ocean_static is a protected name. Do not change this line. +"ocean_model", "deptho", "deptho", "ocean_static", "all", "none", "none", 2 +"ocean_model", "geolon", "geolon", "ocean_static", "all", "none", "none", 2 +"ocean_model", "geolat", "geolat", "ocean_static", "all", "none", "none", 2 +"ocean_model", "wet", "wet", "ocean_static", "all", "none", "none", 2 + +# High-frequency file +"surf_%4yr_%3dy", 1, "hours", 1, "days", "time", 1, "months" +"ocean_model","SSH","SSH","surf_%4yr_%3dy","all","none","none",2 + +# Daily averages +"ocean_daily", 1, "days", 1, "days", "time" +"ocean_model", "tos", "tos", "ocean_daily", "all", "mean", "none",2 + +# Monthly averages +"ocean_month", 1, "months", 1, "days", "time" +"ocean_model", "thetao", "thetao", "ocean_month", "all", "mean", "none",2 + +# Annual averages +"ocean_annual", 12, "months", 1, "days", "time" +"ocean_model", "thetao", "thetao", "ocean_annual", "all", "mean", "none",2 + +# Vertical section +"ocean_Bering_Strait", 5, "days", 1, "days", "time" +"ocean_model", "thetao","thetao", "ocean_Bering_Strait", "all", "mean", "-171.4 -168.7 66.1 66.1 -1 -1",2 +``` + + +\section native_diagnostics Native diagnostics + +The list of available diagnostics is dependent on the particular configuration of the model. +For this reason the model writes a record of the available diagnostic fields at run-time into a file "available_diags.*". +See, for example, [available_diags.000000](https://github.com/NOAA-GFDL/MOM6-examples/blob/dev/master/ocean_only/global_ALE/z/available_diags.000000) for the global_ALE z-coordinate ocean-only test case. + +Diagnostic fields in the module "ocean_model" refer to the native variables or diagnostics in the native grid. +Since the model can be run in arbitrary coordinates, say in hybrid-coordinate mode, then native-space diagnostics can be potentially confusing. +Native diagnostics are useful when examining exactly what the model is doing, +or if the vertical coordinate of the model is configured to be a natural coordinate such as pure isopycnal or z* geopotential. + + +\section remapped_diagnostics Vertically remapped diagnostics + +Alternative vertical coordinates can be configured for diagnostic purposes only. + +The run-time parameter `NUM_DIAG_COORDS` controls how many diagnostic coordinates to use. + +The run-time parameter `DIAG_COORDS` defines the mapping between each coordinate, the name of the module in the diag_table and run-time parameter names that define the coordinate. +A list of string tuples, separated by commas, with each tuple in the form of MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME. +`MODULE_SUFFIX` is the string appended to "ocean_model_" to create a module in the diag_table. +`PARAMETER_SUFFIX` is the string appended to "DiAG_COORD_DEF_", and other paramaters, used to control the generation of the named coordinate. +`COORDINATE_NAME` is a name understood by the MOM6 regridding module. Valid examples are "ZSTAR", "SIGMA", "RHO", etc. + +By default, `NUM_DIAG_COORDS=1` and `DIAG_COORDS="z Z ZSTAR", meaning the module "ocean_model_z" provides diagnostics in "z*" coordinates and uses the parameter `DIAG_COORD_DEF_Z`. + +For example, multiple z*-coordinates could be used for diagnostics with +``` +NUM_DIAG_COORDS = 2 +DIAG_COORDS = "z 01 ZSTAR,abc 02 ZSTAR" +DIAG_COORD_DEF_01 = "WOA09" +DIAG_COORD_DEF_02 = "UNIFORM:10,20." +``` +would create the diag_manager modules "ocean_model_z" and "ocean_model_abc". + +The above is equivalent to +``` +NUM_DIAG_COORDS = 2 +DIAG_COORDS = "z ZA ZSTAR,abc ZB ZSTAR" +DIAG_COORD_DEF_ZA = "WOA09" +DIAG_COORD_DEF_ZB = "UNIFORM:10,20." +``` + +To obtain a diagnostic of monthly-averaged potential temperature in both these coordinate systems the diag_table must include the lines + +``` +"ocean_month_z", 1, "months", 1, "days", "time" +"ocean_month_abc", 1, "months", 1, "days", "time" +"ocean_model", "temp", "temp", "ocean_month_z", "all", "mean", "none",2 +"ocean_model", "temp", "temp", "ocean_month_abc", "all", "mean", "none",2 +``` + + +\subsection diag_table_vertical_coords Diagnostic vertical coordinates + +For each of the `NUM_DIAG_COORDS` vertical coordinates listed in `DIAG_COORDS` the corresponding `DIAG_COORD_DEF_%` parameter must be provided. +It can take the following values: +- PARAM : In this case, a corresponding parameter `DIAG_COORD_RES_%` is read that lists the deltas for each level in the coordinate. + For example, DIAG_COODS="z Z ZTAR", DIAG_COORD_DEF_Z="PARAM", DIAG_COORD_RES_Z=5,10,10,15 creates z*-level with 5,10,10,15 meters thicknesses. +- UNIFORM : Uniform distribution down to the maximum depth of the model using the same number of levels as he model. +- UNIFORM:N : Uniform distribution down to the maximum depth of the model using `N` levels. +- UNIFORM:N,D : Uniform distribution of `N` levels with thickness `D`. +- FILE:filename,varname : Reads vector of coordinate thicknesses from field "varname" from file "filename". +- FILE:filename,interfaces=varname : Reads vector of coordinate positions from field "varname" from file "filename". +- WOA09 : Z-levels that correspond to the World Ocean Atlas, 2009, standard levels down and including the maximum depth of the model. +- WOA09:N : The first `N` levels of the Wordl Ocean Atlas, 2009, standard levels. + + + +\section diagnostics_implementation APIs for diagnostics + +The multiple diagnostic-coordinates are implemented in a layer that sits on top of the FMS diag_manager known as the mom6_diag_mediator. + +A diagnostic is registered with register_diag_field() which is an API that looks similar to the FMS diag_manager routine of the same name: +``` +diag_id = register_diag_field(module_name, diag_name, axes, ...) +``` + +The MOM6 version of this routine optionally allows the specification of CMOR names in addition to the native names which then registers the diagnostic twice, once with each name. + +For each of the native and CMOR names, the diagnostic is registered in the native module "ocean_model", corresponding to the native model coordinate, and a module associated with each of the diagnostic coordinates. + +For each diagnostic coordinate, a horizontally-averaged diagnostic is also registered. + +In all, for each 3D diagnostic, the are 2 + 4N diagnostics registered, where N is the number of diagnostic coordinates. +As a result, the global_ALE examples have of order 900 total diagnostics available in the shipped configuration. + + +The data is made available to the diag_manager via a call to post_data() which is a wrapper that does all the vertical remapping before calling FMS's send_data(): +``` +call post_data(diag_id, data, diag_control) +``` + +*/ diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index e11cce0c74..5b133d49fd 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -83,7 +83,7 @@ module MOM_state_initialization use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid -use MOM_regridding, only : regridding_CS, set_regrid_params +use MOM_regridding, only : regridding_CS, set_regrid_params, getCoordinateResolution use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_tracer_initialization_from_Z, only : horiz_interp_and_extrap_tracer @@ -1854,12 +1854,13 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) call pass_var(tmpS1dIn, G%Domain) ! Build the target grid (and set the model thickness to it) - allocate( hTarget(nz) ) ! This call can be more general but is hard-coded for z* coordinates... ???? - call ALE_initRegridding( GV, G%max_depth, PF, mod, regridCS, hTarget ) ! sets regridCS and hTarget(1:nz) + call ALE_initRegridding( GV, G%max_depth, PF, mod, regridCS ) ! sets regridCS if (.not. remap_general) then ! This is the old way of initializing to z* coordinates only + allocate( hTarget(nz) ) + hTarget = getCoordinateResolution( regridCS ) do j = js, je ; do i = is, ie h(i,j,:) = 0. if (G%mask2dT(i,j)>0.) then diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index fc9668d153..7307188037 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -23,7 +23,7 @@ module MOM_tracer_initialization_from_Z use MOM_verticalGrid, only : setVerticalGridAxes use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type use MOM_EOS, only : int_specific_vol_dp -use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord, ALE_remap_scalar +use MOM_ALE, only : ALE_remap_scalar use MOM_regridding, only : regridding_CS use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 1a6e11d8fb..5edf85af79 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -703,7 +703,7 @@ end subroutine horizontal_viscosity subroutine hor_visc_init(Time, G, param_file, diag, CS) type(time_type), intent(in) :: Time - type(ocean_grid_type), intent(in) :: G + type(ocean_grid_type), intent(inout) :: G type(param_file_type), intent(in) :: param_file type(diag_ctrl), target, intent(inout) :: diag type(hor_visc_CS), pointer :: CS 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 aa47399adf..c03070ace6 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -6,7 +6,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 @@ -321,7 +321,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 d7379f7a0e..093cc19b0b 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 0ffc550558..443fdba524 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_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end @@ -419,7 +419,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 @@ -1055,7 +1055,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. & @@ -1982,10 +1982,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! diagnostics making use of the z-gridding code if (associated(diag_to_Z_CSp)) then - vd = var_desc("Kd", "meter2 second-1", & + vd = var_desc("Kd_interface", "meter2 second-1", & "Diapycnal diffusivity at interfaces, interpolated to z", z_grid='z') CS%id_Kd_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Tflx_dia_dif", "degC meter second-1", & + vd = var_desc("Tflx_dia_diff", "degC meter second-1", & "Diffusive diapycnal temperature flux across interfaces, interpolated to z", & z_grid='z') CS%id_Tdif_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) @@ -1993,7 +1993,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "Advective diapycnal temperature flux across interfaces, interpolated to z",& z_grid='z') CS%id_Tadv_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Sflx_dia_dif", "PSU meter second-1", & + vd = var_desc("Sflx_dia_diff", "PSU meter second-1", & "Diffusive diapycnal salinity flux across interfaces, interpolated to z",& z_grid='z') CS%id_Sdif_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time)