mom_diagnostics Module Reference

Detailed Description

Calculates any requested diagnostic quantities that are not calculated in the various subroutines. Diagnostic quantities are requested by allocating them memory.

Data Types
type	diagnostics_cs
	The control structure for the MOM_diagnostics module. More...
type	surface_diag_ids
	A structure with diagnostic IDs of the surface and integrated variables. More...
type	transport_diag_ids
	A structure with diagnostic IDs of mass transport related diagnostics. More...

Functions/Subroutines
subroutine, public	calculate_diagnostic_fields (u, v, h, uh, vh, tv, ADp, CDp, p_surf, dt, diag_pre_sync, G, GV, US, 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 the 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, p_surf, G, GV, US, CS)
	This 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, GV, US, CS)
	This subroutine calculates terms in the mechanical energy budget. More...
subroutine, public	register_time_deriv (lb, 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	post_surface_dyn_diags (IDs, G, diag, sfc_state, ssh)
	This routine posts diagnostics of various dynamic ocean surface quantities, including velocities, speed and sea surface height, at the time the ocean state is reported back to the caller. More...
subroutine, public	post_surface_thermo_diags (IDs, G, GV, US, diag, dt_int, sfc_state, tv, ssh, ssh_ibc)
	This routine posts diagnostics of various ocean surface and integrated quantities at the time the ocean state is reported back to the caller. More...
subroutine, public	post_transport_diagnostics (G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dyn, diag, dt_trans, Reg)
	This routine posts diagnostics of the transports, including the subgridscale contributions. More...
subroutine, public	mom_diagnostics_init (MIS, ADp, CDp, Time, G, GV, US, param_file, diag, CS, tv)
	This subroutine registers various diagnostics and allocates space for fields that other diagnostis depend upon. More...
subroutine, public	register_surface_diags (Time, G, US, IDs, diag, tv)
	Register diagnostics of the surface state and integrated quantities. More...
subroutine, public	register_transport_diags (Time, G, GV, US, IDs, diag)
	Register certain diagnostics related to transports. More...
subroutine, public	write_static_fields (G, GV, US, tv, diag)
	Offers the static fields in the ocean grid type for output via the diag_manager. More...
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)
	Deallocate memory associated with the diagnostics module. More...

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]	dt	The time interval over which differences occur [T ~> s].
[in,out]	g	The ocean's grid structure.
[in,out]	cs	Control structure returned by previous call to diagnostics_init.

Definition at line 1124 of file MOM_diagnostics.F90.

  real,                  intent(in)    :: dt   !< The time interval over which differences occur [T ~> s].
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
  type(diagnostics_CS),  intent(inout) :: CS   !< Control structure returned by previous call to
                                               !! diagnostics_init.
 
! This subroutine calculates all registered time derivatives.
  real :: Idt  ! The inverse timestep [T-1 ~> s-1]
  integer :: i, j, k, m
 
  if (dt > 0.0) then ; idt = 1.0/dt
  else ; return ; endif
 
  ! Because the field is unknown, its grid index bounds are also unknown.
  ! Additionally, two of the fields (dudt, dvdt) require calculation of spatial
  ! derivatives when computing d(KE)/dt.  This raises issues in non-symmetric
  ! mode, where the symmetric boundaries (west, south) may not be updated.
 
  ! For this reason, we explicitly loop from isc-1:iec and jsc-1:jec, in order
  ! to force boundary value updates, even though it may not be strictly valid
  ! for all fields.  Note this assumes a halo, and that it has been updated.
 
  do m=1,cs%num_time_deriv
    do k=1,cs%nlay(m) ; do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec
      cs%deriv(m)%p(i,j,k) = (cs%var_ptr(m)%p(i,j,k) - cs%prev_val(m)%p(i,j,k)) * idt
      cs%prev_val(m)%p(i,j,k) = cs%var_ptr(m)%p(i,j,k)
    enddo ; enddo ; enddo
  enddo
 

Referenced by calculate_diagnostic_fields().

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,
		real, dimension(:,:), pointer	p_surf,
		real, intent(in)	dt,
		type(diag_grid_storage), intent(in)	diag_pre_sync,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		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]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	u	The zonal velocity [L T-1 ~> m s-1].
[in]	v	The meridional velocity [L T-1 ~> m s-1].
[in]	h	Layer thicknesses [H ~> m or kg m-2].
[in]	uh	Transport through zonal faces = u*h*dy,
[in]	vh	Transport through meridional faces = v*h*dx,
[in]	tv	A structure pointing to various thermodynamic variables.
[in]	adp	structure with pointers to accelerations in momentum equation.
[in]	cdp	structure with pointers to terms in continuity equation.
	p_surf	A pointer to the surface pressure [Pa]. If p_surf is not associated, it is the same as setting the surface pressure to 0.
[in]	dt	The time difference since the last call to this subroutine [T ~> s].
[in]	diag_pre_sync	Target grids from previous timestep
[in,out]	cs	Control structure returned by a previous call to diagnostics_init.
[in]	eta_bt	An optional barotropic

