MOM6
Data Types | Functions/Subroutines
mom_diagnostics Module Reference

Data Types

type  diagnostics_cs
 

Functions/Subroutines

subroutine, public calculate_diagnostic_fields (u, v, h, uh, vh, tv, ADp, CDp, fluxes, dt, G, GV, CS, eta_bt)
 Diagnostics not more naturally calculated elsewhere are computed here. More...
 
subroutine find_weights (Rlist, R_in, k, nz, wt, wt_p)
 This subroutine finds location of R_in in an increasing ordered list, Rlist, returning as k the element such that Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear weights that should be assigned to elements k and k+1. More...
 
subroutine calculate_vertical_integrals (h, tv, fluxes, G, GV, CS)
 Subroutine calculates vertical integrals of several tracers, along with the mass-weight of these tracers, the total column mass, and the carefully calculated column height. More...
 
subroutine calculate_energy_diagnostics (u, v, h, uh, vh, ADp, CDp, G, CS)
 This subroutine calculates terms in the mechanical energy budget. More...
 
subroutine, public register_time_deriv (f_ptr, deriv_ptr, CS)
 This subroutine registers fields to calculate a diagnostic time derivative. More...
 
subroutine calculate_derivs (dt, G, CS)
 This subroutine calculates all registered time derivatives. More...
 
subroutine, public mom_diagnostics_init (MIS, ADp, CDp, Time, G, GV, param_file, diag, CS)
 
subroutine set_dependent_diagnostics (MIS, ADp, CDp, G, CS)
 This subroutine sets up diagnostics upon which other diagnostics depend. More...
 
subroutine, public mom_diagnostics_end (CS, ADp)
 

Function/Subroutine Documentation

◆ calculate_derivs()

subroutine mom_diagnostics::calculate_derivs ( real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(diagnostics_cs), intent(inout)  CS 
)
private

This subroutine calculates all registered time derivatives.

Parameters
[in]dtThe time interval over which differences occur, in s.
[in,out]gThe ocean's grid structure.
[in,out]csControl structure returned by previous call to diagnostics_init.

Definition at line 1086 of file MOM_diagnostics.F90.

Referenced by calculate_diagnostic_fields().

1086  real, intent(in) :: dt !< The time interval over which differences occur,
1087  !! in s.
1088  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
1089  type(diagnostics_cs), intent(inout) :: cs !< Control structure returned by previous call to
1090  !! diagnostics_init.
1091 
1092 ! This subroutine calculates all registered time derivatives.
1093 ! Arguments:
1094 ! (in) dt - time interval in s over which differences occur
1095 ! (in) G - ocean grid structure.
1096 ! (in) CS - control structure returned by previous call to diagnostics_init
1097 
1098  integer i, j, k, m
1099  real idt
1100 
1101  if (dt > 0.0) then ; idt = 1.0/dt
1102  else ; return ; endif
1103 
1104  do m=1,cs%num_time_deriv
1105  do k=1,cs%nlay(m) ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1106  cs%deriv(m)%p(i,j,k) = (cs%var_ptr(m)%p(i,j,k) - cs%prev_val(m)%p(i,j,k)) * idt
1107  cs%prev_val(m)%p(i,j,k) = cs%var_ptr(m)%p(i,j,k)
1108  enddo ; enddo ; enddo
1109  enddo
1110 
Here is the caller graph for this function:

◆ calculate_diagnostic_fields()

subroutine, public mom_diagnostics::calculate_diagnostic_fields ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  uh,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  vh,
type(thermo_var_ptrs), intent(in)  tv,
type(accel_diag_ptrs), intent(in)  ADp,
type(cont_diag_ptrs), intent(in)  CDp,
type(forcing), intent(in)  fluxes,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(diagnostics_cs), intent(inout)  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  eta_bt 
)

Diagnostics not more naturally calculated elsewhere are computed here.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]uThe zonal velocity, in m s-1.
[in]vThe meridional velocity, in m s-1.
[in]hLayer thicknesses, in H (usually m or kg m-2).
[in]uhTransport through zonal faces = u*h*dy, m3/s(Bouss) kg/s(non-Bouss).
[in]vhtransport through meridional faces = v*h*dx, m3/s(Bouss) kg/s(non-Bouss).
[in]tvA structure pointing to various thermodynamic variables.
[in]adpstructure with pointers to accelerations in momentum equation.
[in]cdpstructure with pointers to terms in continuity equation.
[in]fluxesA structure containing the surface fluxes.
[in]dtThe time difference in s since the last call to this subroutine.
[in,out]csControl structure returned by a previous call to diagnostics_init.
[in]eta_btAn optional barotropic variable that gives the "correct" free surface height (Boussinesq) or total water column mass per unit area (non-Boussinesq). This is used to dilate the layer thicknesses when calculating interface heights, in m or kg m-2.

Definition at line 174 of file MOM_diagnostics.F90.

References calculate_derivs(), calculate_energy_diagnostics(), calculate_vertical_integrals(), find_weights(), mom_spatial_means::global_area_mean(), mom_spatial_means::global_layer_mean(), mom_spatial_means::global_volume_mean(), mom_error_handler::mom_error(), and mom_wave_speed::wave_speed().

174  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
175  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
176  !! structure.
177  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1.
178  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity,
179  !! in m s-1.
180  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H
181  !! (usually m or kg m-2).
182  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh !< Transport through zonal faces
183  !! = u*h*dy, m3/s(Bouss)
184  !! kg/s(non-Bouss).
185  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh !< transport through meridional
186  !! faces = v*h*dx, m3/s(Bouss)
187  !! kg/s(non-Bouss).
188  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to
189  !! various thermodynamic
190  !! variables.
191  type(accel_diag_ptrs), intent(in) :: adp !< structure with pointers to
192  !! accelerations in momentum
193  !! equation.
194  type(cont_diag_ptrs), intent(in) :: cdp !< structure with pointers to
195  !! terms in continuity equation.
196  type(forcing), intent(in) :: fluxes !< A structure containing the
197  !! surface fluxes.
198  real, intent(in) :: dt !< The time difference in s since
199  !! the last call to this
200  !! subroutine.
201  type(diagnostics_cs), intent(inout) :: cs !< Control structure returned by
202  !! a previous call to
203  !! diagnostics_init.
204  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< An optional barotropic
205  !! variable that gives the "correct" free surface height (Boussinesq) or total water column
206  !! mass per unit area (non-Boussinesq). This is used to dilate the layer thicknesses when
207  !! calculating interface heights, in m or kg m-2.
208 
209 ! Diagnostics not more naturally calculated elsewhere are computed here.
210 
211 ! Arguments:
212 ! (in) u - zonal velocity component (m/s).
213 ! (in) v - meridional velocity component (m/s).
214 ! (in) h - layer thickness, meter(Bouss) kg/m^2(non-Bouss).
215 ! (in) uh - transport through zonal faces = u*h*dy, m3/s(Bouss) kg/s(non-Bouss).
216 ! (in) vh - transport through meridional faces = v*h*dx, m3/s(Bouss) kg/s(non-Bouss).
217 ! (in) tv - structure pointing to various thermodynamic variables.
218 ! (in) ADp - structure with pointers to accelerations in momentum equation.
219 ! (in) CDp - structure with pointers to terms in continuity equation.
220 ! (in) dt - time difference in s since the last call to this subroutine.
221 ! (in) G - ocean grid structure.
222 ! (in) GV - The ocean's vertical grid structure.
223 ! (in) CS - control structure returned by a previous call to diagnostics_init.
224 ! (in,opt) eta_bt - An optional barotropic variable that gives the "correct"
225 ! free surface height (Boussinesq) or total water column
226 ! mass per unit area (non-Boussinesq). This is used to
227 ! dilate the layer thicknesses when calculating interface
228 ! heights, in m or kg m-2.
229 
230  integer i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
231 
232  ! coordinate variable potential density, in kg m-3.
233  real :: rcv(szi_(g),szj_(g),szk_(g))
234 
235  ! tmp array for surface properties
236  real :: surface_field(szi_(g),szj_(g))
237  real :: pressure_1d(szi_(g)) ! Temporary array for pressure when calling EOS
238  real :: wt, wt_p
239 
240  ! squared Coriolis parameter at to h-points (1/s2)
241  real :: f2_h
242 
243  ! magnitude of the gradient of f (1/(m*s))
244  real :: mag_beta
245 
246  ! frequency squared used to avoid division by 0 (1/s2)
247  ! value is roughly (pi / (the age of the universe) )^2.
248  real, parameter :: absurdly_small_freq2 = 1e-34
249 
250  integer :: k_list
251 
252  real, dimension(SZK_(G)) :: temp_layer_ave, salt_layer_ave
253  real :: thetaoga, soga, masso, tosga, sosga
254 
255  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
256  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
257  nz = g%ke ; nkmb = gv%nk_rho_varies
258 
259  ! smg: is the following robust to ALE? It seems a bit opaque.
260  ! If the model is NOT in isopycnal mode then nkmb=0. But we need all the
261  ! following diagnostics to treat all layers as variable density, so we set
262  ! nkmb = nz, on the expectation that loops nkmb+1,nz will not iterate.
263  ! This behavior is ANSI F77 but some compiler options can force at least
264  ! one iteration that would break the following one-line workaround!
265  if (nkmb==0) nkmb = nz
266 
267  if (loc(cs)==0) call mom_error(fatal, &
268  "calculate_diagnostic_fields: Module must be initialized before used.")
269 
270  call calculate_derivs(dt, g, cs)
271 
272  if (ASSOCIATED(cs%e)) then
273  call find_eta(h, tv, gv%g_Earth, g, gv, cs%e, eta_bt)
274  if (cs%id_e > 0) call post_data(cs%id_e, cs%e, cs%diag)
275  endif
276 
277  if (ASSOCIATED(cs%e_D)) then
278  if (ASSOCIATED(cs%e)) then
279  do k=1,nz+1 ; do j=js,je ; do i=is,ie
280  cs%e_D(i,j,k) = cs%e(i,j,k) + g%bathyT(i,j)
281  enddo ; enddo ; enddo
282  else
283  call find_eta(h, tv, gv%g_Earth, g, gv, cs%e_D, eta_bt)
284  do k=1,nz+1 ; do j=js,je ; do i=is,ie
285  cs%e_D(i,j,k) = cs%e_D(i,j,k) + g%bathyT(i,j)
286  enddo ; enddo ; enddo
287  endif
288 
289  if (cs%id_e_D > 0) call post_data(cs%id_e_D, cs%e_D, cs%diag)
290  endif
291 
292  ! mass per area of grid cell (for Bouss, use Rho0)
293  if (cs%id_masscello > 0) then
294  do k=1,nz; do j=js,je ; do i=is,ie
295  cs%diag_tmp3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)
296  enddo ; enddo ; enddo
297  call post_data(cs%id_masscello, cs%diag_tmp3d, cs%diag)
298  endif
299 
300  ! mass of liquid ocean (for Bouss, use Rho0)
301  if (cs%id_masso > 0) then
302  do k=1,nz; do j=js,je ; do i=is,ie
303  cs%diag_tmp3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)*g%areaT(i,j)
304  enddo ; enddo ; enddo
305  masso = (reproducing_sum(sum(cs%diag_tmp3d,3)))
306  call post_data(cs%id_masso, masso, cs%diag)
307  endif
308 
309  ! diagnose thickness of grid cells (meter)
310  if (cs%id_thkcello > 0) then
311 
312  ! thkcello = h for Boussinesq
313  if (gv%Boussinesq) then
314  call post_data(cs%id_thkcello, gv%H_to_m*h, cs%diag)
315 
316  ! thkcello = dp/(rho*g) for non-Boussinesq
317  else
318  do j=js,je
319 
320  ! pressure loading at top of surface layer (Pa)
321  if(ASSOCIATED(fluxes%p_surf)) then
322  do i=is,ie
323  pressure_1d(i) = fluxes%p_surf(i,j)
324  enddo
325  else
326  do i=is,ie
327  pressure_1d(i) = 0.0
328  enddo
329  endif
330 
331  do k=1,nz
332  ! pressure for EOS at the layer center (Pa)
333  do i=is,ie
334  pressure_1d(i) = pressure_1d(i) + 0.5*(gv%g_Earth*gv%H_to_kg_m2)*h(i,j,k)
335  enddo
336  ! store in-situ density (kg/m3) in diag_tmp3d
337  call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
338  cs%diag_tmp3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
339  ! cell thickness = dz = dp/(g*rho) (meter); store in diag_tmp3d
340  do i=is,ie
341  cs%diag_tmp3d(i,j,k) = (gv%H_to_kg_m2*h(i,j,k))/cs%diag_tmp3d(i,j,k)
342  enddo
343  ! pressure for EOS at the bottom interface (Pa)
344  do i=is,ie
345  pressure_1d(i) = pressure_1d(i) + 0.5*(gv%g_Earth*gv%H_to_kg_m2)*h(i,j,k)
346  enddo
347  enddo ! k
348 
349  enddo ! j
350  call post_data(cs%id_thkcello, cs%diag_tmp3d, cs%diag)
351  endif
352  endif
353 
354  ! volume mean potential temperature
355  if (cs%id_thetaoga>0) then
356  thetaoga = global_volume_mean(tv%T, h, g, gv)
357  call post_data(cs%id_thetaoga, thetaoga, cs%diag)
358  endif
359 
360  ! area mean SST
361  if (cs%id_tosga > 0) then
362  do j=js,je ; do i=is,ie
363  surface_field(i,j) = tv%T(i,j,1)
364  enddo ; enddo
365  tosga = global_area_mean(surface_field, g)
366  call post_data(cs%id_tosga, tosga, cs%diag)
367  endif
368 
369  ! volume mean salinity
370  if (cs%id_soga>0) then
371  soga = global_volume_mean(tv%S, h, g, gv)
372  call post_data(cs%id_soga, soga, cs%diag)
373  endif
374 
375  ! area mean SSS
376  if (cs%id_sosga > 0) then
377  do j=js,je ; do i=is,ie
378  surface_field(i,j) = tv%S(i,j,1)
379  enddo ; enddo
380  sosga = global_area_mean(surface_field, g)
381  call post_data(cs%id_sosga, sosga, cs%diag)
382  endif
383 
384  ! layer mean potential temperature
385  if (cs%id_temp_layer_ave>0) then
386  temp_layer_ave = global_layer_mean(tv%T, h, g, gv)
387  call post_data_1d_k(cs%id_temp_layer_ave, temp_layer_ave, cs%diag)
388  endif
389 
390  ! layer mean salinity
391  if (cs%id_salt_layer_ave>0) then
392  salt_layer_ave = global_layer_mean(tv%S, h, g, gv)
393  call post_data_1d_k(cs%id_salt_layer_ave, salt_layer_ave, cs%diag)
394  endif
395 
396  call calculate_vertical_integrals(h, tv, fluxes, g, gv, cs)
397 
398  if ((cs%id_Rml > 0) .or. (cs%id_Rcv > 0) .or. ASSOCIATED(cs%h_Rlay) .or. &
399  ASSOCIATED(cs%uh_Rlay) .or. ASSOCIATED(cs%vh_Rlay) .or. &
400  ASSOCIATED(cs%uhGM_Rlay) .or. ASSOCIATED(cs%vhGM_Rlay)) then
401 
402  if (associated(tv%eqn_of_state)) then
403  pressure_1d(:) = tv%P_Ref
404 !$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
405  do k=1,nz ; do j=js-1,je+1
406  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, &
407  rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state)
408  enddo ; enddo
409  else ! Rcv should not be used much in this case, so fill in sensible values.
410  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
411  rcv(i,j,k) = gv%Rlay(k)
412  enddo ; enddo ; enddo
413  endif
414  if (cs%id_Rml > 0) call post_data(cs%id_Rml, rcv, cs%diag)
415  if (cs%id_Rcv > 0) call post_data(cs%id_Rcv, rcv, cs%diag)
416 
417  if (ASSOCIATED(cs%h_Rlay)) then
418  k_list = nz/2
419 !$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,CS,Rcv,h,GV) &
420 !$OMP private(wt,wt_p) firstprivate(k_list)
421  do j=js,je
422  do k=1,nkmb ; do i=is,ie
423  cs%h_Rlay(i,j,k) = 0.0
424  enddo ; enddo
425  do k=nkmb+1,nz ; do i=is,ie
426  cs%h_Rlay(i,j,k) = h(i,j,k)
427  enddo ; enddo
428  do k=1,nkmb ; do i=is,ie
429  call find_weights(gv%Rlay, rcv(i,j,k), k_list, nz, wt, wt_p)
430  cs%h_Rlay(i,j,k_list) = cs%h_Rlay(i,j,k_list) + h(i,j,k)*wt
431  cs%h_Rlay(i,j,k_list+1) = cs%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
432  enddo ; enddo
433  enddo
434 
435  if (cs%id_h_Rlay > 0) call post_data(cs%id_h_Rlay, cs%h_Rlay, cs%diag)
436  endif
437 
438  if (ASSOCIATED(cs%uh_Rlay)) then
439  k_list = nz/2
440 !$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CS,GV,uh) &
441 !$OMP private(wt,wt_p) firstprivate(k_list)
442  do j=js,je
443  do k=1,nkmb ; do i=isq,ieq
444  cs%uh_Rlay(i,j,k) = 0.0
445  enddo ; enddo
446  do k=nkmb+1,nz ; do i=isq,ieq
447  cs%uh_Rlay(i,j,k) = uh(i,j,k)
448  enddo ; enddo
449  k_list = nz/2
450  do k=1,nkmb ; do i=isq,ieq
451  call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
452  cs%uh_Rlay(i,j,k_list) = cs%uh_Rlay(i,j,k_list) + uh(i,j,k)*wt
453  cs%uh_Rlay(i,j,k_list+1) = cs%uh_Rlay(i,j,k_list+1) + uh(i,j,k)*wt_p
454  enddo ; enddo
455  enddo
456 
457  if (cs%id_uh_Rlay > 0) call post_data(cs%id_uh_Rlay, cs%uh_Rlay, cs%diag)
458  endif
459 
460  if (ASSOCIATED(cs%vh_Rlay)) then
461  k_list = nz/2
462 !$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,nz,nkmb,Rcv,CS,GV,vh) &
463 !$OMP private(wt,wt_p) firstprivate(k_list)
464  do j=jsq,jeq
465  do k=1,nkmb ; do i=is,ie
466  cs%vh_Rlay(i,j,k) = 0.0
467  enddo ; enddo
468  do k=nkmb+1,nz ; do i=is,ie
469  cs%vh_Rlay(i,j,k) = vh(i,j,k)
470  enddo ; enddo
471  do k=1,nkmb ; do i=is,ie
472  call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
473  cs%vh_Rlay(i,j,k_list) = cs%vh_Rlay(i,j,k_list) + vh(i,j,k)*wt
474  cs%vh_Rlay(i,j,k_list+1) = cs%vh_Rlay(i,j,k_list+1) + vh(i,j,k)*wt_p
475  enddo ; enddo
476  enddo
477 
478  if (cs%id_vh_Rlay > 0) call post_data(cs%id_vh_Rlay, cs%vh_Rlay, cs%diag)
479  endif
480 
481  if (ASSOCIATED(cs%uhGM_Rlay) .and. ASSOCIATED(cdp%uhGM)) then
482  k_list = nz/2
483 !$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CDP,CS,GV) &
484 !$OMP private(wt,wt_p) firstprivate(k_list)
485  do j=js,je
486  do k=1,nkmb ; do i=isq,ieq
487  cs%uhGM_Rlay(i,j,k) = 0.0
488  enddo ; enddo
489  do k=nkmb+1,nz ; do i=isq,ieq
490  cs%uhGM_Rlay(i,j,k) = cdp%uhGM(i,j,k)
491  enddo ; enddo
492  do k=1,nkmb ; do i=isq,ieq
493  call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
494  cs%uhGM_Rlay(i,j,k_list) = cs%uhGM_Rlay(i,j,k_list) + cdp%uhGM(i,j,k)*wt
495  cs%uhGM_Rlay(i,j,k_list+1) = cs%uhGM_Rlay(i,j,k_list+1) + cdp%uhGM(i,j,k)*wt_p
496  enddo ; enddo
497  enddo
498 
499  if (cs%id_uh_Rlay > 0) call post_data(cs%id_uhGM_Rlay, cs%uhGM_Rlay, cs%diag)
500  endif
501 
502  if (ASSOCIATED(cs%vhGM_Rlay) .and. ASSOCIATED(cdp%vhGM)) then
503  k_list = nz/2
504 !$OMP parallel do default(none) shared(is,ie,Jsq,Jeq,nz,nkmb,CS,CDp,Rcv,GV) &
505 !$OMP private(wt,wt_p) firstprivate(k_list)
506  do j=jsq,jeq
507  do k=1,nkmb ; do i=is,ie
508  cs%vhGM_Rlay(i,j,k) = 0.0
509  enddo ; enddo
510  do k=nkmb+1,nz ; do i=is,ie
511  cs%vhGM_Rlay(i,j,k) = cdp%vhGM(i,j,k)
512  enddo ; enddo
513  do k=1,nkmb ; do i=is,ie
514  call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
515  cs%vhGM_Rlay(i,j,k_list) = cs%vhGM_Rlay(i,j,k_list) + cdp%vhGM(i,j,k)*wt
516  cs%vhGM_Rlay(i,j,k_list+1) = cs%vhGM_Rlay(i,j,k_list+1) + cdp%vhGM(i,j,k)*wt_p
517  enddo ; enddo
518  enddo
519 
520  if (cs%id_vhGM_Rlay > 0) call post_data(cs%id_vhGM_Rlay, cs%vhGM_Rlay, cs%diag)
521  endif
522  endif
523 
524  if (associated(tv%eqn_of_state)) then
525  if (cs%id_rhopot0 > 0) then
526  pressure_1d(:) = 0.
527 !$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
528  do k=1,nz ; do j=js,je
529  call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
530  rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
531  enddo ; enddo
532  if (cs%id_rhopot0 > 0) call post_data(cs%id_rhopot0, rcv, cs%diag)
533  endif
534  if (cs%id_rhopot2 > 0) then
535  pressure_1d(:) = 2.e7 ! 2000 dbars
536 !$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
537  do k=1,nz ; do j=js,je
538  call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
539  rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
540  enddo ; enddo
541  if (cs%id_rhopot2 > 0) call post_data(cs%id_rhopot2, rcv, cs%diag)
542  endif
543  if (cs%id_rhoinsitu > 0) then
544 !$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d,h,GV)
545  do j=js,je
546  pressure_1d(:) = 0. ! Start at p=0 Pa at surface
547  do k=1,nz
548  pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure in middle of layer k
549  call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
550  rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
551  pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure at bottom of layer k
552  enddo
553  enddo
554  if (cs%id_rhoinsitu > 0) call post_data(cs%id_rhoinsitu, rcv, cs%diag)
555  endif
556  endif
557 
558  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
559  (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0)) then
560  call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp)
561  if (cs%id_cg1>0) call post_data(cs%id_cg1, cs%cg1, cs%diag)
562  if (cs%id_Rd1>0) then
563 !$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
564 !$OMP private(f2_h,mag_beta)
565  do j=js,je ; do i=is,ie
566  ! Blend the equatorial deformation radius with the standard one.
567  f2_h = absurdly_small_freq2 + 0.25 * &
568  ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
569  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
570  mag_beta = sqrt(0.5 * ( &
571  (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
572  ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
573  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
574  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
575  cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
576 
577  enddo ; enddo
578  call post_data(cs%id_Rd1, cs%Rd1, cs%diag)
579  endif
580  if (cs%id_cfl_cg1>0) then
581  do j=js,je ; do i=is,ie
582  cs%cfl_cg1(i,j) = (dt*cs%cg1(i,j)) * (g%IdxT(i,j) + g%IdyT(i,j))
583  enddo ; enddo
584  call post_data(cs%id_cfl_cg1, cs%cfl_cg1, cs%diag)
585  endif
586  if (cs%id_cfl_cg1_x>0) then
587  do j=js,je ; do i=is,ie
588  cs%cfl_cg1_x(i,j) = (dt*cs%cg1(i,j)) * g%IdxT(i,j)
589  enddo ; enddo
590  call post_data(cs%id_cfl_cg1_x, cs%cfl_cg1_x, cs%diag)
591  endif
592  if (cs%id_cfl_cg1_y>0) then
593  do j=js,je ; do i=is,ie
594  cs%cfl_cg1_y(i,j) = (dt*cs%cg1(i,j)) * g%IdyT(i,j)
595  enddo ; enddo
596  call post_data(cs%id_cfl_cg1_y, cs%cfl_cg1_y, cs%diag)
597  endif
598  endif
599  if ((cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
600  if (cs%id_p_ebt>0) then
601  call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
602  mono_n2_column_fraction=cs%mono_N2_column_fraction, &
603  mono_n2_depth=cs%mono_N2_depth, modal_structure=cs%p_ebt)
604  call post_data(cs%id_p_ebt, cs%p_ebt, cs%diag)
605  else
606  call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
607  mono_n2_column_fraction=cs%mono_N2_column_fraction, &
608  mono_n2_depth=cs%mono_N2_depth)
609  endif
610  if (cs%id_cg_ebt>0) call post_data(cs%id_cg_ebt, cs%cg1, cs%diag)
611  if (cs%id_Rd_ebt>0) then
612 !$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
613 !$OMP private(f2_h,mag_beta)
614  do j=js,je ; do i=is,ie
615  ! Blend the equatorial deformation radius with the standard one.
616  f2_h = absurdly_small_freq2 + 0.25 * &
617  ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
618  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
619  mag_beta = sqrt(0.5 * ( &
620  (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
621  ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
622  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
623  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
624  cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
625 
626  enddo ; enddo
627  call post_data(cs%id_Rd_ebt, cs%Rd1, cs%diag)
628  endif
629  endif
630 
631  if (dt > 0.0) then
632  if (cs%id_du_dt>0) call post_data(cs%id_du_dt, cs%du_dt, cs%diag)
633 
634  if (cs%id_dv_dt>0) call post_data(cs%id_dv_dt, cs%dv_dt, cs%diag)
635 
636  if (cs%id_dh_dt>0) call post_data(cs%id_dh_dt, cs%dh_dt, cs%diag)
637 
638  call calculate_energy_diagnostics(u, v, h, uh, vh, adp, cdp, g, cs)
639  endif
640 
Here is the call graph for this function:

◆ calculate_energy_diagnostics()

subroutine mom_diagnostics::calculate_energy_diagnostics ( 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,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  uh,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  vh,
type(accel_diag_ptrs), intent(in)  ADp,
type(cont_diag_ptrs), intent(in)  CDp,
type(ocean_grid_type), intent(inout)  G,
type(diagnostics_cs), intent(inout)  CS 
)
private

This subroutine calculates terms in the mechanical energy budget.

Parameters
[in,out]gThe ocean's grid structure.
[in]uThe zonal velocity, in m s-1.
[in]vThe meridional velocity, in m s-1.
[in]hLayer thicknesses, in H (usually m or kg m-2).
[in]uhTransport through zonal faces=u*h*dy: m3/s (Bouss) kg/s(non-Bouss).
[in]vhTransport through merid faces=v*h*dx: m3/s (Bouss) kg/s(non-Bouss).
[in]adpStructure pointing to accelerations in momentum equation.
[in]cdpStructure pointing to terms in continuity equations.
[in,out]csControl structure returned by a previous call to diagnostics_init.

Definition at line 839 of file MOM_diagnostics.F90.

Referenced by calculate_diagnostic_fields().

839  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
840  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1.
841  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity,
842  !! in m s-1.
843  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H
844  !! (usually m or kg m-2).
845  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh !< Transport through zonal
846  !! faces=u*h*dy: m3/s (Bouss)
847  !! kg/s(non-Bouss).
848  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh !< Transport through merid
849  !! faces=v*h*dx: m3/s (Bouss)
850  !! kg/s(non-Bouss).
851  type(accel_diag_ptrs), intent(in) :: adp !< Structure pointing to
852  !! accelerations in momentum
853  !! equation.
854  type(cont_diag_ptrs), intent(in) :: cdp !< Structure pointing to terms
855  !! in continuity equations.
856  type(diagnostics_cs), intent(inout) :: cs !< Control structure returned by
857  !! a previous call to
858  !! diagnostics_init.
859 
860 ! This subroutine calculates terms in the mechanical energy budget.
861 
862 ! Arguments:
863 ! (in) u - zonal velocity component (m/s)
864 ! (in) v - meridional velocity componnent (m/s)
865 ! (in) h - layer thickness: metre(Bouss) of kg/m2(non-Bouss)
866 ! (in) uh - transport through zonal faces=u*h*dy: m3/s (Bouss) kg/s(non-Bouss)
867 ! (in) vh - transport through merid faces=v*h*dx: m3/s (Bouss) kg/s(non-Bouss)
868 ! (in) ADp - structure pointing to accelerations in momentum equation
869 ! (in) CDp - structure pointing to terms in continuity equations
870 ! (in) G - ocean grid structure
871 ! (in) CS - control structure returned by a previous call to diagnostics_init
872 
873  real :: ke_u(szib_(g),szj_(g))
874  real :: ke_v(szi_(g),szjb_(g))
875  real :: ke_h(szi_(g),szj_(g))
876 
877  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
878  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
879  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
880 
881  do j=js-1,je ; do i=is-1,ie
882  ke_u(i,j) = 0.0 ; ke_v(i,j) = 0.0
883  enddo ; enddo
884 
885  if (ASSOCIATED(cs%KE)) then
886  do k=1,nz ; do j=js,je ; do i=is,ie
887  cs%KE(i,j,k) = ((u(i,j,k)*u(i,j,k) + u(i-1,j,k)*u(i-1,j,k)) + &
888  (v(i,j,k)*v(i,j,k) + v(i,j-1,k)*v(i,j-1,k)))*0.25
889  ! DELETE THE FOLLOWING... Make this 0 to test the momentum balance,
890  ! or a huge number to test the continuity balance.
891  ! CS%KE(i,j,k) *= 1e20
892  enddo ; enddo ; enddo
893  if (cs%id_KE > 0) call post_data(cs%id_KE, cs%KE, cs%diag)
894  endif
895 
896  if(.not.g%symmetric) then
897  if(ASSOCIATED(cs%dKE_dt) .OR. ASSOCIATED(cs%PE_to_KE) .OR. ASSOCIATED(cs%KE_CorAdv) .OR. &
898  ASSOCIATED(cs%KE_adv) .OR. ASSOCIATED(cs%KE_visc) .OR. ASSOCIATED(cs%KE_horvisc).OR. &
899  ASSOCIATED(cs%KE_dia) ) then
900  call create_group_pass(cs%pass_KE_uv, ke_u, ke_v, g%Domain, to_north+to_east)
901  endif
902  endif
903 
904  if (ASSOCIATED(cs%dKE_dt)) then
905  do k=1,nz
906  do j=js,je ; do i=isq,ieq
907  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*cs%du_dt(i,j,k)
908  enddo ; enddo
909  do j=jsq,jeq ; do i=is,ie
910  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*cs%dv_dt(i,j,k)
911  enddo ; enddo
912  do j=js,je ; do i=is,ie
913  ke_h(i,j) = cs%KE(i,j,k)*cs%dh_dt(i,j,k)
914  enddo ; enddo
915  if (.not.g%symmetric) &
916  call do_group_pass(cs%pass_KE_uv, g%domain)
917  do j=js,je ; do i=is,ie
918  cs%dKE_dt(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
919  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
920  enddo ; enddo
921  enddo
922  if (cs%id_dKEdt > 0) call post_data(cs%id_dKEdt, cs%dKE_dt, cs%diag)
923  endif
924 
925  if (ASSOCIATED(cs%PE_to_KE)) then
926  do k=1,nz
927  do j=js,je ; do i=isq,ieq
928  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%PFu(i,j,k)
929  enddo ; enddo
930  do j=jsq,jeq ; do i=is,ie
931  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%PFv(i,j,k)
932  enddo ; enddo
933  if (.not.g%symmetric) &
934  call do_group_pass(cs%pass_KE_uv, g%domain)
935  do j=js,je ; do i=is,ie
936  cs%PE_to_KE(i,j,k) = 0.5 * g%IareaT(i,j) * &
937  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
938  enddo ; enddo
939  enddo
940  if (cs%id_PE_to_KE > 0) call post_data(cs%id_PE_to_KE, cs%PE_to_KE, cs%diag)
941  endif
942 
943  if (ASSOCIATED(cs%KE_CorAdv)) then
944  do k=1,nz
945  do j=js,je ; do i=isq,ieq
946  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%CAu(i,j,k)
947  enddo ; enddo
948  do j=jsq,jeq ; do i=is,ie
949  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%CAv(i,j,k)
950  enddo ; enddo
951  do j=js,je ; do i=is,ie
952  ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
953  (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
954  enddo ; enddo
955  if (.not.g%symmetric) &
956  call do_group_pass(cs%pass_KE_uv, g%domain)
957  do j=js,je ; do i=is,ie
958  cs%KE_CorAdv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
959  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
960  enddo ; enddo
961  enddo
962  if (cs%id_KE_Coradv > 0) call post_data(cs%id_KE_Coradv, cs%KE_Coradv, cs%diag)
963  endif
964 
965  if (ASSOCIATED(cs%KE_adv)) then
966  do k=1,nz
967  do j=js,je ; do i=isq,ieq
968  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%gradKEu(i,j,k)
969  enddo ; enddo
970  do j=jsq,jeq ; do i=is,ie
971  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%gradKEv(i,j,k)
972  enddo ; enddo
973  do j=js,je ; do i=is,ie
974  ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
975  (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
976  enddo ; enddo
977  if (.not.g%symmetric) &
978  call do_group_pass(cs%pass_KE_uv, g%domain)
979  do j=js,je ; do i=is,ie
980  cs%KE_adv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
981  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
982  enddo ; enddo
983  enddo
984  if (cs%id_KE_adv > 0) call post_data(cs%id_KE_adv, cs%KE_adv, cs%diag)
985  endif
986 
987  if (ASSOCIATED(cs%KE_visc)) then
988  do k=1,nz
989  do j=js,je ; do i=isq,ieq
990  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_visc(i,j,k)
991  enddo ; enddo
992  do j=jsq,jeq ; do i=is,ie
993  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_visc(i,j,k)
994  enddo ; enddo
995  if (.not.g%symmetric) &
996  call do_group_pass(cs%pass_KE_uv, g%domain)
997  do j=js,je ; do i=is,ie
998  cs%KE_visc(i,j,k) = 0.5 * g%IareaT(i,j) * &
999  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
1000  enddo ; enddo
1001  enddo
1002  if (cs%id_KE_visc > 0) call post_data(cs%id_KE_visc, cs%KE_visc, cs%diag)
1003  endif
1004 
1005  if (ASSOCIATED(cs%KE_horvisc)) then
1006  do k=1,nz
1007  do j=js,je ; do i=isq,ieq
1008  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%diffu(i,j,k)
1009  enddo ; enddo
1010  do j=jsq,jeq ; do i=is,ie
1011  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%diffv(i,j,k)
1012  enddo ; enddo
1013  if (.not.g%symmetric) &
1014  call do_group_pass(cs%pass_KE_uv, g%domain)
1015  do j=js,je ; do i=is,ie
1016  cs%KE_horvisc(i,j,k) = 0.5 * g%IareaT(i,j) * &
1017  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
1018  enddo ; enddo
1019  enddo
1020  if (cs%id_KE_horvisc > 0) call post_data(cs%id_KE_horvisc, cs%KE_horvisc, cs%diag)
1021  endif
1022 
1023  if (ASSOCIATED(cs%KE_dia)) then
1024  do k=1,nz
1025  do j=js,je ; do i=isq,ieq
1026  ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_dia(i,j,k)
1027  enddo ; enddo
1028  do j=jsq,jeq ; do i=is,ie
1029  ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_dia(i,j,k)
1030  enddo ; enddo
1031  do j=js,je ; do i=is,ie
1032  ke_h(i,j) = cs%KE(i,j,k) * &
1033  (cdp%diapyc_vel(i,j,k) - cdp%diapyc_vel(i,j,k+1))
1034  enddo ; enddo
1035  if (.not.g%symmetric) &
1036  call do_group_pass(cs%pass_KE_uv, g%domain)
1037  do j=js,je ; do i=is,ie
1038  cs%KE_dia(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
1039  (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
1040  enddo ; enddo
1041  enddo
1042  if (cs%id_KE_dia > 0) call post_data(cs%id_KE_dia, cs%KE_dia, cs%diag)
1043  endif
1044 
Here is the caller graph for this function:

◆ calculate_vertical_integrals()

subroutine mom_diagnostics::calculate_vertical_integrals ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(forcing), intent(in)  fluxes,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(diagnostics_cs), intent(inout)  CS 
)
private

Subroutine calculates vertical integrals of several tracers, along with the mass-weight of these tracers, the total column mass, and the carefully calculated column height.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]hLayer thicknesses, in H (usually m or kg m-2).
[in]tvA structure pointing to various thermodynamic variables.
[in]fluxesA structure containing the surface fluxes.
[in,out]csA control structure returned by a previous call to diagnostics_init.

Definition at line 711 of file MOM_diagnostics.F90.

References mom_eos::int_density_dz().

Referenced by calculate_diagnostic_fields().

711  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
712  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
713  !! structure.
714  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H
715  !! (usually m or kg m-2).
716  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
717  !! thermodynamic variables.
718  type(forcing), intent(in) :: fluxes !< A structure containing the
719  !! surface fluxes.
720  type(diagnostics_cs), intent(inout) :: cs !< A control structure returned
721  !! by a previous call to
722  !! diagnostics_init.
723 
724 ! Subroutine calculates vertical integrals of several tracers, along
725 ! with the mass-weight of these tracers, the total column mass, and the
726 ! carefully calculated column height.
727 
728 ! Arguments:
729 ! (in) h - layer thickness: metre (Bouss) or kg/ m2 (non-Bouss)
730 ! (in) tv - structure pointing to thermodynamic variables
731 ! (in) fluxes - a structure containing the surface fluxes.
732 ! (in) G - ocean grid structure
733 ! (in) GV - The ocean's vertical grid structure.
734 ! (in) CS - control structure returned by a previous call to diagnostics_init
735 
736  real, dimension(SZI_(G), SZJ_(G)) :: &
737  z_top, & ! Height of the top of a layer or the ocean, in m.
738  z_bot, & ! Height of the bottom of a layer (for id_mass) or the
739  ! (positive) depth of the ocean (for id_col_ht), in m.
740  mass, & ! integrated mass of the water column, in kg m-2. For
741  ! non-Boussinesq models this is rho*dz. For Boussiensq
742  ! models, this is either the integral of in-situ density
743  ! (rho*dz for col_mass) or reference dens (Rho_0*dz for mass_wt).
744  btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
745  ! This is the column mass multiplied by gravity plus the pressure
746  ! at the ocean surface.
747  dpress, & ! Change in hydrostatic pressure across a layer, in Pa.
748  tr_int ! vertical integral of a tracer times density,
749  ! (Rho_0 in a Boussinesq model) in TR kg m-2.
750  real :: ig_earth ! Inverse of gravitational acceleration, in s2 m-1.
751 
752  integer :: i, j, k, is, ie, js, je, nz
753  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
754 
755  if (cs%id_mass_wt > 0) then
756  do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
757  do k=1,nz ; do j=js,je ; do i=is,ie
758  mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
759  enddo ; enddo ; enddo
760  call post_data(cs%id_mass_wt, mass, cs%diag)
761  endif
762 
763  if (cs%id_temp_int > 0) then
764  do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
765  do k=1,nz ; do j=js,je ; do i=is,ie
766  tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%T(i,j,k)
767  enddo ; enddo ; enddo
768  call post_data(cs%id_temp_int, tr_int, cs%diag)
769  endif
770 
771  if (cs%id_salt_int > 0) then
772  do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
773  do k=1,nz ; do j=js,je ; do i=is,ie
774  tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%S(i,j,k)
775  enddo ; enddo ; enddo
776  call post_data(cs%id_salt_int, tr_int, cs%diag)
777  endif
778 
779  if (cs%id_col_ht > 0) then
780  call find_eta(h, tv, gv%g_Earth, g, gv, z_top)
781  do j=js,je ; do i=is,ie
782  z_bot(i,j) = z_top(i,j) + g%bathyT(i,j)
783  enddo ; enddo
784  call post_data(cs%id_col_ht, z_bot, cs%diag)
785  endif
786 
787  if (cs%id_col_mass > 0 .or. cs%id_pbo > 0) then
788  do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
789  if (gv%Boussinesq) then
790  if (associated(tv%eqn_of_state)) then
791  ig_earth = 1.0 / gv%g_Earth
792 ! do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo
793  do j=js,je ; do i=is,ie ; z_bot(i,j) = 0.0 ; enddo ; enddo
794  do k=1,nz
795  do j=js,je ; do i=is,ie
796  z_top(i,j) = z_bot(i,j)
797  z_bot(i,j) = z_top(i,j) - gv%H_to_m*h(i,j,k)
798  enddo ; enddo
799  call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), &
800  z_top, z_bot, 0.0, gv%H_to_kg_m2, gv%g_Earth, &
801  g%HI, g%HI, tv%eqn_of_state, dpress)
802  do j=js,je ; do i=is,ie
803  mass(i,j) = mass(i,j) + dpress(i,j) * ig_earth
804  enddo ; enddo
805  enddo
806  else
807  do k=1,nz ; do j=js,je ; do i=is,ie
808  mass(i,j) = mass(i,j) + (gv%H_to_m*gv%Rlay(k))*h(i,j,k)
809  enddo ; enddo ; enddo
810  endif
811  else
812  do k=1,nz ; do j=js,je ; do i=is,ie
813  mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
814  enddo ; enddo ; enddo
815  endif
816  if (cs%id_col_mass > 0) then
817  call post_data(cs%id_col_mass, mass, cs%diag)
818  endif
819  if (cs%id_pbo > 0) then
820  do j=js,je ; do i=is,ie ; btm_pres(i,j) = 0.0 ; enddo ; enddo
821  ! 'pbo' is defined as the sea water pressure at the sea floor
822  ! pbo = (mass * g) + pso
823  ! where pso is the sea water pressure at sea water surface
824  ! note that pso is equivalent to fluxes%p_surf
825  do j=js,je ; do i=is,ie
826  btm_pres(i,j) = mass(i,j) * gv%g_Earth
827  if (ASSOCIATED(fluxes%p_surf)) then
828  btm_pres(i,j) = btm_pres(i,j) + fluxes%p_surf(i,j)
829  endif
830  enddo ; enddo
831  call post_data(cs%id_pbo, btm_pres, cs%diag)
832  endif
833  endif
834 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_weights()

subroutine mom_diagnostics::find_weights ( real, dimension(:), intent(in)  Rlist,
real, intent(in)  R_in,
integer, intent(inout)  k,
integer, intent(in)  nz,
real, intent(out)  wt,
real, intent(out)  wt_p 
)
private

This subroutine finds location of R_in in an increasing ordered list, Rlist, returning as k the element such that Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear weights that should be assigned to elements k and k+1.

Definition at line 648 of file MOM_diagnostics.F90.

Referenced by calculate_diagnostic_fields().

648  real, intent(in) :: rlist(:), r_in
649  integer, intent(inout) :: k
650  integer, intent(in) :: nz
651  real, intent(out) :: wt, wt_p
652 
653  ! This subroutine finds location of R_in in an increasing ordered
654  ! list, Rlist, returning as k the element such that
655  ! Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear
656  ! weights that should be assigned to elements k and k+1.
657 
658  integer :: k_upper, k_lower, k_new, inc
659 
660  ! First, bracket the desired point.
661  if ((k < 1) .or. (k > nz)) k = nz/2
662 
663  k_upper = k ; k_lower = k ; inc = 1
664  if (r_in < rlist(k)) then
665  do
666  k_lower = max(k_lower-inc, 1)
667  if ((k_lower == 1) .or. (r_in >= rlist(k_lower))) exit
668  k_upper = k_lower
669  inc = inc*2
670  end do
671  else
672  do
673  k_upper = min(k_upper+inc, nz)
674  if ((k_upper == nz) .or. (r_in < rlist(k_upper))) exit
675  k_lower = k_upper
676  inc = inc*2
677  end do
678  endif
679 
680  if ((k_lower == 1) .and. (r_in <= rlist(k_lower))) then
681  k = 1 ; wt = 1.0 ; wt_p = 0.0
682  else if ((k_upper == nz) .and. (r_in >= rlist(k_upper))) then
683  k = nz-1 ; wt = 0.0 ; wt_p = 1.0
684  else
685  do
686  if (k_upper <= k_lower+1) exit
687  k_new = (k_upper + k_lower) / 2
688  if (r_in < rlist(k_new)) then
689  k_upper = k_new
690  else
691  k_lower = k_new
692  endif
693  end do
694 
695 ! Uncomment this as a code check
696 ! if ((R_in < Rlist(k_lower)) .or. (R_in >= Rlist(k_upper)) .or. (k_upper-k_lower /= 1)) &
697 ! write (*,*) "Error: ",R_in," is not between R(",k_lower,") = ", &
698 ! Rlist(k_lower)," and R(",k_upper,") = ",Rlist(k_upper),"."
699  k = k_lower
700  wt = (rlist(k_upper) - r_in) / (rlist(k_upper) - rlist(k_lower))
701  wt_p = 1.0 - wt
702 
703  endif
704 
Here is the caller graph for this function:

◆ mom_diagnostics_end()

subroutine, public mom_diagnostics::mom_diagnostics_end ( type(diagnostics_cs), pointer  CS,
type(accel_diag_ptrs), intent(inout)  ADp 
)

Definition at line 1456 of file MOM_diagnostics.F90.

1456  type(diagnostics_cs), pointer :: cs
1457  type(accel_diag_ptrs), intent(inout) :: adp
1458  integer :: m
1459 
1460  if (ASSOCIATED(cs%e)) deallocate(cs%e)
1461  if (ASSOCIATED(cs%e_D)) deallocate(cs%e_D)
1462  if (ASSOCIATED(cs%KE)) deallocate(cs%KE)
1463  if (ASSOCIATED(cs%dKE_dt)) deallocate(cs%dKE_dt)
1464  if (ASSOCIATED(cs%PE_to_KE)) deallocate(cs%PE_to_KE)
1465  if (ASSOCIATED(cs%KE_Coradv)) deallocate(cs%KE_Coradv)
1466  if (ASSOCIATED(cs%KE_adv)) deallocate(cs%KE_adv)
1467  if (ASSOCIATED(cs%KE_visc)) deallocate(cs%KE_visc)
1468  if (ASSOCIATED(cs%KE_horvisc)) deallocate(cs%KE_horvisc)
1469  if (ASSOCIATED(cs%KE_dia)) deallocate(cs%KE_dia)
1470  if (ASSOCIATED(cs%dv_dt)) deallocate(cs%dv_dt)
1471  if (ASSOCIATED(cs%dh_dt)) deallocate(cs%dh_dt)
1472  if (ASSOCIATED(cs%du_dt)) deallocate(cs%du_dt)
1473  if (ASSOCIATED(cs%h_Rlay)) deallocate(cs%h_Rlay)
1474  if (ASSOCIATED(cs%uh_Rlay)) deallocate(cs%uh_Rlay)
1475  if (ASSOCIATED(cs%vh_Rlay)) deallocate(cs%vh_Rlay)
1476  if (ASSOCIATED(cs%uhGM_Rlay)) deallocate(cs%uhGM_Rlay)
1477  if (ASSOCIATED(cs%vhGM_Rlay)) deallocate(cs%vhGM_Rlay)
1478  if (ASSOCIATED(cs%diag_tmp3d)) deallocate(cs%diag_tmp3d)
1479 
1480  if (ASSOCIATED(adp%gradKEu)) deallocate(adp%gradKEu)
1481  if (ASSOCIATED(adp%gradKEu)) deallocate(adp%gradKEu)
1482  if (ASSOCIATED(adp%du_dt_visc)) deallocate(adp%du_dt_visc)
1483  if (ASSOCIATED(adp%dv_dt_visc)) deallocate(adp%dv_dt_visc)
1484  if (ASSOCIATED(adp%du_dt_dia)) deallocate(adp%du_dt_dia)
1485  if (ASSOCIATED(adp%dv_dt_dia)) deallocate(adp%dv_dt_dia)
1486  if (ASSOCIATED(adp%du_other)) deallocate(adp%du_other)
1487  if (ASSOCIATED(adp%dv_other)) deallocate(adp%dv_other)
1488 
1489  do m=1,cs%num_time_deriv ; deallocate(cs%prev_val(m)%p) ; enddo
1490 
1491  deallocate(cs)
1492 

◆ mom_diagnostics_init()

subroutine, public mom_diagnostics::mom_diagnostics_init ( type(ocean_internal_state), intent(in)  MIS,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(diagnostics_cs), pointer  CS 
)
Parameters
[in]misFor "MOM Internal State" a set of pointers to the fields and accelerations that make up the ocean's internal physical state.
[in,out]adpStructure with pointers to momentum equation terms.
[in,out]cdpStructure with pointers to continuity equation terms.
[in]timeCurrent model time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagStructure to regulate diagnostic output.
csPointer set to point to control structure for this module.

Definition at line 1115 of file MOM_diagnostics.F90.

References mom_error_handler::mom_error(), register_time_deriv(), set_dependent_diagnostics(), and mom_wave_speed::wave_speed_init().

1115  type(ocean_internal_state), intent(in) :: mis !< For "MOM Internal State" a set of pointers to
1116  !! the fields and accelerations that make up the
1117  !! ocean's internal physical state.
1118  type(accel_diag_ptrs), intent(inout) :: adp !< Structure with pointers to momentum equation
1119  !! terms.
1120  type(cont_diag_ptrs), intent(inout) :: cdp !< Structure with pointers to continuity
1121  !! equation terms.
1122  type(time_type), intent(in) :: time !< Current model time.
1123  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1124  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1125  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1126  !! parameters.
1127  type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.
1128  type(diagnostics_cs), pointer :: cs !< Pointer set to point to control structure
1129  !! for this module.
1130 
1131 ! Arguments
1132 ! (in) MIS - For "MOM Internal State" a set of pointers to the fields and
1133 ! accelerations that make up the ocean's internal physical
1134 ! state.
1135 ! (inout) ADp - structure with pointers to momentum equation terms
1136 ! (inout) CDp - structure with pointers to continuity equation terms
1137 ! (in) Time - current model time
1138 ! (in) G - ocean grid structure
1139 ! (in) GV - The ocean's vertical grid structure.
1140 ! (in) param_file - structure indicating the open file to parse for
1141 ! model parameter values
1142 ! (in) diag - structure to regulate diagnostic output
1143 ! (in/out) CS - pointer set to point to control structure for this module
1144 
1145 ! This include declares and sets the variable "version".
1146 #include "version_variable.h"
1147 
1148  character(len=40) :: mdl = "MOM_diagnostics" ! This module's name.
1149  real :: omega, f2_min
1150  character(len=48) :: thickness_units, flux_units
1151  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nkml, nkbl
1152  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
1153 
1154  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1155  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1156  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1157  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1158 
1159  if (associated(cs)) then
1160  call mom_error(warning, "MOM_diagnostics_init called with an associated "// &
1161  "control structure.")
1162  return
1163  endif
1164  allocate(cs)
1165 
1166  cs%diag => diag
1167 
1168  ! Read all relevant parameters and write them to the model log.
1169  call log_version(param_file, mdl, version)
1170  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_COLUMN_FRACTION", cs%mono_N2_column_fraction, &
1171  "The lower fraction of water column over which N2 is limited as monotonic\n"// &
1172  "for the purposes of calculating the equivalent barotropic wave speed.", &
1173  units='nondim', default=0.)
1174  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", cs%mono_N2_depth, &
1175  "The depth below which N2 is limited as monotonic for the\n"// &
1176  "purposes of calculating the equivalent barotropic wave speed.", &
1177  units='m', default=-1.)
1178 
1179  if (gv%Boussinesq) then
1180  thickness_units = "meter" ; flux_units = "meter3 second-1"
1181  else
1182  thickness_units = "kilogram meter-2" ; flux_units = "kilogram second-1"
1183  endif
1184 
1185  cs%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', diag%axesZL, time, &
1186  'Layer Average Ocean Temperature', 'Celsius')
1187 
1188  cs%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', diag%axesZL, time, &
1189  'Layer Average Ocean Salinity', 'ppt')
1190 
1191  cs%id_masscello = register_diag_field('ocean_model', 'masscello', diag%axesTL,&
1192  time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2', &
1193  standard_name='sea_water_mass_per_unit_area', v_extensive=.true.)
1194 
1195  cs%id_masso = register_scalar_field('ocean_model', 'masso', time, &
1196  diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass')
1197 
1198  cs%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, time, &
1199  long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.)
1200 
1201  if (((cs%id_masscello>0) .or. (cs%id_masso>0) .or. (cs%id_thkcello>0.and..not.gv%Boussinesq)) &
1202  .and. .not.ASSOCIATED(cs%diag_tmp3d)) then
1203  call safe_alloc_ptr(cs%diag_tmp3d,isd,ied,jsd,jed,nz)
1204  endif
1205 
1206  cs%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga', &
1207  time, diag, 'Global Mean Ocean Potential Temperature', 'Celsius',&
1208  standard_name='sea_water_potential_temperature')
1209 
1210  cs%id_soga = register_scalar_field('ocean_model', 'soga', &
1211  time, diag, 'Global Mean Ocean Salinity', 'ppt', &
1212  standard_name='sea_water_salinity')
1213 
1214  cs%id_tosga = register_scalar_field('ocean_model', 'sst_global', time, diag,&
1215  long_name='Global Area Average Sea Surface Temperature', &
1216  units='degC', standard_name='sea_surface_temperature', &
1217  cmor_field_name='tosga', cmor_standard_name='sea_surface_temperature', &
1218  cmor_units='degC', cmor_long_name='Sea Surface Temperature')
1219 
1220  cs%id_sosga = register_scalar_field('ocean_model', 'sss_global', time, diag,&
1221  long_name='Global Area Average Sea Surface Salinity', &
1222  units='ppt', standard_name='sea_surface_salinity', &
1223  cmor_field_name='sosga', cmor_standard_name='sea_surface_salinity', &
1224  cmor_units='ppt', cmor_long_name='Sea Surface Salinity')
1225 
1226  cs%id_e = register_diag_field('ocean_model', 'e', diag%axesTi, time, &
1227  'Interface Height Relative to Mean Sea Level', 'meter')
1228  if (cs%id_e>0) call safe_alloc_ptr(cs%e,isd,ied,jsd,jed,nz+1)
1229 
1230  cs%id_e_D = register_diag_field('ocean_model', 'e_D', diag%axesTi, time, &
1231  'Interface Height above the Seafloor', 'meter')
1232  if (cs%id_e_D>0) call safe_alloc_ptr(cs%e_D,isd,ied,jsd,jed,nz+1)
1233 
1234  cs%id_Rml = register_diag_field('ocean_model', 'Rml', diag%axesTL, time, &
1235  'Mixed Layer Coordinate Potential Density', 'kg meter-3')
1236 
1237  cs%id_Rcv = register_diag_field('ocean_model', 'Rho_cv', diag%axesTL, time, &
1238  'Coordinate Potential Density', 'kg meter-3')
1239 
1240  cs%id_rhopot0 = register_diag_field('ocean_model', 'rhopot0', diag%axesTL, time, &
1241  'Potential density referenced to surface', 'kg meter-3')
1242  cs%id_rhopot2 = register_diag_field('ocean_model', 'rhopot2', diag%axesTL, time, &
1243  'Potential density referenced to 2000 dbar', 'kg meter-3')
1244  cs%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, time, &
1245  'In situ density', 'kg meter-3')
1246 
1247  cs%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, time, &
1248  'Zonal Acceleration', 'meter second-2')
1249  if ((cs%id_du_dt>0) .and. .not.ASSOCIATED(cs%du_dt)) then
1250  call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
1251  call register_time_deriv(mis%u, cs%du_dt, cs)
1252  endif
1253 
1254  cs%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, time, &
1255  'Meridional Acceleration', 'meter second-2')
1256  if ((cs%id_dv_dt>0) .and. .not.ASSOCIATED(cs%dv_dt)) then
1257  call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
1258  call register_time_deriv(mis%v, cs%dv_dt, cs)
1259  endif
1260 
1261  cs%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, time, &
1262  'Thickness tendency', trim(thickness_units)//" second-1")
1263  if ((cs%id_dh_dt>0) .and. .not.ASSOCIATED(cs%dh_dt)) then
1264  call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
1265  call register_time_deriv(mis%h, cs%dh_dt, cs)
1266  endif
1267 
1268  ! layer thickness variables
1269  !if (GV%nk_rho_varies > 0) then
1270  cs%id_h_Rlay = register_diag_field('ocean_model', 'h_rho', diag%axesTL, time, &
1271  'Layer thicknesses in pure potential density coordinates', thickness_units)
1272  if (cs%id_h_Rlay>0) call safe_alloc_ptr(cs%h_Rlay,isd,ied,jsd,jed,nz)
1273 
1274  cs%id_uh_Rlay = register_diag_field('ocean_model', 'uh_rho', diag%axesCuL, time, &
1275  'Zonal volume transport in pure potential density coordinates', flux_units)
1276  if (cs%id_uh_Rlay>0) call safe_alloc_ptr(cs%uh_Rlay,isdb,iedb,jsd,jed,nz)
1277 
1278  cs%id_vh_Rlay = register_diag_field('ocean_model', 'vh_rho', diag%axesCvL, time, &
1279  'Meridional volume transport in pure potential density coordinates', flux_units)
1280  if (cs%id_vh_Rlay>0) call safe_alloc_ptr(cs%vh_Rlay,isd,ied,jsdb,jedb,nz)
1281 
1282  cs%id_uhGM_Rlay = register_diag_field('ocean_model', 'uhGM_rho', diag%axesCuL, time, &
1283  'Zonal volume transport due to interface height diffusion in pure potential &
1284  &density coordinates', flux_units)
1285  if (cs%id_uhGM_Rlay>0) call safe_alloc_ptr(cs%uhGM_Rlay,isdb,iedb,jsd,jed,nz)
1286 
1287  cs%id_vhGM_Rlay = register_diag_field('ocean_model', 'vhGM_rho', diag%axesCvL, time, &
1288  'Meridional volume transport due to interface height diffusion in pure &
1289  &potential density coordinates', flux_units)
1290  if (cs%id_vhGM_Rlay>0) call safe_alloc_ptr(cs%vhGM_Rlay,isd,ied,jsdb,jedb,nz)
1291  !endif
1292 
1293 
1294  ! terms in the kinetic energy budget
1295  cs%id_KE = register_diag_field('ocean_model', 'KE', diag%axesTL, time, &
1296  'Layer kinetic energy per unit mass', 'meter2 second-2')
1297  if (cs%id_KE>0) call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
1298 
1299  cs%id_dKEdt = register_diag_field('ocean_model', 'dKE_dt', diag%axesTL, time, &
1300  'Kinetic Energy Tendency of Layer', 'meter3 second-3')
1301  if (cs%id_dKEdt>0) call safe_alloc_ptr(cs%dKE_dt,isd,ied,jsd,jed,nz)
1302 
1303  cs%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, time, &
1304  'Potential to Kinetic Energy Conversion of Layer', 'meter3 second-3')
1305  if (cs%id_PE_to_KE>0) call safe_alloc_ptr(cs%PE_to_KE,isd,ied,jsd,jed,nz)
1306 
1307  cs%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, time, &
1308  'Kinetic Energy Source from Coriolis and Advection', 'meter3 second-3')
1309  if (cs%id_KE_Coradv>0) call safe_alloc_ptr(cs%KE_Coradv,isd,ied,jsd,jed,nz)
1310 
1311  cs%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, time, &
1312  'Kinetic Energy Source from Advection', 'meter3 second-3')
1313  if (cs%id_KE_adv>0) call safe_alloc_ptr(cs%KE_adv,isd,ied,jsd,jed,nz)
1314 
1315  cs%id_KE_visc = register_diag_field('ocean_model', 'KE_visc', diag%axesTL, time, &
1316  'Kinetic Energy Source from Vertical Viscosity and Stresses', 'meter3 second-3')
1317  if (cs%id_KE_visc>0) call safe_alloc_ptr(cs%KE_visc,isd,ied,jsd,jed,nz)
1318 
1319  cs%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, time, &
1320  'Kinetic Energy Source from Horizontal Viscosity', 'meter3 second-3')
1321  if (cs%id_KE_horvisc>0) call safe_alloc_ptr(cs%KE_horvisc,isd,ied,jsd,jed,nz)
1322 
1323  cs%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, time, &
1324  'Kinetic Energy Source from Diapycnal Diffusion', 'meter3 second-3')
1325  if (cs%id_KE_dia>0) call safe_alloc_ptr(cs%KE_dia,isd,ied,jsd,jed,nz)
1326 
1327 
1328  ! gravity wave CFLs
1329  cs%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, time, &
1330  'First baroclinic gravity wave speed', 'meter second-1')
1331  cs%id_Rd1 = register_diag_field('ocean_model', 'Rd1', diag%axesT1, time, &
1332  'First baroclinic deformation radius', 'meter')
1333  cs%id_cfl_cg1 = register_diag_field('ocean_model', 'CFL_cg1', diag%axesT1, time, &
1334  'CFL of first baroclinic gravity wave = dt*cg1*(1/dx+1/dy)', 'nondim')
1335  cs%id_cfl_cg1_x = register_diag_field('ocean_model', 'CFL_cg1_x', diag%axesT1, time, &
1336  'i-component of CFL of first baroclinic gravity wave = dt*cg1*/dx', 'nondim')
1337  cs%id_cfl_cg1_y = register_diag_field('ocean_model', 'CFL_cg1_y', diag%axesT1, time, &
1338  'j-component of CFL of first baroclinic gravity wave = dt*cg1*/dy', 'nondim')
1339  cs%id_cg_ebt = register_diag_field('ocean_model', 'cg_ebt', diag%axesT1, time, &
1340  'Equivalent barotropic gravity wave speed', 'meter second-1')
1341  cs%id_Rd_ebt = register_diag_field('ocean_model', 'Rd_ebt', diag%axesT1, time, &
1342  'Equivalent barotropic deformation radius', 'meter')
1343  cs%id_p_ebt = register_diag_field('ocean_model', 'p_ebt', diag%axesTL, time, &
1344  'Equivalent barotropic modal strcuture', 'nondim')
1345 
1346  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
1347  (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0) .or. &
1348  (cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
1349  call wave_speed_init(cs%wave_speed_CSp)
1350  call safe_alloc_ptr(cs%cg1,isd,ied,jsd,jed)
1351  if (cs%id_Rd1>0) call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
1352  if (cs%id_Rd_ebt>0) call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
1353  if (cs%id_cfl_cg1>0) call safe_alloc_ptr(cs%cfl_cg1,isd,ied,jsd,jed)
1354  if (cs%id_cfl_cg1_x>0) call safe_alloc_ptr(cs%cfl_cg1_x,isd,ied,jsd,jed)
1355  if (cs%id_cfl_cg1_y>0) call safe_alloc_ptr(cs%cfl_cg1_y,isd,ied,jsd,jed)
1356  if (cs%id_p_ebt>0) call safe_alloc_ptr(cs%p_ebt,isd,ied,jsd,jed,nz)
1357  endif
1358 
1359  cs%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, time, &
1360  'The column mass for calculating mass-weighted average properties', 'kg m-2')
1361 
1362  cs%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, time, &
1363  'Density weighted column integrated potential temperature', 'degC kg m-2', &
1364  cmor_field_name='opottempmint', &
1365  cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',&
1366  cmor_units='degC kg m-2', cmor_standard_name='Depth integrated density times potential temperature')
1367 
1368  cs%id_salt_int = register_diag_field('ocean_model', 'salt_int', diag%axesT1, time, &
1369  'Density weighted column integrated salinity', 'ppt kg m-2', &
1370  cmor_field_name='somint', &
1371  cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity',&
1372  cmor_units='ppt kg m-2', cmor_standard_name='Depth integrated density times salinity')
1373 
1374  cs%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, time, &
1375  'The column integrated in situ density', 'kg m-2')
1376 
1377  cs%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, time, &
1378  'The height of the water column', 'm')
1379  cs%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, time, &
1380  long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
1381  units='Pa')
1382 
1383  call set_dependent_diagnostics(mis, adp, cdp, g, cs)
1384 
Here is the call graph for this function:

