Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into document_coord_code
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Feb 9, 2024
2 parents 40c370e + 4757bcb commit f6a663b
Show file tree
Hide file tree
Showing 16 changed files with 541 additions and 223 deletions.
36 changes: 24 additions & 12 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1260,16 +1260,17 @@ end subroutine mask_near_bottom_vel
!! h_dst must be dimensioned as a model array with GV%ke layers while h_src can
!! have an arbitrary number of layers specified by nk_src.
subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap, &
answers_2018, answer_date )
answers_2018, answer_date, h_neglect, h_neglect_edge)
type(remapping_CS), intent(in) :: CS !< Remapping control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
integer, intent(in) :: nk_src !< Number of levels on source grid
real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid
!! [H ~> m or kg m-2]
!! [H ~> m or kg m-2] or other units
!! if H_neglect is provided
real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid, in arbitrary units [A]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid
!! [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid in the
!! same units as h_src, often [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid, in the same
!! arbitrary units as s_src [A]
logical, optional, intent(in) :: all_cells !< If false, only reconstruct for
Expand All @@ -1283,10 +1284,16 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
!! use more robust forms of the same expressions.
integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
!! for remapping
real, optional, intent(in) :: h_neglect !< A negligibly small thickness used in
!! remapping cell reconstructions, in the same
!! units as h_src, often [H ~> m or kg m-2]
real, optional, intent(in) :: h_neglect_edge !< A negligibly small thickness used in
!! remapping edge value calculations, in the same
!! units as h_src, often [H ~> m or kg m-2]
! Local variables
integer :: i, j, k, n_points
real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
real :: h_neg, h_neg_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap

ignore_vanished_layers = .false.
Expand All @@ -1297,12 +1304,17 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
use_2018_remap = .true. ; if (present(answers_2018)) use_2018_remap = answers_2018
if (present(answer_date)) use_2018_remap = (answer_date < 20190101)

if (.not.use_2018_remap) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
if (present(h_neglect)) then
h_neg = h_neglect
h_neg_edge = h_neg ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
if (.not.use_2018_remap) then
h_neg = GV%H_subroundoff ; h_neg_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neg = GV%m_to_H*1.0e-30 ; h_neg_edge = GV%m_to_H*1.0e-10
else
h_neg = GV%kg_m2_to_H*1.0e-30 ; h_neg_edge = GV%kg_m2_to_H*1.0e-10
endif
endif

!$OMP parallel do default(shared) firstprivate(n_points,dx)
Expand All @@ -1318,10 +1330,10 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
if (use_remapping_core_w) then
call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx )
call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
GV%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge)
GV%ke, dx, s_dst(i,j,:), h_neg, h_neg_edge)
else
call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge)
GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neg, h_neg_edge)
endif
else
s_dst(i,j,:) = 0.
Expand Down
60 changes: 46 additions & 14 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,7 @@ module MOM_regridding
public getCoordinateResolution, getCoordinateInterfaces
public getCoordinateUnits, getCoordinateShortName, getStaticThickness
public DEFAULT_COORDINATE_MODE
public set_h_neglect, set_dz_neglect
public get_zlike_CS, get_sigma_CS, get_rho_CS

!> Documentation for coordinate options
Expand Down Expand Up @@ -1430,13 +1431,7 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS,
#endif
logical :: ice_shelf

if (CS%remap_answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
endif
h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge)

nz = GV%ke
ice_shelf = present(frac_shelf_h)
Expand Down Expand Up @@ -1592,13 +1587,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface,
logical :: ice_shelf
integer :: i, j, k, nki

if (CS%remap_answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
endif
h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge)

if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_HyCOM1 : "//&
"Target densities must be set before build_grid_HyCOM1 is called.")
Expand Down Expand Up @@ -2133,6 +2122,49 @@ subroutine write_regrid_file( CS, GV, filepath )

end subroutine write_regrid_file

!> Set appropriate values for the negligible thicknesses used for remapping based on an answer date.
function set_h_neglect(GV, remap_answer_date, h_neglect_edge) result(h_neglect)
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
!! for remapping. Values below 20190101 recover the
!! remapping answers from 2018. Higher values use more
!! robust forms of the same remapping algorithms.
real, intent(out) :: h_neglect_edge !< A negligibly small thickness used in
!! remapping edge value calculations [H ~> m or kg m-2]
real :: h_neglect !< A negligibly small thickness used in
!! remapping cell reconstructions [H ~> m or kg m-2]

if (remap_answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
endif
end function set_h_neglect

!> Set appropriate values for the negligible vertical layer extents used for remapping based on an answer date.
function set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) result(dz_neglect)
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
!! for remapping. Values below 20190101 recover the
!! remapping answers from 2018. Higher values use more
!! robust forms of the same remapping algorithms.
real, intent(out) :: dz_neglect_edge !< A negligibly small vertical layer extent
!! used in remapping edge value calculations [Z ~> m]
real :: dz_neglect !< A negligibly small vertical layer extent
!! used in remapping cell reconstructions [Z ~> m]

if (remap_answer_date >= 20190101) then
dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff
elseif (GV%Boussinesq) then
dz_neglect = US%m_to_Z*1.0e-30 ; dz_neglect_edge = US%m_to_Z*1.0e-10
else
dz_neglect = GV%kg_m2_to_H * (GV%H_to_m*US%m_to_Z) * 1.0e-30
dz_neglect_edge = GV%kg_m2_to_H * (GV%H_to_m*US%m_to_Z) * 1.0e-10
endif
end function set_dz_neglect

!------------------------------------------------------------------------------
!> Query the fixed resolution data
Expand Down
20 changes: 11 additions & 9 deletions src/diagnostics/MOM_spatial_means.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,11 +211,13 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale)
! Local variables
! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: tmpForSumming ! An unscaled cell integral [a m3]
real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: weight ! The volume of each cell, used as a weight [m3]
real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: tmpForSumming ! An unscaled cell integral [a m3] or [a kg]
real, dimension(G%isc:G%iec,G%jsc:G%jec,SZK_(GV)) :: weight ! The volume or mass of each cell, depending on
! whether the model is Boussinesq, used as a weight [m3] or [kg]
type(EFP_type), dimension(2*SZK_(GV)) :: laysums
real, dimension(SZK_(GV)) :: global_temp_scalar ! The global integral of the tracer in each layer [a m3]
real, dimension(SZK_(GV)) :: global_weight_scalar ! The global integral of the volume of each layer [m3]
real, dimension(SZK_(GV)) :: global_temp_scalar ! The global integral of the tracer in each layer [a m3] or [a kg]
real, dimension(SZK_(GV)) :: global_weight_scalar ! The global integral of the volume or mass of each
! layer [m3] or [kg]
real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1]
integer :: i, j, k, is, ie, js, je, nz
Expand All @@ -226,7 +228,7 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale)
tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.

do k=1,nz ; do j=js,je ; do i=is,ie
weight(i,j,k) = (GV%H_to_m * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
weight(i,j,k) = (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
tmpForSumming(i,j,k) = scalefac * var(i,j,k) * weight(i,j,k)
enddo ; enddo ; enddo

Expand Down Expand Up @@ -262,9 +264,9 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale)
! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1]
real :: weight_here ! The volume of a grid cell [m3]
real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! The volume integral of the variable in a column [a m3]
real, dimension(SZI_(G),SZJ_(G)) :: sum_weight ! The volume of each column of water [m3]
real :: weight_here ! The volume or mass of a grid cell [m3] or [kg]
real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! The volume integral of the variable in a column [a m3] or [a kg]
real, dimension(SZI_(G),SZJ_(G)) :: sum_weight ! The volume or mass of each column of water [m3] or [kg]
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

Expand All @@ -273,7 +275,7 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale)
tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0.