Definition at line 188 of file MOM_diagnostics.F90.

  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v    !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: uh   !< Transport through zonal faces = u*h*dy,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: vh   !< Transport through meridional faces = v*h*dx,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various
                                                 !! thermodynamic variables.
  type(accel_diag_ptrs),   intent(in)    :: ADp  !< structure with pointers to
                                                 !! accelerations in momentum equation.
  type(cont_diag_ptrs),    intent(in)    :: CDp  !< structure with pointers to
                                                 !! terms in continuity equation.
  real, dimension(:,:),    pointer       :: p_surf !< A pointer to the surface pressure [Pa].
                                                 !! If p_surf is not associated, it is the same
                                                 !! as setting the surface pressure to 0.
  real,                    intent(in)    :: dt   !< The time difference since the last
                                                 !! call to this subroutine [T ~> s].
  type(diag_grid_storage), intent(in)    :: diag_pre_sync !< Target grids from previous timestep
  type(diagnostics_CS),    intent(inout) :: CS   !< Control structure returned by a
                                                 !! previous call to diagnostics_init.
  real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(in)   :: eta_bt !< An 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 [H ~> m or kg m-2].
 
  ! Local variables
  integer i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
 
  real :: Rcv(SZI_(G),SZJ_(G),SZK_(G))  ! Coordinate variable potential density [R ~> kg m-3].
  real :: work_3d(SZI_(G),SZJ_(G),SZK_(G)) ! A 3-d temporary work array.
  real :: work_2d(SZI_(G),SZJ_(G))         ! A 2-d temporary work array.
 
  ! tmp array for surface properties
  real :: surface_field(SZI_(G),SZJ_(G))
  real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS
  real :: wt, wt_p
 
  real :: f2_h     ! Squared Coriolis parameter at to h-points [T-2 ~> s-2]
  real :: mag_beta ! Magnitude of the gradient of f [T-1 L-1 ~> s-1 m-1]
  real :: absurdly_small_freq2 ! Srequency squared used to avoid division by 0 [T-2 ~> s-2]
 
  integer :: k_list
 
  real, dimension(SZK_(G)) :: temp_layer_ave, salt_layer_ave
  real :: thetaoga, soga, masso, tosga, sosga
 
  is  = g%isc  ; ie   = g%iec  ; js  = g%jsc  ; je  = g%jec
  isq = g%IscB ; ieq  = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  nz  = g%ke   ; nkmb = gv%nk_rho_varies
 
  ! This value is roughly (pi / (the age of the universe) )^2.
  absurdly_small_freq2 = 1e-34*us%T_to_s**2
 
  if (loc(cs)==0) call mom_error(fatal, &
         "calculate_diagnostic_fields: Module must be initialized before used.")
 
  call calculate_derivs(dt, g, cs)
 
  if (dt > 0.0) then
    call diag_save_grids(cs%diag)
    call diag_copy_storage_to_diag(cs%diag, diag_pre_sync)
 
    if (cs%id_h_pre_sync > 0) &
        call post_data(cs%id_h_pre_sync, diag_pre_sync%h_state, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_du_dt>0) call post_data(cs%id_du_dt, cs%du_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_dv_dt>0) call post_data(cs%id_dv_dt, cs%dv_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_dh_dt>0) call post_data(cs%id_dh_dt, cs%dh_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    call diag_restore_grids(cs%diag)
 
    call calculate_energy_diagnostics(u, v, h, uh, vh, adp, cdp, g, gv, us, cs)
  endif
 
  ! smg: is the following robust to ALE? It seems a bit opaque.
  ! If the model is NOT in isopycnal mode then nkmb=0. But we need all the
  ! following diagnostics to treat all layers as variable density, so we set
  ! nkmb = nz, on the expectation that loops nkmb+1,nz will not iterate.
  ! This behavior is ANSI F77 but some compiler options can force at least
  ! one iteration that would break the following one-line workaround!
  if (nkmb==0 .and. nz > 1) nkmb = nz
 
  if (cs%id_u > 0) call post_data(cs%id_u, u, cs%diag)
 
  if (cs%id_v > 0) call post_data(cs%id_v, v, cs%diag)
 
  if (cs%id_h > 0) call post_data(cs%id_h, h, cs%diag)
 
  if (associated(cs%e)) then
    call find_eta(h, tv, g, gv, us, cs%e, eta_bt)
    if (cs%id_e > 0) call post_data(cs%id_e, cs%e, cs%diag)
  endif
 
  if (associated(cs%e_D)) then
    if (associated(cs%e)) then
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        cs%e_D(i,j,k) = cs%e(i,j,k) + g%bathyT(i,j)
      enddo ; enddo ; enddo
    else
      call find_eta(h, tv, g, gv, us, cs%e_D, eta_bt)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        cs%e_D(i,j,k) = cs%e_D(i,j,k) + g%bathyT(i,j)
      enddo ; enddo ; enddo
    endif
 
    if (cs%id_e_D > 0) call post_data(cs%id_e_D, cs%e_D, cs%diag)
  endif
 
  ! mass per area of grid cell (for Bouss, use Rho0)
  if (cs%id_masscello > 0) then
    do k=1,nz; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_masscello, work_3d, cs%diag)
  endif
 
  ! mass of liquid ocean (for Bouss, use Rho0). The reproducing sum requires the use of MKS units.
  if (cs%id_masso > 0) then
    work_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_2d(i,j) = work_2d(i,j) + (gv%H_to_kg_m2*h(i,j,k)) * us%L_to_m**2*g%areaT(i,j)
    enddo ; enddo ; enddo
    masso = reproducing_sum(work_2d)
    call post_data(cs%id_masso, masso, cs%diag)
  endif
 
  ! diagnose thickness/volumes of grid cells [m]
  if (cs%id_thkcello>0 .or. cs%id_volcello>0) then
    if (gv%Boussinesq) then ! thkcello = h for Boussinesq
      if (cs%id_thkcello > 0) then ; if (gv%H_to_m == 1.0) then
        call post_data(cs%id_thkcello, h, cs%diag)
      else
        do k=1,nz; do j=js,je ; do i=is,ie
          work_3d(i,j,k) = gv%H_to_m*h(i,j,k)
        enddo ; enddo ; enddo
        call post_data(cs%id_thkcello, work_3d, cs%diag)
      endif ; endif
      if (cs%id_volcello > 0) then ! volcello = h*area for Boussinesq
        do k=1,nz; do j=js,je ; do i=is,ie
          work_3d(i,j,k) = ( gv%H_to_m*h(i,j,k) ) * us%L_to_m**2*g%areaT(i,j)
        enddo ; enddo ; enddo
        call post_data(cs%id_volcello, work_3d, cs%diag)
      endif
    else ! thkcello = dp/(rho*g) for non-Boussinesq
      do j=js,je
        if (associated(p_surf)) then ! Pressure loading at top of surface layer [Pa]
          do i=is,ie
            pressure_1d(i) = p_surf(i,j)
          enddo
        else
          do i=is,ie
            pressure_1d(i) = 0.0
          enddo
        endif
        do k=1,nz ! Integrate vertically downward for pressure
          do i=is,ie ! Pressure for EOS at the layer center [Pa]
            pressure_1d(i) = pressure_1d(i) + 0.5*gv%H_to_Pa*h(i,j,k)
          enddo
          ! Store in-situ density [kg m-3] in work_3d
          call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
                                 work_3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
          do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d
            work_3d(i,j,k) = (gv%H_to_kg_m2*h(i,j,k)) / work_3d(i,j,k)
          enddo
          do i=is,ie ! Pressure for EOS at the bottom interface [Pa]
            pressure_1d(i) = pressure_1d(i) + 0.5*gv%H_to_Pa*h(i,j,k)
          enddo
        enddo ! k
      enddo ! j
      if (cs%id_thkcello > 0) call post_data(cs%id_thkcello, work_3d, cs%diag)
      if (cs%id_volcello > 0) then
        do k=1,nz; do j=js,je ; do i=is,ie ! volcello = dp/(rho*g)*area for non-Boussinesq
          work_3d(i,j,k) = us%L_to_m**2*g%areaT(i,j) * work_3d(i,j,k)
        enddo ; enddo ; enddo
        call post_data(cs%id_volcello, work_3d, cs%diag)
      endif
    endif
  endif
 
  ! Calculate additional, potentially derived temperature diagnostics
  if (tv%T_is_conT) then
    ! Internal T&S variables are conservative temperature & absolute salinity,
    ! so they need to converted to potential temperature and practical salinity
    ! for some diagnostics using TEOS-10 function calls.
    if ((cs%id_Tpot > 0) .or. (cs%id_tob > 0)) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k))
      enddo ; enddo ; enddo
      if (cs%id_Tpot > 0) call post_data(cs%id_Tpot, work_3d, cs%diag)
      if (cs%id_tob > 0) call post_data(cs%id_tob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
    endif
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (cs%id_tob > 0) call post_data(cs%id_tob, tv%T(:,:,nz), cs%diag, mask=g%mask2dT)
  endif
 
  ! Calculate additional, potentially derived salinity diagnostics
  if (tv%S_is_absS) then
    ! Internal T&S variables are conservative temperature & absolute salinity,
    ! so they need to converted to potential temperature and practical salinity
    ! for some diagnostics using TEOS-10 function calls.
    if ((cs%id_Sprac > 0) .or. (cs%id_sob > 0)) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k))
      enddo ; enddo ; enddo
      if (cs%id_Sprac > 0) call post_data(cs%id_Sprac, work_3d, cs%diag)
      if (cs%id_sob > 0) call post_data(cs%id_sob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
    endif
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (cs%id_sob > 0) call post_data(cs%id_sob, tv%S(:,:,nz), cs%diag, mask=g%mask2dT)
  endif
 
  ! volume mean potential temperature
  if (cs%id_thetaoga>0) then
    thetaoga = global_volume_mean(tv%T, h, g, gv)
    call post_data(cs%id_thetaoga, thetaoga, cs%diag)
  endif
 
  ! area mean SST
  if (cs%id_tosga > 0) then
    do j=js,je ; do i=is,ie
       surface_field(i,j) = tv%T(i,j,1)
    enddo ; enddo
    tosga = global_area_mean(surface_field, g)
    call post_data(cs%id_tosga, tosga, cs%diag)
  endif
 
  ! volume mean salinity
  if (cs%id_soga>0) then
    soga = global_volume_mean(tv%S, h, g, gv)
    call post_data(cs%id_soga, soga, cs%diag)
  endif
 
  ! area mean SSS
  if (cs%id_sosga > 0) then
    do j=js,je ; do i=is,ie
       surface_field(i,j) = tv%S(i,j,1)
    enddo ; enddo
    sosga = global_area_mean(surface_field, g)
    call post_data(cs%id_sosga, sosga, cs%diag)
  endif
 
  ! layer mean potential temperature
  if (cs%id_temp_layer_ave>0) then
    temp_layer_ave = global_layer_mean(tv%T, h, g, gv)
    call post_data(cs%id_temp_layer_ave, temp_layer_ave, cs%diag)
  endif
 
  ! layer mean salinity
  if (cs%id_salt_layer_ave>0) then
    salt_layer_ave = global_layer_mean(tv%S, h, g, gv)
    call post_data(cs%id_salt_layer_ave, salt_layer_ave, cs%diag)
  endif
 
  call calculate_vertical_integrals(h, tv, p_surf, g, gv, us, cs)
 
  if ((cs%id_Rml > 0) .or. (cs%id_Rcv > 0) .or. associated(cs%h_Rlay) .or. &
      associated(cs%uh_Rlay) .or. associated(cs%vh_Rlay) .or. &
      associated(cs%uhGM_Rlay) .or. associated(cs%vhGM_Rlay)) then
 
    if (associated(tv%eqn_of_state)) then
      pressure_1d(:) = tv%P_Ref
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=js-1,je+1
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, &
                               rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state, scale=us%kg_m3_to_R)
      enddo ; enddo
    else ! Rcv should not be used much in this case, so fill in sensible values.
      do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
        rcv(i,j,k) = gv%Rlay(k)
      enddo ; enddo ; enddo
    endif
    if (cs%id_Rml > 0) call post_data(cs%id_Rml, rcv, cs%diag)
    if (cs%id_Rcv > 0) call post_data(cs%id_Rcv, rcv, cs%diag)
 
    if (associated(cs%h_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,CS,Rcv,h,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=is,ie
          cs%h_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%h_Rlay(i,j,k) = h(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, rcv(i,j,k), k_list, nz, wt, wt_p)
          cs%h_Rlay(i,j,k_list)   = cs%h_Rlay(i,j,k_list)   + h(i,j,k)*wt
          cs%h_Rlay(i,j,k_list+1) = cs%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
        enddo ; enddo
      enddo
 
      if (cs%id_h_Rlay > 0) call post_data(cs%id_h_Rlay, cs%h_Rlay, cs%diag)
    endif
 
    if (associated(cs%uh_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CS,GV,uh) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=isq,ieq
          cs%uh_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=isq,ieq
          cs%uh_Rlay(i,j,k) = uh(i,j,k)
        enddo ; enddo
        k_list = nz/2
        do k=1,nkmb ; do i=isq,ieq
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
          cs%uh_Rlay(i,j,k_list)   = cs%uh_Rlay(i,j,k_list)   + uh(i,j,k)*wt
          cs%uh_Rlay(i,j,k_list+1) = cs%uh_Rlay(i,j,k_list+1) + uh(i,j,k)*wt_p
        enddo ; enddo
      enddo
 
      if (cs%id_uh_Rlay > 0) call post_data(cs%id_uh_Rlay, cs%uh_Rlay, cs%diag)
    endif
 
    if (associated(cs%vh_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none)  shared(Jsq,Jeq,is,ie,nz,nkmb,Rcv,CS,GV,vh) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=jsq,jeq
        do k=1,nkmb ; do i=is,ie
          cs%vh_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%vh_Rlay(i,j,k) = vh(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
          cs%vh_Rlay(i,j,k_list)   = cs%vh_Rlay(i,j,k_list)   + vh(i,j,k)*wt
          cs%vh_Rlay(i,j,k_list+1) = cs%vh_Rlay(i,j,k_list+1) + vh(i,j,k)*wt_p
        enddo ; enddo
      enddo
 
      if (cs%id_vh_Rlay > 0) call post_data(cs%id_vh_Rlay, cs%vh_Rlay, cs%diag)
    endif
 
    if (associated(cs%uhGM_Rlay) .and. associated(cdp%uhGM)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CDP,CS,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=isq,ieq
          cs%uhGM_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=isq,ieq
          cs%uhGM_Rlay(i,j,k) = cdp%uhGM(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=isq,ieq
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
          cs%uhGM_Rlay(i,j,k_list)   = cs%uhGM_Rlay(i,j,k_list)   + cdp%uhGM(i,j,k)*wt
          cs%uhGM_Rlay(i,j,k_list+1) = cs%uhGM_Rlay(i,j,k_list+1) + cdp%uhGM(i,j,k)*wt_p
        enddo ; enddo
      enddo
 
      if (cs%id_uhGM_Rlay > 0) call post_data(cs%id_uhGM_Rlay, cs%uhGM_Rlay, cs%diag)
    endif
 
    if (associated(cs%vhGM_Rlay) .and. associated(cdp%vhGM)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(is,ie,Jsq,Jeq,nz,nkmb,CS,CDp,Rcv,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=jsq,jeq
        do k=1,nkmb ; do i=is,ie
          cs%vhGM_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%vhGM_Rlay(i,j,k) = cdp%vhGM(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
          cs%vhGM_Rlay(i,j,k_list)   = cs%vhGM_Rlay(i,j,k_list)   + cdp%vhGM(i,j,k)*wt
          cs%vhGM_Rlay(i,j,k_list+1) = cs%vhGM_Rlay(i,j,k_list+1) + cdp%vhGM(i,j,k)*wt_p
        enddo ; enddo
      enddo
 
      if (cs%id_vhGM_Rlay > 0) call post_data(cs%id_vhGM_Rlay, cs%vhGM_Rlay, cs%diag)
    endif
  endif
 
  if (associated(tv%eqn_of_state)) then
    if (cs%id_rhopot0 > 0) then
      pressure_1d(:) = 0.
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                               rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
      enddo ; enddo
      if (cs%id_rhopot0 > 0) call post_data(cs%id_rhopot0, rcv, cs%diag)
    endif
    if (cs%id_rhopot2 > 0) then
      pressure_1d(:) = 2.e7 ! 2000 dbars
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                               rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
      enddo ; enddo
      if (cs%id_rhopot2 > 0) call post_data(cs%id_rhopot2, rcv, cs%diag)
    endif
    if (cs%id_rhoinsitu > 0) then
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d,h,GV)
      do j=js,je
        pressure_1d(:) = 0. ! Start at p=0 Pa at surface
        do k=1,nz
          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure in middle of layer k
          call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                                 rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure at bottom of layer k
        enddo
      enddo
      if (cs%id_rhoinsitu > 0) call post_data(cs%id_rhoinsitu, rcv, cs%diag)
    endif
  endif
 
  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
      (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0)) then
    call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
    if (cs%id_cg1>0) call post_data(cs%id_cg1, cs%cg1, cs%diag)
    if (cs%id_Rd1>0) then
      !$OMP parallel do default(shared) private(f2_h,mag_beta)
      do j=js,je ; do i=is,ie
        ! Blend the equatorial deformation radius with the standard one.
        f2_h = absurdly_small_freq2 + 0.25 * &
            ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
             (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
        mag_beta = sqrt(0.5 * ( &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
        cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
 
      enddo ; enddo
      call post_data(cs%id_Rd1, cs%Rd1, cs%diag)
    endif
    if (cs%id_cfl_cg1>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1(i,j) = (dt*cs%cg1(i,j)) * (g%IdxT(i,j) + g%IdyT(i,j))
      enddo ; enddo
      call post_data(cs%id_cfl_cg1, cs%cfl_cg1, cs%diag)
    endif
    if (cs%id_cfl_cg1_x>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1_x(i,j) = (dt*cs%cg1(i,j)) * g%IdxT(i,j)
      enddo ; enddo
      call post_data(cs%id_cfl_cg1_x, cs%cfl_cg1_x, cs%diag)
    endif
    if (cs%id_cfl_cg1_y>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1_y(i,j) = (dt*cs%cg1(i,j)) * g%IdyT(i,j)
      enddo ; enddo
      call post_data(cs%id_cfl_cg1_y, cs%cfl_cg1_y, cs%diag)
    endif
  endif
  if ((cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
    if (cs%id_p_ebt>0) then
      call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
                      mono_n2_column_fraction=cs%mono_N2_column_fraction, &
                      mono_n2_depth=cs%mono_N2_depth, modal_structure=cs%p_ebt)
      call post_data(cs%id_p_ebt, cs%p_ebt, cs%diag)
    else
      call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
                      mono_n2_column_fraction=cs%mono_N2_column_fraction, &
                      mono_n2_depth=cs%mono_N2_depth)
    endif
    if (cs%id_cg_ebt>0) call post_data(cs%id_cg_ebt, cs%cg1, cs%diag)
    if (cs%id_Rd_ebt>0) then
      !$OMP parallel do default(shared) private(f2_h,mag_beta)
      do j=js,je ; do i=is,ie
        ! Blend the equatorial deformation radius with the standard one.
        f2_h = absurdly_small_freq2 + 0.25 * &
            ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
             (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
        mag_beta = sqrt(0.5 * ( &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
        cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
 
      enddo ; enddo
      call post_data(cs%id_Rd_ebt, cs%Rd1, cs%diag)
    endif
  endif
 

References calculate_derivs(), calculate_energy_diagnostics(), calculate_vertical_integrals(), mom_diag_mediator::diag_copy_storage_to_diag(), mom_diag_mediator::diag_restore_grids(), mom_diag_mediator::diag_save_grids(), find_weights(), mom_spatial_means::global_volume_mean(), mom_error_handler::mom_error(), and mom_wave_speed::wave_speed().

Here is the call graph for this function:

calculate_energy_diagnostics()

subroutine mom_diagnostics::calculate_energy_diagnostics ( 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(accel_diag_ptrs), intent(in)  ADp, type(cont_diag_ptrs), intent(in)  CDp, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(diagnostics_cs), intent(inout)  CS  )

private

This subroutine calculates terms in the mechanical energy budget.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	u	The zonal velocity [L T-1 ~> m s-1].
[in]	v	The meridional velocity [L T-1 ~> m s-1].
[in]	h	Layer thicknesses [H ~> m or kg m-2].
[in]	uh	Transport through zonal faces=u*h*dy,
[in]	vh	Transport through merid faces=v*h*dx,
[in]	adp	Structure pointing to accelerations in momentum equation.
[in]	cdp	Structure pointing to terms in continuity equations.
[in]	us	A dimensional unit scaling type
[in,out]	cs	Control structure returned by a previous call to diagnostics_init.

Definition at line 880 of file MOM_diagnostics.F90.

  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v    !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: uh   !< Transport through zonal faces=u*h*dy,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: vh   !< Transport through merid faces=v*h*dx,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  type(accel_diag_ptrs),   intent(in)    :: ADp  !< Structure pointing to accelerations in momentum equation.
  type(cont_diag_ptrs),    intent(in)    :: CDp  !< Structure pointing to terms in continuity equations.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(diagnostics_CS),    intent(inout) :: CS   !< Control structure returned by a previous call to
                                                 !! diagnostics_init.
 
! This subroutine calculates terms in the mechanical energy budget.
 
  ! Local variables
  real :: KE_u(SZIB_(G),SZJ_(G))
  real :: KE_v(SZI_(G),SZJB_(G))
  real :: KE_h(SZI_(G),SZJ_(G))
 
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
 
  do j=js-1,je ; do i=is-1,ie
    ke_u(i,j) = 0.0 ; ke_v(i,j) = 0.0
  enddo ; enddo
 
  if (associated(cs%KE)) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%KE(i,j,k) = ((u(i,j,k) * u(i,j,k) + u(i-1,j,k) * u(i-1,j,k)) &
          + (v(i,j,k) * v(i,j,k) + v(i,j-1,k) * v(i,j-1,k))) * 0.25
      ! DELETE THE FOLLOWING...  Make this 0 to test the momentum balance,
      ! or a huge number to test the continuity balance.
      ! CS%KE(i,j,k) *= 1e20
    enddo ; enddo ; enddo
    if (cs%id_KE > 0) call post_data(cs%id_KE, cs%KE, cs%diag)
  endif
 
  if (.not.g%symmetric) then
    if (associated(cs%dKE_dt) .OR. associated(cs%PE_to_KE) .OR. associated(cs%KE_CorAdv) .OR. &
       associated(cs%KE_adv) .OR. associated(cs%KE_visc)  .OR. associated(cs%KE_horvisc).OR. &
       associated(cs%KE_dia) ) then
        call create_group_pass(cs%pass_KE_uv, ke_u, ke_v, g%Domain, to_north+to_east)
    endif
  endif
 
  if (associated(cs%dKE_dt)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * cs%du_dt(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * cs%dv_dt(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = cs%KE(i,j,k) * cs%dh_dt(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%dKE_dt(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_dKEdt > 0) call post_data(cs%id_dKEdt, cs%dKE_dt, cs%diag)
  endif
 
  if (associated(cs%PE_to_KE)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%PFu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%PFv(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%PE_to_KE(i,j,k) = 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_PE_to_KE > 0) call post_data(cs%id_PE_to_KE, cs%PE_to_KE, cs%diag)
  endif
 
  if (associated(cs%KE_CorAdv)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%CAu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%CAv(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) &
            * (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_CorAdv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_Coradv > 0) call post_data(cs%id_KE_Coradv, cs%KE_Coradv, cs%diag)
  endif
 
  if (associated(cs%KE_adv)) then
    ! NOTE: All terms in KE_adv are multipled by -1, which can easily produce
    ! negative zeros and may signal a reproducibility issue over land.
    ! We resolve this by re-initializing and only evaluating over water points.
    ke_u(:,:) = 0. ; ke_v(:,:) = 0.
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        if (g%mask2dCu(i,j) /= 0.) &
          ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%gradKEu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        if (g%mask2dCv(i,j) /= 0.) &
          ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%gradKEv(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) &
            * (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_adv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_adv > 0) call post_data(cs%id_KE_adv, cs%KE_adv, cs%diag)
  endif
 
  if (associated(cs%KE_visc)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%du_dt_visc(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%dv_dt_visc(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_visc(i,j,k) = 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_visc > 0) call post_data(cs%id_KE_visc, cs%KE_visc, cs%diag)
  endif
 
  if (associated(cs%KE_horvisc)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%diffu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%diffv(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_horvisc(i,j,k) = 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_horvisc > 0) call post_data(cs%id_KE_horvisc, cs%KE_horvisc, cs%diag)
  endif
 
  if (associated(cs%KE_dia)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k) * g%dxCu(i,j) * adp%du_dt_dia(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k) * g%dyCv(i,j) * adp%dv_dt_dia(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = cs%KE(i,j,k) &
            * (us%T_to_s * (cdp%diapyc_vel(i,j,k) - cdp%diapyc_vel(i,j,k+1)))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_dia(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) &
            * (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_dia > 0) call post_data(cs%id_KE_dia, cs%KE_dia, cs%diag)
  endif
 

Referenced by calculate_diagnostic_fields().

Here is the caller graph for this function:

calculate_vertical_integrals()

subroutine mom_diagnostics::calculate_vertical_integrals ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, real, dimension(:,:), pointer  p_surf, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(diagnostics_cs), intent(inout)  CS  )

private

This 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]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2].
[in]	tv	A structure pointing to various thermodynamic variables.
	p_surf	A pointer to the surface pressure [Pa]. If p_surf is not associated, it is the same as setting the surface pressure to 0.
[in,out]	cs	Control structure returned by a previous call to diagnostics_init.

Definition at line 762 of file MOM_diagnostics.F90.

  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various
                                                 !! thermodynamic variables.
  real, dimension(:,:),    pointer       :: p_surf !< A pointer to the surface pressure [Pa].
                                                 !! If p_surf is not associated, it is the same
                                                 !! as setting the surface pressure to 0.
  type(diagnostics_CS),    intent(inout) :: CS   !< Control structure returned by a
                                                 !! previous call to diagnostics_init.
 
  real, dimension(SZI_(G), SZJ_(G)) :: &
    z_top, &  ! Height of the top of a layer or the ocean [Z ~> m].
    z_bot, &  ! Height of the bottom of a layer (for id_mass) or the
              ! (positive) depth of the ocean (for id_col_ht) [Z ~> m].
    mass, &   ! integrated mass of the water column [kg m-2].  For
              ! non-Boussinesq models this is rho*dz. For Boussinesq
              ! models, this is either the integral of in-situ density
              ! (rho*dz for col_mass) or reference density (Rho_0*dz for mass_wt).
    btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
              ! This is the column mass multiplied by gravity plus the pressure
              ! at the ocean surface [Pa].
    dpress, & ! Change in hydrostatic pressure across a layer [Pa].
    tr_int    ! vertical integral of a tracer times density,
              ! (Rho_0 in a Boussinesq model) [TR kg m-2].
  real    :: IG_Earth  ! Inverse of gravitational acceleration [s2 m-1].
 
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  if (cs%id_mass_wt > 0) then
    do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_mass_wt, mass, cs%diag)
  endif
 
  if (cs%id_temp_int > 0) then
    do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%T(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_temp_int, tr_int, cs%diag)
  endif
 
  if (cs%id_salt_int > 0) then
    do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%S(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_salt_int, tr_int, cs%diag)
  endif
 
  if (cs%id_col_ht > 0) then
    call find_eta(h, tv, g, gv, us, z_top)
    do j=js,je ; do i=is,ie
      z_bot(i,j) = z_top(i,j) + g%bathyT(i,j)
    enddo ; enddo
    call post_data(cs%id_col_ht, z_bot, cs%diag)
  endif
 
  ! NOTE: int_density_z expects z_top and z_btm values from [ij]sq to [ij]eq+1
  if (cs%id_col_mass > 0 .or. cs%id_pbo > 0) then
    do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
    if (gv%Boussinesq) then
      if (associated(tv%eqn_of_state)) then
        ig_earth = 1.0 / gv%mks_g_Earth
!       do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo
        do j=g%jscB,g%jecB+1 ; do i=g%iscB,g%iecB+1
          z_bot(i,j) = 0.0
        enddo ; enddo
        do k=1,nz
          do j=g%jscB,g%jecB+1 ; do i=g%iscB,g%iecB+1
            z_top(i,j) = z_bot(i,j)
            z_bot(i,j) = z_top(i,j) - gv%H_to_Z*h(i,j,k)
          enddo ; enddo
          call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), &
                              z_top, z_bot, 0.0, us%R_to_kg_m3*gv%Rho0, gv%mks_g_Earth*us%Z_to_m, &
                              g%HI, g%HI, tv%eqn_of_state, dpress)
          do j=js,je ; do i=is,ie
            mass(i,j) = mass(i,j) + dpress(i,j) * ig_earth
          enddo ; enddo
        enddo
      else
        do k=1,nz ; do j=js,je ; do i=is,ie
          mass(i,j) = mass(i,j) + (gv%H_to_m*us%R_to_kg_m3*gv%Rlay(k))*h(i,j,k)
        enddo ; enddo ; enddo
      endif
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (cs%id_col_mass > 0) then
      call post_data(cs%id_col_mass, mass, cs%diag)
    endif
    if (cs%id_pbo > 0) then
      do j=js,je ; do i=is,ie ; btm_pres(i,j) = 0.0 ; enddo ; enddo
      ! 'pbo' is defined as the sea water pressure at the sea floor
      !     pbo = (mass * g) + p_surf
      ! where p_surf is the sea water pressure at sea water surface.
      do j=js,je ; do i=is,ie
        btm_pres(i,j) = mass(i,j) * gv%mks_g_Earth
        if (associated(p_surf)) then
          btm_pres(i,j) = btm_pres(i,j) + p_surf(i,j)
        endif
      enddo ; enddo
      call post_data(cs%id_pbo, btm_pres, cs%diag)
    endif
  endif
 

Referenced by calculate_diagnostic_fields().

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 the 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.

Parameters

[in]	rlist	The list of target densities [R ~> kg m-3]
[in]	r_in	The density being inserted into Rlist [R ~> kg m-3]
[in,out]	k	The value of k such that Rlist(k) <= R_in < Rlist(k+1) The input value is a first guess
[in]	nz	The number of layers in Rlist
[out]	wt	The weight of layer k for interpolation, nondim
[out]	wt_p	The weight of layer k+1 for interpolation, nondim

Definition at line 695 of file MOM_diagnostics.F90.

  real, dimension(:), &
            intent(in)    :: Rlist !< The list of target densities [R ~> kg m-3]
  real,     intent(in)    :: R_in !< The density being inserted into Rlist [R ~> kg m-3]
  integer,  intent(inout) :: k    !< The value of k such that Rlist(k) <= R_in < Rlist(k+1)
                                  !! The input value is a first guess
  integer,  intent(in)    :: nz   !< The number of layers in Rlist
  real,     intent(out)   :: wt   !< The weight of layer k for interpolation, nondim
  real,     intent(out)   :: wt_p !< The weight of layer k+1 for interpolation, nondim
 
  ! 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.
 
  integer :: k_upper, k_lower, k_new, inc
 
  ! First, bracket the desired point.
  if ((k < 1) .or. (k > nz)) k = nz/2
 
  k_upper = k ; k_lower = k ; inc = 1
  if (r_in < rlist(k)) then
    do
      k_lower = max(k_lower-inc, 1)
      if ((k_lower == 1) .or. (r_in >= rlist(k_lower))) exit
      k_upper = k_lower
      inc = inc*2
    enddo
  else
    do
      k_upper = min(k_upper+inc, nz)
      if ((k_upper == nz) .or. (r_in < rlist(k_upper))) exit
      k_lower = k_upper
      inc = inc*2
    enddo
  endif
 
  if ((k_lower == 1) .and. (r_in <= rlist(k_lower))) then
    k = 1 ; wt = 1.0 ; wt_p = 0.0
  elseif ((k_upper == nz) .and. (r_in >= rlist(k_upper))) then
    k = nz-1 ; wt = 0.0 ; wt_p = 1.0
  else
    do
      if (k_upper <= k_lower+1) exit
      k_new = (k_upper + k_lower) / 2
      if (r_in < rlist(k_new)) then
        k_upper = k_new
      else
        k_lower = k_new
      endif
    enddo
 
!   Uncomment this as a code check
!    if ((R_in < Rlist(k_lower)) .or. (R_in >= Rlist(k_upper)) .or. (k_upper-k_lower /= 1)) &
!      write (*,*) "Error: ",R_in," is not between R(",k_lower,") = ", &
!        Rlist(k_lower)," and R(",k_upper,") = ",Rlist(k_upper),"."
    k = k_lower
    wt = (rlist(k_upper) - r_in) / (rlist(k_upper) - rlist(k_lower))
    wt_p = 1.0 - wt
 
  endif
 

Referenced by calculate_diagnostic_fields().

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
	)

Deallocate memory associated with the diagnostics module.

Parameters

	cs	Control structure returned by a previous call to diagnostics_init.
[in,out]	adp	structure with pointers to accelerations in momentum equation.

Definition at line 2088 of file MOM_diagnostics.F90.

  type(diagnostics_CS),   pointer       :: CS  !< Control structure returned by a
                                               !! previous call to diagnostics_init.
  type(accel_diag_ptrs),  intent(inout) :: ADp !< structure with pointers to
                                               !! accelerations in momentum equation.
  integer :: m
 
  if (associated(cs%e))          deallocate(cs%e)
  if (associated(cs%e_D))        deallocate(cs%e_D)
  if (associated(cs%KE))         deallocate(cs%KE)
  if (associated(cs%dKE_dt))     deallocate(cs%dKE_dt)
  if (associated(cs%PE_to_KE))   deallocate(cs%PE_to_KE)
  if (associated(cs%KE_Coradv))  deallocate(cs%KE_Coradv)
  if (associated(cs%KE_adv))     deallocate(cs%KE_adv)
  if (associated(cs%KE_visc))    deallocate(cs%KE_visc)
  if (associated(cs%KE_horvisc)) deallocate(cs%KE_horvisc)
  if (associated(cs%KE_dia))     deallocate(cs%KE_dia)
  if (associated(cs%dv_dt))      deallocate(cs%dv_dt)
  if (associated(cs%dh_dt))      deallocate(cs%dh_dt)
  if (associated(cs%du_dt))      deallocate(cs%du_dt)
  if (associated(cs%h_Rlay))     deallocate(cs%h_Rlay)
  if (associated(cs%uh_Rlay))    deallocate(cs%uh_Rlay)
  if (associated(cs%vh_Rlay))    deallocate(cs%vh_Rlay)
  if (associated(cs%uhGM_Rlay))  deallocate(cs%uhGM_Rlay)
  if (associated(cs%vhGM_Rlay))  deallocate(cs%vhGM_Rlay)
 
  if (associated(adp%gradKEu))    deallocate(adp%gradKEu)
  if (associated(adp%gradKEu))    deallocate(adp%gradKEu)
  if (associated(adp%du_dt_visc)) deallocate(adp%du_dt_visc)
  if (associated(adp%dv_dt_visc)) deallocate(adp%dv_dt_visc)
  if (associated(adp%du_dt_dia))  deallocate(adp%du_dt_dia)
  if (associated(adp%dv_dt_dia))  deallocate(adp%dv_dt_dia)
  if (associated(adp%du_other))   deallocate(adp%du_other)
  if (associated(adp%dv_other))   deallocate(adp%dv_other)
 
  do m=1,cs%num_time_deriv ; deallocate(cs%prev_val(m)%p) ; enddo
 
  deallocate(cs)
 

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(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(diagnostics_cs), pointer	CS,
		type(thermo_var_ptrs), intent(in)	tv
	)

This subroutine registers various diagnostics and allocates space for fields that other diagnostis depend upon.

Parameters

[in]	mis	For "MOM Internal State" a set of pointers to the fields and accelerations that make up the ocean's internal physical state.
[in,out]	adp	Structure with pointers to momentum equation terms.
[in,out]	cdp	Structure with pointers to continuity equation terms.
[in]	time	Current model time.
[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	param_file	A structure to parse for run-time parameters.
[in,out]	diag	Structure to regulate diagnostic output.
	cs	Pointer set to point to control structure for this module.
[in]	tv	A structure pointing to various thermodynamic variables.

Definition at line 1423 of file MOM_diagnostics.F90.

  type(ocean_internal_state), intent(in)    :: MIS  !< For "MOM Internal State" a set of pointers to
                                                    !! the fields and accelerations that make up the
                                                    !! ocean's internal physical state.
  type(accel_diag_ptrs),      intent(inout) :: ADp  !< Structure with pointers to momentum equation
                                                    !! terms.
  type(cont_diag_ptrs),       intent(inout) :: CDp  !< Structure with pointers to continuity
                                                    !! equation terms.
  type(time_type),            intent(in)    :: Time !< Current model time.
  type(ocean_grid_type),      intent(in)    :: G    !< The ocean's grid structure.
  type(verticalGrid_type),    intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(unit_scale_type),      intent(in)    :: US   !< A dimensional unit scaling type
  type(param_file_type),      intent(in)    :: param_file !< A structure to parse for run-time
                                                    !! parameters.
  type(diag_ctrl), target,    intent(inout) :: diag !< Structure to regulate diagnostic output.
  type(diagnostics_CS),       pointer       :: CS   !< Pointer set to point to control structure
                                                    !! for this module.
  type(thermo_var_ptrs),      intent(in)    :: tv   !< A structure pointing to various
                                                    !! thermodynamic variables.
 
  ! Local variables
  real :: omega, f2_min, convert_H
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_diagnostics" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  logical :: use_temperature, adiabatic
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nkml, nkbl
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j
 
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  if (associated(cs)) then
    call mom_error(warning, "MOM_diagnostics_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)
 
  cs%diag => diag
  use_temperature = associated(tv%T)
  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
                 do_not_log=.true.)
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_COLUMN_FRACTION", cs%mono_N2_column_fraction, &
                 "The lower fraction of water column over which N2 is limited as monotonic "// &
                 "for the purposes of calculating the equivalent barotropic wave speed.", &
                 units='nondim', default=0.)
  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", cs%mono_N2_depth, &
                 "The depth below which N2 is limited as monotonic for the "// &
                 "purposes of calculating the equivalent barotropic wave speed.", &
                 units='m', scale=us%m_to_Z, default=-1.)
 
  if (gv%Boussinesq) then
    thickness_units = "m" ; flux_units = "m3 s-1" ; convert_h = gv%H_to_m
  else
    thickness_units = "kg m-2" ; flux_units = "kg s-1" ; convert_h = gv%H_to_kg_m2
  endif
 
  cs%id_masscello = register_diag_field('ocean_model', 'masscello', diag%axesTL,&
      time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2',           &
      standard_name='sea_water_mass_per_unit_area', v_extensive=.true.)
 
  cs%id_masso = register_scalar_field('ocean_model', 'masso', time,  &
      diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass')
 
  cs%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, time, &
      long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.)
  cs%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, time, &
      long_name = 'Cell thickness from the previous timestep', units='m', &
      v_extensive=.true., conversion=gv%H_to_m)
 
  ! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and
  ! must be registered first. We earlier stored the handle of volcello but need it here for posting
  ! by this module.
  cs%id_volcello = diag_get_volume_cell_measure_dm_id(diag)
 
  if (use_temperature) then
    if (tv%T_is_conT) then
      cs%id_Tpot = register_diag_field('ocean_model', 'temp', diag%axesTL, &
          time, 'Potential Temperature', 'degC')
    endif
    if (tv%S_is_absS) then
      cs%id_Sprac = register_diag_field('ocean_model', 'salt', diag%axesTL, &
          time, 'Salinity', 'psu')
    endif
 
    cs%id_tob = register_diag_field('ocean_model','tob', diag%axesT1, time, &
        long_name='Sea Water Potential Temperature at Sea Floor', &
        standard_name='sea_water_potential_temperature_at_sea_floor', units='degC')
    cs%id_sob = register_diag_field('ocean_model','sob',diag%axesT1, time, &
        long_name='Sea Water Salinity at Sea Floor', &
        standard_name='sea_water_salinity_at_sea_floor', units='psu')
 
    cs%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', &
        diag%axesZL, time, 'Layer Average Ocean Temperature', 'degC')
    cs%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', &
        diag%axesZL, time, 'Layer Average Ocean Salinity', 'psu')
 
    cs%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga', &
        time, diag, 'Global Mean Ocean Potential Temperature', 'degC',&
        standard_name='sea_water_potential_temperature')
    cs%id_soga = register_scalar_field('ocean_model', 'soga', &
        time, diag, 'Global Mean Ocean Salinity', 'psu', &
        standard_name='sea_water_salinity')
 
    cs%id_tosga = register_scalar_field('ocean_model', 'sst_global', time, diag,&
        long_name='Global Area Average Sea Surface Temperature',                &
        units='degC', standard_name='sea_surface_temperature',                  &
        cmor_field_name='tosga', cmor_standard_name='sea_surface_temperature',  &
        cmor_long_name='Sea Surface Temperature')
    cs%id_sosga = register_scalar_field('ocean_model', 'sss_global', time, diag,&
        long_name='Global Area Average Sea Surface Salinity',                   &
        units='psu', standard_name='sea_surface_salinity',                      &
        cmor_field_name='sosga', cmor_standard_name='sea_surface_salinity',     &
        cmor_long_name='Sea Surface Salinity')
  endif
 
  cs%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, time,              &
      'Zonal velocity', 'm s-1', conversion=us%L_T_to_m_s, cmor_field_name='uo', &
      cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
  cs%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, time,                  &
      'Meridional velocity', 'm s-1', conversion=us%L_T_to_m_s, cmor_field_name='vo', &
      cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
  cs%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, time, &
      'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_h)
 
  cs%id_e = register_diag_field('ocean_model', 'e', diag%axesTi, time, &
      'Interface Height Relative to Mean Sea Level', 'm', conversion=us%Z_to_m)
  if (cs%id_e>0) call safe_alloc_ptr(cs%e,isd,ied,jsd,jed,nz+1)
 
  cs%id_e_D = register_diag_field('ocean_model', 'e_D', diag%axesTi, time, &
      'Interface Height above the Seafloor', 'm', conversion=us%Z_to_m)
  if (cs%id_e_D>0) call safe_alloc_ptr(cs%e_D,isd,ied,jsd,jed,nz+1)
 
  cs%id_Rml = register_diag_field('ocean_model', 'Rml', diag%axesTL, time, &
      'Mixed Layer Coordinate Potential Density', 'kg m-3', conversion=us%R_to_kg_m3)
 
  cs%id_Rcv = register_diag_field('ocean_model', 'Rho_cv', diag%axesTL, time, &
      'Coordinate Potential Density', 'kg m-3', conversion=us%R_to_kg_m3)
 
  cs%id_rhopot0 = register_diag_field('ocean_model', 'rhopot0', diag%axesTL, time, &
      'Potential density referenced to surface', 'kg m-3')
  cs%id_rhopot2 = register_diag_field('ocean_model', 'rhopot2', diag%axesTL, time, &
      'Potential density referenced to 2000 dbar', 'kg m-3')
  cs%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, time, &
      'In situ density', 'kg m-3')
 
  cs%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, time, &
      'Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  if ((cs%id_du_dt>0) .and. .not.associated(cs%du_dt)) then
    call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
    call register_time_deriv(lbound(mis%u), mis%u, cs%du_dt, cs)
  endif
 
  cs%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, time, &
      'Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  if ((cs%id_dv_dt>0) .and. .not.associated(cs%dv_dt)) then
    call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
    call register_time_deriv(lbound(mis%v), mis%v, cs%dv_dt, cs)
  endif
 
  cs%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, time, &
      'Thickness tendency', trim(thickness_units)//" s-1", conversion=convert_h*us%s_to_T, v_extensive=.true.)
  if ((cs%id_dh_dt>0) .and. .not.associated(cs%dh_dt)) then
    call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
    call register_time_deriv(lbound(mis%h), mis%h, cs%dh_dt, cs)
  endif
 
  ! layer thickness variables
  !if (GV%nk_rho_varies > 0) then
    cs%id_h_Rlay = register_diag_field('ocean_model', 'h_rho', diag%axesTL, time, &
        'Layer thicknesses in pure potential density coordinates', thickness_units, &
        conversion=convert_h)
    if (cs%id_h_Rlay>0) call safe_alloc_ptr(cs%h_Rlay,isd,ied,jsd,jed,nz)
 
    cs%id_uh_Rlay = register_diag_field('ocean_model', 'uh_rho', diag%axesCuL, time, &
        'Zonal volume transport in pure potential density coordinates', flux_units, &
        conversion=us%L_to_m**2*us%s_to_T*convert_h)
    if (cs%id_uh_Rlay>0) call safe_alloc_ptr(cs%uh_Rlay,isdb,iedb,jsd,jed,nz)
 
    cs%id_vh_Rlay = register_diag_field('ocean_model', 'vh_rho', diag%axesCvL, time, &
        'Meridional volume transport in pure potential density coordinates', flux_units, &
        conversion=us%L_to_m**2*us%s_to_T*convert_h)
    if (cs%id_vh_Rlay>0) call safe_alloc_ptr(cs%vh_Rlay,isd,ied,jsdb,jedb,nz)
 
    cs%id_uhGM_Rlay = register_diag_field('ocean_model', 'uhGM_rho', diag%axesCuL, time, &
        'Zonal volume transport due to interface height diffusion in pure potential '//&
        'density coordinates', flux_units, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    if (cs%id_uhGM_Rlay>0) call safe_alloc_ptr(cs%uhGM_Rlay,isdb,iedb,jsd,jed,nz)
 
    cs%id_vhGM_Rlay = register_diag_field('ocean_model', 'vhGM_rho', diag%axesCvL, time, &
        'Meridional volume transport due to interface height diffusion in pure potential '//&
        'density coordinates', flux_units, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    if (cs%id_vhGM_Rlay>0) call safe_alloc_ptr(cs%vhGM_Rlay,isd,ied,jsdb,jedb,nz)
  !endif
 
 
  ! terms in the kinetic energy budget
  cs%id_KE = register_diag_field('ocean_model', 'KE', diag%axesTL, time, &
      'Layer kinetic energy per unit mass', 'm2 s-2', &
      conversion=us%L_T_to_m_s**2)
  if (cs%id_KE>0) call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
 
  cs%id_dKEdt = register_diag_field('ocean_model', 'dKE_dt', diag%axesTL, time, &
      'Kinetic Energy Tendency of Layer', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_dKEdt>0) call safe_alloc_ptr(cs%dKE_dt,isd,ied,jsd,jed,nz)
 
  cs%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, time, &
      'Potential to Kinetic Energy Conversion of Layer', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_PE_to_KE>0) call safe_alloc_ptr(cs%PE_to_KE,isd,ied,jsd,jed,nz)
 
  cs%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, time, &
      'Kinetic Energy Source from Coriolis and Advection', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_KE_Coradv>0) call safe_alloc_ptr(cs%KE_Coradv,isd,ied,jsd,jed,nz)
 
  cs%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, time, &
      'Kinetic Energy Source from Advection', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_KE_adv>0) call safe_alloc_ptr(cs%KE_adv,isd,ied,jsd,jed,nz)
 
  cs%id_KE_visc = register_diag_field('ocean_model', 'KE_visc', diag%axesTL, time, &
      'Kinetic Energy Source from Vertical Viscosity and Stresses', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_KE_visc>0) call safe_alloc_ptr(cs%KE_visc,isd,ied,jsd,jed,nz)
 
  cs%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, time, &
      'Kinetic Energy Source from Horizontal Viscosity', 'm3 s-3', &
      conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_KE_horvisc>0) call safe_alloc_ptr(cs%KE_horvisc,isd,ied,jsd,jed,nz)
 
  if (.not. adiabatic) then
    cs%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, time, &
        'Kinetic Energy Source from Diapycnal Diffusion', 'm3 s-3', &
        conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
    if (cs%id_KE_dia>0) call safe_alloc_ptr(cs%KE_dia,isd,ied,jsd,jed,nz)
  endif
 
  ! gravity wave CFLs
  cs%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, time, &
      'First baroclinic gravity wave speed', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_Rd1 = register_diag_field('ocean_model', 'Rd1', diag%axesT1, time, &
      'First baroclinic deformation radius', 'm', conversion=us%L_to_m)
  cs%id_cfl_cg1 = register_diag_field('ocean_model', 'CFL_cg1', diag%axesT1, time, &
      'CFL of first baroclinic gravity wave = dt*cg1*(1/dx+1/dy)', 'nondim')
  cs%id_cfl_cg1_x = register_diag_field('ocean_model', 'CFL_cg1_x', diag%axesT1, time, &
      'i-component of CFL of first baroclinic gravity wave = dt*cg1*/dx', 'nondim')
  cs%id_cfl_cg1_y = register_diag_field('ocean_model', 'CFL_cg1_y', diag%axesT1, time, &
      'j-component of CFL of first baroclinic gravity wave = dt*cg1*/dy', 'nondim')
  cs%id_cg_ebt = register_diag_field('ocean_model', 'cg_ebt', diag%axesT1, time, &
      'Equivalent barotropic gravity wave speed', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_Rd_ebt = register_diag_field('ocean_model', 'Rd_ebt', diag%axesT1, time, &
      'Equivalent barotropic deformation radius', 'm', conversion=us%L_to_m)
  cs%id_p_ebt = register_diag_field('ocean_model', 'p_ebt', diag%axesTL, time, &
      'Equivalent barotropic modal strcuture', 'nondim')
 
  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
      (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0) .or. &
      (cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
    call wave_speed_init(cs%wave_speed_CSp)
    call safe_alloc_ptr(cs%cg1,isd,ied,jsd,jed)
    if (cs%id_Rd1>0)       call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
    if (cs%id_Rd_ebt>0)    call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1>0)   call safe_alloc_ptr(cs%cfl_cg1,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1_x>0) call safe_alloc_ptr(cs%cfl_cg1_x,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1_y>0) call safe_alloc_ptr(cs%cfl_cg1_y,isd,ied,jsd,jed)
    if (cs%id_p_ebt>0) call safe_alloc_ptr(cs%p_ebt,isd,ied,jsd,jed,nz)
  endif
 
  cs%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, time,  &
      'The column mass for calculating mass-weighted average properties', 'kg m-2')
 
  if (use_temperature) then
    cs%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, time,                &
        'Density weighted column integrated potential temperature', 'degC kg m-2',                    &
        cmor_field_name='opottempmint',                                                               &
        cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',&
        cmor_standard_name='Depth integrated density times potential temperature')
 
    cs%id_salt_int = register_diag_field('ocean_model', 'salt_int', diag%axesT1, time,   &
        'Density weighted column integrated salinity', 'psu kg m-2',                     &
        cmor_field_name='somint',                                                        &
        cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity',&
        cmor_standard_name='Depth integrated density times salinity')
  endif
 
  cs%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, time, &
      'The column integrated in situ density', 'kg m-2')
 
  cs%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, time, &
      'The height of the water column', 'm', conversion=us%Z_to_m)
  cs%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, time, &
      long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
      units='Pa')
 
  call set_dependent_diagnostics(mis, adp, cdp, g, cs)
 

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

Here is the call graph for this function:

post_surface_dyn_diags()

subroutine, public mom_diagnostics::post_surface_dyn_diags	(	type(surface_diag_ids), intent(in)	IDs,
		type(ocean_grid_type), intent(in)	G,
		type(diag_ctrl), intent(in)	diag,
		type(surface), intent(in)	sfc_state,
		real, dimension(szi_(g),szj_(g)), intent(in)	ssh
	)

This routine posts diagnostics of various dynamic ocean surface quantities, including velocities, speed and sea surface height, at the time the ocean state is reported back to the caller.

Parameters

[in]	ids	A structure with the diagnostic IDs.
[in]	g	ocean grid structure
[in]	diag	regulates diagnostic output
[in]	sfc_state	structure describing the ocean surface state
[in]	ssh	Time mean surface height without corrections for ice displacement [m]

Definition at line 1158 of file MOM_diagnostics.F90.

  type(surface_diag_IDs),   intent(in) :: IDs !< A structure with the diagnostic IDs.
  type(ocean_grid_type),    intent(in) :: G   !< ocean grid structure
  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
  type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
  real, dimension(SZI_(G),SZJ_(G)), &
                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
 
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
  integer :: i, j, is, ie, js, je
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  if (ids%id_ssh > 0) &
    call post_data(ids%id_ssh, ssh, diag, mask=g%mask2dT)
 
  if (ids%id_ssu > 0) &
    call post_data(ids%id_ssu, sfc_state%u, diag, mask=g%mask2dCu)
 
  if (ids%id_ssv > 0) &
    call post_data(ids%id_ssv, sfc_state%v, diag, mask=g%mask2dCv)
 
  if (ids%id_speed > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = sqrt(0.5*(sfc_state%u(i-1,j)**2 + sfc_state%u(i,j)**2) + &
                            0.5*(sfc_state%v(i,j-1)**2 + sfc_state%v(i,j)**2))
    enddo ; enddo
    call post_data(ids%id_speed, work_2d, diag, mask=g%mask2dT)
  endif
 

post_surface_thermo_diags()

subroutine, public mom_diagnostics::post_surface_thermo_diags	(	type(surface_diag_ids), intent(in)	IDs,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(diag_ctrl), intent(in)	diag,
		real, intent(in)	dt_int,
		type(surface), intent(in)	sfc_state,
		type(thermo_var_ptrs), intent(in)	tv,
		real, dimension(szi_(g),szj_(g)), intent(in)	ssh,
		real, dimension(szi_(g),szj_(g)), intent(in)	ssh_ibc
	)

This routine posts diagnostics of various ocean surface and integrated quantities at the time the ocean state is reported back to the caller.

Parameters

[in]	ids	A structure with the diagnostic IDs.
[in]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	diag	regulates diagnostic output
[in]	dt_int	total time step associated with these diagnostics [T ~> s].
[in]	sfc_state	structure describing the ocean surface state
[in]	tv	A structure pointing to various thermodynamic variables
[in]	ssh	Time mean surface height without corrections for ice displacement [m]
[in]	ssh_ibc	Time mean surface height with corrections for ice displacement and the inverse barometer [m]

Definition at line 1195 of file MOM_diagnostics.F90.

  type(surface_diag_IDs),   intent(in) :: IDs !< A structure with the diagnostic IDs.
  type(ocean_grid_type),    intent(in) :: G   !< ocean grid structure
  type(verticalGrid_type),  intent(in) :: GV  !< ocean vertical grid structure
  type(unit_scale_type),    intent(in) :: US  !< A dimensional unit scaling type
  type(diag_ctrl),          intent(in) :: diag  !< regulates diagnostic output
  real,                     intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s].
  type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
  type(thermo_var_ptrs),    intent(in) :: tv  !< A structure pointing to various thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G)), &
                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
                                                  !! for ice displacement and the inverse barometer [m]
 
  real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
  real, dimension(SZI_(G),SZJ_(G)) :: &
    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m]
  real :: I_time_int    ! The inverse of the time interval [T-1 ~> s-1].
  real :: zos_area_mean, volo, ssh_ga
  integer :: i, j, is, ie, js, je
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  ! area mean SSH
  if (ids%id_ssh_ga > 0) then
    ssh_ga = global_area_mean(ssh, g)
    call post_data(ids%id_ssh_ga, ssh_ga, diag)
  endif
 
  ! post the dynamic sea level, zos, and zossq.
  ! zos is ave_ssh with sea ice inverse barometer removed,
  ! and with zero global area mean.
  if (ids%id_zos > 0 .or. ids%id_zossq > 0) then
    zos(:,:) = 0.0
    do j=js,je ; do i=is,ie
      zos(i,j) = ssh_ibc(i,j)
    enddo ; enddo
    zos_area_mean = global_area_mean(zos, g)
    do j=js,je ; do i=is,ie
      zos(i,j) = zos(i,j) - g%mask2dT(i,j)*zos_area_mean
    enddo ; enddo
    if (ids%id_zos > 0) call post_data(ids%id_zos, zos, diag, mask=g%mask2dT)
    if (ids%id_zossq > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = zos(i,j)*zos(i,j)
      enddo ; enddo
      call post_data(ids%id_zossq, work_2d, diag, mask=g%mask2dT)
    endif
  endif
 
  ! post total volume of the liquid ocean
  if (ids%id_volo > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = g%mask2dT(i,j)*(ssh(i,j) + us%Z_to_m*g%bathyT(i,j))
    enddo ; enddo
    volo = global_area_integral(work_2d, g)
    call post_data(ids%id_volo, volo, diag)
  endif
 
  ! Use Adcroft's rule of reciprocals; it does the right thing here.
  i_time_int = 0.0 ; if (dt_int > 0.0) i_time_int = 1.0 / dt_int
 
  ! post time-averaged rate of frazil formation
  if (associated(tv%frazil) .and. (ids%id_fraz > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%frazil(i,j) * i_time_int
    enddo ; enddo
    call post_data(ids%id_fraz, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post time-averaged salt deficit
  if (associated(tv%salt_deficit) .and. (ids%id_salt_deficit > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%salt_deficit(i,j) * i_time_int
    enddo ; enddo
    call post_data(ids%id_salt_deficit, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post temperature of P-E+R
  if (associated(tv%TempxPmE) .and. (ids%id_Heat_PmE > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%TempxPmE(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(ids%id_Heat_PmE, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post geothermal heating or internal heat source/sinks
  if (associated(tv%internal_heat) .and. (ids%id_intern_heat > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%internal_heat(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(ids%id_intern_heat, work_2d, diag, mask=g%mask2dT)
  endif
 
  if (tv%T_is_conT) then
    ! Internal T&S variables are conservative temperature & absolute salinity
    if (ids%id_sstcon > 0) call post_data(ids%id_sstcon, sfc_state%SST, diag, mask=g%mask2dT)
    ! Use TEOS-10 function calls convert T&S diagnostics from conservative temp
    ! to potential temperature.
    do j=js,je ; do i=is,ie
      work_2d(i,j) = gsw_pt_from_ct(sfc_state%SSS(i,j),sfc_state%SST(i,j))
    enddo ; enddo
    if (ids%id_sst > 0) call post_data(ids%id_sst, work_2d, diag, mask=g%mask2dT)
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (ids%id_sst > 0) call post_data(ids%id_sst, sfc_state%SST, diag, mask=g%mask2dT)
  endif
 
  if (tv%S_is_absS) then
    ! Internal T&S variables are conservative temperature & absolute salinity
    if (ids%id_sssabs > 0) call post_data(ids%id_sssabs, sfc_state%SSS, diag, mask=g%mask2dT)
    ! Use TEOS-10 function calls convert T&S diagnostics from absolute salinity
    ! to practical salinity.
    do j=js,je ; do i=is,ie
      work_2d(i,j) = gsw_sp_from_sr(sfc_state%SSS(i,j))
    enddo ; enddo
    if (ids%id_sss > 0) call post_data(ids%id_sss, work_2d, diag, mask=g%mask2dT)
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (ids%id_sss > 0) call post_data(ids%id_sss, sfc_state%SSS, diag, mask=g%mask2dT)
  endif
 
  if (ids%id_sst_sq > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = sfc_state%SST(i,j)*sfc_state%SST(i,j)
    enddo ; enddo
    call post_data(ids%id_sst_sq, work_2d, diag, mask=g%mask2dT)
  endif
  if (ids%id_sss_sq > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = sfc_state%SSS(i,j)*sfc_state%SSS(i,j)
    enddo ; enddo
    call post_data(ids%id_sss_sq, work_2d, diag, mask=g%mask2dT)
  endif
 
  call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag))
 

References mom_spatial_means::global_area_integral().

Here is the call graph for this function:

post_transport_diagnostics()

subroutine, public mom_diagnostics::post_transport_diagnostics	(	type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)	uhtr,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)	vhtr,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(transport_diag_ids), intent(in)	IDs,
		type(diag_grid_storage), intent(inout)	diag_pre_dyn,
		type(diag_ctrl), intent(inout)	diag,
		real, intent(in)	dt_trans,
		type(tracer_registry_type), pointer	Reg
	)

This routine posts diagnostics of the transports, including the subgridscale contributions.

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	uhtr	Accumulated zonal thickness fluxes used to advect tracers [H L2 ~> m3 or kg]
[in]	vhtr	Accumulated meridional thickness fluxes used to advect tracers [H L2 ~> m3 or kg]
[in]	h	The updated layer thicknesses [H ~> m or kg m-2]
[in]	ids	A structure with the diagnostic IDs.
[in,out]	diag_pre_dyn	Stored grids from before dynamics
[in,out]	diag	regulates diagnostic output
[in]	dt_trans	total time step associated with the transports [T ~> s].
	reg	Pointer to the tracer registry

Definition at line 1338 of file MOM_diagnostics.F90.

  type(ocean_grid_type),    intent(inout) :: G   !< ocean grid structure
  type(verticalGrid_type),  intent(in)    :: GV  !< ocean vertical grid structure
  type(unit_scale_type),    intent(in)    :: US  !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes
                                                 !! used to advect tracers [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes
                                                 !! used to advect tracers [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: h   !< The updated layer thicknesses [H ~> m or kg m-2]
  type(transport_diag_IDs), intent(in)    :: IDs !< A structure with the diagnostic IDs.
  type(diag_grid_storage),  intent(inout) :: diag_pre_dyn !< Stored grids from before dynamics
  type(diag_ctrl),          intent(inout) :: diag !< regulates diagnostic output
  real,                     intent(in)    :: dt_trans !< total time step associated with the transports [T ~> s].
  type(tracer_registry_type), pointer     :: Reg !< Pointer to the tracer registry
 
  ! Local variables
  real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport [kg s-1]
  real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport [kg s-1]
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport [kg s-1]
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport [kg s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))    :: h_tend ! Change in layer thickness due to dynamics
                          ! [H s-1 ~> m s-1 or kg m-2 s-1].
  real :: Idt             ! The inverse of the time interval [T-1 ~> s-1]
  real :: H_to_kg_m2_dt   ! A conversion factor from accumulated transports to fluxes
                          ! [kg L-2 H-1 T-1 ~> kg m-3 s-1 or s-1].
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  idt = 1. / dt_trans
  h_to_kg_m2_dt = gv%H_to_kg_m2 * us%L_to_m**2 * us%s_to_T * idt
 
  call diag_save_grids(diag)
  call diag_copy_storage_to_diag(diag, diag_pre_dyn)
 
  if (ids%id_umo_2d > 0) then
    umo2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo2d(i,j) = umo2d(i,j) + uhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_umo_2d, umo2d, diag)
  endif
  if (ids%id_umo > 0) then
    ! Convert to kg/s.
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo(i,j,k) = uhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_umo, umo, diag, alt_h = diag_pre_dyn%h_state)
  endif
  if (ids%id_vmo_2d > 0) then
    vmo2d(:,:) = 0.0
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo2d(i,j) = vmo2d(i,j) + vhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_vmo_2d, vmo2d, diag)
  endif
  if (ids%id_vmo > 0) then
    ! Convert to kg/s.
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo(i,j,k) = vhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_vmo, vmo, diag, alt_h = diag_pre_dyn%h_state)
  endif
 
  if (ids%id_uhtr > 0) call post_data(ids%id_uhtr, uhtr, diag, alt_h = diag_pre_dyn%h_state)
  if (ids%id_vhtr > 0) call post_data(ids%id_vhtr, vhtr, diag, alt_h = diag_pre_dyn%h_state)
  if (ids%id_dynamics_h > 0) call post_data(ids%id_dynamics_h, diag_pre_dyn%h_state, diag, &
                                            alt_h = diag_pre_dyn%h_state)
  ! Post the change in thicknesses
  if (ids%id_dynamics_h_tendency > 0) then
    h_tend(:,:,:) = 0.
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_tend(i,j,k) = (h(i,j,k) - diag_pre_dyn%h_state(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(ids%id_dynamics_h_tendency, h_tend, diag, alt_h = diag_pre_dyn%h_state)
  endif
 
  call post_tracer_transport_diagnostics(g, gv, reg, diag_pre_dyn%h_state, diag)
 
  call diag_restore_grids(diag)
 

References mom_diag_mediator::diag_copy_storage_to_diag(), mom_diag_mediator::diag_restore_grids(), mom_diag_mediator::diag_save_grids(), and mom_tracer_registry::post_tracer_transport_diagnostics().

Here is the call graph for this function:

register_surface_diags()

subroutine, public mom_diagnostics::register_surface_diags	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(unit_scale_type), intent(in)	US,
		type(surface_diag_ids), intent(inout)	IDs,
		type(diag_ctrl), intent(inout)	diag,
		type(thermo_var_ptrs), intent(in)	tv
	)

Register diagnostics of the surface state and integrated quantities.

Parameters

[in]	time	current model time
[in]	g	ocean grid structure
[in]	us	A dimensional unit scaling type
[in,out]	ids	A structure with the diagnostic IDs.
[in,out]	diag	regulates diagnostic output
[in]	tv	A structure pointing to various thermodynamic variables

Definition at line 1732 of file MOM_diagnostics.F90.

  type(time_type),         intent(in)    :: Time  !< current model time
  type(ocean_grid_type),   intent(in)    :: G     !< ocean grid structure
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(surface_diag_IDs),  intent(inout) :: IDs   !< A structure with the diagnostic IDs.
  type(diag_ctrl),         intent(inout) :: diag  !< regulates diagnostic output
  type(thermo_var_ptrs),   intent(in)    :: tv    !< A structure pointing to various thermodynamic variables
 
  ! Vertically integrated, budget, and surface state diagnostics
  ids%id_volo = register_scalar_field('ocean_model', 'volo', time, diag,&
      long_name='Total volume of liquid ocean', units='m3',            &
      standard_name='sea_water_volume')
  ids%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, time,&
      standard_name = 'sea_surface_height_above_geoid',                   &
      long_name= 'Sea surface height above geoid', units='m')
  ids%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, time,&
      standard_name='square_of_sea_surface_height_above_geoid',             &
      long_name='Square of sea surface height above geoid', units='m2')
  ids%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, time, &
      'Sea Surface Height', 'm')
  ids%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', time, diag,&
      long_name='Area averaged sea surface height', units='m',            &
      standard_name='area_averaged_sea_surface_height')
  ids%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, time, &
      'Sea Surface Zonal Velocity', 'm s-1')
  ids%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, time, &
      'Sea Surface Meridional Velocity', 'm s-1')
  ids%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, time, &
      'Sea Surface Speed', 'm s-1')
 
  if (associated(tv%T)) then
    ids%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, time,     &
        'Sea Surface Temperature', 'degC', cmor_field_name='tos', &
        cmor_long_name='Sea Surface Temperature',                                &
        cmor_standard_name='sea_surface_temperature')
    ids%id_sst_sq = register_diag_field('ocean_model', 'SST_sq', diag%axesT1, time, &
        'Sea Surface Temperature Squared', 'degC2', cmor_field_name='tossq', &
        cmor_long_name='Square of Sea Surface Temperature ',                      &
        cmor_standard_name='square_of_sea_surface_temperature')
    ids%id_sss = register_diag_field('ocean_model', 'SSS', diag%axesT1, time, &
        'Sea Surface Salinity', 'psu', cmor_field_name='sos', &
        cmor_long_name='Sea Surface Salinity',                            &
        cmor_standard_name='sea_surface_salinity')
    ids%id_sss_sq = register_diag_field('ocean_model', 'SSS_sq', diag%axesT1, time, &
        'Sea Surface Salinity Squared', 'psu', cmor_field_name='sossq', &
        cmor_long_name='Square of Sea Surface Salinity ',                     &
        cmor_standard_name='square_of_sea_surface_salinity')
    if (tv%T_is_conT) then
      ids%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, time,     &
          'Sea Surface Conservative Temperature', 'Celsius')
    endif
    if (tv%S_is_absS) then
      ids%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, time,     &
          'Sea Surface Absolute Salinity', 'g kg-1')
    endif
    if (associated(tv%frazil)) then
      ids%id_fraz = register_diag_field('ocean_model', 'frazil', diag%axesT1, time, &
            'Heat from frazil formation', 'W m-2', conversion=us%s_to_T, cmor_field_name='hfsifrazil', &
            cmor_standard_name='heat_flux_into_sea_water_due_to_frazil_ice_formation', &
            cmor_long_name='Heat Flux into Sea Water due to Frazil Ice Formation')
    endif
  endif
 
  ids%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, time, &
         'Salt sink in ocean due to ice flux', &
         'psu m-2 s-1', conversion=g%US%R_to_kg_m3*g%US%Z_to_m*us%s_to_T)
  ids%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, time, &
         'Heat flux into ocean from mass flux into ocean', &
         'W m-2', conversion=g%US%R_to_kg_m3*g%US%Z_to_m*us%s_to_T)
  ids%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, time,&
         'Heat flux into ocean from geothermal or other internal sources', 'W m-2', conversion=us%s_to_T)
 

register_time_deriv()

subroutine, public mom_diagnostics::register_time_deriv	(	integer, dimension(3), intent(in)	lb,
		real, dimension(lb(1):,lb(2):,:), target	f_ptr,
		real, dimension(lb(1):,lb(2):,:), target	deriv_ptr,
		type(diagnostics_cs), pointer	CS
	)

This subroutine registers fields to calculate a diagnostic time derivative.

Parameters

[in]	lb	Lower index bound of f_ptr
	f_ptr	Time derivative operand
	deriv_ptr	Time derivative of f_ptr
	cs	Control structure returned by previous call to diagnostics_init.

Definition at line 1085 of file MOM_diagnostics.F90.

  integer, intent(in), dimension(3) :: lb     !< Lower index bound of f_ptr
  real, dimension(lb(1):,lb(2):,:), target :: f_ptr
                                              !< Time derivative operand
  real, dimension(lb(1):,lb(2):,:), target :: deriv_ptr
                                              !< Time derivative of f_ptr
  type(diagnostics_CS),  pointer :: CS        !< Control structure returned by previous call to
                                              !! diagnostics_init.
 
  ! This subroutine registers fields to calculate a diagnostic time derivative.
  ! NOTE: Lower bound is required for grid indexing in calculate_derivs().
  !       We assume that the vertical axis is 1-indexed.
 
  integer :: m      !< New index of deriv_ptr in CS%deriv
  integer :: ub(3)  !< Upper index bound of f_ptr, based on shape.
 
  if (.not.associated(cs)) call mom_error(fatal, &
         "register_time_deriv: Module must be initialized before it is used.")
 
  if (cs%num_time_deriv >= max_fields_) then
    call mom_error(warning,"MOM_diagnostics:  Attempted to register more than " // &
                   "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
    return
  endif
 
  m = cs%num_time_deriv+1 ; cs%num_time_deriv = m
 
  ub(:) = lb(:) + shape(f_ptr) - 1
 
  cs%nlay(m) = size(f_ptr, 3)
  cs%deriv(m)%p => deriv_ptr
  allocate(cs%prev_val(m)%p(lb(1):ub(1), lb(2):ub(2), cs%nlay(m)))
 
  cs%var_ptr(m)%p => f_ptr
  cs%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
 

References mom_error_handler::mom_error().

Referenced by mom_diagnostics_init(), and set_dependent_diagnostics().

Here is the call graph for this function:

Here is the caller graph for this function:

register_transport_diags()

subroutine, public mom_diagnostics::register_transport_diags	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(transport_diag_ids), intent(inout)	IDs,
		type(diag_ctrl), intent(inout)	diag
	)

Register certain diagnostics related to transports.

Parameters

[in]	time	current model time
[in]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	ids	A structure with the diagnostic IDs.
[in,out]	diag	regulates diagnostic output

Definition at line 1807 of file MOM_diagnostics.F90.

  type(time_type),          intent(in)    :: Time  !< current model time
  type(ocean_grid_type),    intent(in)    :: G     !< ocean grid structure
  type(verticalGrid_type),  intent(in)    :: GV    !< ocean vertical grid structure
  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
  type(transport_diag_IDs), intent(inout) :: IDs   !< A structure with the diagnostic IDs.
  type(diag_ctrl),          intent(inout) :: diag  !< regulates diagnostic output
 
  real :: H_convert
  character(len=48) :: thickness_units
 
  thickness_units = get_thickness_units(gv)
  if (gv%Boussinesq) then
    h_convert = gv%H_to_m
  else
    h_convert = gv%H_to_kg_m2
  endif
 
  ! Diagnostics related to tracer and mass transport
  ids%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, time, &
      'Accumulated zonal thickness fluxes to advect tracers', 'kg', &
      y_cell_method='sum', v_extensive=.true., conversion=h_convert*us%L_to_m**2)
  ids%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, time, &
      'Accumulated meridional thickness fluxes to advect tracers', 'kg', &
      x_cell_method='sum', v_extensive=.true., conversion=h_convert*us%L_to_m**2)
  ids%id_umo = register_diag_field('ocean_model', 'umo', &
      diag%axesCuL, time, 'Ocean Mass X Transport', 'kg s-1', &
      standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.)
  ids%id_vmo = register_diag_field('ocean_model', 'vmo', &
      diag%axesCvL, time, 'Ocean Mass Y Transport', 'kg s-1', &
      standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.)
  ids%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', &
      diag%axesCu1, time, 'Ocean Mass X Transport Vertical Sum', 'kg s-1', &
      standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum')
  ids%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', &
      diag%axesCv1, time, 'Ocean Mass Y Transport Vertical Sum', 'kg s-1', &
      standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')
  ids%id_dynamics_h = register_diag_field('ocean_model','dynamics_h',  &
      diag%axesTl, time, 'Change in layer thicknesses due to horizontal dynamics', &
      'm s-1', v_extensive=.true., conversion=gv%H_to_m)
  ids%id_dynamics_h_tendency = register_diag_field('ocean_model','dynamics_h_tendency',  &
      diag%axesTl, time, 'Change in layer thicknesses due to horizontal dynamics', &
      'm s-1', v_extensive=.true., conversion=gv%H_to_m*us%s_to_T)
 

References mom_verticalgrid::get_thickness_units().

Here is the call 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]	mis	For "MOM Internal State" a set of pointers to the fields and accelerations making up ocean internal physical state.
[in,out]	adp	Structure pointing to accelerations in momentum equation.
[in,out]	cdp	Structure pointing to terms in continuity equation.
[in]	g	The ocean's grid structure.
	cs	Pointer to the control structure for this module.

