MOM6
Data Types | Functions/Subroutines
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 [R L2 T-2 ~> 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 361 of file MOM_PressureForce_Montgomery.F90.

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

◆ 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 890 of file MOM_PressureForce_Montgomery.F90.

890  type(pressureforce_mont_cs), pointer :: cs !< Control structure for Montgomery potential PGF
891  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 822 of file MOM_PressureForce_Montgomery.F90.

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

◆ 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 [R L2 T-2 ~> Pa].
[out]pbceThe baroclinic pressure anomaly in
[out]etaFree surface height [H ~> kg m-1].

Definition at line 64 of file MOM_PressureForce_Montgomery.F90.

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

◆ 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 [R ~> 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 607 of file MOM_PressureForce_Montgomery.F90.

607  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
608  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
609  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface height [Z ~> m].
610  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
611  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
612  real, intent(in) :: rho0 !< The "Boussinesq" ocean density [R ~> kg m-3].
613  real, intent(in) :: gfs_scale !< Ratio between gravity applied to top
614  !! interface and the gravitational acceleration of
615  !! the planet [nondim]. Usually this ratio is 1.
616  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
617  intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
618  !! to free surface height anomalies
619  !! [L2 T-2 H-1 ~> m s-2].
620  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
621  optional, intent(in) :: rho_star !< The layer densities (maybe compressibility
622  !! compensated), times g/rho_0 [L2 Z-1 T-2 ~> m s-2].
623 
624  ! Local variables
625  real :: ihtot(szi_(g)) ! The inverse of the sum of the layer thicknesses [H-1 ~> m-1 or m2 kg-1].
626  real :: press(szi_(g)) ! Interface pressure [R L2 T-2 ~> Pa].
627  real :: t_int(szi_(g)) ! Interface temperature [degC].
628  real :: s_int(szi_(g)) ! Interface salinity [ppt].
629  real :: dr_dt(szi_(g)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
630  real :: dr_ds(szi_(g)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
631  real :: rho_in_situ(szi_(g)) ! In-situ density at the top of a layer [R ~> kg m-3].
632  real :: g_rho0 ! A scaled version of g_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
633  real :: rho0xg ! g_Earth * Rho0 [kg s-2 m-1 Z-1 ~> kg s-2 m-2]
634  logical :: use_eos ! If true, density is calculated from T & S using
635  ! an equation of state.
636  real :: z_neglect ! A thickness that is so small it is usually lost
637  ! in roundoff and can be neglected [Z ~> m].
638  integer, dimension(2) :: eosdom ! The computational domain for the equation of state
639  integer :: isq, ieq, jsq, jeq, nz, i, j, k
640 
641  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
642  eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
643 
644  rho0xg = rho0 * 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  tv%eqn_of_state, eosdom)
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  tv%eqn_of_state, eosdom)
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 

◆ 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 [R L2 T-2 ~> 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 [R L2 T-2 ~> 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
719  !! layer due 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 [H-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 [R L2 T-2 ~> Pa].
739  logical :: use_eos ! If true, density is calculated from T & S using an equation of state.
740  integer, dimension(2) :: eosdom ! The computational domain for the 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  eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
745 
746  use_eos = associated(tv%eqn_of_state)
747 
748  dp_dh = gv%g_Earth * gv%H_to_RZ
749  dp_neglect = gv%g_Earth * gv%H_to_RZ * gv%H_subroundoff
750 
751  if (use_eos) then
752  if (present(alpha_star)) then
753  !$OMP parallel do default(shared)
754  do j=jsq,jeq+1
755  do i=isq,ieq+1
756  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
757  pbce(i,j,nz) = dp_dh * alpha_star(i,j,nz)
758  enddo
759  do k=nz-1,1,-1 ; do i=isq,ieq+1
760  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
761  (alpha_star(i,j,k) - alpha_star(i,j,k+1))
762  enddo ; enddo
763  enddo
764  else
765  !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
766  do j=jsq,jeq+1
767  call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, &
768  tv%eqn_of_state, eosdom)
769  do i=isq,ieq+1
770  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
771  pbce(i,j,nz) = dp_dh / (rho_in_situ(i))
772  enddo
773  do k=nz-1,1,-1
774  do i=isq,ieq+1
775  t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
776  s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
777  enddo
778  call calculate_density(t_int, s_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, eosdom)
779  call calculate_density_derivs(t_int, s_int, p(:,j,k+1), dr_dt, dr_ds, &
780  tv%eqn_of_state, eosdom)
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)) * dalpha_int(k+1)
802  enddo ; enddo
803  enddo
804  endif ! use_EOS
805 
806  if (gfs_scale < 1.0) then
807  ! Adjust the Montgomery potential to make this a reduced gravity model.
808  !$OMP parallel do default(shared)
809  do j=jsq,jeq+1 ; do i=isq,ieq+1
810  dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
811  enddo ; enddo
812  !$OMP parallel do default(shared)
813  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
814  pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
815  enddo ; enddo ; enddo
816  endif
817