MOM6
mom_sum_output Module Reference

Detailed Description

Reports integrated quantities for monitoring the model state.

By Robert Hallberg, April 1994 - June 2002

This file contains the subroutine (write_energy) that writes horizontally integrated quantities, such as energies and layer volumes, and other summary information to an output file. Some of these quantities (APE or resting interface height) are defined relative to the global histogram of topography. The subroutine that compiles that histogram (depth_list_setup) is also included in this file.

In addition, if the number of velocity truncations since the previous call to write_energy exceeds maxtrunc or the total energy exceeds a very large threshold, a fatal termination is triggered.

Data Types

type  depth_list
 A list of depths and corresponding globally integrated ocean area at each depth and the ocean volume below each depth. More...
 
type  sum_output_cs
 The control structure for the MOM_sum_output module. More...
 

Functions/Subroutines

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. More...
 
subroutine mom_sum_output_end (CS)
 MOM_sum_output_end deallocates memory used by the MOM_sum_output module. More...
 
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, and other globally integrated physical quantities. More...
 
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 in diagnosing conservation. More...
 
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 and the volume of fluid deeper than each depth. This might be read from a previously created file or it might be created anew. (For now only new creation occurs. More...
 
subroutine create_depth_list (G, CS)
 create_depth_list makes an ordered list of depths, along with the cross sectional areas at each depth and the volume of fluid deeper than each depth. More...
 
subroutine write_depth_list (G, US, CS, filename, list_size)
 This subroutine writes out the depth list to the specified file. More...
 
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 CSlist_size . More...
 
subroutine get_depth_list_checksums (G, depth_chksum, area_chksum)
 Return the checksums required to verify DEPTH_LIST_FILE contents. More...
 

Variables

integer, parameter num_fields = 17
 Number of diagnostic fields. More...
 
character(*), parameter depth_chksum_attr = "bathyT_checksum"
 Checksum attribute name of GbathyT over the compute domain. More...
 
character(*), parameter area_chksum_attr = "mask2dT_areaT_checksum"
 Checksum attribute of name of Gmask2dT * GareaT over the compute domain. More...
 

Function/Subroutine Documentation

◆ accumulate_net_input()

subroutine, public mom_sum_output::accumulate_net_input ( type(forcing), intent(in)  fluxes,
type(surface), intent(in)  sfc_state,
type(thermo_var_ptrs), intent(in)  tv,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(sum_output_cs), pointer  CS 
)

This subroutine accumates the net input of volume, salt and heat, through the ocean surface for use in diagnosing conservation.

Parameters
[in]fluxesA structure containing pointers to any possible forcing fields. Unused fields are unallocated.
[in]sfc_stateA structure containing fields that describe the surface state of the ocean.
[in]tvA structure pointing to various thermodynamic variables.
[in]dtThe amount of time over which to average [T ~> s].
[in]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to MOM_sum_output_init.

Definition at line 941 of file MOM_sum_output.F90.

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 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ create_depth_list()

subroutine mom_sum_output::create_depth_list ( type(ocean_grid_type), intent(in)  G,
type(sum_output_cs), pointer  CS 
)
private

create_depth_list makes an ordered list of depths, along with the cross sectional areas at each depth and the volume of fluid deeper than each depth.

Parameters
[in]gThe ocean's grid structure.
csThe control structure set up in MOM_sum_output_init, in which the ordered depth list is stored.

Definition at line 1109 of file MOM_sum_output.F90.

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 

Referenced by depth_list_setup(), and read_depth_list().

Here is the caller graph for this function:

◆ depth_list_setup()

subroutine mom_sum_output::depth_list_setup ( type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(sum_output_cs), pointer  CS 
)
private

This subroutine sets up an ordered list of depths, along with the cross sectional areas at each depth and the volume of fluid deeper than each depth. This might be read from a previously created file or it might be created anew. (For now only new creation occurs.

Parameters
[in]gThe ocean's grid structure
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to MOM_sum_output_init.

Definition at line 1079 of file MOM_sum_output.F90.

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 

References create_depth_list(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), read_depth_list(), and write_depth_list().

Referenced by mom_sum_output_init().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_depth_list_checksums()

subroutine mom_sum_output::get_depth_list_checksums ( type(ocean_grid_type), intent(in)  G,
character(len=16), intent(out)  depth_chksum,
character(len=16), intent(out)  area_chksum 
)
private

Return the checksums required to verify DEPTH_LIST_FILE contents.

This function computes checksums for the bathymetry (GbathyT) and masked area (mask2dT * areaT) fields of the model grid G, which are used to compute the depth list. A difference in checksum indicates that a different method was used to compute the grid data, and that any results using the depth list, such as APE, will not be reproducible.

Checksums are saved as hexadecimal strings, in order to avoid potential datatype issues with netCDF attributes.

Parameters
[in]gOcean grid structure
[out]depth_chksumDepth checksum hexstring
[out]area_chksumArea checksum hexstring

Definition at line 1485 of file MOM_sum_output.F90.

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)

Referenced by read_depth_list(), and write_depth_list().

Here is the caller graph for this function:

◆ mom_sum_output_end()

subroutine mom_sum_output::mom_sum_output_end ( type(sum_output_cs), pointer  CS)
private

MOM_sum_output_end deallocates memory used by the MOM_sum_output module.

Parameters
csThe control structure returned by a previous call to MOM_sum_output_init.

Definition at line 290 of file MOM_sum_output.F90.

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

◆ mom_sum_output_init()

subroutine, public mom_sum_output::mom_sum_output_init ( type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
character(len=*), intent(in)  directory,
integer, intent(inout), target  ntrnc,
type(time_type), intent(in)  Input_start_time,
type(sum_output_cs), pointer  CS 
)

MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module.

Parameters
[in]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in]directoryThe directory where the energy file goes.
[in,out]ntrncThe integer that stores the number of times the velocity has been truncated since the last call to write_energy.
[in]input_start_timeThe start time of the simulation.
csA pointer that is set to point to the control structure for this module.

