portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | #include "do_concurrent_compat.h" | |
| 6 | ||
| 7 | !> Reports integrated quantities for monitoring the model state | |
| 8 | module MOM_sum_output | |
| 9 | ||
| 10 | use iso_fortran_env, only : int64 | |
| 11 | use MOM_checksums, only : is_NaN, field_checksum | |
| 12 | use MOM_coms, only : sum_across_PEs, PE_here, root_PE, num_PEs, max_across_PEs | |
| 13 | use MOM_coms, only : reproducing_sum, reproducing_sum_EFP, EFP_to_real, real_to_EFP | |
| 14 | use MOM_coms, only : EFP_type, operator(+), operator(-), assignment(=), EFP_sum_across_PEs | |
| 15 | use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe | |
| 16 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 17 | use MOM_forcing_type, only : forcing | |
| 18 | use MOM_grid, only : ocean_grid_type | |
| 19 | use MOM_interface_heights, only : find_eta | |
| 20 | use MOM_io, only : create_MOM_file, reopen_MOM_file | |
| 21 | use MOM_io, only : MOM_infra_file, MOM_netcdf_file, MOM_field | |
| 22 | use MOM_io, only : file_exists, slasher, vardesc, var_desc, MOM_write_field | |
| 23 | use MOM_io, only : field_size, read_variable, read_attribute, open_ASCII_file, stdout | |
| 24 | use MOM_io, only : axis_info, set_axis_info, delete_axis_info, get_filename_appendix | |
| 25 | use MOM_io, only : attribute_info, set_attribute_info, delete_attribute_info | |
| 26 | use MOM_io, only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE | |
| 27 | use MOM_spatial_means, only : array_global_min_max | |
| 28 | use MOM_time_manager, only : time_type, get_time, get_date, set_time | |
| 29 | use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) | |
| 30 | use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<), operator(>) | |
| 31 | use MOM_time_manager, only : get_calendar_type, time_type_to_real, NO_CALENDAR | |
| 32 | use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_stocks | |
| 33 | use MOM_unit_scaling, only : unit_scale_type | |
| 34 | use MOM_variables, only : surface, thermo_var_ptrs | |
| 35 | use MOM_verticalGrid, only : verticalGrid_type | |
| 36 | ||
| 37 | implicit none ; private | |
| 38 | ||
| 39 | #include <MOM_memory.h> | |
| 40 | ||
| 41 | public write_energy, accumulate_net_input | |
| 42 | public MOM_sum_output_init, MOM_sum_output_end | |
| 43 | ||
| 44 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 45 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 46 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 47 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 48 | ||
| 49 | integer, parameter :: NUM_FIELDS = 17 !< Number of diagnostic fields | |
| 50 | character (*), parameter :: depth_chksum_attr = "bathyT_checksum" | |
| 51 | !< Checksum attribute name of G%bathyT | |
| 52 | !! over the compute domain | |
| 53 | character (*), parameter :: area_chksum_attr = "mask2dT_areaT_checksum" | |
| 54 | !< Checksum attribute of name of | |
| 55 | !! G%mask2dT * G%areaT over the compute | |
| 56 | !! domain | |
| 57 | ||
| 58 | !> A list of depths and corresponding globally integrated ocean area at each | |
| 59 | !! depth and the ocean volume below each depth. | |
| 60 | type :: Depth_List | |
| 61 | integer :: listsize !< length of the list <= niglobal*njglobal + 1 | |
| 62 | real, allocatable, dimension(:) :: depth !< A list of depths [Z ~> m] | |
| 63 | real, allocatable, dimension(:) :: area !< The cross-sectional area of the ocean at that depth [L2 ~> m2] | |
| 64 | real, allocatable, dimension(:) :: vol_below !< The ocean volume below that depth [Z L2 ~> m3] | |
| 65 | end type Depth_List | |
| 66 | ||
| 67 | !> The control structure for the MOM_sum_output module | |
| 68 | type, public :: sum_output_CS ; private | |
| 69 | logical :: initialized = .false. !< True if this control structure has been initialized. | |
| 70 | ||
| 71 | type(Depth_List) :: DL !< The sorted depth list. | |
| 72 | ||
| 73 | integer, allocatable, dimension(:) :: lH | |
| 74 | !< This saves the entry in DL with a volume just | |
| 75 | !! less than the volume of fluid below the interface. | |
| 76 | logical :: do_APE_calc !< If true, calculate the available potential energy of the | |
| 77 | !! interfaces. Disabling this reduces the memory footprint of | |
| 78 | !! high-PE-count models dramatically. | |
| 79 | logical :: read_depth_list !< Read the depth list from a file if it exists | |
| 80 | !! and write it if it doesn't. | |
| 81 | character(len=200) :: depth_list_file !< The name of the depth list file. | |
| 82 | real :: D_list_min_inc !< The minimum increment [Z ~> m], between the depths of the | |
| 83 | !! entries in the depth-list file, 0 by default. | |
| 84 | logical :: require_depth_list_chksum | |
| 85 | !< Require matching checksums in Depth_list.nc when reading | |
| 86 | !! the file. | |
| 87 | logical :: update_depth_list_chksum | |
| 88 | !< Automatically update the Depth_list.nc file if the | |
| 89 | !! checksums are missing or do not match current values. | |
| 90 | logical :: use_temperature !< If true, temperature and salinity are state variables. | |
| 91 | type(EFP_type) :: fresh_water_in_EFP !< The total mass of fresh water added by surface fluxes on | |
| 92 | !! this PE since the last time that write_energy was called [kg]. | |
| 93 | type(EFP_type) :: net_salt_in_EFP !< The total salt added by surface fluxes on this PE since | |
| 94 | !! the last time that write_energy was called [ppt kg]. | |
| 95 | type(EFP_type) :: net_heat_in_EFP !< The total heat added by surface fluxes on this PE since | |
| 96 | !! the last time that write_energy was called [J]. | |
| 97 | type(EFP_type) :: heat_prev_EFP !< The total amount of heat in the ocean the last | |
| 98 | !! time that write_energy was called [J]. | |
| 99 | type(EFP_type) :: salt_prev_EFP !< The total amount of salt in the ocean the last | |
| 100 | !! time that write_energy was called [ppt kg]. | |
| 101 | type(EFP_type) :: mass_prev_EFP !< The total ocean mass the last time that | |
| 102 | !! write_energy was called [kg]. | |
| 103 | real :: dt_in_T !< The baroclinic dynamics time step [T ~> s]. | |
| 104 | ||
| 105 | type(time_type) :: energysavedays !< The interval between writing the energies | |
| 106 | !! and other integral quantities of the run. | |
| 107 | type(time_type) :: energysavedays_geometric !< The starting interval for computing a geometric | |
| 108 | !! progression of time deltas between calls to | |
| 109 | !! write_energy. This interval will increase by a factor of 2. | |
| 110 | !! after each call to write_energy. | |
| 111 | logical :: energysave_geometric !< Logical to control whether calls to write_energy should | |
| 112 | !! follow a geometric progression | |
| 113 | type(time_type) :: write_energy_time !< The next time to write to the energy file. | |
| 114 | type(time_type) :: geometric_end_time !< Time at which to stop the geometric progression | |
| 115 | !! of calls to write_energy and revert to the standard | |
| 116 | !! energysavedays interval | |
| 117 | ||
| 118 | real :: timeunit !< The length of the units for the time axis and certain input parameters | |
| 119 | !! including ENERGYSAVEDAYS [s]. | |
| 120 | ||
| 121 | logical :: date_stamped_output !< If true, use dates (not times) in messages to stdout. | |
| 122 | logical :: ISO_date_stamped_output !< If true, use ISO formatted dates in messages to stdout. | |
| 123 | type(time_type) :: Start_time !< The start time of the simulation. | |
| 124 | ! Start_time is set in MOM_initialization.F90 | |
| 125 | integer, pointer :: ntrunc => NULL() !< The number of times the velocity has been | |
| 126 | !! truncated since the last call to write_energy. | |
| 127 | real :: max_Energy !< The maximum permitted energy per unit mass. If there is | |
| 128 | !! more energy than this, the model should stop [L2 T-2 ~> m2 s-2]. | |
| 129 | integer :: maxtrunc !< The number of truncations per energy save | |
| 130 | !! interval at which the run is stopped. | |
| 131 | logical :: write_stocks !< If true, write the integrated tracer amounts | |
| 132 | !! to stdout when the energy files are written. | |
| 133 | logical :: write_min_max !< If true, write the maximum and minimum values of temperature, | |
| 134 | !! salinity and some tracer concentrations to stdout when the energy | |
| 135 | !! files are written. | |
| 136 | logical :: write_min_max_loc !< If true, write the locations of the maximum and minimum values | |
| 137 | !! of temperature, salinity and some tracer concentrations to stdout | |
| 138 | !! when the energy files are written. | |
| 139 | integer :: previous_calls = 0 !< The number of times write_energy has been called. | |
| 140 | integer :: prev_n = 0 !< The value of n from the last call. | |
| 141 | type(MOM_netcdf_file) :: fileenergy_nc !< The file handle for the netCDF version of the energy file. | |
| 142 | integer :: fileenergy_ascii !< The unit number of the ascii version of the energy file. | |
| 143 | type(MOM_field), dimension(NUM_FIELDS+MAX_FIELDS_) :: & | |
| 144 | fields !< fieldtype variables for the output fields. | |
| 145 | character(len=200) :: energyfile !< The name of the energy file with path. | |
| 146 | end type sum_output_CS | |
| 147 | ||
| 148 | contains | |
| 149 | ||
| 150 | !> MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module. | |
| 151 | 1 | subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, & |
| 152 | Input_start_time, CS) | |
| 153 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 154 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 155 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 156 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time | |
| 157 | !! parameters. | |
| 158 | character(len=*), intent(in) :: directory !< The directory where the energy file goes. | |
| 159 | integer, target, intent(inout) :: ntrnc !< The integer that stores the number of times | |
| 160 | !! the velocity has been truncated since the | |
| 161 | !! last call to write_energy. | |
| 162 | type(time_type), intent(in) :: Input_start_time !< The start time of the simulation. | |
| 163 | type(Sum_output_CS), pointer :: CS !< A pointer that is set to point to the | |
| 164 | !! control structure for this module. | |
| 165 | ! Local variables | |
| 166 | real :: maxvel ! The maximum permitted velocity [L T-1 ~> m s-1] | |
| 167 | ! This include declares and sets the variable "version". | |
| 168 | # include "version_variable.h" | |
| 169 | character(len=40) :: mdl = "MOM_sum_output" ! This module's name. | |
| 170 | character(len=200) :: energyfile ! The name of the energy file. | |
| 171 | character(len=32) :: filename_appendix = '' ! FMS appendix to filename for ensemble runs | |
| 172 | ||
| 173 | 1 | if (associated(CS)) then |
| 174 | 0 | call MOM_error(WARNING, "MOM_sum_output_init called with associated control structure.") |
| 175 | 0 | return |
| 176 | endif | |
| 177 | 68 | allocate(CS) |
| 178 | ||
| 179 | 1 | CS%initialized = .true. |
| 180 | ||
| 181 | ! Read all relevant parameters and write them to the model log. | |
| 182 | 1 | call log_version(param_file, mdl, version, "") |
| 183 | call get_param(param_file, mdl, "CALCULATE_APE", CS%do_APE_calc, & | |
| 184 | "If true, calculate the available potential energy of "//& | |
| 185 | "the interfaces. Setting this to false reduces the "//& | |
| 186 | "memory footprint of high-PE-count models dramatically.", & | |
| 187 | 1 | default=.true.) |
| 188 | call get_param(param_file, mdl, "WRITE_STOCKS", CS%write_stocks, & | |
| 189 | "If true, write the integrated tracer amounts to stdout "//& | |
| 190 | 1 | "when the energy files are written.", default=.true.) |
| 191 | call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & | |
| 192 | "If true, Temperature and salinity are used as state "//& | |
| 193 | 1 | "variables.", default=.true.) |
| 194 | call get_param(param_file, mdl, "WRITE_TRACER_MIN_MAX", CS%write_min_max, & | |
| 195 | "If true, write the maximum and minimum values of temperature, salinity and "//& | |
| 196 | "some tracer concentrations to stdout when the energy files are written.", & | |
| 197 | 1 | default=.false., do_not_log=.not.CS%write_stocks, debuggingParam=.true.) |
| 198 | call get_param(param_file, mdl, "WRITE_TRACER_MIN_MAX_LOC", CS%write_min_max_loc, & | |
| 199 | "If true, write the locations of the maximum and minimum values of "//& | |
| 200 | "temperature, salinity and some tracer concentrations to stdout when the "//& | |
| 201 | "energy files are written.", & | |
| 202 | 1 | default=.false., do_not_log=.not.CS%write_min_max, debuggingParam=.true.) |
| 203 | call get_param(param_file, mdl, "DT", CS%dt_in_T, & | |
| 204 | "The (baroclinic) dynamics time step.", & | |
| 205 | 1 | units="s", scale=US%s_to_T, fail_if_missing=.true.) |
| 206 | call get_param(param_file, mdl, "MAXTRUNC", CS%maxtrunc, & | |
| 207 | "The run will be stopped, and the day set to a very "//& | |
| 208 | "large value if the velocity is truncated more than "//& | |
| 209 | "MAXTRUNC times between energy saves. Set MAXTRUNC to 0 "//& | |
| 210 | "to stop if there is any truncation of velocities.", & | |
| 211 | 1 | units="truncations save_interval-1", default=0) |
| 212 | ||
| 213 | call get_param(param_file, mdl, "MAX_ENERGY", CS%max_Energy, & | |
| 214 | "The maximum permitted average energy per unit mass; the "//& | |
| 215 | "model will be stopped if there is more energy than "//& | |
| 216 | "this. If zero or negative, this is set to 10*MAXVEL^2.", & | |
| 217 | 1 | units="m2 s-2", default=0.0, scale=US%m_s_to_L_T**2) |
| 218 | 1 | if (CS%max_Energy <= 0.0) then |
| 219 | call get_param(param_file, mdl, "MAXVEL", maxvel, & | |
| 220 | "The maximum velocity allowed before the velocity "//& | |
| 221 | 1 | "components are truncated.", units="m s-1", default=3.0e8, scale=US%m_s_to_L_T) |
| 222 | 1 | CS%max_Energy = 10.0 * maxvel**2 |
| 223 | call log_param(param_file, mdl, "MAX_ENERGY as used", CS%max_Energy, & | |
| 224 | 1 | units="m2 s-2", unscale=US%L_T_to_m_s**2) |
| 225 | endif | |
| 226 | ||
| 227 | call get_param(param_file, mdl, "ENERGYFILE", energyfile, & | |
| 228 | "The file to use to write the energies and globally "//& | |
| 229 | 1 | "summed diagnostics.", default="ocean.stats") |
| 230 | ||
| 231 | !query fms_io if there is a filename_appendix (for ensemble runs) | |
| 232 | 1 | call get_filename_appendix(filename_appendix) |
| 233 | 1 | if (len_trim(filename_appendix) > 0) then |
| 234 | 0 | energyfile = trim(energyfile) //'.'//trim(filename_appendix) |
| 235 | endif | |
| 236 | ||
| 237 | 1 | CS%energyfile = trim(slasher(directory))//trim(energyfile) |
| 238 | 1 | call log_param(param_file, mdl, "output_path/ENERGYFILE", CS%energyfile) |
| 239 | #ifdef STATSLABEL | |
| 240 | CS%energyfile = trim(CS%energyfile)//"."//trim(adjustl(STATSLABEL)) | |
| 241 | #endif | |
| 242 | ||
| 243 | call get_param(param_file, mdl, "DATE_STAMPED_STDOUT", CS%date_stamped_output, & | |
| 244 | "If true, use dates (not times) in messages to stdout", & | |
| 245 | 1 | default=.true.) |
| 246 | call get_param(param_file, mdl, "ISO_DATE_STAMPED_STDOUT", CS%ISO_date_stamped_output, & | |
| 247 | "If true, use ISO formatted dates in messages to stdout", & | |
| 248 | 1 | default=.false.) |
| 249 | ! Note that the units of CS%Timeunit are the MKS units of [s]. | |
| 250 | call get_param(param_file, mdl, "TIMEUNIT", CS%Timeunit, & | |
| 251 | "The time unit in seconds a number of input fields", & | |
| 252 | 1 | units="s", default=86400.0) |
| 253 | 1 | if (CS%Timeunit < 0.0) CS%Timeunit = 86400.0 |
| 254 | ||
| 255 | 1 | if (CS%do_APE_calc) then |
| 256 | call get_param(param_file, mdl, "READ_DEPTH_LIST", CS%read_depth_list, & | |
| 257 | "Read the depth list from a file if it exists or "//& | |
| 258 | 1 | "create that file otherwise.", default=.false.) |
| 259 | call get_param(param_file, mdl, "DEPTH_LIST_MIN_INC", CS%D_list_min_inc, & | |
| 260 | "The minimum increment between the depths of the "//& | |
| 261 | "entries in the depth-list file.", & | |
| 262 | 1 | units="m", default=1.0E-10, scale=US%m_to_Z) |
| 263 | 1 | if (CS%read_depth_list) then |
| 264 | call get_param(param_file, mdl, "DEPTH_LIST_FILE", CS%depth_list_file, & | |
| 265 | 0 | "The name of the depth list file.", default="Depth_list.nc") |
| 266 | 0 | if (scan(CS%depth_list_file,'/') == 0) & |
| 267 | 0 | CS%depth_list_file = trim(slasher(directory))//trim(CS%depth_list_file) |
| 268 | ||
| 269 | call get_param(param_file, mdl, "REQUIRE_DEPTH_LIST_CHECKSUMS", & | |
| 270 | CS%require_depth_list_chksum, & | |
| 271 | "Require that matching checksums be in Depth_list.nc "//& | |
| 272 | 0 | "when reading the file.", default=.true.) |
| 273 | 0 | if (.not. CS%require_depth_list_chksum) & |
| 274 | call get_param(param_file, mdl, "UPDATE_DEPTH_LIST_CHECKSUMS", & | |
| 275 | CS%update_depth_list_chksum, & | |
| 276 | "Automatically update the Depth_list.nc file if the "//& | |
| 277 | "checksums are missing or do not match current values.", & | |
| 278 | 0 | default=.false.) |
| 279 | endif | |
| 280 | ||
| 281 | 1 | allocate(CS%lH(GV%ke)) |
| 282 | 1 | call depth_list_setup(G, GV, US, CS%DL, CS) |
| 283 | else | |
| 284 | 0 | CS%DL%listsize = 1 |
| 285 | endif | |
| 286 | ||
| 287 | call get_param(param_file, mdl, "ENERGYSAVEDAYS",CS%energysavedays, & | |
| 288 | "The interval in units of TIMEUNIT between saves of the "//& | |
| 289 | "energies of the run and other globally summed diagnostics.",& | |
| 290 | 1 | default=set_time(0,days=1), timeunit=CS%Timeunit) |
| 291 | call get_param(param_file, mdl, "ENERGYSAVEDAYS_GEOMETRIC",CS%energysavedays_geometric, & | |
| 292 | "The starting interval in units of TIMEUNIT for the first call "//& | |
| 293 | "to save the energies of the run and other globally summed diagnostics. "//& | |
| 294 | "The interval increases by a factor of 2. after each call to write_energy.",& | |
| 295 | 1 | default=set_time(seconds=0), timeunit=CS%Timeunit) |
| 296 | ||
| 297 | 1 | if ((time_type_to_real(CS%energysavedays_geometric) > 0.) .and. & |
| 298 | (CS%energysavedays_geometric < CS%energysavedays)) then | |
| 299 | 0 | CS%energysave_geometric = .true. |
| 300 | else | |
| 301 | 1 | CS%energysave_geometric = .false. |
| 302 | endif | |
| 303 | ||
| 304 | 1 | CS%Start_time = Input_start_time |
| 305 | 1 | CS%ntrunc => ntrnc |
| 306 | ||
| 307 | 1 | end subroutine MOM_sum_output_init |
| 308 | ||
| 309 | !> MOM_sum_output_end deallocates memory used by the MOM_sum_output module. | |
| 310 | 1 | subroutine MOM_sum_output_end(CS) |
| 311 | type(Sum_output_CS), pointer :: CS !< The control structure returned by a | |
| 312 | !! previous call to MOM_sum_output_init. | |
| 313 | 1 | if (associated(CS)) then |
| 314 | 1 | if (CS%do_APE_calc) then |
| 315 | 1 | deallocate(CS%DL%depth, CS%DL%area, CS%DL%vol_below) |
| 316 | 1 | deallocate(CS%lH) |
| 317 | endif | |
| 318 | ||
| 319 | 68 | deallocate(CS) |
| 320 | endif | |
| 321 | 1 | end subroutine MOM_sum_output_end |
| 322 | ||
| 323 | !> This subroutine calculates and writes the total model energy, the energy and | |
| 324 | !! mass of each layer, and other globally integrated physical quantities. | |
| 325 | 13 | subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forcing) |
| 326 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 327 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 328 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 329 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & | |
| 330 | intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. | |
| 331 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & | |
| 332 | intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. | |
| 333 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 334 | intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. | |
| 335 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 336 | !! thermodynamic variables. | |
| 337 | type(time_type), intent(in) :: day !< The current model time. | |
| 338 | integer, intent(in) :: n !< The time step number of the | |
| 339 | !! current execution. | |
| 340 | type(Sum_output_CS), pointer :: CS !< The control structure returned by a | |
| 341 | !! previous call to MOM_sum_output_init. | |
| 342 | type(tracer_flow_control_CS), pointer :: tracer_CSp !< Control structure with the tree of | |
| 343 | !! all registered tracer packages | |
| 344 | type(time_type), optional, intent(in) :: dt_forcing !< The forcing time step | |
| 345 | ||
| 346 | ! Local variables | |
| 347 | 26 | real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The height of interfaces [Z ~> m]. |
| 348 | 3 | real :: areaTm(SZI_(G),SZJ_(G)) ! A masked version of areaT [L2 ~> m2]. |
| 349 | 26 | real :: KE(SZK_(GV)) ! The total kinetic energy of a layer [R Z L4 T-2 ~> J] |
| 350 | 26 | real :: PE(SZK_(GV)+1)! The available potential energy of an interface [R Z L4 T-2 ~> J] |
| 351 | real :: KE_tot ! The total kinetic energy [R Z L4 T-2 ~> J]. | |
| 352 | real :: PE_tot ! The total available potential energy [R Z L4 T-2 ~> J]. | |
| 353 | 26 | real :: Z_0APE(SZK_(GV)+1) ! The uniform depth which overlies the same |
| 354 | ! volume as is below an interface [Z ~> m]. | |
| 355 | real :: toten ! The total kinetic & potential energies of all layers [R Z L4 T-2 ~> J] | |
| 356 | real :: En_mass ! The total kinetic and potential energies divided by | |
| 357 | ! the total mass of the ocean [L2 T-2 ~> m2 s-2]. | |
| 358 | 26 | real :: vol_lay(SZK_(GV)) ! The volume of fluid in a layer [Z L2 ~> m3]. |
| 359 | real :: volbelow ! The volume of all layers beneath an interface [Z L2 ~> m3]. | |
| 360 | 26 | real :: mass_lay(SZK_(GV)) ! The mass of fluid in a layer [R Z L2 ~> kg] |
| 361 | real :: mass_tot ! The total mass of the ocean [R Z L2 ~> kg] | |
| 362 | real :: vol_tot ! The total ocean volume [Z L2 ~> m3] | |
| 363 | real :: mass_chg ! The change in total ocean mass of fresh water since | |
| 364 | ! the last call to this subroutine [R Z L2 ~> kg] | |
| 365 | real :: mass_anom ! The change in fresh water that cannot be accounted for | |
| 366 | ! by the surface fluxes [R Z L2 ~> kg] | |
| 367 | real :: Salt ! The total amount of salt in the ocean [1e-3 R Z L2 ~> g Salt] | |
| 368 | real :: Salt_chg ! The change in total ocean salt since the last call | |
| 369 | ! to this subroutine [1e-3 R Z L2 ~> g Salt] | |
| 370 | real :: Salt_anom ! The change in salt that cannot be accounted for by | |
| 371 | ! the surface fluxes [1e-3 R Z L2 ~> g Salt] | |
| 372 | real :: salin ! The mean salinity of the ocean [ppt]. | |
| 373 | real :: salin_anom ! The change in total salt that cannot be accounted for by | |
| 374 | ! the surface fluxes divided by total mass [ppt]. | |
| 375 | real :: Heat ! The total amount of Heat in the ocean [Q R Z L2 ~> J] | |
| 376 | real :: Heat_chg ! The change in total ocean heat since the last call to this subroutine [Q R Z L2 ~> J] | |
| 377 | real :: Heat_anom ! The change in heat that cannot be accounted for by the surface fluxes [Q R Z L2 ~> J] | |
| 378 | real :: temp ! The mean potential temperature of the ocean [C ~> degC] | |
| 379 | real :: temp_anom ! The change in total heat that cannot be accounted for | |
| 380 | ! by the surface fluxes, divided by the total heat | |
| 381 | ! capacity of the ocean [C ~> degC] | |
| 382 | real :: hint ! The deviation of an interface from H [Z ~> m]. | |
| 383 | real :: hbot ! 0 if the basin is deeper than H, or the | |
| 384 | ! height of the basin depth over H otherwise [Z ~> m]. | |
| 385 | ! This makes PE only include real fluid. | |
| 386 | real :: hbelow ! The depth of fluid in all layers beneath an interface [Z ~> m]. | |
| 387 | type(EFP_type) :: & | |
| 388 | mass_EFP, & ! The total mass of the ocean in extended fixed point form [kg]. | |
| 389 | salt_EFP, & ! The total amount of salt in the ocean in extended fixed point form [ppt kg]. | |
| 390 | heat_EFP, & ! The total amount of heat in the ocean in extended fixed point form [J]. | |
| 391 | salt_chg_EFP, & ! The change in total ocean salt since the last call to this subroutine [ppt kg]. | |
| 392 | heat_chg_EFP, & ! The change in total ocean heat since the last call to this subroutine [J]. | |
| 393 | mass_chg_EFP, & ! The change in total ocean mass of fresh water since | |
| 394 | ! the last call to this subroutine [kg]. | |
| 395 | salt_anom_EFP, & ! The change in salt that cannot be accounted for by the surface | |
| 396 | ! fluxes [ppt kg]. | |
| 397 | heat_anom_EFP, & ! The change in heat that cannot be accounted for by the surface fluxes [J]. | |
| 398 | mass_anom_EFP ! The change in fresh water that cannot be accounted for by the surface | |
| 399 | ! fluxes [kg]. | |
| 400 | type(EFP_type), dimension(5) :: EFP_list ! An array of EFP types for joint global sums. | |
| 401 | real :: CFL_Iarea ! Direction-based inverse area used in CFL test [L-2 ~> m-2]. | |
| 402 | real :: CFL_trans ! A transport-based definition of the CFL number [nondim]. | |
| 403 | real :: CFL_lin ! A simpler definition of the CFL number [nondim]. | |
| 404 | real :: max_CFL(2) ! The maxima of the CFL numbers [nondim]. | |
| 405 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 406 | 26 | tmp1 ! A temporary array used in reproducing sums [various] |
| 407 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & | |
| 408 | 26 | PE_pt ! The potential energy at each point [R Z L4 T-2 ~> J]. |
| 409 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 410 | 13 | Temp_int, Salt_int ! Layer and cell integrated heat and salt [Q R Z L2 ~> J] and [1e-3 R Z L2 ~> g Salt]. |
| 411 | real :: RZL4_T2_to_J ! The combination of unit rescaling factors to convert the spatially integrated | |
| 412 | ! kinetic or potential energies into mks units [T2 kg m2 R-1 Z-1 L-4 s-2 ~> 1] | |
| 413 | real :: QRZL2_to_J ! The combination of unit rescaling factors to convert integrated heat | |
| 414 | ! content into mks units [J Q-1 R-1 Z-1 L-2 ~> 1] | |
| 415 | real :: J_to_QRZL2 ! The combination of unit rescaling factors to rescale integrated heat | |
| 416 | ! content from mks units into the internal units of MOM6 [Q R Z L2 J-1 ~> 1] | |
| 417 | real :: kg_to_RZL2 ! The combination of unit rescaling factors to rescale masses from | |
| 418 | ! mks units into the internal units of MOM6 [R Z L2 kg-1 ~> 1] | |
| 419 | real :: salt_to_kg ! A factor used to rescale salt contents [kg R-1 Z-1 L-2 ~> nondim] | |
| 420 | integer :: num_nc_fields ! The number of fields that will actually go into | |
| 421 | ! the NetCDF file. | |
| 422 | integer :: i, j, k, is, ie, js, je, nz, m, Isq, Ieq, Jsq, Jeq, isr, ier, jsr, jer | |
| 423 | integer :: li, lbelow, labove ! indices of deep_area_vol, used to find Z_0APE. | |
| 424 | ! lbelow & labove are lower & upper limits for li | |
| 425 | ! in the search for the entry in lH to use. | |
| 426 | integer :: start_of_day, num_days | |
| 427 | real :: reday ! Time in units given by CS%Timeunit, but often [days] | |
| 428 | character(len=240) :: energypath_nc | |
| 429 | character(len=200) :: mesg | |
| 430 | character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str, ISO_date_str | |
| 431 | logical :: date_stamped, ISO_date_stamped | |
| 432 | type(time_type) :: dt_force ! A time_type version of the forcing timestep. | |
| 433 | ||
| 434 | real :: S_min ! The global minimum unmasked value of the salinity [ppt] | |
| 435 | real :: S_max ! The global maximum unmasked value of the salinity [ppt] | |
| 436 | real :: S_min_x ! The x-positions of the global salinity minima | |
| 437 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 438 | real :: S_min_y ! The y-positions of the global salinity minima | |
| 439 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 440 | real :: S_min_z ! The z-positions of the global salinity minima [layer] | |
| 441 | real :: S_max_x ! The x-positions of the global salinity maxima | |
| 442 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 443 | real :: S_max_y ! The y-positions of the global salinity maxima | |
| 444 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 445 | real :: S_max_z ! The z-positions of the global salinity maxima [layer] | |
| 446 | ||
| 447 | real :: T_min ! The global minimum unmasked value of the temperature [degC] | |
| 448 | real :: T_max ! The global maximum unmasked value of the temperature [degC] | |
| 449 | real :: T_min_x ! The x-positions of the global temperature minima | |
| 450 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 451 | real :: T_min_y ! The y-positions of the global temperature minima | |
| 452 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 453 | real :: T_min_z ! The z-positions of the global temperature minima [layer] | |
| 454 | real :: T_max_x ! The x-positions of the global temperature maxima | |
| 455 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 456 | real :: T_max_y ! The y-positions of the global temperature maxima | |
| 457 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 458 | real :: T_max_z ! The z-positions of the global temperature maxima [layer] | |
| 459 | ||
| 460 | ||
| 461 | ! The units of the tracer stock vary between tracers, with [conc] given explicitly by Tr_units. | |
| 462 | real :: Tr_stocks(MAX_FIELDS_) ! The total amounts of each of the registered tracers [kg conc] | |
| 463 | real :: Tr_min(MAX_FIELDS_) ! The global minimum unmasked value of the tracers [conc] | |
| 464 | real :: Tr_max(MAX_FIELDS_) ! The global maximum unmasked value of the tracers [conc] | |
| 465 | real :: Tr_min_x(MAX_FIELDS_) ! The x-positions of the global tracer minima | |
| 466 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 467 | real :: Tr_min_y(MAX_FIELDS_) ! The y-positions of the global tracer minima | |
| 468 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 469 | real :: Tr_min_z(MAX_FIELDS_) ! The z-positions of the global tracer minima [layer] | |
| 470 | real :: Tr_max_x(MAX_FIELDS_) ! The x-positions of the global tracer maxima | |
| 471 | ! in the units of G%geoLonT, often [degrees_E] or [km] | |
| 472 | real :: Tr_max_y(MAX_FIELDS_) ! The y-positions of the global tracer maxima | |
| 473 | ! in the units of G%geoLatT, often [degrees_N] or [km] | |
| 474 | real :: Tr_max_z(MAX_FIELDS_) ! The z-positions of the global tracer maxima [layer] | |
| 475 | logical :: Tr_minmax_avail(MAX_FIELDS_) ! A flag indicating whether the global minimum and | |
| 476 | ! maximum information are available for each of the tracers | |
| 477 | character(len=40), dimension(MAX_FIELDS_) :: & | |
| 478 | Tr_names, & ! The short names for each of the tracers | |
| 479 | Tr_units ! The units for each of the tracers | |
| 480 | integer :: nTr_stocks ! The total number of tracers in all registered tracer packages | |
| 481 | integer :: iyear, imonth, iday, ihour, iminute, isecond, itick ! For call to get_date() | |
| 482 | ||
| 483 | ! A description for output of each of the fields. | |
| 484 | 10543 | type(vardesc) :: vars(NUM_FIELDS+MAX_FIELDS_) |
| 485 | ||
| 486 | ! write_energy_time is the next integral multiple of energysavedays. | |
| 487 | 13 | dt_force = set_time(seconds=2) ; if (present(dt_forcing)) dt_force = dt_forcing |
| 488 | 13 | if (CS%previous_calls == 0) then |
| 489 | 1 | if (CS%energysave_geometric) then |
| 490 | 0 | if (CS%energysavedays_geometric < CS%energysavedays) then |
| 491 | 0 | CS%write_energy_time = day + CS%energysavedays_geometric |
| 492 | CS%geometric_end_time = CS%Start_time + CS%energysavedays * & | |
| 493 | 0 | (1 + (day - CS%Start_time) / CS%energysavedays) |
| 494 | else | |
| 495 | CS%write_energy_time = CS%Start_time + CS%energysavedays * & | |
| 496 | 0 | (1 + (day - CS%Start_time) / CS%energysavedays) |
| 497 | endif | |
| 498 | else | |
| 499 | CS%write_energy_time = CS%Start_time + CS%energysavedays * & | |
| 500 | 1 | (1 + (day - CS%Start_time) / CS%energysavedays) |
| 501 | endif | |
| 502 | 12 | elseif (day + (dt_force/2) < CS%write_energy_time) then |
| 503 | 10 | return ! Do not write this step |
| 504 | else ! Determine the next write time before proceeding | |
| 505 | 2 | if (CS%energysave_geometric) then |
| 506 | 0 | if (CS%write_energy_time + CS%energysavedays_geometric >= & |
| 507 | CS%geometric_end_time) then | |
| 508 | 0 | CS%write_energy_time = CS%geometric_end_time |
| 509 | 0 | CS%energysave_geometric = .false. ! stop geometric progression |
| 510 | else | |
| 511 | 0 | CS%write_energy_time = CS%write_energy_time + CS%energysavedays_geometric |
| 512 | endif | |
| 513 | 0 | CS%energysavedays_geometric = CS%energysavedays_geometric*2 |
| 514 | else | |
| 515 | 2 | CS%write_energy_time = CS%write_energy_time + CS%energysavedays |
| 516 | endif | |
| 517 | endif | |
| 518 | ||
| 519 | 3 | RZL4_T2_to_J = US%RZL2_to_kg*US%L_T_to_m_s**2 ! Used to unscale energies |
| 520 | 3 | QRZL2_to_J = US%RZL2_to_kg*US%Q_to_J_kg ! Used to unscale heat contents |
| 521 | 3 | salt_to_kg = 0.001*US%RZL2_to_kg ! Used to unscale salt contents |
| 522 | 3 | kg_to_RZL2 = US%kg_m3_to_R*US%m_to_Z*US%m_to_L**2 ! Used to scale masses |
| 523 | 3 | J_to_QRZL2 = US%J_kg_to_Q*kg_to_RZL2 ! Used to scale heat contents |
| 524 | ||
| 525 | 3 | num_nc_fields = 17 |
| 526 | 3 | if (.not.CS%use_temperature) num_nc_fields = 11 |
| 527 | 18 | vars(1) = var_desc("Ntrunc","Nondim","Number of Velocity Truncations",'1','1') |
| 528 | 18 | vars(2) = var_desc("En","Joules","Total Energy",'1','1', conversion=RZL4_T2_to_J) |
| 529 | 18 | vars(3) = var_desc("APE","Joules","Total Interface APE",'1','i', conversion=RZL4_T2_to_J) |
| 530 | 18 | vars(4) = var_desc("KE","Joules","Total Layer KE",'1','L', conversion=RZL4_T2_to_J) |
| 531 | 18 | vars(5) = var_desc("H0","meter","Zero APE Depth of Interface",'1','i', conversion=US%Z_to_m) |
| 532 | 18 | vars(6) = var_desc("Mass_lay","kg","Total Layer Mass",'1','L', conversion=US%RZL2_to_kg) |
| 533 | 18 | vars(7) = var_desc("Mass","kg","Total Mass",'1','1', conversion=US%RZL2_to_kg) |
| 534 | 18 | vars(8) = var_desc("Mass_chg","kg","Total Mass Change between Entries",'1','1', conversion=US%RZL2_to_kg) |
| 535 | 18 | vars(9) = var_desc("Mass_anom","kg","Anomalous Total Mass Change",'1','1', conversion=US%RZL2_to_kg) |
| 536 | 18 | vars(10) = var_desc("max_CFL_trans","Nondim","Maximum finite-volume CFL",'1','1') |
| 537 | 18 | vars(11) = var_desc("max_CFL_lin","Nondim","Maximum finite-difference CFL",'1','1') |
| 538 | 3 | if (CS%use_temperature) then |
| 539 | 18 | vars(12) = var_desc("Salt","kg","Total Salt",'1','1', conversion=salt_to_kg) |
| 540 | 18 | vars(13) = var_desc("Salt_chg","kg","Total Salt Change between Entries",'1','1', conversion=salt_to_kg) |
| 541 | 18 | vars(14) = var_desc("Salt_anom","kg","Anomalous Total Salt Change",'1','1', conversion=salt_to_kg) |
| 542 | 18 | vars(15) = var_desc("Heat","Joules","Total Heat",'1','1', conversion=QRZL2_to_J) |
| 543 | 18 | vars(16) = var_desc("Heat_chg","Joules","Total Heat Change between Entries",'1','1', conversion=QRZL2_to_J) |
| 544 | 18 | vars(17) = var_desc("Heat_anom","Joules","Anomalous Total Heat Change",'1','1', conversion=QRZL2_to_J) |
| 545 | endif | |
| 546 | ||
| 547 | 3 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 548 | 3 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 549 | 3 | isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) |
| 550 | ||
| 551 | 3 | if (.not.associated(CS)) call MOM_error(FATAL, & |
| 552 | 0 | "write_energy: Module must be initialized before it is used.") |
| 553 | ||
| 554 | 3 | if (.not.CS%initialized) call MOM_error(FATAL, & |
| 555 | 0 | "write_energy: Module must be initialized before it is used.") |
| 556 | ||
| 557 | !$omp target enter data map(alloc: areaTm, tmp1, PE_pt, Salt_int, Temp_int, eta, Z_0APE) | |
| 558 | ||
| 559 | 3 | do concurrent (j=js:je, i=is:ie) |
| 560 | 21963 | areaTm(i,j) = G%mask2dT(i,j)*G%areaT(i,j) |
| 561 | enddo | |
| 562 | ||
| 563 | 3 | do concurrent (k=1:nz, j=js:je, i=is:ie) |
| 564 | 1641963 | tmp1(i,j,k) = h(i,j,k) * (GV%H_to_RZ*areaTm(i,j)) |
| 565 | enddo | |
| 566 | 3 | mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP, unscale=US%RZL2_to_kg) |
| 567 | ||
| 568 | 3 | if (GV%Boussinesq) then |
| 569 | 228 | do k=1,nz ; vol_lay(k) = (1.0 / GV%Rho0) * mass_lay(k) ; enddo |
| 570 | else | |
| 571 | 0 | if (CS%do_APE_calc) then |
| 572 | 0 | call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) |
| 573 | 0 | do concurrent (k=1:nz, j=js:je, i=is:ie) |
| 574 | 0 | tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) |
| 575 | enddo | |
| 576 | 0 | vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay, unscale=US%Z_to_m*US%L_to_m**2) |
| 577 | endif | |
| 578 | endif ! Boussinesq | |
| 579 | ||
| 580 | ! Tracer stock diagnostics (optional; depends on registered tracer packages). | |
| 581 | ! TODO: Verify Tracer data transfers | |
| 582 | 3 | nTr_stocks = 0 |
| 583 | 3 | Tr_minmax_avail(:) = .false. |
| 584 | 3 | if (CS%write_min_max .and. CS%write_min_max_loc) then |
| 585 | call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & | |
| 586 | stock_units=Tr_units, num_stocks=nTr_stocks,& | |
| 587 | got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max, & | |
| 588 | xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,& | |
| 589 | 0 | xgmax=Tr_max_x, ygmax=Tr_max_y, zgmax=Tr_max_z) |
| 590 | 3 | elseif (CS%write_min_max) then |
| 591 | call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & | |
| 592 | stock_units=Tr_units, num_stocks=nTr_stocks,& | |
| 593 | 0 | got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max) |
| 594 | else | |
| 595 | call call_tracer_stocks(h, Tr_stocks, G, GV, US, tracer_CSp, stock_names=Tr_names, & | |
| 596 | 3 | stock_units=Tr_units, num_stocks=nTr_stocks) |
| 597 | endif | |
| 598 | 3 | if (nTr_stocks > 0) then |
| 599 | 6 | do m=1,nTr_stocks |
| 600 | vars(num_nc_fields+m) = var_desc(Tr_names(m), units=Tr_units(m), & | |
| 601 | 21 | longname=Tr_names(m), hor_grid='1', z_grid='1') |
| 602 | enddo | |
| 603 | 3 | num_nc_fields = num_nc_fields + nTr_stocks |
| 604 | endif | |
| 605 | ||
| 606 | 3 | if (CS%use_temperature .and. CS%write_stocks) then |
| 607 | call array_global_min_max(tv%T, G, nz, T_min, T_max, & | |
| 608 | 3 | T_min_x, T_min_y, T_min_z, T_max_x, T_max_y, T_max_z, unscale=US%C_to_degC) |
| 609 | call array_global_min_max(tv%S, G, nz, S_min, S_max, & | |
| 610 | 3 | S_min_x, S_min_y, S_min_z, S_max_x, S_max_y, S_max_z, unscale=US%S_to_ppt) |
| 611 | endif | |
| 612 | ||
| 613 | 3 | if (CS%previous_calls == 0) then |
| 614 | ||
| 615 | 1 | CS%mass_prev_EFP = mass_EFP |
| 616 | 1 | CS%fresh_water_in_EFP = real_to_EFP(0.0) |
| 617 | 1 | if (CS%use_temperature) then |
| 618 | 1 | CS%net_salt_in_EFP = real_to_EFP(0.0) ; CS%net_heat_in_EFP = real_to_EFP(0.0) |
| 619 | endif | |
| 620 | ||
| 621 | ! Reopen or create a text output file, with an explanatory header line. | |
| 622 | 1 | if (is_root_pe()) then |
| 623 | 1 | if (day > CS%Start_time) then |
| 624 | 0 | call open_ASCII_file(CS%fileenergy_ascii, trim(CS%energyfile), action=APPEND_FILE) |
| 625 | else | |
| 626 | 1 | call open_ASCII_file(CS%fileenergy_ascii, trim(CS%energyfile), action=WRITEONLY_FILE) |
| 627 | 1 | if (abs(CS%timeunit - 86400.0) < 1.0) then |
| 628 | 1 | if (CS%use_temperature) then |
| 629 | write(CS%fileenergy_ascii,'(" Step,",7x,"Day, Truncs, & | |
| 630 | &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, & | |
| 631 | 1 | &Mean Temp, Frac Mass Err, Salin Err, Temp Err")') |
| 632 | write(CS%fileenergy_ascii,'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,& | |
| 633 | 1 | &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",8x,"[degC]")') |
| 634 | else | |
| 635 | write(CS%fileenergy_ascii,'(" Step,",7x,"Day, Truncs, & | |
| 636 | 0 | &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")') |
| 637 | write(CS%fileenergy_ascii,'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,& | |
| 638 | 0 | &"[kg]",11x,"[Nondim]")') |
| 639 | endif | |
| 640 | else | |
| 641 | 0 | if ((CS%timeunit >= 0.99) .and. (CS%timeunit < 1.01)) then |
| 642 | 0 | time_units = " [seconds] " |
| 643 | 0 | elseif ((CS%timeunit >= 3599.0) .and. (CS%timeunit < 3601.0)) then |
| 644 | 0 | time_units = " [hours] " |
| 645 | 0 | elseif ((CS%timeunit >= 86399.0) .and. (CS%timeunit < 86401.0)) then |
| 646 | 0 | time_units = " [days] " |
| 647 | 0 | elseif ((CS%timeunit >= 3.0e7) .and. (CS%timeunit < 3.2e7)) then |
| 648 | 0 | time_units = " [years] " |
| 649 | else | |
| 650 | 0 | write(time_units,'(9x,"[",es8.2," s] ")') CS%timeunit |
| 651 | endif | |
| 652 | ||
| 653 | 0 | if (CS%use_temperature) then |
| 654 | write(CS%fileenergy_ascii,'(" Step,",7x,"Time, Truncs, & | |
| 655 | &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, & | |
| 656 | 0 | &Mean Temp, Frac Mass Err, Salin Err, Temp Err")') |
| 657 | write(CS%fileenergy_ascii,'(A25,10x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,& | |
| 658 | &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",6x,& | |
| 659 | 0 | &"[degC]")') time_units |
| 660 | else | |
| 661 | write(CS%fileenergy_ascii,'(" Step,",7x,"Time, Truncs, & | |
| 662 | 0 | &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")') |
| 663 | write(CS%fileenergy_ascii,'(A25,10x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,& | |
| 664 | 0 | &"[kg]",11x,"[Nondim]")') time_units |
| 665 | endif | |
| 666 | endif | |
| 667 | endif | |
| 668 | ||
| 669 | 1 | energypath_nc = trim(CS%energyfile) // ".nc" |
| 670 | 1 | if (day > CS%Start_time) then |
| 671 | call reopen_MOM_file(CS%fileenergy_nc, trim(energypath_nc), vars, & | |
| 672 | 0 | num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, G=G, GV=GV) |
| 673 | else | |
| 674 | call create_MOM_file(CS%fileenergy_nc, trim(energypath_nc), vars, & | |
| 675 | 1 | num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, G=G, GV=GV) |
| 676 | endif | |
| 677 | endif | |
| 678 | endif | |
| 679 | ||
| 680 | 3 | if (CS%do_APE_calc) then |
| 681 | 3 | lbelow = 1 ; volbelow = 0.0 |
| 682 | 228 | do k=nz,1,-1 |
| 683 | 225 | volbelow = volbelow + vol_lay(k) |
| 684 | 225 | if ((volbelow >= CS%DL%vol_below(CS%lH(k))) .and. & |
| 685 | (volbelow < CS%DL%vol_below(CS%lH(k)+1))) then | |
| 686 | 151 | li = CS%lH(k) |
| 687 | else | |
| 688 | 74 | labove=CS%DL%listsize |
| 689 | 74 | li = (labove + lbelow) / 2 |
| 690 | 710 | do while (li > lbelow) |
| 691 | 636 | if (volbelow < CS%DL%vol_below(li)) then ; labove = li |
| 692 | 229 | else ; lbelow = li ; endif |
| 693 | 636 | li = (labove + lbelow) / 2 |
| 694 | enddo | |
| 695 | 74 | CS%lH(k) = li |
| 696 | endif | |
| 697 | 225 | lbelow = li |
| 698 | 228 | Z_0APE(K) = CS%DL%depth(li) - (volbelow - CS%DL%vol_below(li)) / CS%DL%area(li) |
| 699 | enddo | |
| 700 | 3 | Z_0APE(nz+1) = CS%DL%depth(2) |
| 701 | ||
| 702 | ! Calculate the Available Potential Energy integrated over each interface. With a nonlinear | |
| 703 | ! equation of state or with a bulk mixed layer this calculation is only approximate. | |
| 704 | ! With an ALE model this does not make sense and should be revisited. | |
| 705 | !$omp target update to(Z_0APE) | |
| 706 | 26499 | do concurrent (k=1:nz+1, j=1:size(PE_pt,2), i=1:size(PE_pt,1)) |
| 707 | 2011011 | PE_pt(i,j,k) = 0.0 |
| 708 | enddo | |
| 709 | 3 | if (GV%Boussinesq) then |
| 710 | 3 | do concurrent (j=js:je, i=is:ie) DO_LOCALITY(local(hbelow, hint, hbot, k)) |
| 711 | 21600 | hbelow = 0.0 |
| 712 | 1641963 | do K=nz,1,-1 |
| 713 | 1620000 | hbelow = hbelow + h(i,j,k) * GV%H_to_Z |
| 714 | 1620000 | hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref)) |
| 715 | 1620000 | hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref) |
| 716 | 1620000 | hbot = (hbot + ABS(hbot)) * 0.5 |
| 717 | PE_pt(i,j,K) = (0.5 * areaTm(i,j)) * (GV%Rho0*GV%g_prime(K)) * & | |
| 718 | 1641600 | (hint * hint - hbot * hbot) |
| 719 | enddo | |
| 720 | enddo | |
| 721 | 0 | elseif (GV%semi_Boussinesq) then |
| 722 | 0 | do concurrent (j=js:je, i=is:ie) DO_LOCALITY(local(hint, hbot, k)) |
| 723 | 0 | do K=nz,1,-1 |
| 724 | 0 | hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. |
| 725 | 0 | hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) |
| 726 | PE_pt(i,j,K) = (0.5 * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * & | |
| 727 | 0 | (hint * hint - hbot * hbot) |
| 728 | enddo | |
| 729 | enddo | |
| 730 | else | |
| 731 | 0 | do concurrent (j=js:je, i=is:ie) DO_LOCALITY(local(hint, hbot, k)) |
| 732 | 0 | do K=nz,2,-1 |
| 733 | 0 | hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. |
| 734 | 0 | hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) |
| 735 | PE_pt(i,j,K) = (0.25 * areaTm(i,j) * & | |
| 736 | ((GV%Rlay(k)+GV%Rlay(k-1))*GV%g_prime(K))) * & | |
| 737 | 0 | (hint * hint - hbot * hbot) |
| 738 | enddo | |
| 739 | 0 | hint = Z_0APE(1) + eta(i,j,1) ! eta and H_0 have opposite signs. |
| 740 | 0 | hbot = max(Z_0APE(1) - (G%bathyT(i,j) + G%Z_ref), 0.0) |
| 741 | PE_pt(i,j,1) = (0.5 * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * & | |
| 742 | 0 | (hint * hint - hbot * hbot) |
| 743 | enddo | |
| 744 | endif | |
| 745 | ||
| 746 | 3 | PE_tot = reproducing_sum(PE_pt, isr, ier, jsr, jer, sums=PE, unscale=RZL4_T2_to_J) |
| 747 | else | |
| 748 | 0 | PE_tot = 0.0 |
| 749 | 0 | do k=1,nz+1 ; PE(K) = 0.0 ; Z_0APE(K) = 0.0 ; enddo |
| 750 | endif | |
| 751 | ||
| 752 | ! Calculate the Kinetic Energy integrated over each layer. | |
| 753 | 3 | do concurrent (k=1:nz, j=js:je, i=is:ie) |
| 754 | tmp1(i,j,k) = (0.25 * GV%H_to_RZ*(areaTm(i,j) * h(i,j,k))) * & | |
| 755 | 1641963 | (((u(I-1,j,k)**2) + (u(I,j,k)**2)) + ((v(i,J-1,k)**2) + (v(i,J,k)**2))) |
| 756 | enddo | |
| 757 | ||
| 758 | 3 | KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=KE, unscale=RZL4_T2_to_J) |
| 759 | ||
| 760 | ! Use reproducing sums to do global integrals relate to the heat, salinity and water budgets. | |
| 761 | 3 | if (CS%use_temperature) then |
| 762 | !$omp target enter data map(to:tv%S, tv%T) | |
| 763 | 3 | do concurrent (j=js:je, i=is:ie) DO_LOCALITY(local(k)) |
| 764 | 21600 | Salt_int(i,j) = 0.0 ; Temp_int(i,j) = 0.0 |
| 765 | 1641963 | do k=1,nz |
| 766 | 1620000 | Salt_int(i,j) = Salt_int(i,j) + tv%S(i,j,k) * (h(i,j,k)*(GV%H_to_RZ * areaTm(i,j))) |
| 767 | 1641600 | Temp_int(i,j) = Temp_int(i,j) + (tv%C_p * tv%T(i,j,k)) * (h(i,j,k)*(GV%H_to_RZ * areaTm(i,j))) |
| 768 | enddo | |
| 769 | enddo | |
| 770 | salt_EFP = reproducing_sum_EFP(Salt_int, isr, ier, jsr, jer, only_on_PE=.true., & | |
| 771 | 3 | unscale=US%RZL2_to_kg*US%S_to_ppt) |
| 772 | heat_EFP = reproducing_sum_EFP(Temp_int, isr, ier, jsr, jer, only_on_PE=.true., & | |
| 773 | 3 | unscale=US%RZL2_to_kg*US%Q_to_J_kg) |
| 774 | ||
| 775 | ! Combining the sums avoids multiple blocking all-PE updates. | |
| 776 | 3 | EFP_list(1) = salt_EFP ; EFP_list(2) = heat_EFP ; EFP_list(3) = CS%fresh_water_in_EFP |
| 777 | 3 | EFP_list(4) = CS%net_salt_in_EFP ; EFP_list(5) = CS%net_heat_in_EFP |
| 778 | 3 | call EFP_sum_across_PEs(EFP_list, 5) |
| 779 | ! Return the globally summed values to the original variables. | |
| 780 | 3 | salt_EFP = EFP_list(1) ; heat_EFP = EFP_list(2) ; CS%fresh_water_in_EFP = EFP_list(3) |
| 781 | 3 | CS%net_salt_in_EFP = EFP_list(4) ; CS%net_heat_in_EFP = EFP_list(5) |
| 782 | !$omp target exit data map(release: tv%S, tv%T) | |
| 783 | else | |
| 784 | 0 | call EFP_sum_across_PEs(CS%fresh_water_in_EFP) |
| 785 | endif | |
| 786 | ||
| 787 | !$omp target exit data map(delete: areaTm, tmp1, PE_pt, Salt_int, Temp_int, eta, Z_0APE) | |
| 788 | ||
| 789 | ! Calculate the maximum CFL numbers. | |
| 790 | 3 | max_CFL(1:2) = 0.0 |
| 791 | 1647228 | do k=1,nz ; do j=js,je ; do I=Isq,Ieq |
| 792 | 1633500 | CFL_Iarea = G%IareaT(i,j) |
| 793 | 1633500 | if (u(I,j,k) < 0.0) & |
| 794 | 354931 | CFL_Iarea = G%IareaT(i+1,j) |
| 795 | ||
| 796 | 1633500 | CFL_trans = abs(u(I,j,k) * CS%dt_in_T) * (G%dy_Cu(I,j) * CFL_Iarea) |
| 797 | 1633500 | CFL_lin = abs(u(I,j,k) * CS%dt_in_T) * G%IdxCu(I,j) |
| 798 | 1633500 | max_CFL(1) = max(max_CFL(1), CFL_trans) |
| 799 | 1647000 | max_CFL(2) = max(max_CFL(2), CFL_lin) |
| 800 | enddo ; enddo ; enddo | |
| 801 | 1660953 | do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie |
| 802 | 1647000 | CFL_Iarea = G%IareaT(i,j) |
| 803 | 1647000 | if (v(i,J,k) < 0.0) & |
| 804 | 511070 | CFL_Iarea = G%IareaT(i,j+1) |
| 805 | ||
| 806 | 1647000 | CFL_trans = abs(v(i,J,k) * CS%dt_in_T) * (G%dx_Cv(i,J) * CFL_Iarea) |
| 807 | 1647000 | CFL_lin = abs(v(i,J,k) * CS%dt_in_T) * G%IdyCv(i,J) |
| 808 | 1647000 | max_CFL(1) = max(max_CFL(1), CFL_trans) |
| 809 | 1660725 | max_CFL(2) = max(max_CFL(2), CFL_lin) |
| 810 | enddo ; enddo ; enddo | |
| 811 | ||
| 812 | 3 | call sum_across_PEs(CS%ntrunc) |
| 813 | ||
| 814 | 3 | call max_across_PEs(max_CFL, 2) |
| 815 | ||
| 816 | 3 | Salt = 0.0 ; Heat = 0.0 |
| 817 | 3 | if (CS%use_temperature) then |
| 818 | 3 | Salt = kg_to_RZL2 * EFP_to_real(salt_EFP) |
| 819 | 3 | Heat = J_to_QRZL2 * EFP_to_real(heat_EFP) |
| 820 | 3 | if (CS%previous_calls == 0) then |
| 821 | 1 | CS%salt_prev_EFP = salt_EFP ; CS%heat_prev_EFP = heat_EFP |
| 822 | endif | |
| 823 | 3 | Salt_chg_EFP = Salt_EFP - CS%salt_prev_EFP |
| 824 | 3 | Salt_chg = kg_to_RZL2 * EFP_to_real(Salt_chg_EFP) |
| 825 | 3 | Salt_anom_EFP = Salt_chg_EFP - CS%net_salt_in_EFP |
| 826 | 3 | Salt_anom = kg_to_RZL2 * EFP_to_real(Salt_anom_EFP) |
| 827 | 3 | Heat_chg_EFP = Heat_EFP - CS%heat_prev_EFP |
| 828 | 3 | Heat_chg = J_to_QRZL2 * EFP_to_real(Heat_chg_EFP) |
| 829 | 3 | Heat_anom_EFP = Heat_chg_EFP - CS%net_heat_in_EFP |
| 830 | 3 | Heat_anom = J_to_QRZL2 * EFP_to_real(Heat_anom_EFP) |
| 831 | endif | |
| 832 | ||
| 833 | 3 | mass_chg_EFP = mass_EFP - CS%mass_prev_EFP |
| 834 | 3 | mass_anom_EFP = mass_chg_EFP - CS%fresh_water_in_EFP |
| 835 | 3 | mass_anom = kg_to_RZL2 * EFP_to_real(mass_anom_EFP) |
| 836 | 3 | if (CS%use_temperature .and. .not.GV%Boussinesq) then |
| 837 | ! net_salt_input needs to be converted from ppt kg to [R Z L2 ~> kg] | |
| 838 | 0 | mass_anom = mass_anom - 0.001*kg_to_RZL2*EFP_to_real(CS%net_salt_in_EFP) |
| 839 | endif | |
| 840 | 3 | mass_chg = kg_to_RZL2 * EFP_to_real(mass_chg_EFP) |
| 841 | ||
| 842 | 3 | if (CS%use_temperature) then |
| 843 | 3 | salin = Salt / mass_tot |
| 844 | 3 | salin_anom = Salt_anom / mass_tot |
| 845 | ! salin_chg = Salt_chg / mass_tot | |
| 846 | 3 | temp = Heat / (mass_tot*tv%C_p) |
| 847 | 3 | temp_anom = Heat_anom / (mass_tot*tv%C_p) |
| 848 | endif | |
| 849 | 3 | toten = KE_tot + PE_tot |
| 850 | ||
| 851 | 3 | En_mass = toten / mass_tot |
| 852 | ||
| 853 | 3 | call get_time(day, start_of_day, num_days) |
| 854 | 3 | date_stamped = (CS%date_stamped_output .and. (get_calendar_type() /= NO_CALENDAR)) |
| 855 | 3 | ISO_date_stamped = (CS%ISO_date_stamped_output .and. (get_calendar_type() /= NO_CALENDAR)) |
| 856 | 3 | if (date_stamped .or. ISO_date_stamped) & |
| 857 | 3 | call get_date(day, iyear, imonth, iday, ihour, iminute, isecond, itick) |
| 858 | 3 | if (abs(CS%timeunit - 86400.0) < 1.0) then |
| 859 | 3 | reday = REAL(num_days)+ (REAL(start_of_day)/86400.0) |
| 860 | 3 | mesg_intro = "MOM Day " |
| 861 | else | |
| 862 | reday = REAL(num_days)*(86400.0/CS%timeunit) + & | |
| 863 | 0 | REAL(start_of_day)/abs(CS%timeunit) |
| 864 | 0 | mesg_intro = "MOM Time " |
| 865 | endif | |
| 866 | 3 | if (reday < 1.0e8) then ; write(day_str, '(F12.3)') reday |
| 867 | 0 | elseif (reday < 1.0e11) then ; write(day_str, '(F15.3)') reday |
| 868 | 0 | else ; write(day_str, '(ES15.9)') reday ; endif |
| 869 | ||
| 870 | 3 | if (n < 1000000) then ; write(n_str, '(I6)') n |
| 871 | 0 | else ; write(n_str, '(I0)') n ; endif |
| 872 | ||
| 873 | 3 | date_str = trim(mesg_intro)//trim(day_str) |
| 874 | 3 | if (date_stamped) & |
| 875 | write(date_str,'("MOM Date",i7,2("/",i2.2)," ",i2.2,2(":",i2.2))') & | |
| 876 | 3 | iyear, imonth, iday, ihour, iminute, isecond |
| 877 | 3 | if (ISO_date_stamped) & |
| 878 | write(ISO_date_str,'(i7.4,2(i2.2),"T",i2.2,2(i2.2))') & | |
| 879 | 0 | iyear, imonth, iday, ihour, iminute, isecond |
| 880 | ||
| 881 | 3 | if (is_root_pe()) then ! Only the root PE actually writes anything. |
| 882 | 3 | if (CS%use_temperature) then |
| 883 | write(stdout,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", & | |
| 884 | & ES18.12, ", Salt ", F15.11,", Temp ", F15.11)') & | |
| 885 | 3 | trim(date_str), trim(n_str), US%L_T_to_m_s**2*En_mass, max_CFL(1), US%RZL2_to_kg*mass_tot, & |
| 886 | 6 | salin, US%C_to_degC*temp |
| 887 | else | |
| 888 | write(stdout,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", ES18.12)') & | |
| 889 | 0 | trim(date_str), trim(n_str), US%L_T_to_m_s**2*En_mass, max_CFL(1), US%RZL2_to_kg*mass_tot |
| 890 | endif | |
| 891 | ||
| 892 | 3 | if (CS%use_temperature) then |
| 893 | 3 | if (ISO_date_stamped) then |
| 894 | write(CS%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, ", CFL ", F8.5, ", SL ",& | |
| 895 | &es11.4,", M ",ES11.5,", S",f8.4,", T",f8.4,& | |
| 896 | &", Me ",ES9.2,", Se ",ES9.2,", Te ",ES9.2)') & | |
| 897 | 0 | trim(n_str), trim(ISO_date_str), CS%ntrunc, US%L_T_to_m_s**2*En_mass, max_CFL(1), & |
| 898 | 0 | -US%Z_to_m*Z_0APE(1), US%RZL2_to_kg*mass_tot, salin, US%C_to_degC*temp, mass_anom/mass_tot, & |
| 899 | 0 | salin_anom, US%C_to_degC*temp_anom |
| 900 | else | |
| 901 | write(CS%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, ", CFL ", F8.5, ", SL ",& | |
| 902 | &es11.4,", M ",ES11.5,", S",f8.4,", T",f8.4,& | |
| 903 | &", Me ",ES9.2,", Se ",ES9.2,", Te ",ES9.2)') & | |
| 904 | 3 | trim(n_str), trim(day_str), CS%ntrunc, US%L_T_to_m_s**2*En_mass, max_CFL(1), & |
| 905 | 3 | -US%Z_to_m*Z_0APE(1), US%RZL2_to_kg*mass_tot, salin, US%C_to_degC*temp, mass_anom/mass_tot, & |
| 906 | 6 | salin_anom, US%C_to_degC*temp_anom |
| 907 | endif | |
| 908 | else | |
| 909 | 0 | if (ISO_date_stamped) then |
| 910 | write(CS%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, ", CFL ", F8.5, ", SL ",& | |
| 911 | &ES11.4,", Mass ",ES11.5,", Me ",ES9.2)') & | |
| 912 | 0 | trim(n_str), trim(ISO_date_str), CS%ntrunc, US%L_T_to_m_s**2*En_mass, max_CFL(1), & |
| 913 | 0 | -US%Z_to_m*Z_0APE(1), US%RZL2_to_kg*mass_tot, mass_anom/mass_tot |
| 914 | else | |
| 915 | write(CS%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, ", CFL ", F8.5, ", SL ",& | |
| 916 | &ES11.4,", Mass ",ES11.5,", Me ",ES9.2)') & | |
| 917 | 0 | trim(n_str), trim(day_str), CS%ntrunc, US%L_T_to_m_s**2*En_mass, max_CFL(1), & |
| 918 | 0 | -US%Z_to_m*Z_0APE(1), US%RZL2_to_kg*mass_tot, mass_anom/mass_tot |
| 919 | endif | |
| 920 | endif | |
| 921 | ||
| 922 | 3 | if (CS%ntrunc > 0) then |
| 923 | write(stdout,'(A," Energy/Mass:",ES12.5," Truncations ",I0)') & | |
| 924 | 0 | trim(date_str), US%L_T_to_m_s**2*En_mass, CS%ntrunc |
| 925 | endif | |
| 926 | ||
| 927 | 3 | if (CS%write_stocks) then |
| 928 | 3 | write(stdout,'(" Total Energy: ",Z16.16,ES24.16)') RZL4_T2_to_J*toten, RZL4_T2_to_J*toten |
| 929 | write(stdout,'(" Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & | |
| 930 | 3 | US%RZL2_to_kg*mass_tot, US%RZL2_to_kg*mass_chg, US%RZL2_to_kg*mass_anom, mass_anom/mass_tot |
| 931 | 3 | if (CS%use_temperature) then |
| 932 | 3 | if (Salt == 0.) then |
| 933 | write(stdout,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & | |
| 934 | 0 | Salt*salt_to_kg, Salt_chg*salt_to_kg, Salt_anom*salt_to_kg |
| 935 | else | |
| 936 | write(stdout,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & | |
| 937 | 3 | Salt*salt_to_kg, Salt_chg*salt_to_kg, Salt_anom*salt_to_kg, Salt_anom/Salt |
| 938 | endif | |
| 939 | 3 | if (CS%write_min_max .and. CS%write_min_max_loc) then |
| 940 | write(stdout,'(16X,"Salinity Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 941 | 0 | S_min, S_min_x, S_min_y, S_min_z |
| 942 | write(stdout,'(16X,"Salinity Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 943 | 0 | S_max, S_max_x, S_max_y, S_max_z |
| 944 | 3 | elseif (CS%write_min_max) then |
| 945 | 0 | write(stdout,'(16X,"Salinity Global Min & Max:",ES24.16,1X,ES24.16)') S_min, S_max |
| 946 | endif | |
| 947 | ||
| 948 | 3 | if (Heat == 0.) then |
| 949 | write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & | |
| 950 | 0 | QRZL2_to_J*Heat, QRZL2_to_J*Heat_chg, QRZL2_to_J*Heat_anom |
| 951 | else | |
| 952 | write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & | |
| 953 | 3 | QRZL2_to_J*Heat, QRZL2_to_J*Heat_chg, QRZL2_to_J*Heat_anom, Heat_anom/Heat |
| 954 | endif | |
| 955 | 3 | if (CS%write_min_max .and. CS%write_min_max_loc) then |
| 956 | write(stdout,'(16X,"Temperature Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 957 | 0 | T_min, T_min_x, T_min_y, T_min_z |
| 958 | write(stdout,'(16X,"Temperature Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 959 | 0 | T_max, T_max_x, T_max_y, T_max_z |
| 960 | 3 | elseif (CS%write_min_max) then |
| 961 | 0 | write(stdout,'(16X,"Temperature Global Min & Max:",ES24.16,1X,ES24.16)') T_min, T_max |
| 962 | endif | |
| 963 | endif | |
| 964 | 6 | do m=1,nTr_stocks |
| 965 | ||
| 966 | write(stdout,'(" Total ",a,": ",ES24.16,1X,a)') & | |
| 967 | 3 | trim(Tr_names(m)), Tr_stocks(m), trim(Tr_units(m)) |
| 968 | ||
| 969 | 6 | if (CS%write_min_max .and. CS%write_min_max_loc .and. Tr_minmax_avail(m)) then |
| 970 | write(stdout,'(18X,a," Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 971 | 0 | trim(Tr_names(m)), Tr_min(m), Tr_min_x(m), Tr_min_y(m), Tr_min_z(m) |
| 972 | write(stdout,'(18X,a," Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & | |
| 973 | 0 | trim(Tr_names(m)), Tr_max(m), Tr_max_x(m), Tr_max_y(m), Tr_max_z(m) |
| 974 | 3 | elseif (CS%write_min_max .and. Tr_minmax_avail(m)) then |
| 975 | write(stdout,'(18X,a," Global Min & Max:",ES24.16,1X,ES24.16)') & | |
| 976 | 0 | trim(Tr_names(m)), Tr_min(m), Tr_max(m) |
| 977 | endif | |
| 978 | ||
| 979 | enddo | |
| 980 | endif | |
| 981 | ||
| 982 | 3 | call CS%fileenergy_nc%write_field(CS%fields(1), real(CS%ntrunc), reday) |
| 983 | 3 | call CS%fileenergy_nc%write_field(CS%fields(2), toten, reday) |
| 984 | 3 | call CS%fileenergy_nc%write_field(CS%fields(3), PE, reday) |
| 985 | 3 | call CS%fileenergy_nc%write_field(CS%fields(4), KE, reday) |
| 986 | 3 | call CS%fileenergy_nc%write_field(CS%fields(5), Z_0APE, reday) |
| 987 | 3 | call CS%fileenergy_nc%write_field(CS%fields(6), mass_lay, reday) |
| 988 | ||
| 989 | 3 | call CS%fileenergy_nc%write_field(CS%fields(7), mass_tot, reday) |
| 990 | 3 | call CS%fileenergy_nc%write_field(CS%fields(8), mass_chg, reday) |
| 991 | 3 | call CS%fileenergy_nc%write_field(CS%fields(9), mass_anom, reday) |
| 992 | 3 | call CS%fileenergy_nc%write_field(CS%fields(10), max_CFL(1), reday) |
| 993 | 3 | call CS%fileenergy_nc%write_field(CS%fields(11), max_CFL(2), reday) |
| 994 | 3 | if (CS%use_temperature) then |
| 995 | 3 | call CS%fileenergy_nc%write_field(CS%fields(12), Salt, reday) |
| 996 | 3 | call CS%fileenergy_nc%write_field(CS%fields(13), salt_chg, reday) |
| 997 | 3 | call CS%fileenergy_nc%write_field(CS%fields(14), salt_anom, reday) |
| 998 | 3 | call CS%fileenergy_nc%write_field(CS%fields(15), Heat, reday) |
| 999 | 3 | call CS%fileenergy_nc%write_field(CS%fields(16), heat_chg, reday) |
| 1000 | 3 | call CS%fileenergy_nc%write_field(CS%fields(17), heat_anom, reday) |
| 1001 | 6 | do m=1,nTr_stocks |
| 1002 | 6 | call CS%fileenergy_nc%write_field(CS%fields(17+m), Tr_stocks(m), reday) |
| 1003 | enddo | |
| 1004 | else | |
| 1005 | 0 | do m=1,nTr_stocks |
| 1006 | 0 | call CS%fileenergy_nc%write_field(CS%fields(11+m), Tr_stocks(m), reday) |
| 1007 | enddo | |
| 1008 | endif | |
| 1009 | ||
| 1010 | 3 | call CS%fileenergy_nc%flush() |
| 1011 | endif ! Only the root PE actually writes anything. | |
| 1012 | ||
| 1013 | 3 | if (is_NaN(En_mass)) then |
| 1014 | 0 | call MOM_error(FATAL, "write_energy : NaNs in total model energy forced model termination.") |
| 1015 | 3 | elseif (En_mass > CS%max_Energy) then |
| 1016 | write(mesg,'("Energy per unit mass of ",ES11.4," exceeds ",ES11.4)') & | |
| 1017 | 0 | US%L_T_to_m_s**2*En_mass, US%L_T_to_m_s**2*CS%max_Energy |
| 1018 | 0 | call MOM_error(FATAL, "write_energy : Excessive energy per unit mass forced model termination.") |
| 1019 | endif | |
| 1020 | 3 | if (CS%ntrunc>CS%maxtrunc) then |
| 1021 | 0 | call MOM_error(FATAL, "write_energy : Ocean velocity has been truncated too many times.") |
| 1022 | endif | |
| 1023 | 3 | CS%ntrunc = 0 |
| 1024 | 3 | CS%previous_calls = CS%previous_calls + 1 |
| 1025 | ||
| 1026 | 3 | CS%mass_prev_EFP = mass_EFP ; CS%fresh_water_in_EFP = real_to_EFP(0.0) |
| 1027 | 3 | if (CS%use_temperature) then |
| 1028 | 3 | CS%salt_prev_EFP = Salt_EFP ; CS%net_salt_in_EFP = real_to_EFP(0.0) |
| 1029 | 3 | CS%heat_prev_EFP = Heat_EFP ; CS%net_heat_in_EFP = real_to_EFP(0.0) |
| 1030 | endif | |
| 1031 | ||
| 1032 | end subroutine write_energy | |
| 1033 | ||
| 1034 | !> This subroutine accumulates the net input of volume, salt and heat, through | |
| 1035 | !! the ocean surface for use in diagnosing conservation. | |
| 1036 | 12 | subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) |
| 1037 | type(forcing), intent(in) :: fluxes !< A structure containing pointers to any possible | |
| 1038 | !! forcing fields. Unused fields are unallocated. | |
| 1039 | type(surface), intent(in) :: sfc_state !< A structure containing fields that | |
| 1040 | !! describe the surface state of the ocean. | |
| 1041 | type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various | |
| 1042 | !! thermodynamic variables. | |
| 1043 | real, intent(in) :: dt !< The amount of time over which to average [T ~> s]. | |
| 1044 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1045 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1046 | type(Sum_output_CS), pointer :: CS !< The control structure returned by a previous call | |
| 1047 | !! to MOM_sum_output_init. | |
| 1048 | ! Local variables | |
| 1049 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 1050 | 24 | FW_in, & ! The net fresh water input, integrated over a timestep [R Z L2 ~> kg]. |
| 1051 | 24 | salt_in, & ! The total salt added by surface fluxes, integrated |
| 1052 | ! over a time step [1e-3 R Z L2 ~> g Salt]. | |
| 1053 | 24 | heat_in ! The total heat added by surface fluxes, integrated |
| 1054 | ! over a time step [Q R Z L2 ~> J]. | |
| 1055 | ||
| 1056 | type(EFP_type) :: & | |
| 1057 | FW_in_EFP, & ! The net fresh water input, integrated over a timestep | |
| 1058 | ! and summed over space [kg]. | |
| 1059 | salt_in_EFP, & ! The total salt added by surface fluxes, integrated | |
| 1060 | ! over a time step and summed over space [ppt kg]. | |
| 1061 | heat_in_EFP ! The total heat added by boundary fluxes, integrated | |
| 1062 | ! over a time step and summed over space [J]. | |
| 1063 | ||
| 1064 | integer :: i, j, is, ie, js, je, isr, ier, jsr, jer | |
| 1065 | ||
| 1066 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 1067 | ||
| 1068 | 105276 | FW_in(:,:) = 0.0 |
| 1069 | 12 | if (associated(fluxes%evap)) then |
| 1070 | 12 | if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then |
| 1071 | 87132 | do j=js,je ; do i=is,ie |
| 1072 | FW_in(i,j) = dt*G%areaT(i,j)*(fluxes%evap(i,j) + & | |
| 1073 | (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + (fluxes%lrunoff(i,j) + fluxes%lrunoff_glc(i,j))) + & | |
| 1074 | 87120 | (fluxes%fprec(i,j) + (fluxes%frunoff(i,j) + fluxes%frunoff_glc(i,j))))) |
| 1075 | enddo ; enddo | |
| 1076 | else | |
| 1077 | call MOM_error(WARNING, & | |
| 1078 | 0 | "accumulate_net_input called with associated evap field, but no precip field.") |
| 1079 | endif | |
| 1080 | endif | |
| 1081 | ||
| 1082 | 87132 | if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie |
| 1083 | 87120 | FW_in(i,j) = FW_in(i,j) + dt * G%areaT(i,j) * fluxes%seaice_melt(i,j) |
| 1084 | enddo ; enddo ; endif | |
| 1085 | ||
| 1086 | 210540 | salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0 |
| 1087 | 12 | if (CS%use_temperature) then |
| 1088 | ||
| 1089 | 87132 | if (associated(fluxes%sw)) then ; do j=js,je ; do i=is,ie |
| 1090 | heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * (fluxes%sw(i,j) + & | |
| 1091 | 87120 | (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j)))) |
| 1092 | enddo ; enddo ; endif | |
| 1093 | ||
| 1094 | 87132 | if (associated(fluxes%seaice_melt_heat)) then ; do j=js,je ; do i=is,ie |
| 1095 | heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & | |
| 1096 | 87120 | fluxes%seaice_melt_heat(i,j) |
| 1097 | enddo ; enddo ; endif | |
| 1098 | ||
| 1099 | ! smg: new code | |
| 1100 | ! include heat content from water transport across ocean surface | |
| 1101 | ! if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie | |
| 1102 | ! heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & | |
| 1103 | ! (fluxes%heat_content_lprec(i,j) + (fluxes%heat_content_fprec(i,j) & | |
| 1104 | ! + (fluxes%heat_content_lrunoff(i,j) + (fluxes%heat_content_frunoff(i,j) & | |
| 1105 | ! + (fluxes%heat_content_cond(i,j) + (fluxes%heat_content_vprec(i,j) & | |
| 1106 | ! + fluxes%heat_content_massout(i,j))))))) | |
| 1107 | ! enddo ; enddo ; endif | |
| 1108 | ||
| 1109 | ! smg: old code | |
| 1110 | 12 | if (associated(fluxes%heat_content_evap)) then |
| 1111 | 0 | do j=js,je ; do i=is,ie |
| 1112 | heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & | |
| 1113 | (fluxes%heat_content_evap(i,j) + fluxes%heat_content_lprec(i,j) + & | |
| 1114 | fluxes%heat_content_cond(i,j) + fluxes%heat_content_fprec(i,j) + & | |
| 1115 | 0 | fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j)) |
| 1116 | enddo ; enddo | |
| 1117 | 12 | elseif (associated(tv%TempxPmE)) then |
| 1118 | 87132 | do j=js,je ; do i=is,ie |
| 1119 | 87120 | heat_in(i,j) = heat_in(i,j) + (tv%C_p * G%areaT(i,j)) * tv%TempxPmE(i,j) |
| 1120 | enddo ; enddo | |
| 1121 | 0 | elseif (associated(fluxes%evap)) then |
| 1122 | 0 | do j=js,je ; do i=is,ie |
| 1123 | 0 | heat_in(i,j) = heat_in(i,j) + (tv%C_p * sfc_state%SST(i,j)) * FW_in(i,j) |
| 1124 | enddo ; enddo | |
| 1125 | endif | |
| 1126 | ||
| 1127 | ! The following heat sources may or may not be used. | |
| 1128 | 12 | if (associated(tv%internal_heat)) then |
| 1129 | 87132 | do j=js,je ; do i=is,ie |
| 1130 | 87120 | heat_in(i,j) = heat_in(i,j) + (tv%C_p * G%areaT(i,j)) * tv%internal_heat(i,j) |
| 1131 | enddo ; enddo | |
| 1132 | endif | |
| 1133 | 87132 | if (associated(tv%frazil)) then ; do j=js,je ; do i=is,ie |
| 1134 | 87120 | heat_in(i,j) = heat_in(i,j) + G%areaT(i,j) * tv%frazil(i,j) |
| 1135 | enddo ; enddo ; endif | |
| 1136 | 87132 | if (associated(fluxes%heat_added)) then ; do j=js,je ; do i=is,ie |
| 1137 | 87120 | heat_in(i,j) = heat_in(i,j) + dt*G%areaT(i,j) * fluxes%heat_added(i,j) |
| 1138 | enddo ; enddo ; endif | |
| 1139 | ! if (associated(sfc_state%sw_lost)) then ; do j=js,je ; do i=is,ie | |
| 1140 | ! sfc_state%sw_lost must be in units of [Q R Z ~> J m-2] | |
| 1141 | ! heat_in(i,j) = heat_in(i,j) - G%areaT(i,j) * sfc_state%sw_lost(i,j) | |
| 1142 | ! enddo ; enddo ; endif | |
| 1143 | ||
| 1144 | 12 | if (associated(fluxes%salt_flux)) then ; do j=js,je ; do i=is,ie |
| 1145 | ! integrate salt_flux in [R Z T-1 ~> kgSalt m-2 s-1] to give [ppt kg] | |
| 1146 | 0 | salt_in(i,j) = dt * G%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j)) |
| 1147 | enddo ; enddo ; endif | |
| 1148 | endif | |
| 1149 | ||
| 1150 | 12 | if ((CS%use_temperature) .or. associated(fluxes%lprec) .or. & |
| 1151 | associated(fluxes%evap)) then | |
| 1152 | ! The on-PE sums are stored here, but the sums across PEs are deferred to | |
| 1153 | ! the next call to write_energy to avoid extra barriers. | |
| 1154 | 12 | isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) |
| 1155 | FW_in_EFP = reproducing_sum_EFP(FW_in, isr, ier, jsr, jer, only_on_PE=.true., & | |
| 1156 | 12 | unscale=US%RZL2_to_kg) |
| 1157 | heat_in_EFP = reproducing_sum_EFP(heat_in, isr, ier, jsr, jer, only_on_PE=.true., & | |
| 1158 | 12 | unscale=US%RZL2_to_kg*US%Q_to_J_kg) |
| 1159 | salt_in_EFP = reproducing_sum_EFP(salt_in, isr, ier, jsr, jer, only_on_PE=.true., & | |
| 1160 | 12 | unscale=US%RZL2_to_kg) |
| 1161 | ||
| 1162 | 12 | CS%fresh_water_in_EFP = CS%fresh_water_in_EFP + FW_in_EFP |
| 1163 | 12 | CS%net_salt_in_EFP = CS%net_salt_in_EFP + salt_in_EFP |
| 1164 | 12 | CS%net_heat_in_EFP = CS%net_heat_in_EFP + heat_in_EFP |
| 1165 | endif | |
| 1166 | ||
| 1167 | 12 | end subroutine accumulate_net_input |
| 1168 | ||
| 1169 | !> This subroutine sets up an ordered list of depths, along with the | |
| 1170 | !! cross sectional areas at each depth and the volume of fluid deeper | |
| 1171 | !! than each depth. This might be read from a previously created file | |
| 1172 | !! or it might be created anew. (For now only new creation occurs. | |
| 1173 | 1 | subroutine depth_list_setup(G, GV, US, DL, CS) |
| 1174 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 1175 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 1176 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1177 | type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes to set up | |
| 1178 | type(Sum_output_CS), pointer :: CS !< The control structure returned by a | |
| 1179 | !! previous call to MOM_sum_output_init. | |
| 1180 | ! Local variables | |
| 1181 | logical :: valid_DL_read | |
| 1182 | integer :: k | |
| 1183 | ||
| 1184 | 1 | if (CS%read_depth_list) then |
| 1185 | 0 | if (file_exists(CS%depth_list_file)) then |
| 1186 | 0 | if (CS%update_depth_list_chksum) then |
| 1187 | call read_depth_list(G, US, DL, CS%depth_list_file, & | |
| 1188 | 0 | require_chksum=CS%require_depth_list_chksum, file_matches=valid_DL_read) |
| 1189 | else | |
| 1190 | 0 | call read_depth_list(G, US, DL, CS%depth_list_file, require_chksum=CS%require_depth_list_chksum) |
| 1191 | 0 | valid_DL_read = .true. ! Otherwise there would have been a fatal error. |
| 1192 | endif | |
| 1193 | else | |
| 1194 | 0 | if (is_root_pe()) call MOM_error(NOTE, "depth_list_setup: "// & |
| 1195 | 0 | trim(CS%depth_list_file)//" does not exist. Creating a new file.") |
| 1196 | 0 | valid_DL_read = .false. |
| 1197 | endif | |
| 1198 | ||
| 1199 | 0 | if (.not.valid_DL_read) then |
| 1200 | 0 | call create_depth_list(G, DL, CS%D_list_min_inc) |
| 1201 | 0 | call write_depth_list(G, US, DL, CS%depth_list_file) |
| 1202 | endif | |
| 1203 | else | |
| 1204 | 1 | call create_depth_list(G, DL, CS%D_list_min_inc) |
| 1205 | endif | |
| 1206 | ||
| 1207 | 76 | do k=1,GV%ke |
| 1208 | 76 | CS%lH(k) = DL%listsize-1 |
| 1209 | enddo | |
| 1210 | ||
| 1211 | 1 | end subroutine depth_list_setup |
| 1212 | ||
| 1213 | !> create_depth_list makes an ordered list of depths, along with the cross | |
| 1214 | !! sectional areas at each depth and the volume of fluid deeper than each depth. | |
| 1215 | 1 | subroutine create_depth_list(G, DL, min_depth_inc) |
| 1216 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1217 | type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes to create | |
| 1218 | real, intent(in) :: min_depth_inc !< The minimum increment between depths in the list [Z ~> m] | |
| 1219 | ||
| 1220 | ! Local variables | |
| 1221 | real, dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: & | |
| 1222 | 2 | Dlist, & !< The global list of bottom depths [Z ~> m]. |
| 1223 | 2 | AreaList !< The global list of cell areas [L2 ~> m2]. |
| 1224 | integer, dimension(G%Domain%niglobal*G%Domain%njglobal+1) :: & | |
| 1225 | 2 | indx2 !< The position of an element in the original unsorted list. |
| 1226 | real :: Dnow !< The depth now being considered for sorting [Z ~> m]. | |
| 1227 | real :: Dprev !< The most recent depth that was considered [Z ~> m]. | |
| 1228 | real :: vol !< The running sum of open volume below a depth [Z L2 ~> m3]. | |
| 1229 | real :: area !< The open area at the current depth [L2 ~> m2]. | |
| 1230 | real :: D_list_prev !< The most recent depth added to the list [Z ~> m]. | |
| 1231 | logical :: add_to_list !< This depth should be included as an entry on the list. | |
| 1232 | ||
| 1233 | integer :: ir, indxt | |
| 1234 | integer :: mls, list_size | |
| 1235 | integer :: list_pos, i_global, j_global | |
| 1236 | integer :: i, j, k, kl | |
| 1237 | ||
| 1238 | 1 | mls = G%Domain%niglobal*G%Domain%njglobal |
| 1239 | ||
| 1240 | ! Need to collect the global data from compute domains to a 1D array for sorting. | |
| 1241 | 7202 | Dlist(:) = 0.0 |
| 1242 | 7202 | Arealist(:) = 0.0 |
| 1243 | 7261 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 1244 | ! Set global indices that start the global domain at 1 (Fortran convention). | |
| 1245 | 7200 | j_global = j + G%jdg_offset - (G%jsg-1) |
| 1246 | 7200 | i_global = i + G%idg_offset - (G%isg-1) |
| 1247 | ||
| 1248 | 7200 | list_pos = (j_global-1)*G%Domain%niglobal + i_global |
| 1249 | 7200 | Dlist(list_pos) = G%bathyT(i,j) + G%Z_ref |
| 1250 | 7260 | Arealist(list_pos) = G%mask2dT(i,j) * G%areaT(i,j) |
| 1251 | enddo ; enddo | |
| 1252 | ||
| 1253 | ! These sums reproduce across PEs because the arrays are only nonzero on one PE. | |
| 1254 | 1 | call sum_across_PEs(Dlist, mls+1) |
| 1255 | 1 | call sum_across_PEs(Arealist, mls+1) |
| 1256 | ||
| 1257 | 7202 | do j=1,mls+1 ; indx2(j) = j ; enddo |
| 1258 | 1 | k = mls / 2 + 1 ; ir = mls |
| 1259 | 10798 | do |
| 1260 | 10799 | if (k > 1) then |
| 1261 | 3600 | k = k - 1 |
| 1262 | 3600 | indxt = indx2(k) |
| 1263 | 3600 | Dnow = Dlist(indxt) |
| 1264 | else | |
| 1265 | 7199 | indxt = indx2(ir) |
| 1266 | 7199 | Dnow = Dlist(indxt) |
| 1267 | 7199 | indx2(ir) = indx2(1) |
| 1268 | 7199 | ir = ir - 1 |
| 1269 | 7199 | if (ir == 1) then ; indx2(1) = indxt ; exit ; endif |
| 1270 | endif | |
| 1271 | 10798 | i=k ; j=k*2 |
| 1272 | 73713 | do ; if (j > ir) exit |
| 1273 | 62915 | if (j < ir .AND. Dlist(indx2(j)) < Dlist(indx2(j+1))) j = j + 1 |
| 1274 | 62915 | if (Dnow < Dlist(indx2(j))) then ; indx2(i) = indx2(j) ; i = j ; j = j + i |
| 1275 | 7308 | else ; j = ir+1 ; endif |
| 1276 | enddo | |
| 1277 | 10798 | indx2(i) = indxt |
| 1278 | enddo | |
| 1279 | ||
| 1280 | ! At this point, the lists should perhaps be culled to save memory. | |
| 1281 | ! Culling: (1) identical depths (e.g. land) - take the last one. | |
| 1282 | ! (2) the topmost and bottommost depths are always saved. | |
| 1283 | ! (3) very close depths | |
| 1284 | ! (4) equal volume changes. | |
| 1285 | ||
| 1286 | ! Count the unique elements in the list. | |
| 1287 | 1 | D_list_prev = Dlist(indx2(mls)) |
| 1288 | 1 | list_size = 2 |
| 1289 | 7200 | do k=mls-1,1,-1 |
| 1290 | 7200 | if (Dlist(indx2(k)) < D_list_prev-min_depth_inc) then |
| 1291 | 2507 | list_size = list_size + 1 |
| 1292 | 2507 | D_list_prev = Dlist(indx2(k)) |
| 1293 | endif | |
| 1294 | enddo | |
| 1295 | ||
| 1296 | 1 | DL%listsize = list_size+1 |
| 1297 | 1 | allocate(DL%depth(DL%listsize), DL%area(DL%listsize), DL%vol_below(DL%listsize)) |
| 1298 | ||
| 1299 | 1 | vol = 0.0 ; area = 0.0 |
| 1300 | 1 | Dprev = Dlist(indx2(mls)) |
| 1301 | 1 | D_list_prev = Dprev |
| 1302 | ||
| 1303 | 1 | kl = 0 |
| 1304 | 7201 | do k=mls,1,-1 |
| 1305 | 7200 | i = indx2(k) |
| 1306 | 7200 | vol = vol + area * (Dprev - Dlist(i)) |
| 1307 | 7200 | area = area + AreaList(i) |
| 1308 | ||
| 1309 | 7200 | add_to_list = .false. |
| 1310 | 7200 | if ((kl == 0) .or. (k==1)) then |
| 1311 | 2 | add_to_list = .true. |
| 1312 | 7198 | elseif (Dlist(indx2(k-1)) < D_list_prev-min_depth_inc) then |
| 1313 | 2507 | add_to_list = .true. |
| 1314 | 2507 | D_list_prev = Dlist(indx2(k-1)) |
| 1315 | endif | |
| 1316 | ||
| 1317 | 7200 | if (add_to_list) then |
| 1318 | 2509 | kl = kl+1 |
| 1319 | 2509 | DL%depth(kl) = Dlist(i) |
| 1320 | 2509 | DL%area(kl) = area |
| 1321 | 2509 | DL%vol_below(kl) = vol |
| 1322 | endif | |
| 1323 | 7201 | Dprev = Dlist(i) |
| 1324 | enddo | |
| 1325 | ||
| 1326 | 1 | do while (kl+1 < DL%listsize) |
| 1327 | ! I don't understand why this is needed... RWH | |
| 1328 | 0 | kl = kl+1 |
| 1329 | 0 | DL%vol_below(kl) = DL%vol_below(kl-1) * 1.000001 |
| 1330 | 0 | DL%area(kl) = DL%area(kl-1) |
| 1331 | 0 | DL%depth(kl) = DL%depth(kl-1) |
| 1332 | enddo | |
| 1333 | ||
| 1334 | 1 | DL%vol_below(DL%listsize) = DL%vol_below(DL%listsize-1) * 1000.0 |
| 1335 | 1 | DL%area(DL%listsize) = DL%area(DL%listsize-1) |
| 1336 | 1 | DL%depth(DL%listsize) = DL%depth(DL%listsize-1) |
| 1337 | ||
| 1338 | 1 | end subroutine create_depth_list |
| 1339 | ||
| 1340 | !> This subroutine writes out the depth list to the specified file. | |
| 1341 | 0 | subroutine write_depth_list(G, US, DL, filename) |
| 1342 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 1343 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1344 | type(Depth_List), intent(in) :: DL !< The list of depths, areas and volumes to write | |
| 1345 | character(len=*), intent(in) :: filename !< The path to the depth list file to write. | |
| 1346 | ||
| 1347 | ! Local variables | |
| 1348 | type(vardesc), dimension(:), allocatable :: & | |
| 1349 | 0 | vars ! Types that described the staggering and metadata for the fields |
| 1350 | type(MOM_field), dimension(:), allocatable :: & | |
| 1351 | 0 | fields ! Types with metadata about the variables that will be written |
| 1352 | type(axis_info), dimension(:), allocatable :: & | |
| 1353 | 0 | extra_axes ! Descriptors for extra axes that might be used |
| 1354 | type(attribute_info), dimension(:), allocatable :: & | |
| 1355 | 0 | global_atts ! Global attributes and their values |
| 1356 | 0 | type(MOM_netcdf_file) :: IO_handle ! The I/O handle of the fileset |
| 1357 | character(len=16) :: depth_chksum, area_chksum | |
| 1358 | ||
| 1359 | ! All ranks are required to compute the global checksum | |
| 1360 | 0 | call get_depth_list_checksums(G, US, depth_chksum, area_chksum) |
| 1361 | ||
| 1362 | 0 | if (.not.is_root_pe()) return |
| 1363 | ||
| 1364 | 0 | allocate(vars(3)) |
| 1365 | 0 | allocate(fields(3)) |
| 1366 | 0 | allocate(extra_axes(1)) |
| 1367 | 0 | allocate(global_atts(2)) |
| 1368 | ||
| 1369 | 0 | call set_axis_info(extra_axes(1), "list", ax_size=DL%listsize) |
| 1370 | 0 | vars(1) = var_desc("depth", "m", "Sorted depth", '1', dim_names=(/"list"/), fixed=.true.) |
| 1371 | 0 | vars(2) = var_desc("area", "m2", "Open area at depth", '1', dim_names=(/"list"/), fixed=.true.) |
| 1372 | 0 | vars(3) = var_desc("vol_below", "m3", "Open volume below depth", '1', dim_names=(/"list"/), fixed=.true.) |
| 1373 | 0 | call set_attribute_info(global_atts(1), depth_chksum_attr, depth_chksum) |
| 1374 | 0 | call set_attribute_info(global_atts(2), area_chksum_attr, area_chksum) |
| 1375 | ||
| 1376 | call create_MOM_file(IO_handle, filename, vars, 3, fields, SINGLE_FILE, & | |
| 1377 | 0 | extra_axes=extra_axes, global_atts=global_atts) |
| 1378 | 0 | call MOM_write_field(IO_handle, fields(1), DL%depth, unscale=US%Z_to_m) |
| 1379 | 0 | call MOM_write_field(IO_handle, fields(2), DL%area, unscale=US%L_to_m**2) |
| 1380 | 0 | call MOM_write_field(IO_handle, fields(3), DL%vol_below, unscale=US%Z_to_m*US%L_to_m**2) |
| 1381 | ||
| 1382 | 0 | call delete_axis_info(extra_axes) |
| 1383 | 0 | call delete_attribute_info(global_atts) |
| 1384 | 0 | deallocate(vars, extra_axes, fields, global_atts) |
| 1385 | 0 | call IO_handle%close() |
| 1386 | 0 | end subroutine write_depth_list |
| 1387 | ||
| 1388 | !> This subroutine reads in the depth list from the specified file | |
| 1389 | !! and allocates the memory within and sets up DL. | |
| 1390 | 0 | subroutine read_depth_list(G, US, DL, filename, require_chksum, file_matches) |
| 1391 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 1392 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1393 | type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes | |
| 1394 | character(len=*), intent(in) :: filename !< The path to the depth list file to read. | |
| 1395 | logical, intent(in) :: require_chksum !< If true, missing or mismatched depth | |
| 1396 | !! and area checksums result in a fatal error. | |
| 1397 | logical, optional, intent(out) :: file_matches !< If present, this indicates whether the file | |
| 1398 | !! has been read with matching depth and area checksums | |
| 1399 | ||
| 1400 | ! Local variables | |
| 1401 | character(len=240) :: var_msg | |
| 1402 | integer :: list_size, ndim, sizes(4) | |
| 1403 | 0 | character(len=:), allocatable :: depth_file_chksum, area_file_chksum |
| 1404 | character(len=16) :: depth_grid_chksum, area_grid_chksum | |
| 1405 | logical :: depth_att_found, area_att_found | |
| 1406 | ||
| 1407 | ! Check bathymetric consistency between this configuration and the depth list file. | |
| 1408 | 0 | call read_attribute(filename, depth_chksum_attr, depth_file_chksum, found=depth_att_found) |
| 1409 | 0 | call read_attribute(filename, area_chksum_attr, area_file_chksum, found=area_att_found) |
| 1410 | ||
| 1411 | 0 | if ((.not.depth_att_found) .or. (.not.area_att_found)) then |
| 1412 | 0 | var_msg = trim(filename) // " checksums are missing;" |
| 1413 | 0 | if (require_chksum) then |
| 1414 | 0 | call MOM_error(FATAL, trim(var_msg) // " aborting.") |
| 1415 | 0 | elseif (present(file_matches)) then |
| 1416 | 0 | call MOM_error(WARNING, trim(var_msg) // " updating file.") |
| 1417 | 0 | file_matches = .false. |
| 1418 | 0 | return |
| 1419 | else | |
| 1420 | 0 | call MOM_error(WARNING, trim(var_msg) // " some diagnostics may not be reproducible.") |
| 1421 | endif | |
| 1422 | else | |
| 1423 | 0 | call get_depth_list_checksums(G, US, depth_grid_chksum, area_grid_chksum) |
| 1424 | ||
| 1425 | 0 | if ((trim(depth_grid_chksum) /= trim(depth_file_chksum)) .or. & |
| 1426 | (trim(area_grid_chksum) /= trim(area_file_chksum)) ) then | |
| 1427 | 0 | var_msg = trim(filename) // " checksums do not match;" |
| 1428 | 0 | if (require_chksum) then |
| 1429 | 0 | call MOM_error(FATAL, trim(var_msg) // " aborting.") |
| 1430 | 0 | elseif (present(file_matches)) then |
| 1431 | 0 | call MOM_error(WARNING, trim(var_msg) // " updating file.") |
| 1432 | 0 | file_matches = .false. |
| 1433 | 0 | return |
| 1434 | else | |
| 1435 | 0 | call MOM_error(WARNING, trim(var_msg) // " some diagnostics may not be reproducible.") |
| 1436 | endif | |
| 1437 | endif | |
| 1438 | endif | |
| 1439 | 0 | if (allocated(area_file_chksum)) deallocate(area_file_chksum) |
| 1440 | 0 | if (allocated(depth_file_chksum)) deallocate(depth_file_chksum) |
| 1441 | ||
| 1442 | ! Get the length of the list. | |
| 1443 | 0 | call field_size(filename, "depth", sizes, ndims=ndim) |
| 1444 | 0 | if (ndim /= 1) call MOM_ERROR(FATAL, "MOM_sum_output read_depth_list: depth in "//& |
| 1445 | 0 | trim(filename)//" has too many or too few dimensions.") |
| 1446 | 0 | list_size = sizes(1) |
| 1447 | ||
| 1448 | 0 | DL%listsize = list_size |
| 1449 | 0 | allocate(DL%depth(list_size), DL%area(list_size), DL%vol_below(list_size)) |
| 1450 | ||
| 1451 | 0 | call read_variable(filename, "depth", DL%depth, scale=US%m_to_Z) |
| 1452 | 0 | call read_variable(filename, "area", DL%area, scale=US%m_to_L**2) |
| 1453 | 0 | call read_variable(filename, "vol_below", DL%vol_below, scale=US%m_to_Z*US%m_to_L**2) |
| 1454 | ||
| 1455 | 0 | if (present(file_matches)) file_matches = .true. |
| 1456 | ||
| 1457 | 0 | end subroutine read_depth_list |
| 1458 | ||
| 1459 | ||
| 1460 | !> Return the checksums required to verify DEPTH_LIST_FILE contents. | |
| 1461 | !! | |
| 1462 | !! This function computes checksums for the bathymetry (G%bathyT) and masked | |
| 1463 | !! area (mask2dT * areaT) fields of the model grid G, which are used to compute | |
| 1464 | !! the depth list. A difference in checksum indicates that a different method | |
| 1465 | !! was used to compute the grid data, and that any results using the depth | |
| 1466 | !! list, such as APE, will not be reproducible. | |
| 1467 | !! | |
| 1468 | !! Checksums are saved as hexadecimal strings, in order to avoid potential | |
| 1469 | !! datatype issues with netCDF attributes. | |
| 1470 | 0 | subroutine get_depth_list_checksums(G, US, depth_chksum, area_chksum) |
| 1471 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure | |
| 1472 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1473 | character(len=16), intent(out) :: depth_chksum !< Depth checksum hexstring | |
| 1474 | character(len=16), intent(out) :: area_chksum !< Area checksum hexstring | |
| 1475 | ||
| 1476 | ! Local variables | |
| 1477 | 0 | real, allocatable :: field(:,:) ! A temporary array with no halos [Z ~> m] or [L2 ~> m2] |
| 1478 | integer :: i, j | |
| 1479 | ||
| 1480 | 0 | allocate(field(G%isc:G%iec, G%jsc:G%jec)) |
| 1481 | ||
| 1482 | ! Depth checksum | |
| 1483 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 1484 | 0 | field(i,j) = G%bathyT(i,j) + G%Z_ref |
| 1485 | enddo ; enddo | |
| 1486 | 0 | write(depth_chksum, '(Z16)') field_checksum(field(:,:), unscale=US%Z_to_m) |
| 1487 | ||
| 1488 | ! Area checksum | |
| 1489 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 1490 | 0 | field(i,j) = G%mask2dT(i,j) * G%areaT(i,j) |
| 1491 | enddo ; enddo | |
| 1492 | 0 | write(area_chksum, '(Z16)') field_checksum(field(:,:), unscale=US%L_to_m**2) |
| 1493 | ||
| 1494 | 0 | deallocate(field) |
| 1495 | 0 | end subroutine get_depth_list_checksums |
| 1496 | ||
| 1497 | !> \namespace mom_sum_output | |
| 1498 | !! | |
| 1499 | !! By Robert Hallberg, April 1994 - June 2002 | |
| 1500 | !! | |
| 1501 | !! This file contains the subroutine (write_energy) that writes | |
| 1502 | !! horizontally integrated quantities, such as energies and layer | |
| 1503 | !! volumes, and other summary information to an output file. Some | |
| 1504 | !! of these quantities (APE or resting interface height) are defined | |
| 1505 | !! relative to the global histogram of topography. The subroutine | |
| 1506 | !! that compiles that histogram (depth_list_setup) is also included | |
| 1507 | !! in this file. | |
| 1508 | !! | |
| 1509 | !! In addition, if the number of velocity truncations since the | |
| 1510 | !! previous call to write_energy exceeds maxtrunc or the total energy | |
| 1511 | !! exceeds a very large threshold, a fatal termination is triggered. | |
| 1512 | ||
| 1513 | 0 | end module MOM_sum_output |