MOM_sum_output.F90

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