Definition at line 142 of file MOM_sum_output.F90.

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 

References depth_list_setup(), and mom_error_handler::mom_error().

Referenced by mom::initialize_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ read_depth_list()

subroutine mom_sum_output::read_depth_list ( type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(sum_output_cs), pointer  CS,
character(len=*), intent(in)  filename 
)
private

This subroutine reads in the depth list to the specified file and allocates and sets up CSDL and CSlist_size .

Parameters
[in]gThe ocean's grid structure
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to MOM_sum_output_init.
[in]filenameThe path to the depth list file to read.

Definition at line 1329 of file MOM_sum_output.F90.

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) // ":" &
1375  // depth_chksum_attr
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 

References area_chksum_attr, create_depth_list(), depth_chksum_attr, get_depth_list_checksums(), mom_error_handler::mom_error(), and write_depth_list().

Referenced by depth_list_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ write_depth_list()

subroutine mom_sum_output::write_depth_list ( type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(sum_output_cs), pointer  CS,
character(len=*), intent(in)  filename,
integer, intent(in)  list_size 
)
private

This subroutine writes out the depth list to the specified file.

Parameters
[in]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to MOM_sum_output_init.
[in]filenameThe path to the depth list file to write.
[in]list_sizeThe size of the depth list.

Definition at line 1234 of file MOM_sum_output.F90.

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 

References area_chksum_attr, depth_chksum_attr, get_depth_list_checksums(), mom_error_handler::is_root_pe(), and mom_error_handler::mom_error().

Referenced by depth_list_setup(), and read_depth_list().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ write_energy()

subroutine, public mom_sum_output::write_energy ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(time_type), intent(in)  day,
integer, intent(in)  n,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(sum_output_cs), pointer  CS,
type(tracer_flow_control_cs), optional, pointer  tracer_CSp,
type(ocean_obc_type), optional, pointer  OBC,
type(time_type), intent(in), optional  dt_forcing 
)

This subroutine calculates and writes the total model energy, the energy and mass of each layer, and other globally integrated physical quantities.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1].
[in]vThe meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvA structure pointing to various thermodynamic variables.
[in]dayThe current model time.
[in]nThe time step number of the current execution.
csThe control structure returned by a previous call to MOM_sum_output_init.
tracer_csptracer control structure.
obcOpen boundaries control structure.
[in]dt_forcingThe forcing time step

Definition at line 304 of file MOM_sum_output.F90.

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

References mom_tracer_flow_control::call_tracer_stocks(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), mom_open_boundary::obc_direction_e, mom_open_boundary::obc_direction_n, mom_open_boundary::obc_direction_s, and mom_open_boundary::obc_direction_w.

Here is the call graph for this function:

Variable Documentation

◆ area_chksum_attr

character (*), parameter mom_sum_output::area_chksum_attr = "mask2dT_areaT_checksum"
private

Checksum attribute of name of Gmask2dT * GareaT over the compute domain.

Definition at line 47 of file MOM_sum_output.F90.

47 character (*), parameter :: area_chksum_attr = "mask2dT_areaT_checksum"

Referenced by read_depth_list(), and write_depth_list().

◆ depth_chksum_attr

character (*), parameter mom_sum_output::depth_chksum_attr = "bathyT_checksum"
private

Checksum attribute name of GbathyT over the compute domain.

Definition at line 44 of file MOM_sum_output.F90.

44 character (*), parameter :: depth_chksum_attr = "bathyT_checksum"

Referenced by read_depth_list(), and write_depth_list().

◆ num_fields

integer, parameter mom_sum_output::num_fields = 17
private

Number of diagnostic fields.

Definition at line 43 of file MOM_sum_output.F90.

43 integer, parameter :: NUM_FIELDS = 17 !< Number of diagnostic fields
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54