Skip to content

Commit

Permalink
Rescaled density units in MOM_regridding.F90
Browse files Browse the repository at this point in the history 
  Rescaled density units in MOM_regridding.F90, including using the new optional
arguments to init_coord_hycom, init_coord_rho, and init_coord_slight to rescale
densities in those modules as well.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Oct 4, 2019
1 parent e9b36cc commit 1175786
Showing 1 changed file with 33 additions and 27 deletions.
60 changes: 33 additions & 27 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ module MOM_regridding

!> This array is set by function setCoordinateResolution()
!! It contains the "resolution" or delta coordinate of the target
!! coorindate. It has the units of the target coordinate, e.g.
!! coordinate. It has the units of the target coordinate, e.g.
!! [Z ~> m] for z*, non-dimensional for sigma, etc.
real, dimension(:), allocatable :: coordinateResolution

Expand All @@ -56,9 +56,9 @@ module MOM_regridding

!> This array is set by function set_target_densities()
!! This array is the nominal coordinate of interfaces and is the
!! running sum of coordinateResolution. i.e.
!! running sum of coordinateResolution, in [R ~> kg m-3]. i.e.
!! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
!! It is only used in "rho" mode.
!! It is only used in "rho", "SLight" or "Hycom" mode.
real, dimension(:), allocatable :: target_density

!> A flag to indicate that the target_density arrays has been filled with data.
Expand Down Expand Up @@ -199,8 +199,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
integer :: nz_fixed_sfc, k, nzf(4)
real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be
! [m] or [Z ~> m] or [H ~> m or kg m-2] or [kg m-3] or other units.
real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
! units depending on the coordinate
Expand Down Expand Up @@ -310,13 +310,9 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
'Unable to interpret "'//trim(string)//'".')
endif
allocate(dz(ke))
if (ke==1) then
dz(:) = uniformResolution(ke, coord_mode, tmpReal, US%R_to_kg_m3*GV%Rlay(1), US%R_to_kg_m3*GV%Rlay(1))
else
dz(:) = uniformResolution(ke, coord_mode, tmpReal, &
US%R_to_kg_m3*(GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2))), &
US%R_to_kg_m3*(GV%Rlay(ke)+0.5*(GV%Rlay(ke)-GV%Rlay(ke-1))) )
endif
dz(:) = uniformResolution(ke, coord_mode, tmpReal, &
US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), &
US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) )
if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
trim(message), units=trim(coord_units))
elseif (trim(string)=='PARAM') then
Expand Down Expand Up @@ -469,13 +465,15 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
allocate( CS%coordinateResolution(CS%nk) ); CS%coordinateResolution(:) = -1.E30
if (state_dependent(CS%regridding_scheme)) then
! Target values
allocate( CS%target_density(CS%nk+1) ); CS%target_density(:) = -1.E30
allocate( CS%target_density(CS%nk+1) ); CS%target_density(:) = -1.E30*US%kg_m3_to_R
endif

