From f03495fb65280dbed9a8eb9c89436298c53af54c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 1 Oct 2019 07:42:07 -0400 Subject: [PATCH] +Rescaled density units in MOM_T_S_init_from_Z Rescaled density units in MOM_temp_salt_initialize_from_Z for dimensional consistency testing, and added a new optional argument, eps_rho, to find_interfaces. All answers are bitwise identical, but there is a new optional argument to a public interface. --- src/initialization/MOM_state_initialization.F90 | 17 +++++++++++------ src/initialization/midas_vertmap.F90 | 17 +++++++++-------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index edd29d426e..961f965bde 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1969,6 +1969,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param integer :: kd, inconsistent integer :: nkd ! number of levels to use for regridding input arrays real :: eps_Z ! A negligibly thin layer thickness [Z ~> m]. + real :: eps_rho ! A negligibly small density difference [R ~> kg m-3]. real :: PI_180 ! for conversion from degrees to radians real, dimension(:,:), pointer :: shelf_area => NULL() @@ -1988,9 +1989,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param logical :: debug = .false. ! manually set this to true for verbose output ! data arrays - real, dimension(:), allocatable :: z_edges_in, z_in, Rb + real, dimension(:), allocatable :: z_edges_in, z_in + real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3] real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z - real, dimension(:,:,:), allocatable :: rho_z + real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi ! Interface heights [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)) :: nlevs real, dimension(SZI_(G)) :: press ! Pressures [Pa]. @@ -2115,6 +2117,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param !### Change this to GV%Angstrom_Z eps_z = 1.0e-10*US%m_to_Z + eps_rho = 1.0e-10*US%kg_m3_to_R ! Read input grid coordinates for temperature and salinity field ! in z-coordinate dataset. The file is REQUIRED to contain the @@ -2154,7 +2157,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call convert_temp_salt_for_TEOS10(temp_z, salt_z, press, G, kd, mask_z, eos) do k=1,kd ; do j=js,je - call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos) + call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, & + eos, scale=US%kg_m3_to_R) enddo ; enddo call pass_var(temp_z,G%Domain) @@ -2286,11 +2290,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! Rb contains the layer interface densities allocate(Rb(nz+1)) - do k=2,nz ; Rb(k) = 0.5*US%R_to_kg_m3*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo - Rb(1) = 0.0 ; Rb(nz+1) = US%R_to_kg_m3*( 2.0*GV%Rlay(nz) - GV%Rlay(nz-1) ) + do k=2,nz ; Rb(k) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo + Rb(1) = 0.0 ; Rb(nz+1) = 2.0*GV%Rlay(nz) - GV%Rlay(nz-1) zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, Rb, G%bathyT(is:ie,js:je), & - nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z) + nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z, & + eps_rho=eps_rho) if (correct_thickness) then call adjustEtaToFitBathymetry(G, GV, US, zi, h) diff --git a/src/initialization/midas_vertmap.F90 b/src/initialization/midas_vertmap.F90 index 9869877b68..f33d476cf0 100644 --- a/src/initialization/midas_vertmap.F90 +++ b/src/initialization/midas_vertmap.F90 @@ -559,12 +559,12 @@ function find_limited_slope(val, e, k) result(slope) end function find_limited_slope !> Find interface positions corresponding to density profile -function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z) result(zi) +function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z, eps_rho) result(zi) real, dimension(:,:,:), & - intent(in) :: rho !< potential density in z-space [R ~> kg m-3] + intent(in) :: rho !< potential density in z-space [kg m-3 or R ~> kg m-3] real, dimension(size(rho,3)), & - intent(in) :: zin !< Input data levels [Z ~> m or m]. - real, dimension(:), intent(in) :: Rb !< target interface densities [R ~> kg m-3] + intent(in) :: zin !< Input data levels [m or Z ~> m]. + real, dimension(:), intent(in) :: Rb !< target interface densities [kg m-3 or R ~> kg m-3] real, dimension(size(rho,1),size(rho,2)), & intent(in) :: depth !< ocean depth [Z ~> m]. real, dimension(size(rho,1),size(rho,2)), & @@ -573,7 +573,8 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps integer, optional, intent(in) :: nkml !< number of mixed layer pieces integer, optional, intent(in) :: nkbl !< number of buffer layer pieces real, optional, intent(in) :: hml !< mixed layer depth [Z ~> m]. - real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m or m]. + real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [m or Z ~> m]. + real, optional, intent(in) :: eps_rho !< A negligibly small density difference [kg m-3 or R ~> kg m-3]. real, dimension(size(rho,1),size(rho,2),size(Rb,1)) :: zi !< The returned interface, in the same units az zin. ! Local variables @@ -589,8 +590,8 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps integer :: n,i,j,k,l,nx,ny,nz,nt integer :: nlay,kk,nkml_,nkbl_ logical :: debug_ = .false. - real :: epsln_Z ! A negligibly thin layer thickness [Z ~> m]. - real :: epsln_rho ! A negligibly small density change [kg m-3]. + real :: epsln_Z ! A negligibly thin layer thickness [m or Z ~> m]. + real :: epsln_rho ! A negligibly small density change [kg m-3 or R ~> kg m-3]. real, parameter :: zoff=0.999 nlay=size(Rb)-1 @@ -606,7 +607,7 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps nkbl_ = 0 ; if (PRESENT(nkbl)) nkbl_ = max(0, nkbl) hml_ = 0.0 ; if (PRESENT(hml)) hml_ = hml epsln_Z = 1.0e-10 ; if (PRESENT(eps_z)) epsln_Z = eps_z - epsln_rho = 1.0e-10 + epsln_rho = 1.0e-10 ; if (PRESENT(eps_rho)) epsln_rho = eps_rho if (PRESENT(nlevs)) then nlevs_data(:,:) = nlevs(:,:)