do k=1,nz ; do j=js,je ; do i=is,ie
weight_here = (GV%H_to_m * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
weight_here = (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
tmpForSumming(i,j) = tmpForSumming(i,j) + scalefac * var(i,j,k) * weight_here
sum_weight(i,j) = sum_weight(i,j) + weight_here
enddo ; enddo ; enddo
Expand Down
9 changes: 4 additions & 5 deletions src/diagnostics/MOM_wave_speed.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,10 +197,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N
enddo ; enddo ; enddo
endif

g_Rho0 = GV%g_Earth*GV%H_to_Z / GV%Rho0
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
H_to_pres = GV%H_to_RZ * GV%g_Earth
! Note that g_Rho0 = H_to_pres / GV%Rho0**2
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
if (.not.nonBous) g_Rho0 = GV%g_Earth*GV%H_to_Z / GV%Rho0
use_EOS = associated(tv%eqn_of_state)

better_est = CS%better_cg1_est
Expand Down Expand Up @@ -900,9 +900,9 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo
endif

g_Rho0 = GV%g_Earth * GV%H_to_Z / GV%Rho0
H_to_pres = GV%H_to_RZ * GV%g_Earth
nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
H_to_pres = GV%H_to_RZ * GV%g_Earth
if (.not.nonBous) g_Rho0 = GV%g_Earth * GV%H_to_Z / GV%Rho0
use_EOS = associated(tv%eqn_of_state)

if (CS%c1_thresh < 0.0) &
Expand Down Expand Up @@ -1057,7 +1057,6 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
enddo
endif
endif
cg1_est = g_Rho0 * drxh_sum
else ! Not use_EOS
drxh_sum = 0.0 ; dSpVxh_sum = 0.0
if (better_est) then
Expand Down
2 changes: 1 addition & 1 deletion src/equation_of_state/MOM_EOS.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1703,7 +1703,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
do k=1,kd ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec
if (mask_z(i,j,k) >= 1.0) then
S(i,j,k) = Sref_Sprac * S(i,j,k)
T(i,j,k) = EOS%degC_to_C*poTemp_to_consTemp(EOS%S_to_ppt*S(i,j,k), EOS%S_to_ppt*T(i,j,k))
T(i,j,k) = EOS%degC_to_C*poTemp_to_consTemp(EOS%C_to_degC*T(i,j,k), EOS%S_to_ppt*S(i,j,k))
endif
enddo ; enddo ; enddo
end subroutine convert_temp_salt_for_TEOS10
Expand Down
17 changes: 12 additions & 5 deletions src/framework/MOM_diag_remap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -658,8 +658,15 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field [nondim]

! Local variables
real, dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff
real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages
real :: volume(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area [m2], volume [m3] or mass [kg] of each cell.
real :: stuff(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area, volume or mass-weighted integral of the
! field being averaged in each cell, in [m2 A], [m3 A] or [kg A],
! depending on the weighting for the averages and whether the
! model makes the Boussinesq approximation.
real, dimension(size(field, 3)) :: vol_sum ! The global sum of the areas [m2], volumes [m3] or mass [kg]
! in the cells that used in the weighted averages.
real, dimension(size(field, 3)) :: stuff_sum ! The global sum of the weighted field in all cells, in
! [A m2], [A m3] or [A kg]
type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer
real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2]
integer :: i, j, k, nz
Expand Down Expand Up @@ -688,7 +695,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
I1 = i - G%isdB + 1
height = 0.5 * (h(i,j,k) + h(i+1,j,k))
volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) &
* (GV%H_to_m * height) * G%mask2dCu(I,j)
* (GV%H_to_MKS * height) * G%mask2dCu(I,j)
stuff(I,j,k) = volume(I,j,k) * field(I1,j,k)
enddo ; enddo
endif
Expand Down Expand Up @@ -717,7 +724,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
J1 = J - G%jsdB + 1
height = 0.5 * (h(i,j,k) + h(i,j+1,k))
volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) &
* (GV%H_to_m * height) * G%mask2dCv(i,J)
* (GV%H_to_MKS * height) * G%mask2dCv(i,J)
stuff(i,J,k) = volume(i,J,k) * field(i,J1,k)
enddo ; enddo
endif
Expand Down Expand Up @@ -748,7 +755,7 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i
else ! Intensive
do j=G%jsc, G%jec ; do i=G%isc, G%iec
volume(i,j,k) = (G%US%L_to_m**2 * G%areaT(i,j)) &
* (GV%H_to_m * h(i,j,k)) * G%mask2dT(i,j)
* (GV%H_to_MKS * h(i,j,k)) * G%mask2dT(i,j)
stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
enddo ; enddo
endif
Expand Down
Loading

0 comments on commit f6a663b

Please sign in to comment.