portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Provides the Montgomery potential form of pressure gradient | |
| 6 | module MOM_PressureForce_Mont | |
| 7 | ||
| 8 | use MOM_density_integrals, only : int_specific_vol_dp | |
| 9 | use MOM_diag_mediator, only : post_data, register_diag_field | |
| 10 | use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type | |
| 11 | use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe | |
| 12 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 13 | use MOM_grid, only : ocean_grid_type | |
| 14 | use MOM_self_attr_load, only : calc_SAL, SAL_CS | |
| 15 | use MOM_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_CS | |
| 16 | use MOM_unit_scaling, only : unit_scale_type | |
| 17 | use MOM_variables, only : thermo_var_ptrs | |
| 18 | use MOM_verticalGrid, only : verticalGrid_type | |
| 19 | use MOM_EOS, only : calculate_density, calculate_density_derivs | |
| 20 | use MOM_EOS, only : query_compressible | |
| 21 | ||
| 22 | implicit none ; private | |
| 23 | ||
| 24 | #include <MOM_memory.h> | |
| 25 | ||
| 26 | public PressureForce_Mont_Bouss, PressureForce_Mont_nonBouss, Set_pbce_Bouss | |
| 27 | public Set_pbce_nonBouss, PressureForce_Mont_init | |
| 28 | ||
| 29 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 30 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 31 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 32 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 33 | ||
| 34 | !> Control structure for the Montgomery potential form of pressure gradient | |
| 35 | type, public :: PressureForce_Mont_CS ; private | |
| 36 | logical :: initialized = .false. !< True if this control structure has been initialized. | |
| 37 | logical :: calculate_SAL !< If true, calculate self-attraction and loading. | |
| 38 | logical :: tides !< If true, apply tidal momentum forcing. | |
| 39 | real :: Rho0 !< The density used in the Boussinesq | |
| 40 | !! approximation [R ~> kg m-3]. | |
| 41 | real :: GFS_scale !< Ratio between gravity applied to top interface and the | |
| 42 | !! gravitational acceleration of the planet [nondim]. | |
| 43 | !! Usually this ratio is 1. | |
| 44 | type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. | |
| 45 | type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate | |
| 46 | !! the timing of diagnostic output. | |
| 47 | real, allocatable :: PFu_bc(:,:,:) !< Zonal accelerations due to pressure gradients | |
| 48 | !! deriving from density gradients within layers [L T-2 ~> m s-2]. | |
| 49 | real, allocatable :: PFv_bc(:,:,:) !< Meridional accelerations due to pressure gradients | |
| 50 | !! deriving from density gradients within layers [L T-2 ~> m s-2]. | |
| 51 | !>@{ Diagnostic IDs | |
| 52 | integer :: id_PFu_bc = -1, id_PFv_bc = -1, id_e_sal = -1 | |
| 53 | integer :: id_e_tide = -1, id_e_tide_eq = -1, id_e_tide_sal = -1 | |
| 54 | !>@} | |
| 55 | type(SAL_CS), pointer :: SAL_CSp => NULL() !< SAL control structure | |
| 56 | type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< The tidal forcing control structure | |
| 57 | end type PressureForce_Mont_CS | |
| 58 | ||
| 59 | contains | |
| 60 | ||
| 61 | !> \brief Non-Boussinesq Montgomery-potential form of pressure gradient | |
| 62 | !! | |
| 63 | !! Determines the acceleration due to pressure forces in a | |
| 64 | !! non-Boussinesq fluid using the compressibility compensated (if appropriate) | |
| 65 | !! Montgomery-potential form described in Hallberg (Ocean Mod., 2005). | |
| 66 | !! | |
| 67 | !! To work, the following fields must be set outside of the usual (is:ie,js:je) | |
| 68 | !! range before this subroutine is called: | |
| 69 | !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1). | |
| 70 | 0 | subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta) |
| 71 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. | |
| 72 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 73 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 74 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness, [H ~> kg m-2]. | |
| 75 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables. | |
| 76 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: PFu !< Zonal acceleration due to pressure gradients | |
| 77 | !! (equal to -dM/dx) [L T-2 ~> m s-2]. | |
| 78 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients | |
| 79 | !! (equal to -dM/dy) [L T-2 ~> m s-2]. | |
| 80 | type(PressureForce_Mont_CS), intent(inout) :: CS !< Control structure for Montgomery potential PGF | |
| 81 | real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean or | |
| 82 | !! atmosphere-ocean [R L2 T-2 ~> Pa]. | |
| 83 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 84 | optional, intent(out) :: pbce !< The baroclinic pressure anomaly in | |
| 85 | !! each layer due to free surface height anomalies, | |
| 86 | !! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2]. | |
| 87 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to calculate | |
| 88 | !! PFu and PFv [H ~> kg m-2]. | |
| 89 | ! Local variables | |
| 90 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 91 | 0 | M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2]. |
| 92 | 0 | alpha_star, & ! Compression adjusted specific volume [R-1 ~> m3 kg-1]. |
| 93 | 0 | dz_geo ! The change in geopotential across a layer [L2 T-2 ~> m2 s-2]. |
| 94 | 0 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa]. |
| 95 | ! p may be adjusted (with a nonlinear equation of state) so that | |
| 96 | ! its derivative compensates for the adiabatic compressibility | |
| 97 | ! in seawater, but p will still be close to the pressure. | |
| 98 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: & | |
| 99 | 0 | T_tmp, & ! Temporary array of temperatures where layers that are lighter |
| 100 | ! than the mixed layer have the mixed layer's properties [C ~> degC]. | |
| 101 | 0 | S_tmp ! Temporary array of salinities where layers that are lighter |
| 102 | ! than the mixed layer have the mixed layer's properties [S ~> ppt]. | |
| 103 | ||
| 104 | 0 | real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the |
| 105 | ! deepest variable density near-surface layer [R ~> kg m-3]. | |
| 106 | ||
| 107 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 108 | 0 | dM, & ! A barotropic correction to the Montgomery potentials to enable the use |
| 109 | ! of a reduced gravity form of the equations [L2 T-2 ~> m2 s-2]. | |
| 110 | 0 | dp_star, & ! Layer thickness after compensation for compressibility [R L2 T-2 ~> Pa]. |
| 111 | 0 | SSH, & ! The sea surface height anomaly, in depth units [Z ~> m]. |
| 112 | 0 | e_sal, & ! Bottom geopotential anomaly due to self-attraction and loading [Z ~> m]. |
| 113 | 0 | e_tide_eq, & ! Bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m]. |
| 114 | 0 | e_tide_sal, & ! Bottom geopotential anomaly due to harmonic self-attraction and loading |
| 115 | ! specific to tides [Z ~> m]. | |
| 116 | 0 | geopot_bot ! Bottom geopotential relative to a temporally fixed reference value, |
| 117 | ! including any tidal contributions [L2 T-2 ~> m2 s-2]. | |
| 118 | 0 | real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate |
| 119 | ! density [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). | |
| 120 | 0 | real :: rho_in_situ(SZI_(G)) !In-situ density of a layer [R ~> kg m-3]. |
| 121 | real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer | |
| 122 | ! compensated density gradients [L T-2 ~> m s-2] | |
| 123 | real :: dp_neglect ! A thickness that is so small it is usually lost | |
| 124 | ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. | |
| 125 | logical :: use_p_atm ! If true, use the atmospheric pressure. | |
| 126 | logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. | |
| 127 | logical :: is_split ! A flag indicating whether the pressure gradient terms are to be | |
| 128 | ! split into barotropic and baroclinic pieces. | |
| 129 | 0 | type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. |
| 130 | ||
| 131 | real :: I_gEarth ! The inverse of g_Earth [T2 Z L-2 ~> s2 m-1] | |
| 132 | ! real :: dalpha | |
| 133 | real :: Pa_to_H ! A factor to convert from R L2 T-2 to the thickness units (H) | |
| 134 | ! [H T2 R-1 L-2 ~> m2 s2 kg-1 or s2 m-1]. | |
| 135 | 0 | real :: alpha_Lay(SZK_(GV)) ! The specific volume of each layer [R-1 ~> m3 kg-1]. |
| 136 | 0 | real :: dalpha_int(SZK_(GV)+1) ! The change in specific volume across each |
| 137 | ! interface [R-1 ~> m3 kg-1]. | |
| 138 | integer, dimension(2) :: EOSdom ! The computational domain for the equation of state | |
| 139 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb | |
| 140 | integer :: i, j, k | |
| 141 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 142 | 0 | nkmb=GV%nk_rho_varies |
| 143 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 144 | 0 | EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) |
| 145 | ||
| 146 | 0 | use_p_atm = associated(p_atm) |
| 147 | 0 | is_split = present(pbce) |
| 148 | 0 | use_EOS = associated(tv%eqn_of_state) |
| 149 | ||
| 150 | 0 | if (.not.CS%initialized) call MOM_error(FATAL, & |
| 151 | 0 | "MOM_PressureForce_Mont: Module must be initialized before it is used.") |
| 152 | ||
| 153 | 0 | if (use_EOS) then |
| 154 | 0 | if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, & |
| 155 | "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//& | |
| 156 | 0 | "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.") |
| 157 | endif | |
| 158 | ||
| 159 | 0 | I_gEarth = 1.0 / GV%g_Earth |
| 160 | 0 | dp_neglect = GV%g_Earth * GV%H_to_RZ * GV%H_subroundoff |
| 161 | 0 | do k=1,nz ; alpha_Lay(k) = 1.0 / (GV%Rlay(k)) ; enddo |
| 162 | 0 | do k=2,nz ; dalpha_int(K) = alpha_Lay(k-1) - alpha_Lay(k) ; enddo |
| 163 | ||
| 164 | 0 | if (use_p_atm) then |
| 165 | !$OMP parallel do default(shared) | |
| 166 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; p(i,j,1) = p_atm(i,j) ; enddo ; enddo |
| 167 | else | |
| 168 | !$OMP parallel do default(shared) | |
| 169 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; p(i,j,1) = 0.0 ; enddo ; enddo |
| 170 | endif | |
| 171 | !$OMP parallel do default(shared) | |
| 172 | 0 | do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 |
| 173 | 0 | p(i,j,K+1) = p(i,j,K) + GV%g_Earth * GV%H_to_RZ * h(i,j,k) |
| 174 | enddo ; enddo ; enddo | |
| 175 | ||
| 176 | 0 | if (present(eta)) then |
| 177 | 0 | Pa_to_H = 1.0 / (GV%g_Earth * GV%H_to_RZ) |
| 178 | 0 | if (use_p_atm) then |
| 179 | !$OMP parallel do default(shared) | |
| 180 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 181 | 0 | eta(i,j) = (p(i,j,nz+1) - p_atm(i,j)) * Pa_to_H ! eta has the same units as h. |
| 182 | enddo ; enddo | |
| 183 | else | |
| 184 | !$OMP parallel do default(shared) | |
| 185 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 186 | 0 | eta(i,j) = p(i,j,nz+1) * Pa_to_H ! eta has the same units as h. |
| 187 | enddo ; enddo | |
| 188 | endif | |
| 189 | endif | |
| 190 | ||
| 191 | !$OMP parallel do default(shared) | |
| 192 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 193 | 0 | geopot_bot(i,j) = -GV%g_Earth * G%bathyT(i,j) |
| 194 | enddo ; enddo | |
| 195 | ||
| 196 | ! Calculate and add the self-attraction and loading geopotential anomaly. | |
| 197 | 0 | if (CS%calculate_SAL) then |
| 198 | ! Determine the sea surface height anomalies, to enable the calculation | |
| 199 | ! of self-attraction and loading. | |
| 200 | !$OMP parallel do default(shared) | |
| 201 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 202 | 0 | SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0) |
| 203 | enddo ; enddo | |
| 204 | 0 | if (use_EOS) then |
| 205 | !$OMP parallel do default(shared) | |
| 206 | 0 | do k=1,nz |
| 207 | call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), & | |
| 208 | 0 | 0.0, G%HI, tv%eqn_of_state, US, dz_geo(:,:,k), halo_size=1) |
| 209 | enddo | |
| 210 | !$OMP parallel do default(shared) | |
| 211 | 0 | do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 |
| 212 | 0 | SSH(i,j) = SSH(i,j) + I_gEarth * dz_geo(i,j,k) |
| 213 | enddo ; enddo ; enddo | |
| 214 | else | |
| 215 | !$OMP parallel do default(shared) | |
| 216 | 0 | do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 |
| 217 | 0 | SSH(i,j) = SSH(i,j) + GV%H_to_RZ * h(i,j,k) * alpha_Lay(k) |
| 218 | enddo ; enddo ; enddo | |
| 219 | endif | |
| 220 | ||
| 221 | 0 | call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) |
| 222 | !$OMP parallel do default(shared) | |
| 223 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 224 | 0 | geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*e_sal(i,j) |
| 225 | enddo ; enddo | |
| 226 | endif | |
| 227 | ||
| 228 | ! Calculate and add the tidal geopotential anomaly. | |
| 229 | 0 | if (CS%tides) then |
| 230 | 0 | call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) |
| 231 | !$OMP parallel do default(shared) | |
| 232 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 233 | 0 | geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*(e_tide_eq(i,j) + e_tide_sal(i,j)) |
| 234 | enddo ; enddo | |
| 235 | endif | |
| 236 | ||
| 237 | 0 | if (use_EOS) then |
| 238 | ! Calculate in-situ specific volumes (alpha_star). | |
| 239 | ||
| 240 | ! With a bulk mixed layer, replace the T & S of any layers that are | |
| 241 | ! lighter than the buffer layer with the properties of the buffer | |
| 242 | ! layer. These layers will be massless anyway, and it avoids any | |
| 243 | ! formal calculations with hydrostatically unstable profiles. | |
| 244 | 0 | if (nkmb>0) then |
| 245 | 0 | tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp |
| 246 | 0 | tv_tmp%eqn_of_state => tv%eqn_of_state |
| 247 | 0 | do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo |
| 248 | !$OMP parallel do default(shared) private(Rho_cv_BL) | |
| 249 | 0 | do j=Jsq,Jeq+1 |
| 250 | 0 | do k=1,nkmb ; do i=Isq,Ieq+1 |
| 251 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) |
| 252 | enddo ; enddo | |
| 253 | call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & | |
| 254 | 0 | tv%eqn_of_state, EOSdom) |
| 255 | 0 | do k=nkmb+1,nz ; do i=Isq,Ieq+1 |
| 256 | 0 | if (GV%Rlay(k) < Rho_cv_BL(i)) then |
| 257 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) |
| 258 | else | |
| 259 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) |
| 260 | endif | |
| 261 | enddo ; enddo | |
| 262 | enddo | |
| 263 | else | |
| 264 | 0 | tv_tmp%T => tv%T ; tv_tmp%S => tv%S |
| 265 | 0 | tv_tmp%eqn_of_state => tv%eqn_of_state |
| 266 | 0 | do i=Isq,Ieq+1 ; p_ref(i) = 0 ; enddo |
| 267 | endif | |
| 268 | !$OMP parallel do default(shared) private(rho_in_situ) | |
| 269 | 0 | do k=1,nz ; do j=Jsq,Jeq+1 |
| 270 | call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, & | |
| 271 | 0 | tv%eqn_of_state, EOSdom) |
| 272 | 0 | do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo |
| 273 | enddo ; enddo | |
| 274 | endif ! use_EOS | |
| 275 | ||
| 276 | 0 | if (use_EOS) then |
| 277 | !$OMP parallel do default(shared) | |
| 278 | 0 | do j=Jsq,Jeq+1 |
| 279 | 0 | do i=Isq,Ieq+1 |
| 280 | 0 | M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_star(i,j,nz) |
| 281 | enddo | |
| 282 | 0 | do k=nz-1,1,-1 ; do i=Isq,Ieq+1 |
| 283 | 0 | 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)) |
| 284 | enddo ; enddo | |
| 285 | enddo | |
| 286 | else ! not use_EOS | |
| 287 | !$OMP parallel do default(shared) | |
| 288 | 0 | do j=Jsq,Jeq+1 |
| 289 | 0 | do i=Isq,Ieq+1 |
| 290 | 0 | M(i,j,nz) = geopot_bot(i,j) + p(i,j,nz+1) * alpha_Lay(nz) |
| 291 | enddo | |
| 292 | 0 | do k=nz-1,1,-1 ; do i=Isq,Ieq+1 |
| 293 | 0 | M(i,j,k) = M(i,j,k+1) + p(i,j,K+1) * dalpha_int(K+1) |
| 294 | enddo ; enddo | |
| 295 | enddo | |
| 296 | endif ! use_EOS | |
| 297 | ||
| 298 | 0 | if (CS%GFS_scale < 1.0) then |
| 299 | ! Adjust the Montgomery potential to make this a reduced gravity model. | |
| 300 | !$OMP parallel do default(shared) | |
| 301 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 302 | 0 | dM(i,j) = (CS%GFS_scale - 1.0) * M(i,j,1) |
| 303 | enddo ; enddo | |
| 304 | !$OMP parallel do default(shared) | |
| 305 | 0 | do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 306 | 0 | M(i,j,k) = M(i,j,k) + dM(i,j) |
| 307 | enddo ; enddo ; enddo | |
| 308 | ||
| 309 | ! Could instead do the following, to avoid taking small differences | |
| 310 | ! of large numbers... | |
| 311 | ! do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 | |
| 312 | ! M(i,j,1) = CS%GFS_scale * M(i,j,1) | |
| 313 | ! enddo ; enddo | |
| 314 | ! if (use_EOS) then | |
| 315 | ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 | |
| 316 | ! M(i,j,k) = M(i,j,k-1) - p(i,j,K) * (alpha_star(i,j,k-1) - alpha_star(i,j,k)) | |
| 317 | ! enddo ; enddo ; enddo | |
| 318 | ! else ! not use_EOS | |
| 319 | ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 | |
| 320 | ! M(i,j,k) = M(i,j,k-1) - p(i,j,K) * dalpha_int(K) | |
| 321 | ! enddo ; enddo ; enddo | |
| 322 | ! endif ! use_EOS | |
| 323 | ||
| 324 | endif | |
| 325 | ||
| 326 | ! Note that ddM/dPb = alpha_star(i,j,1) | |
| 327 | 0 | if (present(pbce)) then |
| 328 | 0 | call Set_pbce_nonBouss(p, tv_tmp, G, GV, US, CS%GFS_scale, pbce, alpha_star) |
| 329 | endif | |
| 330 | ||
| 331 | ! Calculate the pressure force. On a Cartesian grid, | |
| 332 | ! PFu = - dM/dx and PFv = - dM/dy. | |
| 333 | 0 | if (use_EOS) then |
| 334 | !$OMP parallel do default(shared) private(dp_star,PFu_bc,PFv_bc) | |
| 335 | 0 | do k=1,nz |
| 336 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 337 | 0 | dp_star(i,j) = (p(i,j,K+1) - p(i,j,K)) + dp_neglect |
| 338 | enddo ; enddo | |
| 339 | 0 | do j=js,je ; do I=Isq,Ieq |
| 340 | ! PFu_bc = p* grad alpha* | |
| 341 | PFu_bc = (alpha_star(i+1,j,k) - alpha_star(i,j,k)) * (G%IdxCu(I,j) * & | |
| 342 | ((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)))) / & | |
| 343 | 0 | (dp_star(i,j) + dp_star(i+1,j)))) |
| 344 | 0 | PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) + PFu_bc |
| 345 | 0 | if (allocated(CS%PFu_bc)) CS%PFu_bc(i,j,k) = PFu_bc |
| 346 | enddo ; enddo | |
| 347 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 348 | PFv_bc = (alpha_star(i,j+1,k) - alpha_star(i,j,k)) * (G%IdyCv(i,J) * & | |
| 349 | ((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)))) / & | |
| 350 | 0 | (dp_star(i,j) + dp_star(i,j+1)))) |
| 351 | 0 | PFv(i,J,k) = -(M(i,j+1,k) - M(i,j,k)) * G%IdyCv(i,J) + PFv_bc |
| 352 | 0 | if (allocated(CS%PFv_bc)) CS%PFv_bc(i,j,k) = PFv_bc |
| 353 | enddo ; enddo | |
| 354 | enddo ! k-loop | |
| 355 | else ! .not. use_EOS | |
| 356 | !$OMP parallel do default(shared) | |
| 357 | 0 | do k=1,nz |
| 358 | 0 | do j=js,je ; do I=Isq,Ieq |
| 359 | 0 | PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) |
| 360 | enddo ; enddo | |
| 361 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 362 | 0 | PFv(i,J,k) = -(M(i,j+1,k) - M(i,j,k)) * G%IdyCv(i,J) |
| 363 | enddo ; enddo | |
| 364 | enddo | |
| 365 | endif ! use_EOS | |
| 366 | ||
| 367 | 0 | if (CS%id_PFu_bc>0) call post_data(CS%id_PFu_bc, CS%PFu_bc, CS%diag) |
| 368 | 0 | if (CS%id_PFv_bc>0) call post_data(CS%id_PFv_bc, CS%PFv_bc, CS%diag) |
| 369 | ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. | |
| 370 | ! New diagnostics are given for each individual field. | |
| 371 | 0 | if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tide_eq+e_tide_sal, CS%diag) |
| 372 | 0 | if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) |
| 373 | 0 | if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) |
| 374 | 0 | if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) |
| 375 | ||
| 376 | 0 | end subroutine PressureForce_Mont_nonBouss |
| 377 | ||
| 378 | !> \brief Boussinesq Montgomery-potential form of pressure gradient | |
| 379 | !! | |
| 380 | !! Determines the acceleration due to pressure forces. | |
| 381 | !! | |
| 382 | !! To work, the following fields must be set outside of the usual (is:ie,js:je) | |
| 383 | !! range before this subroutine is called: | |
| 384 | !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1). | |
| 385 | 0 | subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta) |
| 386 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. | |
| 387 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 388 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 389 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m]. | |
| 390 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables. | |
| 391 | real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: PFu !< Zonal acceleration due to pressure gradients | |
| 392 | !! (equal to -dM/dx) [L T-2 ~> m s-2]. | |
| 393 | real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients | |
| 394 | !! (equal to -dM/dy) [L T-2 ~> m s-2]. | |
| 395 | type(PressureForce_Mont_CS), intent(inout) :: CS !< Control structure for Montgomery potential PGF | |
| 396 | real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean or | |
| 397 | !! atmosphere-ocean [R L2 T-2 ~> Pa]. | |
| 398 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure anomaly in | |
| 399 | !! each layer due to free surface height anomalies | |
| 400 | !! [L2 T-2 H-1 ~> m s-2]. | |
| 401 | real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> m]. | |
| 402 | ! Local variables | |
| 403 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & | |
| 404 | 0 | M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2]. |
| 405 | 0 | rho_star ! In-situ density divided by the derivative with depth of the |
| 406 | ! corrected e times (G_Earth/Rho0) [L2 Z-1 T-2 ~> m s-2]. | |
| 407 | 0 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height [Z ~> m]. |
| 408 | ! e may be adjusted (with a nonlinear equation of state) so that | |
| 409 | ! its derivative compensates for the adiabatic compressibility | |
| 410 | ! in seawater, but e will still be close to the interface depth. | |
| 411 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: & | |
| 412 | 0 | T_tmp, & ! Temporary array of temperatures where layers that are lighter |
| 413 | ! than the mixed layer have the mixed layer's properties [C ~> degC]. | |
| 414 | 0 | S_tmp ! Temporary array of salinities where layers that are lighter |
| 415 | ! than the mixed layer have the mixed layer's properties [S ~> ppt]. | |
| 416 | ||
| 417 | 0 | real :: Rho_cv_BL(SZI_(G)) ! The coordinate potential density in |
| 418 | ! the deepest variable density near-surface layer [R ~> kg m-3]. | |
| 419 | 0 | real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation |
| 420 | ! for compressibility [Z ~> m]. | |
| 421 | 0 | real :: SSH(SZI_(G),SZJ_(G)) ! The sea surface height anomaly, in depth units [Z ~> m]. |
| 422 | 0 | real :: e_sal(SZI_(G),SZJ_(G)) ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m]. |
| 423 | 0 | real :: e_tide_eq(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal forces from astronomical sources |
| 424 | ! [Z ~> m]. | |
| 425 | 0 | real :: e_tide_sal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to harmonic self-attraction and loading |
| 426 | ! specific to tides, in depth units [Z ~> m]. | |
| 427 | 0 | real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate |
| 428 | ! density [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). | |
| 429 | real :: I_Rho0 ! 1/Rho0 [R-1 ~> m3 kg-1]. | |
| 430 | real :: G_Rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. | |
| 431 | real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer | |
| 432 | ! compensated density gradients [L T-2 ~> m s-2] | |
| 433 | real :: dz_neglect ! A vertical distance that is so small it is usually lost | |
| 434 | ! in roundoff and can be neglected [Z ~> m]. | |
| 435 | logical :: use_p_atm ! If true, use the atmospheric pressure. | |
| 436 | logical :: use_EOS ! If true, density is calculated from T & S using | |
| 437 | ! an equation of state. | |
| 438 | logical :: is_split ! A flag indicating whether the pressure | |
| 439 | ! gradient terms are to be split into | |
| 440 | ! barotropic and baroclinic pieces. | |
| 441 | 0 | type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. |
| 442 | integer, dimension(2) :: EOSdom ! The computational domain for the equation of state | |
| 443 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb | |
| 444 | integer :: i, j, k | |
| 445 | ||
| 446 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 447 | 0 | nkmb=GV%nk_rho_varies |
| 448 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 449 | 0 | EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) |
| 450 | ||
| 451 | 0 | use_p_atm = associated(p_atm) |
| 452 | 0 | is_split = present(pbce) |
| 453 | 0 | use_EOS = associated(tv%eqn_of_state) |
| 454 | ||
| 455 | 0 | if (.not.CS%initialized) call MOM_error(FATAL, & |
| 456 | 0 | "MOM_PressureForce_Mont: Module must be initialized before it is used.") |
| 457 | ||
| 458 | 0 | if (use_EOS) then |
| 459 | 0 | if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, & |
| 460 | "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//& | |
| 461 | 0 | "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.") |
| 462 | endif | |
| 463 | ||
| 464 | 0 | dz_neglect = GV%dZ_subroundoff |
| 465 | 0 | I_Rho0 = 1.0/CS%Rho0 |
| 466 | 0 | G_Rho0 = GV%g_Earth / GV%Rho0 |
| 467 | ||
| 468 | !$OMP parallel do default(shared) | |
| 469 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 470 | 0 | e(i,j,nz+1) = -G%bathyT(i,j) |
| 471 | enddo ; enddo | |
| 472 | ||
| 473 | ! Calculate and add the self-attraction and loading geopotential anomaly. | |
| 474 | 0 | if (CS%calculate_SAL) then |
| 475 | ! Determine the surface height anomaly for calculating self attraction | |
| 476 | ! and loading. This should really be based on bottom pressure anomalies, | |
| 477 | ! but that is not yet implemented, and the current form is correct for | |
| 478 | ! barotropic tides. | |
| 479 | !$OMP parallel do default(shared) | |
| 480 | 0 | do j=Jsq,Jeq+1 |
| 481 | 0 | do i=Isq,Ieq+1 ; SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0) ; enddo |
| 482 | 0 | do k=1,nz ; do i=Isq,Ieq+1 |
| 483 | 0 | SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z |
| 484 | enddo ; enddo | |
| 485 | enddo | |
| 486 | 0 | call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) |
| 487 | !$OMP parallel do default(shared) | |
| 488 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 489 | 0 | e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j) |
| 490 | enddo ; enddo | |
| 491 | endif | |
| 492 | ||
| 493 | ! Calculate and add the tidal geopotential anomaly. | |
| 494 | 0 | if (CS%tides) then |
| 495 | 0 | call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) |
| 496 | !$OMP parallel do default(shared) | |
| 497 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 498 | 0 | e(i,j,nz+1) = e(i,j,nz+1) - (e_tide_eq(i,j) + e_tide_sal(i,j)) |
| 499 | enddo ; enddo | |
| 500 | endif | |
| 501 | ||
| 502 | !$OMP parallel do default(shared) | |
| 503 | 0 | do j=Jsq,Jeq+1 ; do k=nz,1,-1 ; do i=Isq,Ieq+1 |
| 504 | 0 | e(i,j,K) = e(i,j,K+1) + h(i,j,k)*GV%H_to_Z |
| 505 | enddo ; enddo ; enddo | |
| 506 | ||
| 507 | 0 | if (use_EOS) then |
| 508 | ! Calculate in-situ densities (rho_star). | |
| 509 | ||
| 510 | ! With a bulk mixed layer, replace the T & S of any layers that are | |
| 511 | ! lighter than the buffer layer with the properties of the buffer | |
| 512 | ! layer. These layers will be massless anyway, and it avoids any | |
| 513 | ! formal calculations with hydrostatically unstable profiles. | |
| 514 | ||
| 515 | 0 | if (nkmb>0) then |
| 516 | 0 | tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp |
| 517 | 0 | tv_tmp%eqn_of_state => tv%eqn_of_state |
| 518 | ||
| 519 | 0 | do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo |
| 520 | !$OMP parallel do default(shared) private(Rho_cv_BL) | |
| 521 | 0 | do j=Jsq,Jeq+1 |
| 522 | 0 | do k=1,nkmb ; do i=Isq,Ieq+1 |
| 523 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) |
| 524 | enddo ; enddo | |
| 525 | call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), & | |
| 526 | 0 | tv%eqn_of_state, EOSdom) |
| 527 | ||
| 528 | 0 | do k=nkmb+1,nz ; do i=Isq,Ieq+1 |
| 529 | 0 | if (GV%Rlay(k) < Rho_cv_BL(i)) then |
| 530 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) |
| 531 | else | |
| 532 | 0 | tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) |
| 533 | endif | |
| 534 | enddo ; enddo | |
| 535 | enddo | |
| 536 | else | |
| 537 | 0 | tv_tmp%T => tv%T ; tv_tmp%S => tv%S |
| 538 | 0 | tv_tmp%eqn_of_state => tv%eqn_of_state |
| 539 | 0 | do i=Isq,Ieq+1 ; p_ref(i) = 0.0 ; enddo |
| 540 | endif | |
| 541 | ||
| 542 | ! This no longer includes any pressure dependency, since this routine | |
| 543 | ! will come down with a fatal error if there is any compressibility. | |
| 544 | !$OMP parallel do default(shared) | |
| 545 | 0 | do k=1,nz ; do j=Jsq,Jeq+1 |
| 546 | call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), & | |
| 547 | 0 | tv%eqn_of_state, EOSdom) |
| 548 | 0 | do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo |
| 549 | enddo ; enddo | |
| 550 | endif ! use_EOS | |
| 551 | ||
| 552 | ! Here the layer Montgomery potentials, M, are calculated. | |
| 553 | 0 | if (use_EOS) then |
| 554 | !$OMP parallel do default(shared) | |
| 555 | 0 | do j=Jsq,Jeq+1 |
| 556 | 0 | do i=Isq,Ieq+1 |
| 557 | 0 | M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1)) |
| 558 | 0 | if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 |
| 559 | enddo | |
| 560 | 0 | do k=2,nz ; do i=Isq,Ieq+1 |
| 561 | 0 | M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K) |
| 562 | enddo ; enddo | |
| 563 | enddo | |
| 564 | else ! not use_EOS | |
| 565 | !$OMP parallel do default(shared) | |
| 566 | 0 | do j=Jsq,Jeq+1 |
| 567 | 0 | do i=Isq,Ieq+1 |
| 568 | 0 | M(i,j,1) = GV%g_prime(1) * e(i,j,1) |
| 569 | 0 | if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 |
| 570 | enddo | |
| 571 | 0 | do k=2,nz ; do i=Isq,Ieq+1 |
| 572 | 0 | M(i,j,k) = M(i,j,k-1) + GV%g_prime(K) * e(i,j,K) |
| 573 | enddo ; enddo | |
| 574 | enddo | |
| 575 | endif ! use_EOS | |
| 576 | ||
| 577 | 0 | if (present(pbce)) then |
| 578 | 0 | call Set_pbce_Bouss(e, tv_tmp, G, GV, US, CS%Rho0, CS%GFS_scale, pbce, rho_star, 1) |
| 579 | endif | |
| 580 | ||
| 581 | ! Calculate the pressure force. On a Cartesian grid, | |
| 582 | ! PFu = - dM/dx and PFv = - dM/dy. | |
| 583 | 0 | if (use_EOS) then |
| 584 | !$OMP parallel do default(shared) private(h_star,PFu_bc,PFv_bc) | |
| 585 | 0 | do k=1,nz |
| 586 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 587 | 0 | h_star(i,j) = (e(i,j,K) - e(i,j,K+1)) + dz_neglect |
| 588 | enddo ; enddo | |
| 589 | 0 | do j=js,je ; do I=Isq,Ieq |
| 590 | PFu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (G%IdxCu(I,j) * & | |
| 591 | ((h_star(i,j) * h_star(i+1,j) - ((e(i,j,K) * h_star(i+1,j)) + & | |
| 592 | 0 | (e(i+1,j,K) * h_star(i,j)))) / (h_star(i,j) + h_star(i+1,j)))) |
| 593 | 0 | PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) + PFu_bc |
| 594 | 0 | if (allocated(CS%PFu_bc)) CS%PFu_bc(i,j,k) = PFu_bc |
| 595 | enddo ; enddo | |
| 596 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 597 | PFv_bc = -1.0*(rho_star(i,j+1,k) - rho_star(i,j,k)) * (G%IdyCv(i,J) * & | |
| 598 | ((h_star(i,j) * h_star(i,j+1) - ((e(i,j,K) * h_star(i,j+1)) + & | |
| 599 | 0 | (e(i,j+1,K) * h_star(i,j)))) / (h_star(i,j) + h_star(i,j+1)))) |
| 600 | 0 | PFv(i,J,k) = -(M(i,j+1,k) - M(i,j,k)) * G%IdyCv(i,J) + PFv_bc |
| 601 | 0 | if (allocated(CS%PFv_bc)) CS%PFv_bc(i,j,k) = PFv_bc |
| 602 | enddo ; enddo | |
| 603 | enddo ! k-loop | |
| 604 | else ! .not. use_EOS | |
| 605 | !$OMP parallel do default(shared) | |
| 606 | 0 | do k=1,nz |
| 607 | 0 | do j=js,je ; do I=Isq,Ieq |
| 608 | 0 | PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) |
| 609 | enddo ; enddo | |
| 610 | 0 | do J=Jsq,Jeq ; do i=is,ie |
| 611 | 0 | PFv(i,J,k) = -(M(i,j+1,k) - M(i,j,k)) * G%IdyCv(i,J) |
| 612 | enddo ; enddo | |
| 613 | enddo | |
| 614 | endif ! use_EOS | |
| 615 | ||
| 616 | 0 | if (present(eta)) then |
| 617 | ! eta is the sea surface height relative to a time-invariant geoid, for | |
| 618 | ! comparison with what is used for eta in btstep. See how e was calculated | |
| 619 | ! about 200 lines above. | |
| 620 | !$OMP parallel do default(shared) | |
| 621 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 622 | 0 | eta(i,j) = e(i,j,1)*GV%Z_to_H |
| 623 | enddo ; enddo | |
| 624 | 0 | if (CS%tides) then |
| 625 | !$OMP parallel do default(shared) | |
| 626 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 627 | 0 | eta(i,j) = eta(i,j) + (e_tide_eq(i,j)+e_tide_sal(i,j))*GV%Z_to_H |
| 628 | enddo ; enddo | |
| 629 | endif | |
| 630 | 0 | if (CS%calculate_SAL) then |
| 631 | !$OMP parallel do default(shared) | |
| 632 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 633 | 0 | eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H |
| 634 | enddo ; enddo | |
| 635 | endif | |
| 636 | endif | |
| 637 | ||
| 638 | 0 | if (CS%id_PFu_bc>0) call post_data(CS%id_PFu_bc, CS%PFu_bc, CS%diag) |
| 639 | 0 | if (CS%id_PFv_bc>0) call post_data(CS%id_PFv_bc, CS%PFv_bc, CS%diag) |
| 640 | ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. | |
| 641 | ! New diagnostics are given for each individual field. | |
| 642 | 0 | if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tide_eq+e_tide_sal, CS%diag) |
| 643 | 0 | if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) |
| 644 | 0 | if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) |
| 645 | 0 | if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) |
| 646 | ||
| 647 | 0 | end subroutine PressureForce_Mont_Bouss |
| 648 | ||
| 649 | !> Determines the partial derivative of the acceleration due | |
| 650 | !! to pressure forces with the free surface height. | |
| 651 | 0 | subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star, nkblock) |
| 652 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure | |
| 653 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 654 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m]. | |
| 655 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables | |
| 656 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 657 | real, intent(in) :: Rho0 !< The "Boussinesq" ocean density [R ~> kg m-3]. | |
| 658 | real, intent(in) :: GFS_scale !< Ratio between gravity applied to top | |
| 659 | !! interface and the gravitational acceleration of | |
| 660 | !! the planet [nondim]. Usually this ratio is 1. | |
| 661 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 662 | intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due | |
| 663 | !! to free surface height anomalies | |
| 664 | !! [L2 T-2 H-1 ~> m s-2]. | |
| 665 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 666 | optional, intent(in) :: rho_star !< The layer densities (maybe compressibility | |
| 667 | !! compensated), times g/rho_0 [L2 Z-1 T-2 ~> m s-2]. | |
| 668 | integer, intent(in) :: nkblock !< The number of vertical levels in each block for the variable density | |
| 669 | !! near-surface layer. This is needed to determine whether rho_star is | |
| 670 | !! compressibility compensated or not. | |
| 671 | 48 | real :: Ihtot(SZI_(G),SZJ_(G)) ! The inverse of the sum of the layer thicknesses [H-1 ~> m-1 or m2 kg-1]. |
| 672 | 48 | real :: press(SZI_(G),SZJ_(G),nkblock) ! Interface pressure [R L2 T-2 ~> Pa]. |
| 673 | 48 | real :: T_int(SZI_(G),SZJ_(G),nkblock) ! Interface temperature [C ~> degC] |
| 674 | 48 | real :: S_int(SZI_(G),SZJ_(G),nkblock) ! Interface salinity [S ~> ppt] |
| 675 | 48 | real :: dR_dT(SZI_(G),SZJ_(G),nkblock) ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] |
| 676 | 48 | real :: dR_dS(SZI_(G),SZJ_(G),nkblock) ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. |
| 677 | 48 | real :: rho_in_situ(SZI_(G),SZJ_(G)) ! In-situ density at the top of a layer [R ~> kg m-3]. |
| 678 | real :: G_Rho0 ! A scaled version of g_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] | |
| 679 | real :: Rho0xG ! g_Earth * Rho0 [R L2 Z-1 T-2 ~> kg s-2 m-2] | |
| 680 | logical :: use_EOS ! If true, density is calculated from T & S using | |
| 681 | ! an equation of state. | |
| 682 | real :: dz_neglect ! A vertical distance that is so small it is usually lost | |
| 683 | ! in roundoff and can be neglected [Z ~> m]. | |
| 684 | integer :: EOSdom(3,2) ! The computational domain for the equation of state | |
| 685 | integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, kstart, kend, kk | |
| 686 | ||
| 687 | !$omp target data map(alloc: Ihtot) | |
| 688 | 24 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke |
| 689 | ||
| 690 | 24 | use_EOS = associated(tv%eqn_of_state) |
| 691 | 24 | dz_neglect = GV%dZ_subroundoff |
| 692 | ||
| 693 | 24 | if (use_EOS) then |
| 694 | 24 | if (present(rho_star)) then |
| 695 | 0 | do concurrent (j=Jsq:Jeq+1) |
| 696 | 0 | do concurrent (i=Isq:Ieq+1) |
| 697 | 0 | Ihtot(i,j) = GV%H_to_Z / ((e(i,j,1) - e(i,j,nz+1)) + dz_neglect) |
| 698 | 0 | pbce(i,j,1) = GFS_scale * rho_star(i,j,1) * GV%H_to_Z |
| 699 | enddo | |
| 700 | ||
| 701 | 0 | do k=2,nz |
| 702 | 0 | do concurrent (i=Isq:Ieq+1) |
| 703 | pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) & | |
| 704 | 0 | * ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i,j)) |
| 705 | enddo | |
| 706 | enddo | |
| 707 | enddo | |
| 708 | else | |
| 709 | !$omp target data & | |
| 710 | !$omp map(alloc: EOSdom, press, T_int, S_int, rho_in_situ) & | |
| 711 | !$omp map(alloc: dR_dT, dR_dS) | |
| 712 | ||
| 713 | 24 | Rho0xG = Rho0 * GV%g_Earth |
| 714 | 24 | G_Rho0 = GV%g_Earth / GV%Rho0 |
| 715 | ||
| 716 | 72 | EOSdom(1,:) = [Isq - (G%isd-1), G%iec+1 - (G%isd-1)] |
| 717 | 72 | EOSdom(2,:) = [Jsq - (G%jsd-1), G%jec+1 - (G%jsd-1)] |
| 718 | 2952 | do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1) |
| 719 | 181536 | Ihtot(i,j) = GV%H_to_Z / ((e(i,j,1) - e(i,j,nz+1)) + dz_neglect) |
| 720 | 184488 | press(i,j,1) = -Rho0xG * (e(i,j,1) - G%meanSL(i,j)) |
| 721 | enddo | |
| 722 | ||
| 723 | call calculate_density(tv%T(:,:,1), tv%S(:,:,1), press(:,:,1), rho_in_situ, & | |
| 724 | 24 | tv%eqn_of_state, EOSdom(1:2,:)) |
| 725 | ||
| 726 | 2952 | do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1) |
| 727 | 184488 | pbce(i,j,1) = G_Rho0 * (GFS_scale * rho_in_situ(i,j)) * GV%H_to_Z |
| 728 | enddo | |
| 729 | ||
| 730 | 24 | EOSdom(3,1) = 1 |
| 731 | 1800 | do kstart=2,nz,nkblock |
| 732 | 1776 | kend = min(kstart+nkblock-1, nz) |
| 733 | 218448 | do concurrent (k=kstart:kend, j=Jsq:Jeq+1, i=Isq:Ieq+1) |
| 734 | 13433664 | kk=k-kstart+1 |
| 735 | 13433664 | press(i,j,kk) = -Rho0xG * (e(i,j,K) - G%meanSL(i,j)) |
| 736 | 13433664 | T_int(i,j,kk) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k)) |
| 737 | 27085776 | S_int(i,j,kk) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k)) |
| 738 | enddo | |
| 739 | 1776 | EOSdom(3,2) = kend-kstart+1 |
| 740 | call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, & | |
| 741 | 1776 | tv%eqn_of_state, EOSdom) |
| 742 | ||
| 743 | 24 | do concurrent(j=Jsq:Jeq+1) |
| 744 | 222000 | do k=kstart,kend ; do concurrent (i=Isq:Ieq+1) |
| 745 | 13433664 | kk=k-kstart+1 |
| 746 | pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * & | |
| 747 | ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i,j)) * & | |
| 748 | (dR_dT(i,j,kk) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + & | |
| 749 | 13543776 | dR_dS(i,j,kk) * (tv%S(i,j,k) - tv%S(i,j,k-1))) |
| 750 | enddo ; enddo | |
| 751 | enddo | |
| 752 | enddo | |
| 753 | !$omp end target data | |
| 754 | endif | |
| 755 | else | |
| 756 | 0 | do concurrent (j=Jsq:Jeq+1) |
| 757 | 0 | do concurrent (i=Isq:Ieq+1) |
| 758 | 0 | Ihtot(i,j) = 1.0 / ((e(i,j,1) - e(i,j,nz+1)) + dz_neglect) |
| 759 | 0 | pbce(i,j,1) = GV%g_prime(1) * GV%H_to_Z |
| 760 | enddo | |
| 761 | 0 | do k=2,nz |
| 762 | 0 | do concurrent (i=Isq:Ieq+1) |
| 763 | pbce(i,j,k) = pbce(i,j,k-1) + (GV%g_prime(K) * GV%H_to_Z) & | |
| 764 | 0 | * ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i,j)) |
| 765 | enddo | |
| 766 | enddo | |
| 767 | enddo | |
| 768 | endif | |
| 769 | !$omp end target data | |
| 770 | 24 | end subroutine Set_pbce_Bouss |
| 771 | ||
| 772 | !> Determines the partial derivative of the acceleration due | |
| 773 | !! to pressure forces with the column mass. | |
| 774 | 0 | subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) |
| 775 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure | |
| 776 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 777 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: p !< Interface pressures [R L2 T-2 ~> Pa]. | |
| 778 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables | |
| 779 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 780 | real, intent(in) :: GFS_scale !< Ratio between gravity applied to top | |
| 781 | !! interface and the gravitational acceleration of | |
| 782 | !! the planet [nondim]. Usually this ratio is 1. | |
| 783 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: pbce !< The baroclinic pressure anomaly in each | |
| 784 | !! layer due to free surface height anomalies | |
| 785 | !! [L2 H-1 T-2 ~> m4 kg-1 s-2]. | |
| 786 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(in) :: alpha_star !< The layer specific volumes | |
| 787 | !! (maybe compressibility compensated) [R-1 ~> m3 kg-1]. | |
| 788 | ! Local variables | |
| 789 | real, dimension(SZI_(G),SZJ_(G)) :: & | |
| 790 | 0 | dpbce, & ! A barotropic correction to the pbce to enable the use of |
| 791 | ! a reduced gravity form of the equations [L2 H-1 T-2 ~> m4 kg-1 s-2]. | |
| 792 | 0 | C_htot ! dP_dH divided by the total ocean pressure [H-1 ~> m2 kg-1]. |
| 793 | 0 | real :: T_int(SZI_(G)) ! Interface temperature [C ~> degC] |
| 794 | 0 | real :: S_int(SZI_(G)) ! Interface salinity [S ~> ppt] |
| 795 | 0 | real :: dR_dT(SZI_(G)) ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] |
| 796 | 0 | real :: dR_dS(SZI_(G)) ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. |
| 797 | 0 | real :: rho_in_situ(SZI_(G)) ! In-situ density at an interface [R ~> kg m-3]. |
| 798 | 0 | real :: alpha_Lay(SZK_(GV)) ! The specific volume of each layer [R-1 ~> m3 kg-1]. |
| 799 | 0 | real :: dalpha_int(SZK_(GV)+1) ! The change in specific volume across each interface [R-1 ~> m3 kg-1]. |
| 800 | real :: dP_dH ! A factor that converts from thickness to pressure times other dimensional | |
| 801 | ! conversion factors [R L2 T-2 H-1 ~> Pa m2 kg-1]. | |
| 802 | real :: dp_neglect ! A thickness that is so small it is usually lost | |
| 803 | ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. | |
| 804 | logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. | |
| 805 | integer, dimension(2) :: EOSdom ! The computational domain for the equation of state | |
| 806 | integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k | |
| 807 | ||
| 808 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke |
| 809 | 0 | EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) |
| 810 | ||
| 811 | 0 | use_EOS = associated(tv%eqn_of_state) |
| 812 | ||
| 813 | 0 | dP_dH = GV%g_Earth * GV%H_to_RZ |
| 814 | 0 | dp_neglect = GV%g_Earth * GV%H_to_RZ * GV%H_subroundoff |
| 815 | ||
| 816 | 0 | if (use_EOS) then |
| 817 | 0 | if (present(alpha_star)) then |
| 818 | !$OMP parallel do default(shared) | |
| 819 | 0 | do j=Jsq,Jeq+1 |
| 820 | 0 | do i=Isq,Ieq+1 |
| 821 | 0 | C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) |
| 822 | 0 | pbce(i,j,nz) = dP_dH * alpha_star(i,j,nz) |
| 823 | enddo | |
| 824 | 0 | do k=nz-1,1,-1 ; do i=Isq,Ieq+1 |
| 825 | pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1)) * C_htot(i,j)) * & | |
| 826 | 0 | (alpha_star(i,j,k) - alpha_star(i,j,k+1)) |
| 827 | enddo ; enddo | |
| 828 | enddo | |
| 829 | else | |
| 830 | !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ) | |
| 831 | 0 | do j=Jsq,Jeq+1 |
| 832 | call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, & | |
| 833 | 0 | tv%eqn_of_state, EOSdom) |
| 834 | 0 | do i=Isq,Ieq+1 |
| 835 | 0 | C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) |
| 836 | 0 | pbce(i,j,nz) = dP_dH / (rho_in_situ(i)) |
| 837 | enddo | |
| 838 | 0 | do k=nz-1,1,-1 |
| 839 | 0 | do i=Isq,Ieq+1 |
| 840 | 0 | T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1)) |
| 841 | 0 | S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1)) |
| 842 | enddo | |
| 843 | 0 | call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, EOSdom) |
| 844 | call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & | |
| 845 | 0 | tv%eqn_of_state, EOSdom) |
| 846 | 0 | do i=Isq,Ieq+1 |
| 847 | pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * & | |
| 848 | ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + & | |
| 849 | 0 | dR_dS(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / (rho_in_situ(i)**2)) |
| 850 | enddo | |
| 851 | enddo | |
| 852 | enddo | |
| 853 | endif | |
| 854 | else ! not use_EOS | |
| 855 | ||
| 856 | 0 | do k=1,nz ; alpha_Lay(k) = 1.0 / (GV%Rlay(k)) ; enddo |
| 857 | 0 | do k=2,nz ; dalpha_int(K) = alpha_Lay(k-1) - alpha_Lay(k) ; enddo |
| 858 | ||
| 859 | !$OMP parallel do default(shared) | |
| 860 | 0 | do j=Jsq,Jeq+1 |
| 861 | 0 | do i=Isq,Ieq+1 |
| 862 | 0 | C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) |
| 863 | 0 | pbce(i,j,nz) = dP_dH * alpha_Lay(nz) |
| 864 | enddo | |
| 865 | 0 | do k=nz-1,1,-1 ; do i=Isq,Ieq+1 |
| 866 | 0 | 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) |
| 867 | enddo ; enddo | |
| 868 | enddo | |
| 869 | endif ! use_EOS | |
| 870 | ||
| 871 | 0 | if (GFS_scale < 1.0) then |
| 872 | ! Adjust the Montgomery potential to make this a reduced gravity model. | |
| 873 | !$OMP parallel do default(shared) | |
| 874 | 0 | do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 875 | 0 | dpbce(i,j) = (GFS_scale - 1.0) * pbce(i,j,1) |
| 876 | enddo ; enddo | |
| 877 | !$OMP parallel do default(shared) | |
| 878 | 0 | do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 |
| 879 | 0 | pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j) |
| 880 | enddo ; enddo ; enddo | |
| 881 | endif | |
| 882 | ||
| 883 | 0 | end subroutine Set_pbce_nonBouss |
| 884 | ||
| 885 | !> Initialize the Montgomery-potential form of PGF control structure | |
| 886 | 0 | subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, tides_CSp) |
| 887 | type(time_type), target, intent(in) :: Time !< Current model time | |
| 888 | type(ocean_grid_type), intent(in) :: G !< ocean grid structure | |
| 889 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 890 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 891 | type(param_file_type), intent(in) :: param_file !< Parameter file handles | |
| 892 | type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure | |
| 893 | type(PressureForce_Mont_CS), intent(inout) :: CS !< Montgomery PGF control structure | |
| 894 | type(SAL_CS), intent(in), target, optional :: SAL_CSp !< SAL control structure | |
| 895 | type(tidal_forcing_CS), intent(in), target, optional :: tides_CSp !< Tides control structure | |
| 896 | ||
| 897 | ! Local variables | |
| 898 | logical :: use_EOS | |
| 899 | ! This include declares and sets the variable "version". | |
| 900 | # include "version_variable.h" | |
| 901 | character(len=40) :: mdl ! This module's name. | |
| 902 | ||
| 903 | 0 | CS%initialized = .true. |
| 904 | 0 | CS%diag => diag ; CS%Time => Time |
| 905 | 0 | if (present(tides_CSp)) & |
| 906 | 0 | CS%tides_CSp => tides_CSp |
| 907 | 0 | if (present(SAL_CSp)) & |
| 908 | 0 | CS%SAL_CSp => SAL_CSp |
| 909 | ||
| 910 | 0 | mdl = "MOM_PressureForce_Mont" |
| 911 | 0 | call log_version(param_file, mdl, version, "") |
| 912 | call get_param(param_file, mdl, "RHO_0", CS%Rho0, & | |
| 913 | "The mean ocean density used with BOUSSINESQ true to "//& | |
| 914 | "calculate accelerations and the mass for conservation "//& | |
| 915 | "properties, or with BOUSSINESQ false to convert some "//& | |
| 916 | "parameters from vertical units of m to kg m-2.", & | |
| 917 | 0 | units="kg m-3", default=1035.0, scale=US%R_to_kg_m3) |
| 918 | call get_param(param_file, mdl, "TIDES", CS%tides, & | |
| 919 | 0 | "If true, apply tidal momentum forcing.", default=.false.) |
| 920 | call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & | |
| 921 | 0 | "If true, calculate self-attraction and loading.", default=CS%tides) |
| 922 | call get_param(param_file, mdl, "USE_EOS", use_EOS, default=.true., & | |
| 923 | 0 | do_not_log=.true.) ! Input for diagnostic use only. |
| 924 | ||
| 925 | 0 | if (use_EOS) then |
| 926 | CS%id_PFu_bc = register_diag_field('ocean_model', 'PFu_bc', diag%axesCuL, Time, & | |
| 927 | 0 | 'Density Gradient Zonal Pressure Force Accel.', "meter second-2", conversion=US%L_T2_to_m_s2) |
| 928 | CS%id_PFv_bc = register_diag_field('ocean_model', 'PFv_bc', diag%axesCvL, Time, & | |
| 929 | 0 | 'Density Gradient Meridional Pressure Force Accel.', "meter second-2", conversion=US%L_T2_to_m_s2) |
| 930 | 0 | if (CS%id_PFu_bc > 0) & |
| 931 | 0 | allocate(CS%PFu_bc(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.) |
| 932 | 0 | if (CS%id_PFv_bc > 0) & |
| 933 | 0 | allocate(CS%PFv_bc(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.) |
| 934 | endif | |
| 935 | ||
| 936 | 0 | if (CS%calculate_SAL) then |
| 937 | CS%id_e_sal = register_diag_field('ocean_model', 'e_sal', diag%axesT1, Time, & | |
| 938 | 0 | 'Self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m) |
| 939 | endif | |
| 940 | 0 | if (CS%tides) then |
| 941 | CS%id_e_tide = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, Time, & | |
| 942 | 0 | 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=US%Z_to_m) |
| 943 | CS%id_e_tide_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, & | |
| 944 | 0 | 'Equilibrium tides height anomaly', 'meter', conversion=US%Z_to_m) |
| 945 | CS%id_e_tide_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, & | |
| 946 | 0 | 'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m) |
| 947 | endif | |
| 948 | ||
| 949 | 0 | CS%GFS_scale = 1.0 |
| 950 | 0 | if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth |
| 951 | ||
| 952 | 0 | call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale, units="nondim") |
| 953 | ||
| 954 | 0 | end subroutine PressureForce_Mont_init |
| 955 | ||
| 956 | !>\namespace mom_pressureforce_mont | |
| 957 | !! | |
| 958 | !! Provides the Boussunesq and non-Boussinesq forms of the horizontal | |
| 959 | !! accelerations due to pressure gradients using the Montgomery potential. A | |
| 960 | !! second-order accurate, centered scheme is used. If a split time stepping | |
| 961 | !! scheme is used, the vertical decomposition into barotropic and baroclinic | |
| 962 | !! contributions described by Hallberg (J Comp Phys 1997) is used. With a | |
| 963 | !! nonlinear equation of state, compressibility is added along the lines proposed | |
| 964 | !! by Sun et al. (JPO 1999), but with compressibility coefficients based on a fit | |
| 965 | !! to a user-provided reference profile. | |
| 966 | ||
| 967 | 0 | end module MOM_PressureForce_Mont |