Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into DOME_initialization_cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Nov 28, 2021
2 parents 14f807a + a65fbed commit 459ecc0
Show file tree
Hide file tree
Showing 20 changed files with 410 additions and 403 deletions.
30 changes: 14 additions & 16 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,9 +338,9 @@ module MOM_barotropic
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_u_EE. uBT_EE must be non-positive.
real :: uh_crvW !< The curvature of face area with velocity for flow from the west [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
!! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: uh_crvE !< The curvature of face area with velocity for flow from the east [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
!! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: uh_WW !< The zonal transport when ubt=ubt_WW [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
real :: uh_EE !< The zonal transport when ubt=ubt_EE [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
Expand All @@ -364,9 +364,9 @@ module MOM_barotropic
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_v_NN. vBT_NN must be non-positive.
real :: vh_crvS !< The curvature of face area with velocity for flow from the south [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
!! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: vh_crvN !< The curvature of face area with velocity for flow from the north [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
!! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: vh_SS !< The meridional transport when vbt=vbt_SS [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
real :: vh_NN !< The meridional transport when vbt=vbt_NN [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
Expand Down Expand Up @@ -622,24 +622,23 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
vhbt_prev, vhbt_sum_prev, & ! Previous transports stored for OBCs [L2 H T-1 ~> m3 s-1]
vbt_int_prev, & ! Previous value of time-integrated velocity stored for OBCs [L ~> m]
vhbt_int_prev ! Previous value of time-integrated transport stored for OBCs [L2 H ~> m3]
real :: mass_to_Z ! The depth unit conversion divided by the mean density (Rho0) [Z m-1 R-1 ~> m3 kg-1] !### R-1
real :: mass_accel_to_Z ! The inverse of the mean density (Rho0) [R-1 ~> m3 kg-1].
real :: visc_rem ! A work variable that may equal visc_rem_[uv]. Nondim.
real :: mass_to_Z ! The inverse of the the mean density (Rho0) [R-1 ~> m3 kg-1]
real :: visc_rem ! A work variable that may equal visc_rem_[uv] [nondim]
real :: vel_prev ! The previous velocity [L T-1 ~> m s-1].
real :: dtbt ! The barotropic time step [T ~> s].
real :: bebt ! A copy of CS%bebt [nondim].
real :: be_proj ! The fractional amount by which velocities are projected
! when project_velocity is true. For now be_proj is set
! when project_velocity is true [nondim]. For now be_proj is set
! to equal bebt, as they have similar roles and meanings.
real :: Idt ! The inverse of dt [T-1 ~> s-1].
real :: det_de ! The partial derivative due to self-attraction and loading
! of the reference geopotential with the sea surface height.
! of the reference geopotential with the sea surface height [nondim].
! This is typically ~0.09 or less.
real :: dgeo_de ! The constant of proportionality between geopotential and
! sea surface height. It is a nondimensional number of
! sea surface height [nondim]. It is a nondimensional number of
! order 1. For stability, this may be made larger
! than the physical problem would suggest.
real :: Instep ! The inverse of the number of barotropic time steps to take.
real :: Instep ! The inverse of the number of barotropic time steps to take [nondim].
real :: wt_end ! The weighting of the final value of eta_PF [nondim]
integer :: nstep ! The number of barotropic time steps to take.
type(time_type) :: &
Expand Down Expand Up @@ -772,8 +771,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
Idtbt = 1.0 / dtbt
bebt = CS%bebt
be_proj = CS%bebt
mass_accel_to_Z = 1.0 / GV%Rho0
mass_to_Z = US%m_to_Z / GV%Rho0 !### THis should be the same as mass_accel_to_Z.
mass_to_Z = 1.0 / GV%Rho0

!--- setup the weight when computing vbt_trans and ubt_trans
if (project_velocity) then
Expand Down Expand Up @@ -1243,7 +1241,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
endif
endif

BT_force_u(I,j) = forces%taux(I,j) * mass_accel_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1)
BT_force_u(I,j) = forces%taux(I,j) * mass_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1)
else
BT_force_u(I,j) = 0.0
endif ; enddo ; enddo
Expand All @@ -1269,11 +1267,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
endif
endif

BT_force_v(i,J) = forces%tauy(i,J) * mass_accel_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1)
BT_force_v(i,J) = forces%tauy(i,J) * mass_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1)
else
BT_force_v(i,J) = 0.0
endif ; enddo ; enddo
if (associated(taux_bot) .and. associated(tauy_bot)) then
if (associated(taux_bot) .and. associated(tauy_bot)) then
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j)
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_checksum_packages.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric)
integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully
!! symmetric computational domain.
real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> nondim] or [nondim]
real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> 1] or [1]
integer :: hs
logical :: sym

Expand Down
52 changes: 26 additions & 26 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -925,34 +925,32 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< salinity [ppt]
type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type
integer, intent(in) :: j !< j-row to work on
real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G)), intent(inout) :: netHeatMinusSW !< Surface heat flux excluding shortwave
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G)), intent(inout) :: netSalt !< surface salt flux
!! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1]
real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G)), intent(out) :: netHeatMinusSW !< Surface heat flux excluding shortwave
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G)), intent(out) :: netSalt !< surface salt flux
!! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

! local variables
integer :: k
real, parameter :: dt = 1. ! to return a rate from extractFluxes1d
real, dimension(SZI_(G)) :: netH ! net FW flux [H s-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netH ! net FW flux [H T-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netEvap ! net FW flux leaving ocean via evaporation
! [H s-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netHeat ! net temp flux [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! [H T-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netHeat ! net temp flux [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(max(nsw,1), SZI_(G)) :: penSWbnd ! penetrating SW radiation by band
! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G)) :: pressure ! pressure at the surface [R L2 T-2 ~> Pa]
real, dimension(SZI_(G)) :: dRhodT ! density partial derivative wrt temp [R degC-1 ~> kg m-3 degC-1]
real, dimension(SZI_(G)) :: dRhodS ! density partial derivative wrt saln [R ppt-1 ~> kg m-3 ppt-1]
real, dimension(SZI_(G),SZK_(GV)+1) :: netPen ! The net penetrating shortwave radiation at each level
! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]

logical :: useRiverHeatContent
logical :: useCalvingHeatContent
real :: depthBeforeScalingFluxes ! A depth scale [H ~> m or kg m-2]
real :: GoRho ! The gravitational acceleration divided by mean density times some
! unit conversion factors [L2 H-1 s R-1 T-3 ~> m4 kg-1 s-2 or m7 kg-2 s-2]
real :: GoRho ! The gravitational acceleration divided by mean density times a
! unit conversion factor [L2 H-1 R-1 T-2 ~> m4 kg-1 s-2 or m7 kg-2 s-2]
real :: H_limit_fluxes ! Another depth scale [H ~> m or kg m-2]
integer :: i
integer :: i, k

! smg: what do we do when have heat fluxes from calving and river?
useRiverHeatContent = .False.
Expand All @@ -961,38 +959,40 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
depthBeforeScalingFluxes = max( GV%Angstrom_H, 1.e-30*GV%m_to_H )
pressure(:) = 0.
if (associated(tv%p_surf)) then ; do i=G%isc,G%iec ; pressure(i) = tv%p_surf(i,j) ; enddo ; endif
GoRho = (GV%g_Earth * GV%H_to_Z*US%T_to_s) / GV%Rho0
GoRho = (GV%g_Earth * GV%H_to_Z) / GV%Rho0

H_limit_fluxes = depthBeforeScalingFluxes

! The surface forcing is contained in the fluxes type.
! We aggregate the thermodynamic forcing for a time step into the following:
! netH = water added/removed via surface fluxes [H s-1 ~> m s-1 or kg m-2 s-1]
! netHeat = heat via surface fluxes [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! netSalt = salt via surface fluxes [ppt H s-1 ~> ppt m s-1 or gSalt m-2 s-1]
! netH = water added/removed via surface fluxes [H T-1 ~> m s-1 or kg m-2 s-1]
! netHeat = heat via surface fluxes [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
! netSalt = salt via surface fluxes [ppt H T-1 ~> ppt m s-1 or gSalt m-2 s-1]
! Note that unlike other calls to extractFLuxes1d() that return the time-integrated flux
! this call returns the rate because dt=1
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt*US%s_to_T, &
! this call returns the rate because dt=1 (in arbitrary time units)
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, 1.0, &
depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW, &
netSalt, penSWbnd, tv, .false.)

! Sum over bands and attenuate as a function of depth
! netPen is the netSW as a function of depth
call sumSWoverBands(G, GV, US, h(:,j,:), optics_nbands(optics), optics, j, dt*US%s_to_T, &
call sumSWoverBands(G, GV, US, h(:,j,:), optics_nbands(optics), optics, j, 1.0, &
H_limit_fluxes, .true., penSWbnd, netPen)

! Density derivatives
call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, &
tv%eqn_of_state, EOS_domain(G%HI))

! Adjust netSalt to reflect dilution effect of FW flux
netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec) ! ppt H/s
! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec)

! Add in the SW heating for purposes of calculating the net
! surface buoyancy flux affecting the top layer.
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
!netHeat(:) = netHeatMinusSW(:) + sum( penSWbnd, dim=1 )
netHeat(G%isc:G%iec) = netHeatMinusSW(G%isc:G%iec) + netPen(G%isc:G%iec,1) ! K H/s
netHeat(G%isc:G%iec) = netHeatMinusSW(G%isc:G%iec) + netPen(G%isc:G%iec,1)

! Convert to a buoyancy flux, excluding penetrating SW heating
buoyancyFlux(G%isc:G%iec,1) = - GoRho * ( dRhodS(G%isc:G%iec) * netSalt(G%isc:G%iec) + &
Expand Down Expand Up @@ -1020,9 +1020,9 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv,
type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: netHeatMinusSW !< surface heat flux excluding shortwave
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: netSalt !< Net surface salt flux
!! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1]
!! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
! local variables
integer :: j

Expand Down
9 changes: 5 additions & 4 deletions src/core/MOM_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,10 +171,11 @@ module MOM_grid
! initialization routines (but not all)
real :: south_lat !< The latitude (or y-coordinate) of the first v-line
real :: west_lon !< The longitude (or x-coordinate) of the first u-line
real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain
real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain
real :: Rad_Earth = 6.378e6 !< The radius of the planet [m].
real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m].
real :: len_lat !< The latitudinal (or y-coord) extent of physical domain
real :: len_lon !< The longitudinal (or x-coord) extent of physical domain
real :: Rad_Earth !< The radius of the planet [m]
real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m]
real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m]
end type ocean_grid_type

contains
Expand Down
6 changes: 4 additions & 2 deletions src/core/MOM_transcribe_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,8 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global
oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon
oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon
oG%Rad_Earth = dG%Rad_Earth ; oG%max_depth = dG%max_depth
oG%Rad_Earth = dG%Rad_Earth ; oG%Rad_Earth_L = dG%Rad_Earth_L
oG%max_depth = dG%max_depth

! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
call pass_var(oG%areaT, oG%Domain)
Expand Down Expand Up @@ -272,7 +273,8 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global
dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon
dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon
dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth
dG%Rad_Earth = oG%Rad_Earth ; dG%Rad_Earth_L = oG%Rad_Earth_L
dG%max_depth = oG%max_depth

! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
call pass_var(dG%areaT, dG%Domain)
Expand Down
Loading

0 comments on commit 459ecc0

Please sign in to comment.