if (allocated(dz)) then
if ((coordinateMode(coord_mode) == REGRIDDING_SIGMA) .or. &
(coordinateMode(coord_mode) == REGRIDDING_RHO)) then
if (coordinateMode(coord_mode) == REGRIDDING_SIGMA) then
call setCoordinateResolution(dz, CS, scale=1.0)
elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
call setCoordinateResolution(dz, CS, scale=US%kg_m3_to_R)
CS%coord_scale = US%R_to_kg_m3
elseif (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
call setCoordinateResolution(dz, CS, scale=GV%m_to_H)
CS%coord_scale = GV%H_to_m
Expand All @@ -486,18 +484,18 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
endif

if (allocated(rho_target)) then
call set_target_densities(CS, rho_target)
call set_target_densities(CS, US%kg_m3_to_R*rho_target)
deallocate(rho_target)

! \todo This line looks like it would overwrite the target densities set just above?
elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
call set_target_densities_from_GV(GV, US, CS)
call log_param(param_file, mdl, "!TARGET_DENSITIES", CS%target_density, &
call log_param(param_file, mdl, "!TARGET_DENSITIES", US%R_to_kg_m3*CS%target_density(:), &
'RHO target densities for interfaces', units=coordinateUnits(coord_mode))
endif

! initialise coordinate-specific control structure
call initCoord(CS, GV, coord_mode)
call initCoord(CS, GV, US, coord_mode)

if (main_parameters .and. coord_is_state_dependent) then
call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpReal, &
Expand Down Expand Up @@ -1947,12 +1945,13 @@ end function uniformResolution

!> Initialize the coordinate resolutions by calling the appropriate initialization
!! routine for the specified coordinate mode.
subroutine initCoord(CS, GV, coord_mode)
subroutine initCoord(CS, GV, US, coord_mode)
type(regridding_CS), intent(inout) :: CS !< Regridding control structure
character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode.
!! See the documenttion for regrid_consts
!! for the recognized values.
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

select case (coordinateMode(coord_mode))
case (REGRIDDING_ZSTAR)
Expand All @@ -1962,11 +1961,14 @@ subroutine initCoord(CS, GV, coord_mode)
case (REGRIDDING_SIGMA)
call init_coord_sigma(CS%sigma_CS, CS%nk, CS%coordinateResolution)
case (REGRIDDING_RHO)
call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS)
call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS, &
rho_scale=US%kg_m3_to_R)
case (REGRIDDING_HYCOM1)
call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, CS%interp_CS)
call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, &
CS%interp_CS, rho_scale=US%kg_m3_to_R)
case (REGRIDDING_SLIGHT)
call init_coord_slight(CS%slight_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS, GV%m_to_H)
call init_coord_slight(CS%slight_CS, CS%nk, CS%ref_pressure, CS%target_density, &
CS%interp_CS, GV%m_to_H, rho_scale=US%kg_m3_to_R)
case (REGRIDDING_ADAPTIVE)
call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H)
end select
Expand Down Expand Up @@ -1999,8 +2001,8 @@ subroutine set_target_densities_from_GV( GV, US, CS )
integer :: k, nz

nz = CS%nk
CS%target_density(1) = US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(2)))
CS%target_density(nz+1) = US%R_to_kg_m3*(GV%Rlay(nz) + 0.5*(GV%Rlay(nz)-GV%Rlay(nz-1)))
CS%target_density(1) = (GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(2)))
CS%target_density(nz+1) = (GV%Rlay(nz) + 0.5*(GV%Rlay(nz)-GV%Rlay(nz-1)))
do k = 2,nz
CS%target_density(k) = CS%target_density(k-1) + CS%coordinateResolution(k)
enddo
Expand All @@ -2011,7 +2013,7 @@ end subroutine set_target_densities_from_GV
!> Set target densities based on vector of interface values
subroutine set_target_densities( CS, rho_int )
type(regridding_CS), intent(inout) :: CS !< Regridding control structure
real, dimension(CS%nk+1), intent(in) :: rho_int !< Interface densities
real, dimension(CS%nk+1), intent(in) :: rho_int !< Interface densities [R ~> kg m-3]

if (size(CS%target_density)/=size(rho_int)) then
call MOM_error(FATAL, "set_target_densities inconsistent args!")
Expand Down Expand Up @@ -2124,7 +2126,11 @@ function getCoordinateInterfaces( CS, undo_scaling )
call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//&
'target densities not set!')

getCoordinateInterfaces(:) = CS%target_density(:)
if (unscale) then
getCoordinateInterfaces(:) = CS%coord_scale * CS%target_density(:)
else
getCoordinateInterfaces(:) = CS%target_density(:)
endif
else
if (unscale) then
getCoordinateInterfaces(1) = 0.
Expand Down Expand Up @@ -2402,7 +2408,7 @@ end subroutine dz_function1
integer function rho_function1( string, rho_target )
character(len=*), intent(in) :: string !< String with list of parameters in form
!! dz_min, H_total, power, precision
real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities
real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
! Local variables
integer :: nki, k, nk
real :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
Expand Down

0 comments on commit 1175786

Please sign in to comment.