Definition at line 2029 of file MOM_diagnostics.F90.

  type(ocean_internal_state), intent(in)    :: MIS !< For "MOM Internal State" a set of pointers to
                                                   !! the fields and accelerations making up ocean
                                                   !! internal physical state.
  type(accel_diag_ptrs),      intent(inout) :: ADp !< Structure pointing to accelerations in
                                                   !! momentum equation.
  type(cont_diag_ptrs),       intent(inout) :: CDp !< Structure pointing to terms in continuity
                                                   !! equation.
  type(ocean_grid_type),      intent(in)    :: G   !< The ocean's grid structure.
  type(diagnostics_CS),       pointer       :: CS  !< Pointer to the control structure for this
                                                   !! module.
 
! This subroutine sets up diagnostics upon which other diagnostics depend.
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  if (associated(cs%dKE_dt) .or. associated(cs%PE_to_KE) .or. &
      associated(cs%KE_CorAdv) .or. associated(cs%KE_adv) .or. &
      associated(cs%KE_visc) .or. associated(cs%KE_horvisc) .or. &
      associated(cs%KE_dia)) &
    call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
 
  if (associated(cs%dKE_dt)) then
    if (.not.associated(cs%du_dt)) then
      call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
      call register_time_deriv(lbound(mis%u), mis%u, cs%du_dt, cs)
    endif
    if (.not.associated(cs%dv_dt)) then
      call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
      call register_time_deriv(lbound(mis%v), mis%v, cs%dv_dt, cs)
    endif
    if (.not.associated(cs%dh_dt)) then
      call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
      call register_time_deriv(lbound(mis%h), mis%h, cs%dh_dt, cs)
    endif
  endif
 
  if (associated(cs%KE_adv)) then
    call safe_alloc_ptr(adp%gradKEu,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%gradKEv,isd,ied,jsdb,jedb,nz)
  endif
 
  if (associated(cs%KE_visc)) then
    call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
  endif
 
  if (associated(cs%KE_dia)) then
    call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
  endif
 
  if (associated(cs%uhGM_Rlay)) call safe_alloc_ptr(cdp%uhGM,isdb,iedb,jsd,jed,nz)
  if (associated(cs%vhGM_Rlay)) call safe_alloc_ptr(cdp%vhGM,isd,ied,jsdb,jedb,nz)
 

References register_time_deriv().

Referenced by mom_diagnostics_init().

Here is the call graph for this function:

Here is the caller graph for this function:

write_static_fields()

subroutine, public mom_diagnostics::write_static_fields	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(thermo_var_ptrs), intent(in)	tv,
		type(diag_ctrl), intent(inout), target	diag
	)

Offers the static fields in the ocean grid type for output via the diag_manager.

Parameters

[in]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	tv	A structure pointing to various thermodynamic variables
[in,out]	diag	regulates diagnostic output

Definition at line 1854 of file MOM_diagnostics.F90.

  type(ocean_grid_type),   intent(in)    :: G    !< ocean grid structure
  type(verticalGrid_type), intent(in)    :: GV   !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
 
  ! Local variables
  integer :: id
  logical :: use_temperature
 
  id = register_static_field('ocean_model', 'geolat', diag%axesT1, &
        'Latitude of tracer (T) points', 'degrees_north')
  if (id > 0) call post_data(id, g%geoLatT, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon', diag%axesT1, &
        'Longitude of tracer (T) points', 'degrees_east')
  if (id > 0) call post_data(id, g%geoLonT, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_c', diag%axesB1, &
        'Latitude of corner (Bu) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_c', diag%axesB1, &
        'Longitude of corner (Bu) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_v', diag%axesCv1, &
        'Latitude of meridional velocity (Cv) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_v', diag%axesCv1, &
        'Longitude of meridional velocity (Cv) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_u', diag%axesCu1, &
        'Latitude of zonal velocity (Cu) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_u', diag%axesCu1, &
        'Longitude of zonal velocity (Cu) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_t', diag%axesT1, &
        'Surface area of tracer (T) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) then
    call post_data(id, g%areaT, diag, .true.)
    call diag_register_area_ids(diag, id_area_t=id)
  endif
 
  id = register_static_field('ocean_model', 'area_u', diag%axesCu1, &
        'Surface area of x-direction flow (U) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_cu', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_v', diag%axesCv1, &
        'Surface area of y-direction flow (V) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_cv', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_q', diag%axesB1, &
        'Surface area of B-grid flow (Q) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_bu', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'depth_ocean', diag%axesT1,  &
        'Depth of the ocean at tracer points', 'm',                      &
        standard_name='sea_floor_depth_below_geoid',                     &
        cmor_field_name='deptho', cmor_long_name='Sea Floor Depth',      &
        cmor_standard_name='sea_floor_depth_below_geoid', area=diag%axesT1%id_area, &
        x_cell_method='mean', y_cell_method='mean', area_cell_method='mean', &
        conversion=us%Z_to_m)
  if (id > 0) call post_data(id, g%bathyT, diag, .true., mask=g%mask2dT)
 
  id = register_static_field('ocean_model', 'wet', diag%axesT1, &
        '0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area)
  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_c', diag%axesB1, &
        '0 if land, 1 if ocean at corner (Bu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_u', diag%axesCu1, &
        '0 if land, 1 if ocean at zonal velocity (Cu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_v', diag%axesCv1, &
        '0 if land, 1 if ocean at meridional velocity (Cv) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'Coriolis', diag%axesB1, &
        'Coriolis parameter at corner (Bu) points', 's-1', interp_method='none', conversion=us%s_to_T)
  if (id > 0) call post_data(id, g%CoriolisBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxt', diag%axesT1, &
        'Delta(x) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxT, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyt', diag%axesT1, &
        'Delta(y) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyT, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCu', diag%axesCu1, &
        'Delta(x) at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCu', diag%axesCu1, &
        'Delta(y) at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCv', diag%axesCv1, &
        'Delta(x) at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCv', diag%axesCv1, &
        'Delta(y) at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCuo', diag%axesCu1, &
        'Open meridional grid spacing at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dy_Cu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCvo', diag%axesCv1, &
        'Open zonal grid spacing at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dx_Cv, diag, .true.)
 
  id = register_static_field('ocean_model', 'sin_rot', diag%axesT1, &
        'sine of the clockwise angle of the ocean grid north to true north', 'none')
  if (id > 0) call post_data(id, g%sin_rot, diag, .true.)
 
  id = register_static_field('ocean_model', 'cos_rot', diag%axesT1, &
        'cosine of the clockwise angle of the ocean grid north to true north', 'none')
  if (id > 0) call post_data(id, g%cos_rot, diag, .true.)
 
 
  ! This static diagnostic is from CF 1.8, and is the fraction of a cell
  ! covered by ocean, given as a percentage (poorly named).
  id = register_static_field('ocean_model', 'area_t_percent', diag%axesT1, &
        'Percentage of cell area covered by ocean', '%', &
        cmor_field_name='sftof', cmor_standard_name='SeaAreaFraction', &
        cmor_long_name='Sea Area Fraction', &
        x_cell_method='mean', y_cell_method='mean', area_cell_method='mean', &
        conversion=100.0)
  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)
 
 
  id = register_static_field('ocean_model','Rho_0', diag%axesNull, &
       'mean ocean density used with the Boussinesq approximation', &
       'kg m-3', cmor_field_name='rhozero', conversion=us%R_to_kg_m3, &
       cmor_standard_name='reference_sea_water_density_for_boussinesq_approximation', &
       cmor_long_name='reference sea water density for boussinesq approximation')
  if (id > 0) call post_data(id, gv%Rho0, diag, .true.)
 
  use_temperature = associated(tv%T)
  if (use_temperature) then
     id = register_static_field('ocean_model','C_p', diag%axesNull, &
          'heat capacity of sea water', 'J kg-1 K-1', cmor_field_name='cpocean', &
          cmor_standard_name='specific_heat_capacity_of_sea_water', &
          cmor_long_name='specific_heat_capacity_of_sea_water')
     if (id > 0) call post_data(id, tv%C_p, diag, .true.)
  endif