◆ register_time_deriv()

subroutine, public mom_diagnostics::register_time_deriv ( real, dimension(:,:,:), target  f_ptr,
real, dimension(:,:,:), target  deriv_ptr,
type(diagnostics_cs), pointer  CS 
)

This subroutine registers fields to calculate a diagnostic time derivative.

Parameters
f_ptrField whose derivative is taken.
deriv_ptrField in which the calculated time derivatives placed.
csControl structure returned by previous call to diagnostics_init.

Definition at line 1049 of file MOM_diagnostics.F90.

References mom_error_handler::mom_error().

Referenced by mom_diagnostics_init(), and set_dependent_diagnostics().

1049  real, dimension(:,:,:), target :: f_ptr !< Field whose derivative is taken.
1050  real, dimension(:,:,:), target :: deriv_ptr !< Field in which the calculated time derivatives
1051  !! placed.
1052  type(diagnostics_cs), pointer :: cs !< Control structure returned by previous call to
1053  !! diagnostics_init.
1054 
1055 ! This subroutine registers fields to calculate a diagnostic time derivative.
1056 ! Arguments:
1057 ! (target) f_ptr - field whose derivative is taken
1058 ! (in) deriv_ptr - field in which the calculated time derivatives placed
1059 ! (in) num_lay - number of layers in this field
1060 ! (in) CS - control structure returned by previous call to diagnostics_init
1061 
1062  integer :: m
1063 
1064  if (.not.associated(cs)) call mom_error(fatal, &
1065  "register_time_deriv: Module must be initialized before it is used.")
1066 
1067  if (cs%num_time_deriv >= max_fields_) then
1068  call mom_error(warning,"MOM_diagnostics: Attempted to register more than " // &
1069  "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
1070  return
1071  endif
1072 
1073  m = cs%num_time_deriv+1 ; cs%num_time_deriv = m
1074 
1075  cs%nlay(m) = size(f_ptr(:,:,:),3)
1076  cs%deriv(m)%p => deriv_ptr
1077  allocate(cs%prev_val(m)%p(size(f_ptr(:,:,:),1), size(f_ptr(:,:,:),2), cs%nlay(m)) )
1078 
1079  cs%var_ptr(m)%p => f_ptr
1080  cs%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
1081 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_dependent_diagnostics()

subroutine mom_diagnostics::set_dependent_diagnostics ( type(ocean_internal_state), intent(in)  MIS,
type(accel_diag_ptrs), intent(inout)  ADp,
type(cont_diag_ptrs), intent(inout)  CDp,
type(ocean_grid_type), intent(in)  G,
type(diagnostics_cs), pointer  CS 
)
private

This subroutine sets up diagnostics upon which other diagnostics depend.

Parameters
[in]misFor "MOM Internal State" a set of pointers to the fields and accelerations making up ocean internal physical state.
[in,out]adpStructure pointing to accelerations in momentum equation.
[in,out]cdpStructure pointing to terms in continuity equation.
[in]gThe ocean's grid structure.
csPointer to the control structure for this module.

Definition at line 1389 of file MOM_diagnostics.F90.

References register_time_deriv().

Referenced by mom_diagnostics_init().

1389  type(ocean_internal_state), intent(in) :: mis !< For "MOM Internal State" a set of pointers to
1390  !! the fields and accelerations making up ocean
1391  !! internal physical state.
1392  type(accel_diag_ptrs), intent(inout) :: adp !< Structure pointing to accelerations in
1393  !! momentum equation.
1394  type(cont_diag_ptrs), intent(inout) :: cdp !< Structure pointing to terms in continuity
1395  !! equation.
1396  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1397  type(diagnostics_cs), pointer :: cs !< Pointer to the control structure for this
1398  !! module.
1399 
1400 ! This subroutine sets up diagnostics upon which other diagnostics depend.
1401 ! Arguments:
1402 ! (in) MIS - For "MOM Internal State" a set of pointers to the fields and
1403 ! accelerations making up ocean internal physical state.
1404 ! (inout) ADp - structure pointing to accelerations in momentum equation
1405 ! (inout) CDp - structure pointing to terms in continuity equation
1406 ! (in) G - ocean grid structure
1407 ! (in) CS - pointer to the control structure for this module
1408 
1409  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
1410  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1411  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1412 
1413  if (ASSOCIATED(cs%dKE_dt) .or. ASSOCIATED(cs%PE_to_KE) .or. &
1414  ASSOCIATED(cs%KE_CorAdv) .or. ASSOCIATED(cs%KE_adv) .or. &
1415  ASSOCIATED(cs%KE_visc) .or. ASSOCIATED(cs%KE_horvisc) .or. &
1416  ASSOCIATED(cs%KE_dia)) &
1417  call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
1418 
1419  if (ASSOCIATED(cs%dKE_dt)) then
1420  if (.not.ASSOCIATED(cs%du_dt)) then
1421  call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
1422  call register_time_deriv(mis%u, cs%du_dt, cs)
1423  endif
1424  if (.not.ASSOCIATED(cs%dv_dt)) then
1425  call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
1426  call register_time_deriv(mis%v, cs%dv_dt, cs)
1427  endif
1428  if (.not.ASSOCIATED(cs%dh_dt)) then
1429  call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
1430  call register_time_deriv(mis%h, cs%dh_dt, cs)
1431  endif
1432  endif
1433 
1434  if (ASSOCIATED(cs%KE_adv)) then
1435  call safe_alloc_ptr(adp%gradKEu,isdb,iedb,jsd,jed,nz)
1436  call safe_alloc_ptr(adp%gradKEv,isd,ied,jsdb,jedb,nz)
1437  endif
1438 
1439  if (ASSOCIATED(cs%KE_visc)) then
1440  call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
1441  call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
1442  endif
1443 
1444  if (ASSOCIATED(cs%KE_dia)) then
1445  call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
1446  call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
1447  endif
1448 
1449  if (ASSOCIATED(cs%uhGM_Rlay)) call safe_alloc_ptr(cdp%uhGM,isdb,iedb,jsd,jed,nz)
1450  if (ASSOCIATED(cs%vhGM_Rlay)) call safe_alloc_ptr(cdp%vhGM,isd,ied,jsdb,jedb,nz)
1451 
Here is the call graph for this function:
Here is the caller graph for this function: