From a1653f6b5e6178882e9cfec0ceba1d6b599b35c0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 15:36:45 -0500 Subject: [PATCH 01/10] +Rescale depth in DOME_initialize_topography Added an optional unit_scale_type argument to DOME_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/DOME_initialization.F90 | 35 ++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 40e6e422df..101e52eb30 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -32,16 +32,19 @@ module DOME_initialization ! ----------------------------------------------------------------------------- !> This subroutine sets up the DOME topography -subroutine DOME_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine DOME_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m - - real :: min_depth ! The minimum and maximum depths in m. -! This include declares and sets the variable "version". -#include "version_variable.h" + 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 + real :: m_to_Z ! A dimensional rescaling factor. + real :: min_depth ! The minimum and maximum depths in Z. + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "DOME_initialize_topography" ! This subroutine's name. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -49,22 +52,24 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth) call MOM_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5) + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call log_version(param_file, mdl, version, "") 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) do j=js,je ; do i=is,ie if (G%geoLatT(i,j) < 600.0) then if (G%geoLatT(i,j) < 300.0) then - D(i,j)=max_depth + D(i,j) = max_depth else - D(i,j)=max_depth-10.0*(G%geoLatT(i,j)-300.0) + D(i,j) = max_depth - 10.0*m_to_Z * (G%geoLatT(i,j)-300.0) endif else - if ((G%geoLonT(i,j) > 1000.0).AND.(G%geoLonT(i,j) < 1100.0)) then - D(i,j)=600.0 + if ((G%geoLonT(i,j) > 1000.0) .AND. (G%geoLonT(i,j) < 1100.0)) then + D(i,j) = 600.0*m_to_Z else - D(i,j)=0.5*min_depth + D(i,j) = 0.5*min_depth endif endif From 771e44f218d8ad709599f38852362298741865f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 15:37:10 -0500 Subject: [PATCH 02/10] +Rescale depth in ISOMIP_initialize_topography Added an optional unit_scale_type argument to ISOMIP_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/ISOMIP_initialization.F90 | 49 +++++++++++++++++------------- 1 file changed, 28 insertions(+), 21 deletions(-) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index c5d55640a6..6a3bf8007e 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -36,23 +36,27 @@ module ISOMIP_initialization contains !> Initialization of topography for the ISOMIP configuration -subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine ISOMIP_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + 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 - real :: min_depth ! The minimum and maximum depths in m. + real :: min_depth ! The minimum and maximum depths in Z. + real :: m_to_Z ! A dimensional rescaling factor. ! The following variables are used to set up the bathymetry in the ISOMIP example. real :: bmax ! max depth of bedrock topography real :: b0,b2,b4,b6 ! first, second, third and fourth bedrock topography coeff - real :: xbar ! characteristic along-flow lenght scale of the bedrock - real :: dc ! depth of the trough compared with side walls + real :: xbar ! characteristic along-flow lenght scale of the bedrock + real :: dc ! depth of the trough compared with side walls in Z real :: fc ! characteristic width of the side walls of the channel real :: wc ! half-width of the trough real :: ly ! domain width (across ice flow) - real :: bx, by, xtil ! dummy vatiables + real :: bx, by ! dummy vatiables in Z + real :: xtil ! dummy vatiable logical :: is_2D ! If true, use 2D setup ! This include declares and sets the variable "version". #include "version_variable.h" @@ -63,15 +67,18 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) call MOM_mesg(" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5) + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call log_version(param_file, mdl, version, "") 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, "ISOMIP_2D",is_2D,'If true, use a 2D setup.', default=.false.) ! The following variables should be transformed into runtime parameters? - bmax=720.0; b0=-150.0; b2=-728.8; b4=343.91; b6=-50.57 - xbar=300.0E3; dc=500.0; fc=4.0E3; wc=24.0E3; ly=80.0E3 - bx = 0.0; by = 0.0; xtil = 0.0 + bmax = 720.0*m_to_Z ; dc = 500.0*m_to_Z + b0 = -150.0*m_to_Z ; b2 = -728.8*m_to_Z ; b4 = 343.91*m_to_Z ; b6 = -50.57*m_to_Z + xbar = 300.0e3 ; fc = 4.0e3 ; wc = 24.0e3 ; ly = 80.0e3 + bx = 0.0 ; by = 0.0 ; xtil = 0.0 if (is_2D) then @@ -79,15 +86,15 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) ! 2D setup xtil = G%geoLonT(i,j)*1.0e3/xbar !xtil = 450*1.0e3/xbar - bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) ! slice at y = 40 km - by = (dc/(1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc))) + by = (dc / (1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + & + (dc / (1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc))) - D(i,j) = -max(bx+by,-bmax) + D(i,j) = -max(bx+by, -bmax) if (D(i,j) > max_depth) D(i,j) = max_depth if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth enddo ; enddo @@ -105,11 +112,11 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) xtil = G%geoLonT(i,j)*1.0e3/xbar - bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + by = (dc / (1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + (dc / (1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) - D(i,j) = -max(bx+by,-bmax) + D(i,j) = -max(bx+by, -bmax) if (D(i,j) > max_depth) D(i,j) = max_depth if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth enddo ; enddo From f272edc30bda516ddd86b2fff3736bbf2837636a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 15:37:30 -0500 Subject: [PATCH 03/10] +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 From e543439096cb7c8bdd88e3c55710bd5363dbfca8 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:50:57 -0500 Subject: [PATCH 04/10] +Rescale depth in Phillips_initialize_topography Added an optional unit_scale_type argument to Phillips_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/Phillips_initialization.F90 | 54 +++++++++++++++------------- 1 file changed, 29 insertions(+), 25 deletions(-) diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index e5eea1d2d2..1e513460b5 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -283,48 +283,52 @@ function sech(x) end function sech !> Initialize topography. -subroutine Phillips_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine Phillips_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + 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 - real :: PI, Htop, Wtop, Ltop, offset, dist, & - x1, x2, x3, x4, y1, y2 + ! Local variables + real :: m_to_Z ! A dimensional rescaling factor. + real :: PI, Htop, Wtop, Ltop, offset, dist + real :: x1, x2, x3, x4, y1, y2 integer :: i,j,is,ie,js,je character(len=40) :: mdl = "Phillips_initialize_topography" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec PI = 4.0*atan(1.0) + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl, "PHILLIPS_HTOP", Htop, & - "The maximum height of the topography.", units="m", & + call get_param(param_file, mdl, "PHILLIPS_HTOP", Htop, & + "The maximum height of the topography.", units="m", scale=m_to_Z, & fail_if_missing=.true.) ! Htop=0.375*max_depth ! max height of topog. above max_depth - Wtop=0.5*G%len_lat ! meridional width of drake and mount - Ltop=0.25*G%len_lon ! zonal width of topographic features - offset=0.1*G%len_lat ! meridional offset from center - dist=0.333*G%len_lon ! distance between drake and mount + Wtop = 0.5*G%len_lat ! meridional width of drake and mount + Ltop = 0.25*G%len_lon ! zonal width of topographic features + offset = 0.1*G%len_lat ! meridional offset from center + dist = 0.333*G%len_lon ! distance between drake and mount ! should be longer than Ltop/2 y1=G%south_lat+0.5*G%len_lat+offset-0.5*Wtop; y2=y1+Wtop x1=G%west_lon+0.1*G%len_lon; x2=x1+Ltop; x3=x1+dist; x4=x3+3.0/2.0*Ltop do i=is,ie ; do j=js,je - D(i,j)=0.0 - if (G%geoLonT(i,j)>x1 .and. G%geoLonT(i,j)y1 .and. G%geoLatT(i,j)x3 .and. G%geoLonT(i,j)y1 .and. G%geoLatT(i,j)x1 .and. G%geoLonT(i,j)y1 .and. G%geoLatT(i,j)x3 .and. G%geoLonT(i,j)y1 .and. G%geoLatT(i,j) Date: Fri, 16 Nov 2018 16:51:46 -0500 Subject: [PATCH 05/10] +Rescale depth in benchmark_initialize_topography Added an optional unit_scale_type argument to benchmark_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/benchmark_initialization.F90 | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index d50fe3fa05..0a5589a0b6 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -26,17 +26,20 @@ module benchmark_initialization contains !> This subroutine sets up the benchmark test case topography. -subroutine benchmark_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine benchmark_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + 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 - real :: min_depth ! The minimum and maximum depths in m. + real :: min_depth ! The minimum and maximum depths in Z. real :: PI ! 3.1415926... calculated as 4*atan(1) real :: D0 ! A constant to make the maximum ! ! basin depth MAXIMUM_DEPTH. ! + real :: m_to_Z ! A dimensional rescaling factor. real :: x, y ! This include declares and sets the variable "version". # include "version_variable.h" @@ -47,17 +50,19 @@ subroutine benchmark_initialize_topography(D, G, param_file, max_depth) call MOM_mesg(" benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5) + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call log_version(param_file, mdl, version, "") 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) PI = 4.0*atan(1.0) D0 = max_depth / 0.5 ! Calculate the depth of the bottom. - do i=is,ie ; do j=js,je - x=(G%geoLonT(i,j)-G%west_lon)/G%len_lon - y=(G%geoLatT(i,j)-G%south_lat)/G%len_lat + do j=js,je ; do i=is,ie + x = (G%geoLonT(i,j)-G%west_lon) / G%len_lon + y = (G%geoLatT(i,j)-G%south_lat) / G%len_lat ! This sets topography that has a reentrant channel to the south. D(i,j) = -D0 * ( y*(1.0 + 0.6*cos(4.0*PI*x)) & + 0.75*exp(-6.0*y) & From 8560c62f39d323dd1007788d9fcf68e38e2bb8b0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:52:29 -0500 Subject: [PATCH 06/10] +Rescale depth in shelfwave_initialize_topography Added an optional unit_scale_type argument to shelfwave_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/shelfwave_initialization.F90 | 33 +++++++++++++++------------ 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 9207830032..1a52519122 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -12,6 +12,7 @@ module shelfwave_initialization use MOM_open_boundary, only : OBC_segment_type, register_OBC use MOM_open_boundary, only : OBC_registry_type use MOM_time_manager, only : time_type, time_type_to_real +use MOM_unit_scaling, only : unit_scale_type implicit none ; private @@ -93,31 +94,33 @@ subroutine shelfwave_OBC_end(CS) end subroutine shelfwave_OBC_end !> Initialization of topography. -subroutine shelfwave_initialize_topography ( D, G, param_file, max_depth ) - ! Arguments - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine shelfwave_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + 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 + real :: m_to_Z ! A dimensional rescaling factor. integer :: i, j real :: y, rLy, Ly, H0 + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE",Ly, & default=50., do_not_log=.true.) - call get_param(param_file, mdl,"MINIMUM_DEPTH",H0, & - default=10., do_not_log=.true.) + call get_param(param_file, mdl,"MINIMUM_DEPTH", H0, & + default=10., units="m", scale=m_to_Z, do_not_log=.true.) rLy = 0. ; if (Ly>0.) rLy = 1. / Ly - do i=G%isc,G%iec - do j=G%jsc,G%jec - ! Compute normalized zonal coordinates (x,y=0 at center of domain) - y = ( G%geoLatT(i,j) - G%south_lat ) - D(i,j) = H0 * exp(2 * rLy * y) - enddo - enddo + + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! Compute normalized zonal coordinates (x,y=0 at center of domain) + y = ( G%geoLatT(i,j) - G%south_lat ) + D(i,j) = H0 * exp(2 * rLy * y) + enddo ; enddo end subroutine shelfwave_initialize_topography From c9517fedbb27e085c6742c253625247802260a2f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:55:48 -0500 Subject: [PATCH 07/10] +Rescale depth in USER_initialize_topography Added an optional unit_scale_type argument to USER_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/user_initialization.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 0a4ce7ccaa..c94613117a 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -13,6 +13,7 @@ module user_initialization use MOM_open_boundary, only : OBC_DIRECTION_S use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS use MOM_tracer_registry, only : tracer_registry_type +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type @@ -54,12 +55,13 @@ subroutine USER_set_coord(Rlay, g_prime, GV, param_file, eqn_of_state) end subroutine USER_set_coord !> Initialize topography. -subroutine USER_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine USER_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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + 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 call MOM_error(FATAL, & "USER_initialization.F90, USER_initialize_topography: " // & From 74014e7b77fe77decf43ca56a05a0202c89ca5f6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:56:37 -0500 Subject: [PATCH 08/10] Better comments in initialize_topography routines Clarified comments in several initialize_topography routines to explain that they will set the topography using the same units as the input value of mad_depth. All answers are bitwise identical. --- src/user/DOME2d_initialization.F90 | 16 ++++++------- src/user/Neverland_initialization.F90 | 23 +++++++++--------- src/user/dense_water_initialization.F90 | 16 ++++++++----- src/user/dumbbell_initialization.F90 | 32 ++++++++++++------------- src/user/seamount_initialization.F90 | 26 ++++++++++---------- src/user/sloshing_initialization.F90 | 10 ++++---- 6 files changed, 61 insertions(+), 62 deletions(-) diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index bd4f652dec..f20e4466bd 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -33,19 +33,19 @@ module DOME2d_initialization contains !> Initialize topography with a shelf and slope in a 2D domain -subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth ) - ! Arguments - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine DOME2d_initialize_topography( D, G, param_file, max_depth ) + 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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + intent(out) :: D !< Ocean bottom depth in the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + ! Local variables integer :: i, j real :: x, bay_depth, l1, l2 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "DOME2D_SHELF_WIDTH", dome2d_width_bay, & diff --git a/src/user/Neverland_initialization.F90 b/src/user/Neverland_initialization.F90 index 5f5b081f1c..67845a522f 100644 --- a/src/user/Neverland_initialization.F90 +++ b/src/user/Neverland_initialization.F90 @@ -26,18 +26,19 @@ module Neverland_initialization !> This subroutine sets up the Neverland test case topography. subroutine Neverland_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + 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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + intent(out) :: D !< Ocean bottom depth in the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + ! Local variables real :: PI ! 3.1415926... calculated as 4*atan(1) real :: D0 ! A constant to make the maximum ! ! basin depth MAXIMUM_DEPTH. ! real :: x, y -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "Neverland_initialize_topography" ! This subroutine's name. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed real :: nl_roughness_amp, nl_top_amp @@ -55,10 +56,9 @@ subroutine Neverland_initialize_topography(D, G, param_file, max_depth) PI = 4.0*atan(1.0) ! Calculate the depth of the bottom. - do i=is,ie - do j=js,je - x=(G%geoLonT(i,j)-G%west_lon)/G%len_lon - y=(G%geoLatT(i,j)-G%south_lat)/G%len_lat + do j=js,je ; do i=is,ie + x = (G%geoLonT(i,j)-G%west_lon) / G%len_lon + y =( G%geoLatT(i,j)-G%south_lat) / G%len_lat ! This sets topography that has a reentrant channel to the south. D(i,j) = 1.0 - 1.1 * spike(y-1,0.12) - 1.1 * spike(y,0.12) - & !< The great northern wall and Antarctica nl_top_amp*( & @@ -73,8 +73,7 @@ subroutine Neverland_initialize_topography(D, G, param_file, max_depth) - nl_roughness_amp * cos(20*PI*x) * cos(20*PI*y) !< roughness if (D(i,j) < 0.0) D(i,j) = 0.0 D(i,j) = D(i,j) * max_depth - enddo - enddo + enddo ; enddo end subroutine Neverland_initialize_topography ! ----------------------------------------------------------------------------- diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index f1a7bd6492..f857978c8e 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -11,6 +11,7 @@ module dense_water_initialization use MOM_file_parser, only : get_param, param_file_type use MOM_grid, only : ocean_grid_type use MOM_sponge, only : sponge_CS +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -32,10 +33,12 @@ module dense_water_initialization !> Initialize the topography field for the dense water experiment subroutine dense_water_initialize_topography(D, G, param_file, max_depth) - type(dyn_horgrid_type), intent(in) :: G !< Grid control structure - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Output topography field - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of the model + 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 the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + ! Local variables real, dimension(5) :: domain_params ! nondimensional widths of all domain sections real :: sill_frac, shelf_frac @@ -63,8 +66,8 @@ subroutine dense_water_initialize_topography(D, G, param_file, max_depth) domain_params(i) = domain_params(i-1) + domain_params(i) enddo - do i = G%isc,G%iec - do j = G%jsc,G%jec + do j = G%jsc,G%jec + do i = G%isc,G%iec ! compute normalised zonal coordinate x = (G%geoLonT(i,j) - G%west_lon) / G%len_lon @@ -88,6 +91,7 @@ subroutine dense_water_initialize_topography(D, G, param_file, max_depth) endif enddo enddo + end subroutine dense_water_initialize_topography !> Initialize the temperature and salinity for the dense water experiment diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index 5a5a845449..02222c9865 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -34,13 +34,13 @@ module dumbbell_initialization contains !> Initialization of topography. -subroutine dumbbell_initialize_topography ( D, G, param_file, max_depth ) - ! Arguments - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine dumbbell_initialize_topography( D, G, param_file, max_depth ) + 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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + intent(out) :: D !< Ocean bottom depth in the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + ! Local variables integer :: i, j real :: x, y, delta, dblen, dbfrac @@ -56,17 +56,15 @@ subroutine dumbbell_initialize_topography ( D, G, param_file, max_depth ) dblen=dblen*1.e3 endif - do i=G%isc,G%iec - do j=G%jsc,G%jec - ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) ) / dblen - y = ( G%geoLatT(i,j) ) / G%len_lat - D(i,j) = G%max_depth - if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then - D(i,j) = 0.0 - endif - enddo - enddo + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! Compute normalized zonal coordinates (x,y=0 at center of domain) + x = ( G%geoLonT(i,j) ) / dblen + y = ( G%geoLatT(i,j) ) / G%len_lat + D(i,j) = G%max_depth + if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then + D(i,j) = 0.0 + endif + enddo ; enddo end subroutine dumbbell_initialize_topography diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 3a5d7d186f..8f1ed97b06 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -33,13 +33,12 @@ module seamount_initialization contains !> Initialization of topography. -subroutine seamount_initialize_topography ( D, G, param_file, max_depth ) - ! Arguments - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine seamount_initialize_topography( D, G, param_file, max_depth ) + 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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + intent(out) :: D !< Ocean bottom depth in the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units ! Local variables integer :: i, j @@ -61,14 +60,13 @@ subroutine seamount_initialize_topography ( D, G, param_file, max_depth ) Ly = Ly / G%len_lat rLx = 0. ; if (Lx>0.) rLx = 1. / Lx rLy = 0. ; if (Ly>0.) rLy = 1. / Ly - do i=G%isc,G%iec - do j=G%jsc,G%jec - ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 - y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 - D(i,j) = G%max_depth * ( 1.0 - delta * exp(-(rLx*x)**2 -(rLy*y)**2) ) - enddo - enddo + + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! Compute normalized zonal coordinates (x,y=0 at center of domain) + x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 + y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 + D(i,j) = G%max_depth * ( 1.0 - delta * exp(-(rLx*x)**2 -(rLy*y)**2) ) + enddo ; enddo end subroutine seamount_initialize_topography diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index 5c4c5ec5b6..f192af804d 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -28,12 +28,12 @@ module sloshing_initialization contains !> Initialization of topography. -subroutine sloshing_initialize_topography ( D, G, param_file, max_depth ) - type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type +subroutine sloshing_initialize_topography( D, G, param_file, max_depth ) + 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 - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in m + intent(out) :: D !< Ocean bottom depth in the units of depth_max + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units ! Local variables integer :: i, j From bc808be8e3cd8a50a4fbef4f9aabbb38366b71f9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:57:55 -0500 Subject: [PATCH 09/10] +Rescale topography in MOM_shared_initialization Added an optional unit_scale_type argument to initialize_topography_from_file, apply_topography_edits_from_file, initialize_topography_named and limit_topography, which if present causes the topography to be scaled to units of Z when it is first set up. In addition, optional unit_scale_type arguments were added to all of the routines in MOM_shared_initialization and MOM_grid_initialize where they will eventually be needed so that the various components using this code (i.e., MOM6 and SIS2) can have a graceful transition to the new interfaces. All answers are bitwise identical. --- src/initialization/MOM_grid_initialize.F90 | 13 +- .../MOM_shared_initialization.F90 | 173 +++++++++--------- 2 files changed, 93 insertions(+), 93 deletions(-) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 7019564586..a0a354858f 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -54,9 +54,10 @@ module MOM_grid_initialize !> set_grid_metrics is used to set the primary values in the model's horizontal !! grid. The bathymetry, land-sea mask and any restricted channel widths are !! not known yet, so these are set later. -subroutine set_grid_metrics(G, param_file) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: param_file !< Parameter file structure +subroutine set_grid_metrics(G, param_file, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1216,9 +1217,9 @@ end function Adcroft_reciprocal !! any land or boundary point. For points in the interior, mask2dCu, !! mask2dCv, and mask2dBu are all 1.0. subroutine initialize_masks(G, PF, US) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: PF !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: PF !< Parameter file structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables real :: m_to_Z_scale ! A unit conversion factor from m to Z. real :: Dmin ! The depth for masking in the same units as G%bathyT (Z). diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 9ca352c3ed..113c3b3a85 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -49,10 +49,11 @@ end subroutine MOM_shared_init_init ! ----------------------------------------------------------------------------- !> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter. -subroutine MOM_initialize_rotation(f, G, PF) +subroutine MOM_initialize_rotation(f, G, PF, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f !< The Coriolis parameter in s-1 type(param_file_type), intent(in) :: PF !< Parameter file structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine makes the appropriate call to set up the Coriolis parameter. ! This is a separate subroutine so that it can be made public and shared with @@ -81,12 +82,13 @@ subroutine MOM_initialize_rotation(f, G, PF) end subroutine MOM_initialize_rotation !> Calculates the components of grad f (Coriolis parameter) -subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G) +subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(out) :: dF_dx !< x-component of grad f real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(out) :: dF_dy !< y-component of grad f + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables integer :: i,j real :: f1, f2 @@ -109,37 +111,39 @@ subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G) call pass_vector(dF_dx, dF_dy, G%Domain, stagger=AGRID) end subroutine MOM_calculate_grad_Coriolis -!> Return the global maximum ocean bottom depth in m. -function diagnoseMaximumDepth(D,G) +!> Return the global maximum ocean bottom depth in the same units as the input depth. +function diagnoseMaximumDepth(D, G) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(in) :: D !< Ocean bottom depth in m - real :: diagnoseMaximumDepth !< The global maximum ocean bottom depth in m + intent(in) :: D !< Ocean bottom depth in m or Z + real :: diagnoseMaximumDepth !< The global maximum ocean bottom depth in m or Z ! Local variables integer :: i,j - diagnoseMaximumDepth=D(G%isc,G%jsc) - do j=G%jsc, G%jec - do i=G%isc, G%iec - diagnoseMaximumDepth=max(diagnoseMaximumDepth,D(i,j)) - enddo - enddo + diagnoseMaximumDepth = D(G%isc,G%jsc) + do j=G%jsc, G%jec ; do i=G%isc, G%iec + diagnoseMaximumDepth = max(diagnoseMaximumDepth,D(i,j)) + enddo ; enddo call max_across_PEs(diagnoseMaximumDepth) end function diagnoseMaximumDepth !> Read gridded depths from file -subroutine initialize_topography_from_file(D, G, param_file) +subroutine initialize_topography_from_file(D, G, param_file, 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 + 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 + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables + real :: m_to_Z ! A dimensional rescaling factor. character(len=200) :: filename, topo_file, inputdir ! Strings for file/path character(len=200) :: topo_varname ! Variable name in file character(len=40) :: mdl = "initialize_topography_from_file" ! This subroutine's name. call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mdl, "TOPO_FILE", topo_file, & @@ -155,28 +159,29 @@ subroutine initialize_topography_from_file(D, G, param_file) if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_topography_from_file: Unable to open "//trim(filename)) - D(:,:) = -9.E30 ! Initializing to a very large negative depth (tall mountains) - ! everywhere before reading from a file should do nothing. - ! However, in the instance of masked-out PEs, halo regions - ! are not updated when a processor does not exist. We need to - ! ensure the depth in masked-out PEs appears to be that of land - ! so this line does that in the halo regions. For non-masked PEs - ! the halo region is filled properly with a later pass_var(). - call MOM_read_data(filename, trim(topo_varname), D, G%Domain) + D(:,:) = -9.e30*m_to_Z ! Initializing to a very large negative depth (tall mountains) everywhere + ! before reading from a file should do nothing. However, in the instance of + ! masked-out PEs, halo regions are not updated when a processor does not + ! exist. We need to ensure the depth in masked-out PEs appears to be that + ! of land so this line does that in the halo regions. For non-masked PEs + ! the halo region is filled properly with a later pass_var(). + call MOM_read_data(filename, trim(topo_varname), D, G%Domain, scale=m_to_Z) - call apply_topography_edits_from_file(D, G, param_file) + call apply_topography_edits_from_file(D, G, param_file, US) call callTree_leave(trim(mdl)//'()') end subroutine initialize_topography_from_file !> Applies a list of topography overrides read from a netcdf file -subroutine apply_topography_edits_from_file(D, G, param_file) +subroutine apply_topography_edits_from_file(D, G, param_file, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(inout) :: D !< Ocean bottom depth in m + intent(inout) :: D !< Ocean bottom depth in m or Z if US is present type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables + real :: m_to_Z ! A dimensional rescaling factor. character(len=200) :: topo_edits_file, inputdir ! Strings for file/path character(len=40) :: mdl = "apply_topography_edits_from_file" ! This subroutine's name. integer :: n_edits, n, ashape(5), i, j, ncid, id, ncstatus, iid, jid, zid @@ -185,6 +190,8 @@ subroutine apply_topography_edits_from_file(D, G, param_file) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, & @@ -267,8 +274,8 @@ subroutine apply_topography_edits_from_file(D, G, param_file) if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then if (new_depth(n)/=0.) then write(*,'(a,3i5,f8.2,a,f8.2,2i4)') & - 'Ocean topography edit: ',n,ig(n),jg(n),D(i,j),'->',abs(new_depth(n)),i,j - D(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) + 'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)/m_to_Z,'->',abs(new_depth(n)),i,j + D(i,j) = abs(m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) else call MOM_error(FATAL, ' apply_topography_edits_from_file: '//& "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file)) @@ -282,30 +289,25 @@ subroutine apply_topography_edits_from_file(D, G, param_file) end subroutine apply_topography_edits_from_file !> initialize the bathymetry based on one of several named idealized configurations -subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth) +subroutine initialize_topography_named(D, G, param_file, topog_config, 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 + 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 character(len=*), intent(in) :: topog_config !< The name of an idealized !! topographic configuration - real, intent(in) :: max_depth !< Maximum depth of model in m + real, intent(in) :: max_depth !< Maximum depth of model in the units of D + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type -! Arguments: D - the bottom depth in m. Intent out. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) topog_config - The name of an idealized topographic configuration. -! (in) max_depth - The maximum depth in m. + ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config. -! This subroutine places the bottom depth in m into D(:,:), shaped in a spoon - real :: min_depth ! The minimum depth in m. + ! Local variables + real :: m_to_Z ! A dimensional rescaling factor. + real :: min_depth ! The minimum depth in Z. real :: PI ! 3.1415926... calculated as 4*atan(1) - real :: D0 ! A constant to make the maximum ! - ! basin depth MAXIMUM_DEPTH. ! - real :: expdecay ! A decay scale of associated with ! - ! the sloping boundaries, in m. ! - real :: Dedge ! The depth in m at the basin edge. ! + real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH. + real :: expdecay ! A decay scale of associated with the sloping boundaries, in m. + real :: Dedge ! The depth in Z at the basin edge ! real :: south_lat, west_lon, len_lon, len_lat, Rad_earth integer :: i, j, is, ie, js, je, isd, ied, jsd, jed character(len=40) :: mdl = "initialize_topography_named" ! This subroutine's name. @@ -316,15 +318,17 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth call MOM_mesg(" MOM_shared_initialization.F90, initialize_topography_named: "//& "TOPO_CONFIG = "//trim(topog_config), 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) if (max_depth<=0.) call MOM_error(FATAL,"initialize_topography_named: "// & "MAXIMUM_DEPTH has a non-sensical value! Was it set?") if (trim(topog_config) /= "flat") then call get_param(param_file, mdl, "EDGE_DEPTH", Dedge, & "The depth at the edge of one of the named topographies.", & - units="m", default=100.0) + units="m", default=100.0, scale=m_to_Z) ! call get_param(param_file, mdl, "SOUTHLAT", south_lat, & ! "The southern latitude of the domain.", units="degrees", & ! fail_if_missing=.true.) @@ -398,37 +402,36 @@ end subroutine initialize_topography_named ! ----------------------------------------------------------------------------- !> limit_topography ensures that min_depth < D(x,y) < max_depth -subroutine limit_topography(D, G, param_file, max_depth) +subroutine limit_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(inout) :: D !< Ocean bottom depth in m + intent(inout) :: 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 depth of model in m -! Arguments: D - the bottom depth in m. Intent in/out. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) max_depth - The maximum depth in m. - -! This subroutine ensures that min_depth < D(x,y) < max_depth + real, intent(in) :: max_depth !< Maximum depth of model in the units of D + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + + ! Local variables + real :: m_to_Z ! A dimensional rescaling factor. integer :: i, j character(len=40) :: mdl = "limit_topography" ! This subroutine's name. real :: min_depth, mask_depth call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") + 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, & "If MASKING_DEPTH is unspecified, then anything shallower than\n"//& "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out.\n"//& "If MASKING_DEPTH is specified, then all depths shallower than\n"//& "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & - units="m", default=0.0) + units="m", default=0.0, scale=m_to_Z) call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & - "The depth below which to mask the ocean as land.", units="m", & - default=-9999.0, do_not_log=.true.) + "The depth below which to mask the ocean as land.", & + units="m", default=-9999.0, scale=m_to_Z, do_not_log=.true.) ! Make sure that min_depth < D(x,y) < max_depth - if (mask_depth<-9990.) then + if (mask_depth < -9990.*m_to_Z) then do j=G%jsd,G%jed ; do i=G%isd,G%ied D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth ) enddo ; enddo @@ -448,11 +451,12 @@ end subroutine limit_topography ! ----------------------------------------------------------------------------- !> This subroutine sets up the Coriolis parameter for a sphere -subroutine set_rotation_planetary(f, G, param_file) +subroutine set_rotation_planetary(f, G, param_file, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(out) :: f !< Coriolis parameter (vertical component) in s^-1 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine sets up the Coriolis parameter for a sphere character(len=30) :: mdl = "set_rotation_planetary" ! This subroutine's name. @@ -476,11 +480,12 @@ end subroutine set_rotation_planetary ! ----------------------------------------------------------------------------- !> This subroutine sets up the Coriolis parameter for a beta-plane or f-plane -subroutine set_rotation_beta_plane(f, G, param_file) +subroutine set_rotation_beta_plane(f, G, param_file, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(out) :: f !< Coriolis parameter (vertical component) in s^-1 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine sets up the Coriolis parameter for a beta-plane integer :: I, J @@ -594,18 +599,14 @@ end function modulo_around_point ! ----------------------------------------------------------------------------- !> This subroutine sets the open face lengths at selected points to restrict !! passages to their observed widths based on a named set of sizes. -subroutine reset_face_lengths_named(G, param_file, name) +subroutine reset_face_lengths_named(G, param_file, name, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters character(len=*), intent(in) :: name !< The name for the set of face lengths. Only "global_1deg" !! is currently implemented. -! This subroutine sets the open face lengths at selected points to restrict -! passages to their observed widths. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type -! Arguments: G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) name - The name for the set of face lengths. + ! Local variables character(len=256) :: mesg ! Message for error messages. real :: dx_2 = -1.0, dy_2 = -1.0 real :: pi_180 @@ -722,15 +723,12 @@ end subroutine reset_face_lengths_named ! ----------------------------------------------------------------------------- !> This subroutine sets the open face lengths at selected points to restrict !! passages to their observed widths from a arrays read from a file. -subroutine reset_face_lengths_file(G, param_file) +subroutine reset_face_lengths_file(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters -! This subroutine sets the open face lengths at selected points to restrict -! passages to their observed widths. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type -! Arguments: G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. + ! Local variables character(len=40) :: mdl = "reset_face_lengths_file" ! This subroutine's name. character(len=256) :: mesg ! Message for error messages. character(len=200) :: filename, chan_file, inputdir ! Strings for file/path @@ -791,15 +789,12 @@ end subroutine reset_face_lengths_file ! ----------------------------------------------------------------------------- !> This subroutine sets the open face lengths at selected points to restrict !! passages to their observed widths from a list read from a file. -subroutine reset_face_lengths_list(G, param_file) +subroutine reset_face_lengths_list(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters -! This subroutine sets the open face lengths at selected points to restrict -! passages to their observed widths. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type -! Arguments: G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. + ! Local variables character(len=120), pointer, dimension(:) :: lines => NULL() character(len=120) :: line character(len=200) :: filename, chan_file, inputdir, mesg ! Strings for file/path @@ -1121,8 +1116,8 @@ end subroutine set_velocity_depth_min !! later use in reporting diagnostics subroutine compute_global_grid_integrals(G) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid - ! Subroutine to pre-compute global integrals of grid quantities for - ! later use in reporting diagnostics + + ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming integer :: i,j @@ -1282,10 +1277,14 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) call write_field(unit, fields(19), G%Domain%mpp_domain, G%mask2dT) if (G%bathymetry_at_vel) then - call write_field(unit, fields(20), G%Domain%mpp_domain, G%Dblock_u) - call write_field(unit, fields(21), G%Domain%mpp_domain, G%Dopen_u) - call write_field(unit, fields(22), G%Domain%mpp_domain, G%Dblock_v) - call write_field(unit, fields(23), G%Domain%mpp_domain, G%Dopen_v) + do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dblock_u(I,j) ; enddo ; enddo + call write_field(unit, fields(20), G%Domain%mpp_domain, out_u) + do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dopen_u(I,j) ; enddo ; enddo + call write_field(unit, fields(21), G%Domain%mpp_domain, out_u) + do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dblock_v(i,J) ; enddo ; enddo + call write_field(unit, fields(22), G%Domain%mpp_domain, out_v) + do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dopen_v(i,J) ; enddo ; enddo + call write_field(unit, fields(23), G%Domain%mpp_domain, out_v) endif call close_file(unit) From 604a716ca95ee63975fd0f221399ad0159b00856 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 16 Nov 2018 16:58:39 -0500 Subject: [PATCH 10/10] Rescale topography during initialization Use new use_scale_type arguments in the call to MOM_initialize_topography and eliminated the call to rescale_dyn_horgrid_bathymetry in MOM_initialize_fixed, so the rescaling of the topography to units of Z happens when it is being set up via the various initialize_topography routines. All answers in the GFDL MOM6-examples test cases are bitwise identical, including rescaling Z over a large range. --- .../MOM_fixed_initialization.F90 | 45 ++++++++++--------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index c7ebb6a251..63774b98f8 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -83,8 +83,8 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) ! Set up the bottom depth, G%bathyT either analytically or from file ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH, ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) ) - call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF) - call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) + call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF, US) +! call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) ! To initialize masks, the bathymetry in halo regions must be filled in call pass_var(G%bathyT, G%Domain) @@ -170,19 +170,24 @@ end subroutine MOM_initialize_fixed !> MOM_initialize_topography makes the appropriate call to set up the bathymetry. At this !! point the topography is in units of m, but this can be changed later. -subroutine MOM_initialize_topography(D, max_depth, G, PF) +subroutine MOM_initialize_topography(D, max_depth, G, PF, 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 type(param_file_type), intent(in) :: PF !< Parameter file structure real, intent(out) :: max_depth !< Maximum depth of model in m + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine makes the appropriate call to set up the bottom depth. ! This is a separate subroutine so that it can be made public and shared with ! the ice-sheet code or other components. + real :: m_to_Z, Z_to_m ! Dimensional rescaling factors character(len=40) :: mdl = "MOM_initialize_topography" ! This subroutine's name. character(len=200) :: config + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + Z_to_m = 1.0 ; if (present(US)) Z_to_m = US%Z_to_m + call get_param(PF, mdl, "TOPO_CONFIG", config, & "This specifies how bathymetry is specified: \n"//& " \t file - read bathymetric information from the file \n"//& @@ -210,39 +215,39 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) " \t dense - Denmark Strait-like dense water formation and overflow.\n"//& " \t USER - call a user modified routine.", & fail_if_missing=.true.) - max_depth = -1.e9; call read_param(PF, "MAXIMUM_DEPTH", max_depth) + max_depth = -1.e9*m_to_Z ; call read_param(PF, "MAXIMUM_DEPTH", max_depth, scale=m_to_Z) select case ( trim(config) ) - case ("file"); call initialize_topography_from_file(D, G, PF) - case ("flat"); call initialize_topography_named(D, G, PF, config, max_depth) - case ("spoon"); call initialize_topography_named(D, G, PF, config, max_depth) - case ("bowl"); call initialize_topography_named(D, G, PF, config, max_depth) - case ("halfpipe"); call initialize_topography_named(D, G, PF, config, max_depth) - case ("DOME"); call DOME_initialize_topography(D, G, PF, max_depth) - case ("ISOMIP"); call ISOMIP_initialize_topography(D, G, PF, max_depth) - case ("benchmark"); call benchmark_initialize_topography(D, G, PF, max_depth) + case ("file"); call initialize_topography_from_file(D, G, PF, US) + case ("flat"); call initialize_topography_named(D, G, PF, config, max_depth, US) + case ("spoon"); call initialize_topography_named(D, G, PF, config, max_depth, US) + case ("bowl"); call initialize_topography_named(D, G, PF, config, max_depth, US) + case ("halfpipe"); call initialize_topography_named(D, G, PF, config, max_depth, US) + case ("DOME"); call DOME_initialize_topography(D, G, PF, max_depth, US) + case ("ISOMIP"); call ISOMIP_initialize_topography(D, G, PF, max_depth, US) + case ("benchmark"); call benchmark_initialize_topography(D, G, PF, max_depth, US) case ("Neverland"); call Neverland_initialize_topography(D, G, PF, max_depth) case ("DOME2D"); call DOME2d_initialize_topography(D, G, PF, max_depth) - case ("Kelvin"); call Kelvin_initialize_topography(D, G, PF, max_depth) + case ("Kelvin"); call Kelvin_initialize_topography(D, G, PF, max_depth, US) case ("sloshing"); call sloshing_initialize_topography(D, G, PF, max_depth) case ("seamount"); call seamount_initialize_topography(D, G, PF, max_depth) - case ("dumbbell"); call dumbbell_initialize_topography(D, G, PF, max_depth) - case ("shelfwave"); call shelfwave_initialize_topography(D, G, PF, max_depth) - case ("Phillips"); call Phillips_initialize_topography(D, G, PF, max_depth) + case ("dumbbell"); call dumbbell_initialize_topography(D, G, PF, max_depth) + case ("shelfwave"); call shelfwave_initialize_topography(D, G, PF, max_depth, US) + case ("Phillips"); call Phillips_initialize_topography(D, G, PF, max_depth, US) case ("dense"); call dense_water_initialize_topography(D, G, PF, max_depth) - case ("USER"); call user_initialize_topography(D, G, PF, max_depth) + case ("USER"); call user_initialize_topography(D, G, PF, max_depth, US) case default ; call MOM_error(FATAL,"MOM_initialize_topography: "// & "Unrecognized topography setup '"//trim(config)//"'") end select if (max_depth>0.) then - call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, & + call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth*Z_to_m, & "The maximum depth of the ocean.", units="m") else max_depth = diagnoseMaximumDepth(D,G) - call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth, & + call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth*Z_to_m, & "The (diagnosed) maximum depth of the ocean.", units="m") endif if (trim(config) /= "DOME") then - call limit_topography(D, G, PF, max_depth) + call limit_topography(D, G, PF, max_depth, US) endif end subroutine MOM_initialize_topography