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, 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, CS, p_atm, pbce, eta)
 Boussinesq Montgomery-potential form of pressure gradient. More...
 
subroutine, public set_pbce_bouss (e, tv, G, GV, g_Earth, 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, g_Earth, 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, 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(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 ie to ie, je to je range before this subroutine is called: h[ie+1] and h[je+1] and (if tvform_of_EOS is set) T[ie+1], S[ie+1], T[je+1], and S[je+1].

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

Definition at line 368 of file MOM_PressureForce_Montgomery.F90.

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

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

Definition at line 905 of file MOM_PressureForce_Montgomery.F90.

905  type(pressureforce_mont_cs), pointer :: cs
906  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(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]param_fileParameter file handles
[in,out]diagDiagnostics control structure
csMontgomery PGF control structure
tides_cspTides control structure

Definition at line 839 of file MOM_PressureForce_Montgomery.F90.

References mom_error_handler::mom_error().

Referenced by mom_pressureforce::pressureforce_init().

839  type(time_type), target, intent(in) :: time !< Current model time
840  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
841  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
842  type(param_file_type), intent(in) :: param_file !< Parameter file handles
843  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
844  type(pressureforce_mont_cs), pointer :: cs !< Montgomery PGF control structure
845  type(tidal_forcing_cs), optional, pointer :: tides_csp !< Tides control structure
846  ! Local variables
847  logical :: use_temperature, use_eos
848 ! This include declares and sets the variable "version".
849 #include "version_variable.h"
850  character(len=40) :: mdl ! This module's name.
851 
852  if (associated(cs)) then
853  call mom_error(warning, "PressureForce_init called with an associated "// &
854  "control structure.")
855  return
856  else ; allocate(cs) ; endif
857 
858  cs%diag => diag ; cs%Time => time
859  if (present(tides_csp)) then
860  if (associated(tides_csp)) cs%tides_CSp => tides_csp
861  endif
862 
863  mdl = "MOM_PressureForce_Mont"
864  call log_version(param_file, mdl, version, "")
865  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
866  "The mean ocean density used with BOUSSINESQ true to \n"//&
867  "calculate accelerations and the mass for conservation \n"//&
868  "properties, or with BOUSSINSEQ false to convert some \n"//&
869  "parameters from vertical units of m to kg m-2.", &
870  units="kg m-3", default=1035.0)
871  call get_param(param_file, mdl, "TIDES", cs%tides, &
872  "If true, apply tidal momentum forcing.", default=.false.)
873  call get_param(param_file, mdl, "USE_EOS", use_eos, default=.true., &
874  do_not_log=.true.) ! Input for diagnostic use only.
875 
876  if (use_eos) then
877  cs%id_PFu_bc = register_diag_field('ocean_model', 'PFu_bc', diag%axesCuL, time, &
878  'Density Gradient Zonal Pressure Force Accel.', "meter second-2")
879  cs%id_PFv_bc = register_diag_field('ocean_model', 'PFv_bc', diag%axesCvL, time, &
880  'Density Gradient Meridional Pressure Force Accel.', "meter second-2")
881  if (cs%id_PFu_bc > 0) then
882  call safe_alloc_ptr(cs%PFu_bc,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
883  cs%PFu_bc(:,:,:) = 0.0
884  endif
885  if (cs%id_PFv_bc > 0) then
886  call safe_alloc_ptr(cs%PFv_bc,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
887  cs%PFv_bc(:,:,:) = 0.0
888  endif
889  endif
890 
891  if (cs%tides) then
892  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
893  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter')
894  endif
895 
896  cs%GFS_scale = 1.0
897  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
898 
899  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
900 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller 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(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 ie to ie, je to je range before this subroutine is called: h[ie+1] and h[je+1] and and (if tvform_of_EOS is set) T[ie+1], S[ie+1], T[je+1], and S[je+1].

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

Definition at line 57 of file MOM_PressureForce_Montgomery.F90.

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().

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

Definition at line 619 of file MOM_PressureForce_Montgomery.F90.

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

619  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
620  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
621  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface height, in H.
622  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
623  real, intent(in) :: g_earth !< The gravitational acceleration, in m s-2.
624  real, intent(in) :: rho0 !< The "Boussinesq" ocean density, in kg m-3.
625  real, intent(in) :: gfs_scale !< Ratio between gravity applied to top interface
626  !! and the gravitational acceleration of the planet.
627  !! Usually this ratio is 1.
628  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
629  !! to free surface height anomalies, in m2 H-1 s-2.
630  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: rho_star !< The layer densities (maybe
631  !! compressibility compensated), times g/rho_0, in m s-2.
632  ! Local variables
633  real :: ihtot(szi_(g)) ! The inverse of the sum of the layer
634  ! thicknesses, in m-1.
635  real :: press(szi_(g)) ! Interface pressure, in Pa.
636  real :: t_int(szi_(g)) ! Interface temperature in C.
637  real :: s_int(szi_(g)) ! Interface salinity in PSU.
638  real :: dr_dt(szi_(g)) ! Partial derivatives of density with temperature
639  real :: dr_ds(szi_(g)) ! and salinity in kg m-3 K-1 and kg m-3 PSU-1.
640  real :: rho_in_situ(szi_(g)) !In-situ density at the top of a layer.
641  real :: g_rho0 ! g_Earth / Rho0 in m4 s-2 kg-1.
642  real :: rho0xg ! g_Earth * Rho0 in kg s-2 m-2.
643  logical :: use_eos ! If true, density is calculated from T & S using
644  ! an equation of state.
645  real :: h_neglect ! A thickness that is so small it is usually lost
646  ! in roundoff and can be neglected, in m.
647  integer :: isq, ieq, jsq, jeq, nz, i, j, k
648 
649  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
650 
651  rho0xg = rho0*g_earth
652  g_rho0 = g_earth/rho0
653  use_eos = associated(tv%eqn_of_state)
654  h_neglect = gv%H_subroundoff*gv%H_to_m
655 
656  if (use_eos) then
657  if (present(rho_star)) then
658 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,h_neglect,pbce,rho_star,&
659 !$OMP GFS_scale,GV) &
660 !$OMP private(Ihtot)
661  do j=jsq,jeq+1
662  do i=isq,ieq+1
663  ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * gv%m_to_H)
664  pbce(i,j,1) = gfs_scale * rho_star(i,j,1) * gv%H_to_m
665  enddo
666  do k=2,nz ; do i=isq,ieq+1
667  pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * &
668  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
669  enddo ; enddo
670  enddo ! end of j loop
671  else
672 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,tv,h_neglect,G_Rho0,Rho0xG,&
673 !$OMP pbce,GFS_scale,GV) &
674 !$OMP private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
675  do j=jsq,jeq+1
676  do i=isq,ieq+1
677  ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * gv%m_to_H)
678  press(i) = -rho0xg*e(i,j,1)
679  enddo
680  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
681  isq, ieq-isq+2, tv%eqn_of_state)
682  do i=isq,ieq+1
683  pbce(i,j,1) = g_rho0*(gfs_scale * rho_in_situ(i)) * gv%H_to_m
684  enddo
685  do k=2,nz
686  do i=isq,ieq+1
687  press(i) = -rho0xg*e(i,j,k)
688  t_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
689  s_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
690  enddo
691  call calculate_density_derivs(t_int, s_int, press, dr_dt, dr_ds, &
692  isq, ieq-isq+2, tv%eqn_of_state)
693  do i=isq,ieq+1
694  pbce(i,j,k) = pbce(i,j,k-1) + g_rho0 * &
695  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i)) * &
696  (dr_dt(i)*(tv%T(i,j,k)-tv%T(i,j,k-1)) + &
697  dr_ds(i)*(tv%S(i,j,k)-tv%S(i,j,k-1)))
698  enddo
699  enddo
700  enddo ! end of j loop
701  endif
702  else ! not use_EOS
703 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,GV,h_neglect,pbce) private(Ihtot)
704  do j=jsq,jeq+1
705  do i=isq,ieq+1
706  ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * gv%m_to_H)
707  pbce(i,j,1) = gv%g_prime(1) * gv%H_to_m
708  enddo
709  do k=2,nz ; do i=isq,ieq+1
710  pbce(i,j,k) = pbce(i,j,k-1) + &
711  gv%g_prime(k) * ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
712  enddo ; enddo
713  enddo ! end of j loop
714  endif ! use_EOS
715 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
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,
real, intent(in)  g_Earth,
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, in Pa.
[in]tvThermodynamic variables
[in]g_earthThe gravitational acceleration, in m s-2.
[in]gfs_scaleRatio between gravity applied to top interface and the gravitational acceleration of the planet. Usually this ratio is 1.
[out]pbceThe baroclinic pressure anomaly in each layer due to free surface height anomalies, in m2 H-1 s-2.
[in]alpha_starThe layer specific volumes (maybe compressibility compensated), in m3 kg-1.

Definition at line 721 of file MOM_PressureForce_Montgomery.F90.

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

721  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
722  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
723  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: p !< Interface pressures, in Pa.
724  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
725  real, intent(in) :: g_earth !< The gravitational acceleration, in m s-2.
726  real, intent(in) :: gfs_scale !< Ratio between gravity applied to top interface
727  !! and the gravitational acceleration of the planet.
728  !! Usually this ratio is 1.
729  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
730  !! to free surface height anomalies, in m2 H-1 s-2.
731  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: alpha_star !< The layer specific volumes
732  !! (maybe compressibility compensated), in m3 kg-1.
733  ! Local variables
734  real, dimension(SZI_(G),SZJ_(G)) :: &
735  dpbce, & ! A barotropic correction to the pbce to enable the use of
736  ! a reduced gravity form of the equations, in m4 s-2 kg-1.
737  c_htot ! dP_dH divided by the total ocean pressure, m2 kg-1.
738  real :: t_int(szi_(g)) ! Interface temperature in C.
739  real :: s_int(szi_(g)) ! Interface salinity in PSU.
740  real :: dr_dt(szi_(g)) ! Partial derivatives of density with temperature
741  real :: dr_ds(szi_(g)) ! and salinity in kg m-3 K-1 and kg m-3 PSU-1.
742  real :: rho_in_situ(szi_(g)) !In-situ density at an interface, in kg m-3.
743  real :: alpha_lay(szk_(g)) ! The specific volume of each layer, in kg m-3.
744  real :: dalpha_int(szk_(g)+1) ! The change in specific volume across each
745  ! interface, in kg m-3.
746  real :: dp_dh ! A factor that converts from thickness to pressure,
747  ! usually in Pa m2 kg-1.
748  real :: dp_neglect ! A thickness that is so small it is usually lost
749  ! in roundoff and can be neglected, in Pa.
750  logical :: use_eos ! If true, density is calculated from T & S using
751  ! an equation of state.
752  integer :: isq, ieq, jsq, jeq, nz, i, j, k
753 
754  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
755 
756  use_eos = associated(tv%eqn_of_state)
757 
758  dp_dh = g_earth * gv%H_to_kg_m2
759  dp_neglect = dp_dh * gv%H_subroundoff
760 
761  do k=1,nz ; alpha_lay(k) = 1.0 / gv%Rlay(k) ; enddo
762  do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ; enddo
763 
764  if (use_eos) then
765  if (present(alpha_star)) then
766 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,C_htot,dP_dH,p,dp_neglect, &
767 !$OMP pbce,alpha_star)
768  do j=jsq,jeq+1
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 * alpha_star(i,j,nz)
772  enddo
773  do k=nz-1,1,-1 ; do i=isq,ieq+1
774  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
775  (alpha_star(i,j,k) - alpha_star(i,j,k+1))
776  enddo ; enddo
777  enddo
778  else
779 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,tv,p,C_htot, &
780 !$OMP dP_dH,dp_neglect,pbce) &
781 !$OMP private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
782  do j=jsq,jeq+1
783  call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), &
784  rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
785  do i=isq,ieq+1
786  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
787  pbce(i,j,nz) = dp_dh / rho_in_situ(i)
788  enddo
789  do k=nz-1,1,-1
790  do i=isq,ieq+1
791  t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
792  s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
793  enddo
794  call calculate_density(t_int, s_int, p(:,j,k+1), rho_in_situ, &
795  isq, ieq-isq+2, tv%eqn_of_state)
796  call calculate_density_derivs(t_int, s_int, p(:,j,k+1), dr_dt, dr_ds, &
797  isq, ieq-isq+2, tv%eqn_of_state)
798  do i=isq,ieq+1
799  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
800  ((dr_dt(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
801  dr_ds(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / rho_in_situ(i)**2)
802  enddo
803  enddo
804  enddo
805  endif
806  else ! not use_EOS
807 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,C_htot,dP_dH,p,dp_neglect, &
808 !$OMP pbce,alpha_Lay,dalpha_int)
809  do j=jsq,jeq+1
810  do i=isq,ieq+1
811  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
812  pbce(i,j,nz) = dp_dh * alpha_lay(nz)
813  enddo
814  do k=nz-1,1,-1 ; do i=isq,ieq+1
815  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
816  dalpha_int(k+1)
817  enddo ; enddo
818  enddo
819  endif ! use_EOS
820 
821  if (gfs_scale < 1.0) then
822  ! Adjust the Montgomery potential to make this a reduced gravity model.
823 !$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,dpbce,GFS_scale,pbce,nz)
824 !$OMP do
825  do j=jsq,jeq+1 ; do i=isq,ieq+1
826  dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
827  enddo ; enddo
828 !$OMP do
829  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
830  pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
831  enddo ; enddo ; enddo
832 !$OMP end parallel
833  endif
834 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the caller graph for this function: