Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into mono_N2_depth_units
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Jul 25, 2023
2 parents 60089b3 + d5ba107 commit 7012d26
Show file tree
Hide file tree
Showing 12 changed files with 472 additions and 177 deletions.
13 changes: 12 additions & 1 deletion config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ module MOM_surface_forcing_gfdl
!! from MOM_domains) to indicate the staggering of
!! the winds that are being provided in calls to
!! update_ocean_model.
logical :: use_temperature !< If true, temp and saln used as state variables
logical :: use_temperature !< If true, temp and saln used as state variables.
real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim].

real :: Rho0 !< Boussinesq reference density [R ~> kg m-3]
Expand Down Expand Up @@ -175,6 +175,7 @@ module MOM_surface_forcing_gfdl
real, pointer, dimension(:,:) :: t_flux =>NULL() !< sensible heat flux [W m-2]
real, pointer, dimension(:,:) :: q_flux =>NULL() !< specific humidity flux [kg m-2 s-1]
real, pointer, dimension(:,:) :: salt_flux =>NULL() !< salt flux [kg m-2 s-1]
real, pointer, dimension(:,:) :: excess_salt =>NULL() !< salt left behind by brine rejection [kg m-2 s-1]
real, pointer, dimension(:,:) :: lw_flux =>NULL() !< long wave radiation [W m-2]
real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation [W m-2]
real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation [W m-2]
Expand Down Expand Up @@ -304,6 +305,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)

if (associated(IOB%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed)

do j=js-2,je+2 ; do i=is-2,ie+2
fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j)
fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j)
Expand Down Expand Up @@ -576,6 +579,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
call check_mask_val_consistency(IOB%salt_flux(i-i0,j-j0), G%mask2dT(i,j), i, j, 'salt_flux', G)
enddo ; enddo
endif
if (associated(IOB%excess_salt)) then
do j=js,je ; do i=is,ie
fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*(kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0))
enddo ; enddo
endif

!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
!#CTRL# do j=js,je ; do i=is,ie
Expand Down Expand Up @@ -1729,6 +1737,9 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
if (associated(iobt%mass_berg)) then
chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks
endif
if (associated(iobt%excess_salt)) then
chks = field_chksum( iobt%excess_salt ) ; if (root) write(outunit,100) 'iobt%excess_salt ', chks
endif
100 FORMAT(" CHECKSUM::",A20," = ",Z20)

call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')
Expand Down
44 changes: 31 additions & 13 deletions config_src/external/drifters/MOM_particles.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,28 +11,29 @@ module MOM_particles_mod
implicit none ; private

public particles, particles_run, particles_init, particles_save_restart, particles_end
public particles_to_k_space, particles_to_z_space

contains

!> Initializes particles container "parts"
subroutine particles_init(parts, Grid, Time, dt, u, v)
subroutine particles_init(parts, Grid, Time, dt, u, v, h)
! Arguments
type(particles), pointer, intent(out) :: parts !< Container for all types and memory
type(ocean_grid_type), target, intent(in) :: Grid !< Grid type from parent model
type(time_type), intent(in) :: Time !< Time type from parent model
real, intent(in) :: dt !< particle timestep [s]
real, dimension(:,:,:), intent(in) :: u !< Zonal velocity field [m s-1]
real, dimension(:,:,:), intent(in) :: v !< Meridional velocity field [m s-1]

real, intent(in) :: dt !< particle timestep in seconds [T ~> s]
real, dimension(:,:,:),intent(in) :: u !< Zonal velocity field [L T-1 ~> m s-1]
real, dimension(:,:,:),intent(in) :: v !< Meridional velocity field [L T-1 ~> m s-1]
real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2]
end subroutine particles_init

!> The main driver the steps updates particles
subroutine particles_run(parts, time, uo, vo, ho, tv, stagger)
! Arguments
type(particles), pointer :: parts !< Container for all types and memory
type(time_type), intent(in) :: time !< Model time
real, dimension(:,:,:), intent(in) :: uo !< Ocean zonal velocity [m s-1]
real, dimension(:,:,:), intent(in) :: vo !< Ocean meridional velocity [m s-1]
real, dimension(:,:,:), intent(in) :: uo !< Ocean zonal velocity [L T-1 ~>m s-1]
real, dimension(:,:,:), intent(in) :: vo !< Ocean meridional velocity [L T-1~> m s-1]
real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< structure containing pointers to available thermodynamic fields
integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered
Expand All @@ -41,21 +42,38 @@ end subroutine particles_run


