MOM6
mom_pressureforce_mont Module Reference

Detailed Description

Provides the Montgomery potential form of pressure gradient.

Provides the Boussunesq and non-Boussinesq forms of the horizontal accelerations due to pressure gradients using the Montgomery potential. A second-order accurate, centered scheme is used. If a split time stepping scheme is used, the vertical decomposition into barotropic and baroclinic contributions described by Hallberg (J Comp Phys 1997) is used. With a nonlinear equation of state, compressibility is added along the lines proposed by Sun et al. (JPO 1999), but with compressibility coefficients based on a fit to a user-provided reference profile.

Data Types

type  pressureforce_mont_cs
 Control structure for the Montgomery potential form of pressure gradient. More...
 

Functions/Subroutines

subroutine, public pressureforce_mont_nonbouss (h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
 Non-Boussinesq Montgomery-potential form of pressure gradient. More...
 
subroutine, public pressureforce_mont_bouss (h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
 Boussinesq Montgomery-potential form of pressure gradient. More...
 
subroutine, public set_pbce_bouss (e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
 Determines the partial derivative of the acceleration due to pressure forces with the free surface height. More...
 
subroutine, public set_pbce_nonbouss (p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
 Determines the partial derivative of the acceleration due to pressure forces with the column mass. More...
 
subroutine, public pressureforce_mont_init (Time, G, GV, US, param_file, diag, CS, tides_CSp)
 Initialize the Montgomery-potential form of PGF control structure. More...
 
subroutine, public pressureforce_mont_end (CS)
 Deallocates the Montgomery-potential form of PGF control structure. More...
 

Function/Subroutine Documentation

◆ pressureforce_mont_bouss()

subroutine, public mom_pressureforce_mont::pressureforce_mont_bouss ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)  PFu,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  PFv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(pressureforce_mont_cs), pointer  CS,
real, dimension(:,:), optional, pointer  p_atm,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional  pbce,
real, dimension(szi_(g),szj_(g)), intent(out), optional  eta 
)

Boussinesq Montgomery-potential form of pressure gradient.

Determines the acceleration due to pressure forces.

To work, the following fields must be set outside of the usual (is:ie,js:je) range before this subroutine is called: h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).

Parameters
[in]gOcean grid structure.
[in]gvVertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m].
[in]tvThermodynamic variables.
[out]pfuZonal acceleration due to pressure gradients (equal to -dM/dx) [L T-2 ~> m s-2].
[out]pfvMeridional acceleration due to pressure gradients (equal to -dM/dy) [L T-2 ~> m s2].
csControl structure for Montgomery potential PGF
p_atmThe pressure at the ice-ocean or atmosphere-ocean [Pa].
[out]pbceThe baroclinic pressure anomaly in each layer due to free surface height anomalies [L2 T-2 H-1 ~> m s-2].
[out]etaFree surface height [H ~> m].

Definition at line 366 of file MOM_PressureForce_Montgomery.F90.

366  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure.
367  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
368  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
369  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m].
370  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables.
371  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: PFu !< Zonal acceleration due to pressure gradients
372  !! (equal to -dM/dx) [L T-2 ~> m s-2].
373  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients
374  !! (equal to -dM/dy) [L T-2 ~> m s2].
375  type(PressureForce_Mont_CS), pointer :: CS !< Control structure for Montgomery potential PGF
376  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or
377  !! atmosphere-ocean [Pa].
378  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure anomaly in
379  !! each layer due to free surface height anomalies
380  !! [L2 T-2 H-1 ~> m s-2].
381  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> m].
382  ! Local variables
383  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
384  M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2].
385  rho_star ! In-situ density divided by the derivative with depth of the
386  ! corrected e times (G_Earth/Rho0) [m2 Z-1 s-2 ~> m s-2].
387  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in m.
388  ! e may be adjusted (with a nonlinear equation of state) so that
389  ! its derivative compensates for the adiabatic compressibility
390  ! in seawater, but e will still be close to the interface depth.
391  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
392  T_tmp, & ! Temporary array of temperatures where layers that are lighter
393  ! than the mixed layer have the mixed layer's properties [degC].
394  s_tmp ! Temporary array of salinities where layers that are lighter
395  ! than the mixed layer have the mixed layer's properties [ppt].
396 
397  real :: Rho_cv_BL(SZI_(G)) ! The coordinate potential density in
398  ! the deepest variable density near-surface layer [R ~> kg m-3].
399  real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation
400  ! for compressibility [Z ~> m].
401  real :: e_tidal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal
402  ! forces from astronomical sources and self-
403  ! attraction and loading, in depth units [Z ~> m].
404  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
405  ! density [Pa] (usually 2e7 Pa = 2000 dbar).
406  real :: I_Rho0 ! 1/Rho0 [m3 kg-1].
407  real :: G_Rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
408  real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer
409  ! compensated density gradients [L T-2 ~> m s-2]
410  real :: h_neglect ! A thickness that is so small it is usually lost
411  ! in roundoff and can be neglected [Z ~> m].
412  logical :: use_p_atm ! If true, use the atmospheric pressure.
413  logical :: use_EOS ! If true, density is calculated from T & S using
414  ! an equation of state.
415  logical :: is_split ! A flag indicating whether the pressure
416  ! gradient terms are to be split into
417  ! barotropic and baroclinic pieces.
418  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
419  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
420  integer :: i, j, k
421 
422  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
423  nkmb=gv%nk_rho_varies
424  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
425 
426  use_p_atm = .false.
427  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
428  is_split = .false. ; if (present(pbce)) is_split = .true.
429  use_eos = associated(tv%eqn_of_state)
430 
431  if (.not.associated(cs)) call mom_error(fatal, &
432  "MOM_PressureForce_Mont: Module must be initialized before it is used.")
433  if (use_eos) then
434  if (query_compressible(tv%eqn_of_state)) call mom_error(fatal, &
435  "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
436  "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
437  endif
438 
439  h_neglect = gv%H_subroundoff * gv%H_to_Z
440  i_rho0 = 1.0/cs%Rho0
441  g_rho0 = gv%g_Earth / gv%Rho0
442 
443  if (cs%tides) then
444  ! Determine the surface height anomaly for calculating self attraction
445  ! and loading. This should really be based on bottom pressure anomalies,
446  ! but that is not yet implemented, and the current form is correct for
447  ! barotropic tides.
448  !$OMP parallel do default(shared)
449  do j=jsq,jeq+1
450  do i=isq,ieq+1 ; e(i,j,1) = -g%bathyT(i,j) ; enddo
451  do k=1,nz ; do i=isq,ieq+1
452  e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
453  enddo ; enddo
454  enddo
455  call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
456  endif
457 
458 ! Here layer interface heights, e, are calculated.
459  if (cs%tides) then
460  !$OMP parallel do default(shared)
461  do j=jsq,jeq+1 ; do i=isq,ieq+1
462  e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
463  enddo ; enddo
464  else
465  !$OMP parallel do default(shared)
466  do j=jsq,jeq+1 ; do i=isq,ieq+1
467  e(i,j,nz+1) = -g%bathyT(i,j)
468  enddo ; enddo
469  endif
470  !$OMP parallel do default(shared)
471  do j=jsq,jeq+1 ; do k=nz,1,-1 ; do i=isq,ieq+1
472  e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
473  enddo ; enddo ; enddo
474 
475  if (use_eos) then
476 ! Calculate in-situ densities (rho_star).
477 
478 ! With a bulk mixed layer, replace the T & S of any layers that are
479 ! lighter than the the buffer layer with the properties of the buffer
480 ! layer. These layers will be massless anyway, and it avoids any
481 ! formal calculations with hydrostatically unstable profiles.
482 
483  if (nkmb>0) then
484  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
485  tv_tmp%eqn_of_state => tv%eqn_of_state
486 
487  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
488  !$OMP parallel do default(shared) private(Rho_cv_BL)
489  do j=jsq,jeq+1
490  do k=1,nkmb ; do i=isq,ieq+1
491  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
492  enddo ; enddo
493  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
494  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
495 
496  do k=nkmb+1,nz ; do i=isq,ieq+1
497  if (gv%Rlay(k) < rho_cv_bl(i)) then
498  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
499  else
500  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
501  endif
502  enddo ; enddo
503  enddo
504  else
505  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
506  tv_tmp%eqn_of_state => tv%eqn_of_state
507  do i=isq,ieq+1 ; p_ref(i) = 0.0 ; enddo
508  endif
509 
510  ! This no longer includes any pressure dependency, since this routine
511  ! will come down with a fatal error if there is any compressibility.
512  !$OMP parallel do default(shared)
513  do k=1,nz+1 ; do j=jsq,jeq+1
514  call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
515  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R*g_rho0)
516  enddo ; enddo
517  endif ! use_EOS
518 
519 ! Here the layer Montgomery potentials, M, are calculated.
520  if (use_eos) then
521  !$OMP parallel do default(shared)
522  do j=jsq,jeq+1
523  do i=isq,ieq+1
524  m(i,j,1) = cs%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
525  if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
526  enddo
527  do k=2,nz ; do i=isq,ieq+1
528  m(i,j,k) = m(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,k)
529  enddo ; enddo
530  enddo
531  else ! not use_EOS
532  !$OMP parallel do default(shared)
533  do j=jsq,jeq+1
534  do i=isq,ieq+1
535  m(i,j,1) = gv%g_prime(1) * e(i,j,1)
536  if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
537  enddo
538  do k=2,nz ; do i=isq,ieq+1
539  m(i,j,k) = m(i,j,k-1) + gv%g_prime(k) * e(i,j,k)
540  enddo ; enddo
541  enddo
542  endif ! use_EOS
543 
544  if (present(pbce)) then
545  call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce, rho_star)
546  endif
547 
548 ! Calculate the pressure force. On a Cartesian grid,
549 ! PFu = - dM/dx and PFv = - dM/dy.
550  if (use_eos) then
551  !$OMP parallel do default(shared) private(h_star,PFu_bc,PFv_bc)
552  do k=1,nz
553  do j=jsq,jeq+1 ; do i=isq,ieq+1
554  h_star(i,j) = (e(i,j,k) - e(i,j,k+1)) + h_neglect
555  enddo ; enddo
556  do j=js,je ; do i=isq,ieq
557  pfu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (g%IdxCu(i,j) * &
558  ((h_star(i,j) * h_star(i+1,j) - (e(i,j,k) * h_star(i+1,j) + &
559  e(i+1,j,k) * h_star(i,j))) / (h_star(i,j) + h_star(i+1,j))))
560  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
561  if (associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
562  enddo ; enddo
563  do j=jsq,jeq ; do i=is,ie
564  pfv_bc = -1.0*(rho_star(i,j+1,k) - rho_star(i,j,k)) * (g%IdyCv(i,j) * &
565  ((h_star(i,j) * h_star(i,j+1) - (e(i,j,k) * h_star(i,j+1) + &
566  e(i,j+1,k) * h_star(i,j))) / (h_star(i,j) + h_star(i,j+1))))
567  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
568  if (associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
569  enddo ; enddo
570  enddo ! k-loop
571  else ! .not. use_EOS
572  !$OMP parallel do default(shared)
573  do k=1,nz
574  do j=js,je ; do i=isq,ieq
575  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
576  enddo ; enddo
577  do j=jsq,jeq ; do i=is,ie
578  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
579  enddo ; enddo
580  enddo
581  endif ! use_EOS
582 
583  if (present(eta)) then
584  if (cs%tides) then
585  ! eta is the sea surface height relative to a time-invariant geoid, for
586  ! comparison with what is used for eta in btstep. See how e was calculated
587  ! about 200 lines above.
588  !$OMP parallel do default(shared)
589  do j=jsq,jeq+1 ; do i=isq,ieq+1
590  eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
591  enddo ; enddo
592  else
593  !$OMP parallel do default(shared)
594  do j=jsq,jeq+1 ; do i=isq,ieq+1
595  eta(i,j) = e(i,j,1)*gv%Z_to_H
596  enddo ; enddo
597  endif
598  endif
599 
600  if (cs%id_PFu_bc>0) call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
601  if (cs%id_PFv_bc>0) call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
602  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
603 

References mom_tidal_forcing::calc_tidal_forcing(), mom_error_handler::mom_error(), mom_eos::query_compressible(), and set_pbce_bouss().

Here is the call graph for this function:

◆ pressureforce_mont_end()

subroutine, public mom_pressureforce_mont::pressureforce_mont_end ( type(pressureforce_mont_cs), pointer  CS)

Deallocates the Montgomery-potential form of PGF control structure.

Parameters
csControl structure for Montgomery potential PGF

Definition at line 891 of file MOM_PressureForce_Montgomery.F90.

891  type(PressureForce_Mont_CS), pointer :: CS !< Control structure for Montgomery potential PGF
892  if (associated(cs)) deallocate(cs)

◆ pressureforce_mont_init()

subroutine, public mom_pressureforce_mont::pressureforce_mont_init ( type(time_type), intent(in), target  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(pressureforce_mont_cs), pointer  CS,
type(tidal_forcing_cs), optional, pointer  tides_CSp 
)

Initialize the Montgomery-potential form of PGF control structure.

Parameters
[in]timeCurrent model time
[in]gocean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileParameter file handles
[in,out]diagDiagnostics control structure
csMontgomery PGF control structure
tides_cspTides control structure

Definition at line 823 of file MOM_PressureForce_Montgomery.F90.

823  type(time_type), target, intent(in) :: Time !< Current model time
824  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
825  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
826  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
827  type(param_file_type), intent(in) :: param_file !< Parameter file handles
828  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
829  type(PressureForce_Mont_CS), pointer :: CS !< Montgomery PGF control structure
830  type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tides control structure
831 
832  ! Local variables
833  logical :: use_temperature, use_EOS
834  ! This include declares and sets the variable "version".
835 # include "version_variable.h"
836  character(len=40) :: mdl ! This module's name.
837 
838  if (associated(cs)) then
839  call mom_error(warning, "PressureForce_init called with an associated "// &
840  "control structure.")
841  return
842  else ; allocate(cs) ; endif
843 
844  cs%diag => diag ; cs%Time => time
845  if (present(tides_csp)) then
846  if (associated(tides_csp)) cs%tides_CSp => tides_csp
847  endif
848 
849  mdl = "MOM_PressureForce_Mont"
850  call log_version(param_file, mdl, version, "")
851  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
852  "The mean ocean density used with BOUSSINESQ true to "//&
853  "calculate accelerations and the mass for conservation "//&
854  "properties, or with BOUSSINSEQ false to convert some "//&
855  "parameters from vertical units of m to kg m-2.", &
856  units="kg m-3", default=1035.0)
857  call get_param(param_file, mdl, "TIDES", cs%tides, &
858  "If true, apply tidal momentum forcing.", default=.false.)
859  call get_param(param_file, mdl, "USE_EOS", use_eos, default=.true., &
860  do_not_log=.true.) ! Input for diagnostic use only.
861 
862  if (use_eos) then
863  cs%id_PFu_bc = register_diag_field('ocean_model', 'PFu_bc', diag%axesCuL, time, &
864  'Density Gradient Zonal Pressure Force Accel.', "meter second-2", conversion=us%L_T2_to_m_s2)
865  cs%id_PFv_bc = register_diag_field('ocean_model', 'PFv_bc', diag%axesCvL, time, &
866  'Density Gradient Meridional Pressure Force Accel.', "meter second-2", conversion=us%L_T2_to_m_s2)
867  if (cs%id_PFu_bc > 0) then
868  call safe_alloc_ptr(cs%PFu_bc,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
869  cs%PFu_bc(:,:,:) = 0.0
870  endif
871  if (cs%id_PFv_bc > 0) then
872  call safe_alloc_ptr(cs%PFv_bc,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
873  cs%PFv_bc(:,:,:) = 0.0
874  endif
875  endif
876 
877  if (cs%tides) then
878  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
879  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
880  endif
881 
882  cs%GFS_scale = 1.0
883  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
884 
885  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
886 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ pressureforce_mont_nonbouss()

subroutine, public mom_pressureforce_mont::pressureforce_mont_nonbouss ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)  PFu,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  PFv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(pressureforce_mont_cs), pointer  CS,
real, dimension(:,:), optional, pointer  p_atm,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional  pbce,
real, dimension(szi_(g),szj_(g)), intent(out), optional  eta 
)

Non-Boussinesq Montgomery-potential form of pressure gradient.

Determines the acceleration due to pressure forces in a non-Boussinesq fluid using the compressibility compensated (if appropriate) Montgomery-potential form described in Hallberg (Ocean Mod., 2005).

To work, the following fields must be set outside of the usual (is:ie,js:je) range before this subroutine is called: h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).

Parameters
[in]gOcean grid structure.
[in]gvVertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thickness, [H ~> kg m-2].
[in]tvThermodynamic variables.
[out]pfuZonal acceleration due to pressure gradients (equal to -dM/dx) [L T-2 ~> m s-2].
[out]pfvMeridional acceleration due to pressure gradients (equal to -dM/dy) [L T-2 ~> m s-2].
csControl structure for Montgomery potential PGF
p_atmThe pressure at the ice-ocean or atmosphere-ocean [Pa].
[out]pbceThe baroclinic pressure anomaly in
[out]etaFree surface height [H ~> kg m-1].

Definition at line 65 of file MOM_PressureForce_Montgomery.F90.

65  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure.
66  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
67  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
68  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, [H ~> kg m-2].
69  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables.
70  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: PFu !< Zonal acceleration due to pressure gradients
71  !! (equal to -dM/dx) [L T-2 ~> m s-2].
72  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients
73  !! (equal to -dM/dy) [L T-2 ~> m s-2].
74  type(PressureForce_Mont_CS), pointer :: CS !< Control structure for Montgomery potential PGF
75  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or
76  !! atmosphere-ocean [Pa].
77  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78  optional, intent(out) :: pbce !< The baroclinic pressure anomaly in
79  !! each layer due to free surface height anomalies,
80  !! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2].
81  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> kg m-1].
82 
83  ! Local variables
84  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
85  M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2].
86  alpha_star, & ! Compression adjusted specific volume [R-1 ~> m3 kg-1].
87  dz_geo ! The change in geopotential across a layer [m2 s-2].
88  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure [Pa].
89  ! p may be adjusted (with a nonlinear equation of state) so that
90  ! its derivative compensates for the adiabatic compressibility
91  ! in seawater, but p will still be close to the pressure.
92  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
93  T_tmp, & ! Temporary array of temperatures where layers that are lighter
94  ! than the mixed layer have the mixed layer's properties [degC].
95  s_tmp ! Temporary array of salinities where layers that are lighter
96  ! than the mixed layer have the mixed layer's properties [ppt].
97 
98  real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the
99  ! deepest variable density near-surface layer [R ~> kg m-3].
100 
101  real, dimension(SZI_(G),SZJ_(G)) :: &
102  dM, & ! A barotropic correction to the Montgomery potentials to
103  ! enable the use of a reduced gravity form of the equations
104  ! [m2 s-2].
105  dp_star, & ! Layer thickness after compensation for compressibility [Pa].
106  ssh, & ! The sea surface height anomaly, in depth units [Z ~> m].
107  e_tidal, & ! Bottom geopotential anomaly due to tidal forces from
108  ! astronomical sources and self-attraction and loading [Z ~> m].
109  geopot_bot ! Bottom geopotential relative to time-mean sea level,
110  ! including any tidal contributions [L2 T-2 ~> m2 s-2].
111  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
112  ! density [Pa] (usually 2e7 Pa = 2000 dbar).
113  real :: rho_in_situ(SZI_(G)) !In-situ density of a layer [R ~> kg m-3].
114  real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer
115  ! compensated density gradients [L T-2 ~> m s-2]
116  real :: dp_neglect ! A thickness that is so small it is usually lost
117  ! in roundoff and can be neglected [Pa].
118  logical :: use_p_atm ! If true, use the atmospheric pressure.
119  logical :: use_EOS ! If true, density is calculated from T & S using
120  ! an equation of state.
121  logical :: is_split ! A flag indicating whether the pressure
122  ! gradient terms are to be split into
123  ! barotropic and baroclinic pieces.
124  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
125 
126  real :: I_gEarth ! The inverse of g_Earth [s2 Z m-2 ~> s2 m-1]
127 ! real :: dalpha
128  real :: Pa_to_p_dyn ! A conversion factor from Pa (= kg m-1 s-2) to the units of
129  ! dynamic pressure (R L2 T-2) [ R L2 T-2 m s2 kg-1 ~> nondim]
130  real :: Pa_to_H ! A factor to convert from Pa to the thicknesss units (H).
131  real :: alpha_Lay(SZK_(G)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
132  real :: dalpha_int(SZK_(G)+1) ! The change in specific volume across each
133  ! interface [R-1 ~> m3 kg-1].
134  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
135  integer :: i, j, k
136  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
137  nkmb=gv%nk_rho_varies
138  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
139 
140  use_p_atm = .false.
141  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
142  is_split = .false. ; if (present(pbce)) is_split = .true.
143  use_eos = associated(tv%eqn_of_state)
144 
145  if (.not.associated(cs)) call mom_error(fatal, &
146  "MOM_PressureForce_Mont: Module must be initialized before it is used.")
147  if (use_eos) then
148  if (query_compressible(tv%eqn_of_state)) call mom_error(fatal, &
149  "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
150  "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
151  endif
152 
153  pa_to_p_dyn = us%kg_m3_to_R * us%m_s_to_L_T**2
154  i_gearth = 1.0 / (us%L_T_to_m_s**2 * gv%g_Earth)
155  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
156  do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ; enddo
157  do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ; enddo
158 
159  if (use_p_atm) then
160  !$OMP parallel do default(shared)
161  do j=jsq,jeq+1 ; do i=isq,ieq+1 ; p(i,j,1) = p_atm(i,j) ; enddo ; enddo
162  else
163  !$OMP parallel do default(shared)
164  do j=jsq,jeq+1 ; do i=isq,ieq+1 ; p(i,j,1) = 0.0 ; enddo ; enddo
165  endif
166  !$OMP parallel do default(shared)
167  do j=jsq,jeq+1 ; do k=1,nz ; do i=isq,ieq+1
168  p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa * h(i,j,k)
169  enddo ; enddo ; enddo
170 
171  if (present(eta)) then
172  pa_to_h = 1.0 / gv%H_to_Pa
173  if (use_p_atm) then
174  !$OMP parallel do default(shared)
175  do j=jsq,jeq+1 ; do i=isq,ieq+1
176  eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.
177  enddo ; enddo
178  else
179  !$OMP parallel do default(shared)
180  do j=jsq,jeq+1 ; do i=isq,ieq+1
181  eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.
182  enddo ; enddo
183  endif
184  endif
185 
186  if (cs%tides) then
187  ! Determine the sea surface height anomalies, to enable the calculation
188  ! of self-attraction and loading.
189  !$OMP parallel do default(shared)
190  do j=jsq,jeq+1 ; do i=isq,ieq+1
191  ssh(i,j) = -g%bathyT(i,j)
192  enddo ; enddo
193  if (use_eos) then
194  !$OMP parallel do default(shared)
195  do k=1,nz
196  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
197  0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1)
198  enddo
199  !$OMP parallel do default(shared)
200  do j=jsq,jeq+1 ; do k=1,nz; do i=isq,ieq+1
201  ssh(i,j) = ssh(i,j) + i_gearth * dz_geo(i,j,k)
202  enddo ; enddo ; enddo
203  else
204  !$OMP parallel do default(shared)
205  do j=jsq,jeq+1 ; do k=1,nz ; do i=isq,ieq+1
206  ssh(i,j) = ssh(i,j) + gv%H_to_RZ * h(i,j,k) * alpha_lay(k)
207  enddo ; enddo ; enddo
208  endif
209 
210  call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
211  !$OMP parallel do default(shared)
212  do j=jsq,jeq+1 ; do i=isq,ieq+1
213  geopot_bot(i,j) = -gv%g_Earth*(e_tidal(i,j) + g%bathyT(i,j))
214  enddo ; enddo
215  else
216  !$OMP parallel do default(shared)
217  do j=jsq,jeq+1 ; do i=isq,ieq+1
218  geopot_bot(i,j) = -gv%g_Earth*g%bathyT(i,j)
219  enddo ; enddo
220  endif
221 
222  if (use_eos) then
223  ! Calculate in-situ specific volumes (alpha_star).
224 
225  ! With a bulk mixed layer, replace the T & S of any layers that are
226  ! lighter than the the buffer layer with the properties of the buffer
227  ! layer. These layers will be massless anyway, and it avoids any
228  ! formal calculations with hydrostatically unstable profiles.
229  if (nkmb>0) then
230  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
231  tv_tmp%eqn_of_state => tv%eqn_of_state
232  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
233  !$OMP parallel do default(shared) private(Rho_cv_BL)
234  do j=jsq,jeq+1
235  do k=1,nkmb ; do i=isq,ieq+1
236  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
237  enddo ; enddo
238  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
239  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
240  do k=nkmb+1,nz ; do i=isq,ieq+1
241  if (gv%Rlay(k) < rho_cv_bl(i)) then
242  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
243  else
244  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
245  endif
246  enddo ; enddo
247  enddo
248  else
249  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
250  tv_tmp%eqn_of_state => tv%eqn_of_state
251  do i=isq,ieq+1 ; p_ref(i) = 0 ; enddo
252  endif
253  !$OMP parallel do default(shared) private(rho_in_situ)
254  do k=1,nz ; do j=jsq,jeq+1
255  call calculate_density(tv_tmp%T(:,j,k),tv_tmp%S(:,j,k),p_ref, &
256  rho_in_situ,isq,ieq-isq+2,tv%eqn_of_state, scale=us%kg_m3_to_R)
257  do i=isq,ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo
258  enddo ; enddo
259  endif ! use_EOS
260 
261  if (use_eos) then
262  !$OMP parallel do default(shared)
263  do j=jsq,jeq+1
264  do i=isq,ieq+1
265  m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_star(i,j,nz)
266  enddo
267  do k=nz-1,1,-1 ; do i=isq,ieq+1
268  m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1))
269  enddo ; enddo
270  enddo
271  else ! not use_EOS
272  !$OMP parallel do default(shared)
273  do j=jsq,jeq+1
274  do i=isq,ieq+1
275  m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_lay(nz)
276  enddo
277  do k=nz-1,1,-1 ; do i=isq,ieq+1
278  m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * dalpha_int(k+1)
279  enddo ; enddo
280  enddo
281  endif ! use_EOS
282 
283  if (cs%GFS_scale < 1.0) then
284  ! Adjust the Montgomery potential to make this a reduced gravity model.
285  !$OMP parallel do default(shared)
286  do j=jsq,jeq+1 ; do i=isq,ieq+1
287  dm(i,j) = (cs%GFS_scale - 1.0) * m(i,j,1)
288  enddo ; enddo
289  !$OMP parallel do default(shared)
290  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
291  m(i,j,k) = m(i,j,k) + dm(i,j)
292  enddo ; enddo ; enddo
293 
294  ! Could instead do the following, to avoid taking small differences
295  ! of large numbers...
296 ! do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
297 ! M(i,j,1) = CS%GFS_scale * M(i,j,1)
298 ! enddo ; enddo
299 ! if (use_EOS) then
300 ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
301 ! M(i,j,k) = M(i,j,k-1) - Pa_to_p_dyn*p(i,j,K) * (alpha_star(i,j,k-1) - alpha_star(i,j,k))
302 ! enddo ; enddo ; enddo
303 ! else ! not use_EOS
304 ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
305 ! M(i,j,k) = M(i,j,k-1) - Pa_to_p_dyn*p(i,j,K) * dalpha_int(K)
306 ! enddo ; enddo ; enddo
307 ! endif ! use_EOS
308 
309  endif
310 
311  ! Note that ddM/dPb = alpha_star(i,j,1)
312  if (present(pbce)) then
313  call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce, alpha_star)
314  endif
315 
316 ! Calculate the pressure force. On a Cartesian grid,
317 ! PFu = - dM/dx and PFv = - dM/dy.
318  if (use_eos) then
319  !$OMP parallel do default(shared) private(dp_star,PFu_bc,PFv_bc)
320  do k=1,nz
321  do j=jsq,jeq+1 ; do i=isq,ieq+1
322  dp_star(i,j) = (p(i,j,k+1) - p(i,j,k)) + dp_neglect
323  enddo ; enddo
324  do j=js,je ; do i=isq,ieq
325  ! PFu_bc = p* grad alpha*
326  pfu_bc = (alpha_star(i+1,j,k) - alpha_star(i,j,k)) * (g%IdxCu(i,j) * pa_to_p_dyn * &
327  ((dp_star(i,j) * dp_star(i+1,j) + (p(i,j,k) * dp_star(i+1,j) + &
328  p(i+1,j,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i+1,j))))
329  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
330  if (associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
331  enddo ; enddo
332  do j=jsq,jeq ; do i=is,ie
333  pfv_bc = (alpha_star(i,j+1,k) - alpha_star(i,j,k)) * (g%IdyCv(i,j) * pa_to_p_dyn * &
334  ((dp_star(i,j) * dp_star(i,j+1) + (p(i,j,k) * dp_star(i,j+1) + &
335  p(i,j+1,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i,j+1))))
336  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
337  if (associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
338  enddo ; enddo
339  enddo ! k-loop
340  else ! .not. use_EOS
341  !$OMP parallel do default(shared)
342  do k=1,nz
343  do j=js,je ; do i=isq,ieq
344  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
345  enddo ; enddo
346  do j=jsq,jeq ; do i=is,ie
347  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
348  enddo ; enddo
349  enddo
350  endif ! use_EOS
351 
352  if (cs%id_PFu_bc>0) call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
353  if (cs%id_PFv_bc>0) call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
354  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
355 

References mom_tidal_forcing::calc_tidal_forcing(), mom_eos::int_specific_vol_dp(), mom_error_handler::mom_error(), mom_eos::query_compressible(), and set_pbce_nonbouss().

Here is the call graph for this function:

◆ set_pbce_bouss()

subroutine, public mom_pressureforce_mont::set_pbce_bouss ( real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  e,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, intent(in)  Rho0,
real, intent(in)  GFS_scale,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  pbce,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional  rho_star 
)

Determines the partial derivative of the acceleration due to pressure forces with the free surface height.

Parameters
[in]gOcean grid structure
[in]gvVertical grid structure
[in]eInterface height [Z ~> m].
[in]tvThermodynamic variables
[in]usA dimensional unit scaling type
[in]rho0The "Boussinesq" ocean density [kg m-3].
[in]gfs_scaleRatio between gravity applied to top interface and the gravitational acceleration of the planet [nondim]. Usually this ratio is 1.
[out]pbceThe baroclinic pressure anomaly in each layer due
[in]rho_starThe layer densities (maybe compressibility

Definition at line 609 of file MOM_PressureForce_Montgomery.F90.

609  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
610  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
611  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface height [Z ~> m].
612  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
613  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
614  real, intent(in) :: Rho0 !< The "Boussinesq" ocean density [kg m-3].
615  real, intent(in) :: GFS_scale !< Ratio between gravity applied to top
616  !! interface and the gravitational acceleration of
617  !! the planet [nondim]. Usually this ratio is 1.
618  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
619  intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
620  !! to free surface height anomalies
621  !! [L2 T-2 H-1 ~> m s-2].
622  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
623  optional, intent(in) :: rho_star !< The layer densities (maybe compressibility
624  !! compensated), times g/rho_0 [L2 Z-1 T-2 ~> m s-2].
625 
626  ! Local variables
627  real :: Ihtot(SZI_(G)) ! The inverse of the sum of the layer thicknesses [H-1 ~> m-1 or m2 kg-1].
628  real :: press(SZI_(G)) ! Interface pressure [Pa].
629  real :: T_int(SZI_(G)) ! Interface temperature [degC].
630  real :: S_int(SZI_(G)) ! Interface salinity [ppt].
631  real :: dR_dT(SZI_(G)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
632  real :: dR_dS(SZI_(G)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
633  real :: rho_in_situ(SZI_(G)) ! In-situ density at the top of a layer [R ~> kg m-3].
634  real :: G_Rho0 ! A scaled version of g_Earth / Rho0 [L2 m3 Z-1 T-2 kg-1 ~> m4 s-2 kg-1]
635  real :: Rho0xG ! g_Earth * Rho0 [kg s-2 m-1 Z-1 ~> kg s-2 m-2]
636  logical :: use_EOS ! If true, density is calculated from T & S using
637  ! an equation of state.
638  real :: z_neglect ! A thickness that is so small it is usually lost
639  ! in roundoff and can be neglected [Z ~> m].
640  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
641 
642  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
643 
644  rho0xg = rho0*us%L_T_to_m_s**2 * gv%g_Earth
645  g_rho0 = gv%g_Earth / gv%Rho0
646  use_eos = associated(tv%eqn_of_state)
647  z_neglect = gv%H_subroundoff*gv%H_to_Z
648 
649  if (use_eos) then
650  if (present(rho_star)) then
651  !$OMP parallel do default(shared) private(Ihtot)
652  do j=jsq,jeq+1
653  do i=isq,ieq+1
654  ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
655  pbce(i,j,1) = gfs_scale * rho_star(i,j,1) * gv%H_to_Z
656  enddo
657  do k=2,nz ; do i=isq,ieq+1
658  pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * &
659  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
660  enddo ; enddo
661  enddo ! end of j loop
662  else
663  !$OMP parallel do default(shared) private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
664  do j=jsq,jeq+1
665  do i=isq,ieq+1
666  ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
667  press(i) = -rho0xg*e(i,j,1)
668  enddo
669  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
670  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
671  do i=isq,ieq+1
672  pbce(i,j,1) = g_rho0*(gfs_scale * rho_in_situ(i)) * gv%H_to_Z
673  enddo
674  do k=2,nz
675  do i=isq,ieq+1
676  press(i) = -rho0xg*e(i,j,k)
677  t_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
678  s_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
679  enddo
680  call calculate_density_derivs(t_int, s_int, press, dr_dt, dr_ds, &
681  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
682  do i=isq,ieq+1
683  pbce(i,j,k) = pbce(i,j,k-1) + g_rho0 * &
684  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i)) * &
685  (dr_dt(i)*(tv%T(i,j,k)-tv%T(i,j,k-1)) + &
686  dr_ds(i)*(tv%S(i,j,k)-tv%S(i,j,k-1)))
687  enddo
688  enddo
689  enddo ! end of j loop
690  endif
691  else ! not use_EOS
692  !$OMP parallel do default(shared) private(Ihtot)
693  do j=jsq,jeq+1
694  do i=isq,ieq+1
695  ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
696  pbce(i,j,1) = gv%g_prime(1) * gv%H_to_Z
697  enddo
698  do k=2,nz ; do i=isq,ieq+1
699  pbce(i,j,k) = pbce(i,j,k-1) + &
700  (gv%g_prime(k)*gv%H_to_Z) * ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
701  enddo ; enddo
702  enddo ! end of j loop
703  endif ! use_EOS
704 

Referenced by mom_pressureforce_afv::pressureforce_afv_bouss(), mom_pressureforce_blk_afv::pressureforce_blk_afv_bouss(), and pressureforce_mont_bouss().

Here is the caller graph for this function:

◆ set_pbce_nonbouss()

subroutine, public mom_pressureforce_mont::set_pbce_nonbouss ( real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  p,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, intent(in)  GFS_scale,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  pbce,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional  alpha_star 
)

Determines the partial derivative of the acceleration due to pressure forces with the column mass.

Parameters
[in]gOcean grid structure
[in]gvVertical grid structure
[in]pInterface pressures [Pa].
[in]tvThermodynamic variables
[in]usA dimensional unit scaling type
[in]gfs_scaleRatio between gravity applied to top interface and the gravitational acceleration of the planet [nondim]. Usually this ratio is 1.
[out]pbceThe baroclinic pressure anomaly in each layer due to free surface height anomalies [L2 H-1 T-2 ~> m4 kg-1 s-2].
[in]alpha_starThe layer specific volumes (maybe compressibility compensated) [R-1 ~> m3 kg-1].

Definition at line 710 of file MOM_PressureForce_Montgomery.F90.

710  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
711  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
712  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: p !< Interface pressures [Pa].
713  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
714  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
715  real, intent(in) :: GFS_scale !< Ratio between gravity applied to top
716  !! interface and the gravitational acceleration of
717  !! the planet [nondim]. Usually this ratio is 1.
718  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
719  !! to free surface height anomalies
720  !! [L2 H-1 T-2 ~> m4 kg-1 s-2].
721  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: alpha_star !< The layer specific volumes
722  !! (maybe compressibility compensated) [R-1 ~> m3 kg-1].
723  ! Local variables
724  real, dimension(SZI_(G),SZJ_(G)) :: &
725  dpbce, & ! A barotropic correction to the pbce to enable the use of
726  ! a reduced gravity form of the equations [L2 H-1 T-2 ~> m4 kg-1 s-2].
727  c_htot ! dP_dH divided by the total ocean pressure [R L2 T-2 H-1 Pa-1 ~> m2 kg-1].
728  real :: T_int(SZI_(G)) ! Interface temperature [degC].
729  real :: S_int(SZI_(G)) ! Interface salinity [ppt].
730  real :: dR_dT(SZI_(G)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
731  real :: dR_dS(SZI_(G)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
732  real :: rho_in_situ(SZI_(G)) ! In-situ density at an interface [R ~> kg m-3].
733  real :: alpha_Lay(SZK_(G)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
734  real :: dalpha_int(SZK_(G)+1) ! The change in specific volume across each interface [R-1 ~> m3 kg-1].
735  real :: dP_dH ! A factor that converts from thickness to pressure times other dimensional
736  ! conversion factors [R L2 T-2 H-1 ~> Pa m2 kg-1].
737  real :: dp_neglect ! A thickness that is so small it is usually lost
738  ! in roundoff and can be neglected [Pa].
739  logical :: use_EOS ! If true, density is calculated from T & S using
740  ! an equation of state.
741  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
742 
743  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
744 
745  use_eos = associated(tv%eqn_of_state)
746 
747  dp_dh = gv%g_Earth * gv%H_to_RZ
748  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
749 
750  if (use_eos) then
751  if (present(alpha_star)) then
752  !$OMP parallel do default(shared)
753  do j=jsq,jeq+1
754  do i=isq,ieq+1
755  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
756  pbce(i,j,nz) = dp_dh * alpha_star(i,j,nz)
757  enddo
758  do k=nz-1,1,-1 ; do i=isq,ieq+1
759  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
760  (alpha_star(i,j,k) - alpha_star(i,j,k+1))
761  enddo ; enddo
762  enddo
763  else
764  !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
765  do j=jsq,jeq+1
766  call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), &
767  rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
768  do i=isq,ieq+1
769  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
770  pbce(i,j,nz) = dp_dh / (rho_in_situ(i))
771  enddo
772  do k=nz-1,1,-1
773  do i=isq,ieq+1
774  t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
775  s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
776  enddo
777  call calculate_density(t_int, s_int, p(:,j,k+1), rho_in_situ, &
778  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
779  call calculate_density_derivs(t_int, s_int, p(:,j,k+1), dr_dt, dr_ds, &
780  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
781  do i=isq,ieq+1
782  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
783  ((dr_dt(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
784  dr_ds(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / (rho_in_situ(i)**2))
785  enddo
786  enddo
787  enddo
788  endif
789  else ! not use_EOS
790 
791  do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ; enddo
792  do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ; enddo
793 
794  !$OMP parallel do default(shared)
795  do j=jsq,jeq+1
796  do i=isq,ieq+1
797  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
798  pbce(i,j,nz) = dp_dh * alpha_lay(nz)
799  enddo
800  do k=nz-1,1,-1 ; do i=isq,ieq+1
801  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
802  dalpha_int(k+1)
803  enddo ; enddo
804  enddo
805  endif ! use_EOS
806 
807  if (gfs_scale < 1.0) then
808  ! Adjust the Montgomery potential to make this a reduced gravity model.
809  !$OMP parallel do default(shared)
810  do j=jsq,jeq+1 ; do i=isq,ieq+1
811  dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
812  enddo ; enddo
813  !$OMP parallel do default(shared)
814  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
815  pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
816  enddo ; enddo ; enddo
817  endif
818 

Referenced by mom_pressureforce_afv::pressureforce_afv_nonbouss(), mom_pressureforce_blk_afv::pressureforce_blk_afv_nonbouss(), and pressureforce_mont_nonbouss().

Here is the caller graph for this function: