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