!>Save particle locations (and sometimes other vars) to restart file
subroutine particles_save_restart(parts, temp, salt)
subroutine particles_save_restart(parts, h, temp, salt)
! Arguments
type(particles), pointer :: parts !< Container for all types and memory
real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature
real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity
real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2]
real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature [C ~> degC]
real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity [S ~> ppt]

end subroutine particles_save_restart

!> Deallocate all memory and disassociated pointer
subroutine particles_end(parts, temp, salt)
subroutine particles_end(parts, h, temp, salt)
! Arguments
type(particles), pointer :: parts !< Container for all types and memory
real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature
real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity
real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2]
real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature [C ~> degC]
real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity [S ~> ppt]

end subroutine particles_end

subroutine particles_to_k_space(parts, h)
! Arguments
type(particles), pointer :: parts !< Container for all types and memory
real, dimension(:,:,:),intent(in) :: h !< Thickness of layers [H ~> m or kg m-2]

end subroutine particles_to_k_space


subroutine particles_to_z_space(parts, h)
! Arguments
type(particles), pointer :: parts !< Container for all types and memory
real, dimension(:,:,:),intent(in) :: h !< Thickness of layers [H ~> m or kg m-2]

end subroutine particles_to_z_space

end module MOM_particles_mod
9 changes: 9 additions & 0 deletions docs/zotero.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2738,3 +2738,12 @@ @article{kraus1967
journal = {Tellus}
}

@article{Nguyen2009,
doi = {10.1029/2008JC005121},
year = {2009},
journal = {JGR Oceans},
volume = {114},
author = {A. T. Nguyen and D. Menemenlis and R. Kwok},
title = {Improved modeling of the Arctic halocline with a subgrid-scale brine rejection parameterization},
pages = {C11014}
}
17 changes: 13 additions & 4 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ module MOM
use MOM_offline_main, only : offline_advection_layer, offline_transport_end
use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf
use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end

use MOM_particles_mod, only : particles_to_k_space, particles_to_z_space
implicit none ; private

#include <MOM_memory.h>
Expand Down Expand Up @@ -1541,6 +1541,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &

call preAle_tracer_diagnostics(CS%tracer_Reg, G, GV)

if (CS%use_particles) then
call particles_to_z_space(CS%particles, h)
endif

if (CS%debug) then
call MOM_state_chksum("Pre-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US, omit_corners=.true.)
call hchksum(tv%T,"Pre-ALE T", G%HI, haloshift=1, omit_corners=.true., scale=US%C_to_degC)
Expand Down Expand Up @@ -1588,6 +1592,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
call cpu_clock_end(id_clock_ALE)
endif ! endif for the block "if ( CS%use_ALE_algorithm )"


if (CS%use_particles) then
call particles_to_k_space(CS%particles, h)
endif

dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo)
call create_group_pass(pass_uv_T_S_h, u, v, G%Domain, halo=dynamics_stencil)
if (associated(tv%T)) &
Expand Down Expand Up @@ -3232,7 +3241,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
G => CS%G ; GV => CS%GV ; US => CS%US

if (CS%use_particles) then
call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v)
call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v, CS%h)
endif

! Write initial conditions
Expand Down Expand Up @@ -3935,7 +3944,7 @@ subroutine save_MOM6_internal_state(CS, dirs, time, stamp_time)

! Could call save_restart(CS%restart_CSp) here

if (CS%use_particles) call particles_save_restart(CS%particles)
if (CS%use_particles) call particles_save_restart(CS%particles, CS%h)

end subroutine save_MOM6_internal_state

Expand Down Expand Up @@ -3978,7 +3987,7 @@ subroutine MOM_end(CS)
endif

if (CS%use_particles) then
call particles_end(CS%particles)
call particles_end(CS%particles, CS%h)
deallocate(CS%particles)
endif

Expand Down
8 changes: 4 additions & 4 deletions src/core/MOM_continuity_PPM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,9 +378,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are
dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif

if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)) &
if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
du_max_CFL(I) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(I,k)
if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)) &
if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
du_min_CFL(I) = -(dx_E*CFL_dt + u(I,j,k)) / visc_rem(I,k)
enddo ; enddo
endif
Expand Down Expand Up @@ -1201,9 +1201,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac
dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)) &
if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
dv_max_CFL(i) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(i,k)
if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)) &
if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
dv_min_CFL(i) = -(dy_N*CFL_dt + v(i,J,k)) / visc_rem(i,k)
enddo ; enddo
endif
Expand Down
16 changes: 11 additions & 5 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,10 @@ module MOM_forcing_type
real, pointer, dimension(:,:) :: &
salt_flux => NULL(), & !< net salt flux into the ocean [R Z T-1 ~> kgSalt m-2 s-1]
salt_flux_in => NULL(), & !< salt flux provided to the ocean from coupler [R Z T-1 ~> kgSalt m-2 s-1]
salt_flux_added => NULL() !< additional salt flux from restoring or flux adjustment before adjustment
salt_flux_added => NULL(), & !< additional salt flux from restoring or flux adjustment before adjustment
!! to net zero [R Z T-1 ~> kgSalt m-2 s-1]
salt_left_behind => NULL() !< salt left in ocean at the surface from brine rejection
!! [R Z T-1 ~> kgSalt m-2 s-1]

! applied surface pressure from other component models (e.g., atmos, sea ice, land ice)
real, pointer, dimension(:,:) :: p_surf_full => NULL()
Expand Down Expand Up @@ -746,15 +748,15 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
endif

! Salt fluxes
Net_salt(i) = 0.0
if (do_NSR) Net_salt_rate(i) = 0.0
net_salt(i) = 0.0
if (do_NSR) net_salt_rate(i) = 0.0
! Convert salt_flux from kg (salt)/(m^2 * s) to
! Boussinesq: (ppt * m)
! non-Bouss: (g/m^2)
if (associated(fluxes%salt_flux)) then
Net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
!Repeat above code for 'rate' term
if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
if (do_NSR) net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
endif

! Diagnostics follow...
Expand Down Expand Up @@ -1947,6 +1949,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
diag%axesT1,Time,'Salt flux into ocean at surface due to restoring or flux adjustment', &
units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)

handles%id_saltFluxAdded = register_diag_field('ocean_model', 'salt_left_behind', &
diag%axesT1,Time,'Salt left in ocean at surface due to ice formation', &
units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)

handles%id_saltFluxGlobalAdj = register_scalar_field('ocean_model', &
'salt_flux_global_restoring_adjustment', Time, diag, &
'Adjustment needed to balance net global salt flux into ocean at surface', &
Expand Down
38 changes: 23 additions & 15 deletions src/diagnostics/MOM_sum_output.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -510,24 +510,18 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
do k=1,nz ; vol_lay(k) = (US%m_to_L**2*GV%H_to_Z/GV%H_to_kg_m2)*mass_lay(k) ; enddo
else
tmp1(:,:,:) = 0.0
if (CS%do_APE_calc) then
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)

if (CS%do_APE_calc) then
call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = US%Z_to_m*US%L_to_m**2*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
enddo ; enddo ; enddo
vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay)
do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2 * vol_lay(k) ; enddo
else
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)
do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2*US%kg_m3_to_R * (mass_lay(k) / GV%Rho0) ; enddo
endif
endif ! Boussinesq

Expand Down Expand Up @@ -643,7 +637,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
if (GV%Boussinesq) then
do j=js,je ; do i=is,ie
hbelow = 0.0
do k=nz,1,-1
do K=nz,1,-1
hbelow = hbelow + h(i,j,k) * GV%H_to_Z
hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref))
hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref)
Expand All @@ -652,14 +646,28 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
(hint * hint - hbot * hbot)
enddo
enddo ; enddo
else
elseif (GV%semi_Boussinesq) then
do j=js,je ; do i=is,ie
do k=nz,1,-1
do K=nz,1,-1
hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * &
(hint * hint - hbot * hbot)
(hint * hint - hbot * hbot)
enddo
enddo ; enddo
else
do j=js,je ; do i=is,ie
do K=nz,2,-1
hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,K) = (0.25 * PE_scale_factor * areaTm(i,j) * &
((GV%Rlay(k)+GV%Rlay(k-1))*GV%g_prime(K))) * &
(hint * hint - hbot * hbot)
enddo
hint = Z_0APE(1) + eta(i,j,1) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(1) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,1) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * &
(hint * hint - hbot * hbot)
enddo ; enddo
endif

Expand Down
Loading

0 comments on commit 7012d26

Please sign in to comment.