From f272edc30bda516ddd86b2fff3736bbf2837636a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 15:37:30 -0500 Subject: [PATCH] +Rescale depth in Kelvin_initialize_topography Added an optional unit_scale_type argument to Kelvin_initialize_topography, and if it is present the unit conversion of the topography from m to Z occurs within this initialization routine. All answers are bitwise identical. --- src/user/Kelvin_initialization.F90 | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 341f6e99f6..1fd9edb3ea 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -110,22 +110,28 @@ end subroutine Kelvin_OBC_end ! ----------------------------------------------------------------------------- !> This subroutine sets up the Kelvin topography and land mask -subroutine Kelvin_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Ocean bottom depth in m - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m +subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m or Z if US is present + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum model depth in the units of D + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + ! Local variables character(len=40) :: mdl = "Kelvin_initialize_topography" ! This subroutine's name. - real :: min_depth ! The minimum and maximum depths in m. + real :: m_to_Z ! A dimensional rescaling factor. + real :: min_depth ! The minimum and maximum depths in Z. real :: PI ! 3.1415... real :: coast_offset1, coast_offset2, coast_angle, right_angle integer :: i, j call MOM_mesg(" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5) + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0) + "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", coast_offset1, & default=100.0, do_not_log=.true.) call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, & @@ -137,17 +143,17 @@ subroutine Kelvin_initialize_topography(D, G, param_file, max_depth) right_angle = 2 * atan(1.0) do j=G%jsc,G%jec ; do i=G%isc,G%iec - D(i,j)=max_depth + D(i,j) = max_depth ! Southern side if ((G%geoLonT(i,j) - G%west_lon > coast_offset1) .AND. & (atan2(G%geoLatT(i,j) - G%south_lat + coast_offset2, & G%geoLonT(i,j) - G%west_lon - coast_offset1) < coast_angle)) & - D(i,j)=0.5*min_depth + D(i,j) = 0.5*min_depth ! Northern side if ((G%geoLonT(i,j) - G%west_lon < G%len_lon - coast_offset1) .AND. & (atan2(G%len_lat + G%south_lat + coast_offset2 - G%geoLatT(i,j), & G%len_lon + G%west_lon - coast_offset1 - G%geoLonT(i,j)) < coast_angle)) & - D(i,j)=0.5*min_depth + D(i,j) = 0.5*min_depth if (D(i,j) > max_depth) D(i,j) = max_depth if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth