diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 6fb8426395..5e5b307103 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -190,6 +190,10 @@ module MOM_grid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index f8ae58d9e1..4c1b87f37e 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -135,6 +135,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) ! Copy various scalar variables and strings. oG%x_axis_units = dG%x_axis_units ; oG%y_axis_units = dG%y_axis_units oG%x_ax_unit_short = dG%x_ax_unit_short ; oG%y_ax_unit_short = dG%y_ax_unit_short + oG%grid_unit_to_L = dG%grid_unit_to_L oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon @@ -296,6 +297,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) ! Copy various scalar variables and strings. dG%x_axis_units = oG%x_axis_units ; dG%y_axis_units = oG%y_axis_units dG%x_ax_unit_short = oG%x_ax_unit_short ; dG%y_ax_unit_short = oG%y_ax_unit_short + dG%grid_unit_to_L = oG%grid_unit_to_L dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 987d5bf502..440ea351b9 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -181,6 +181,10 @@ module MOM_dyn_horgrid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] @@ -400,6 +404,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%len_lon = G_in%len_lon ! Rotation-invariant fields + G%grid_unit_to_L = G_in%grid_unit_to_L G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global G%Rad_Earth = G_in%Rad_Earth diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 5ecfb9e788..fe54dd6533 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -132,6 +132,7 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) G%x_axis_units = "degrees_E" ; G%y_axis_units = "degrees_N" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" + G%grid_unit_to_L = 0.0 if (index(lowercase(trim(grid_config)),"cartesian") > 0) then ! This is a cartesian grid, and may have different axis units. @@ -145,9 +146,11 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" G%x_ax_unit_short = "km" ; G%y_ax_unit_short = "km" + G%grid_unit_to_L = 1000.0*diag_cs%US%m_to_L elseif (units_temp(1:1) == 'm') then G%x_axis_units = "meters" ; G%y_axis_units = "meters" G%x_ax_unit_short = "m" ; G%y_ax_unit_short = "m" + G%grid_unit_to_L = diag_cs%US%m_to_L endif call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units) else diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index ef78a896c3..429ce24d79 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -84,7 +84,7 @@ subroutine set_grid_metrics(G, param_file, US) ! These are defaults that may be changed in the next select block. G%x_axis_units = "degrees_east" ; G%y_axis_units = "degrees_north" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" - + G%grid_unit_to_L = 0.0 G%Rad_Earth_L = -1.0*US%m_to_L ; G%len_lat = 0.0 ; G%len_lon = 0.0 select case (trim(config)) case ("mosaic"); call set_grid_metrics_from_mosaic(G, param_file, US) @@ -251,6 +251,11 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) G%geoLatCv(i,J) = tmpZ(i2-1,j2) enddo ; enddo + ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates, + ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be + ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates + ! are assumed to be degrees of longitude and latitude. + ! Read DX,DY from the supergrid tmpU(:,:) = 0. ; tmpV(:,:) = 0. call MOM_read_data(filename, 'dx', tmpV, SGdom, position=NORTH_FACE, scale=US%m_to_L) @@ -440,9 +445,11 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) enddo if (units_temp(1:1) == 'k') then ! Axes are measured in km. + G%grid_unit_to_L = 1000.0*US%m_to_L dx_everywhere = 1000.0*US%m_to_L * G%len_lon / (REAL(niglobal)) dy_everywhere = 1000.0*US%m_to_L * G%len_lat / (REAL(njglobal)) elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. + G%grid_unit_to_L = US%m_to_L dx_everywhere = US%m_to_L*G%len_lon / (REAL(niglobal)) dy_everywhere = US%m_to_L*G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 0a257b6ca1..1f7a519d9e 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -485,7 +485,6 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1] real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1] real :: beta_lat_ref ! The reference latitude for the beta plane [degrees_N] or [km] or [m] - real :: Rad_Earth_L ! The radius of the planet in rescaled units [L ~> m] real :: y_scl ! A scaling factor from the units of latitude [L lat-1 ~> m lat-1] real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name. @@ -503,18 +502,16 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees") PI = 4.0*atan(1.0) + y_scl = G%grid_unit_to_L + if (G%grid_unit_to_L <= 0.0) y_scl = PI * G%Rad_Earth_L / 180. + select case (axis_units(1:1)) case ("d") - call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) beta_lat_ref_units = "degrees" - y_scl = PI * Rad_Earth_L / 180. case ("k") beta_lat_ref_units = "kilometers" - y_scl = 1.0e3 * US%m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1.0 * US%m_to_L case default ; call MOM_error(FATAL, & " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units)) end select diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7909bddb80..0f5d06aa4e 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1636,7 +1636,10 @@ subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read) if (just_read) return ! All run-time parameters have been read, so return. - dpi=acos(0.0)*2.0 ! pi + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "MOM_state_initialization.F90: "//& + "initialize_velocity_circular() is only set to work with Cartesian axis units.") + + dpi = acos(0.0)*2.0 ! pi do k=1,nz ; do j=js,je ; do I=Isq,Ieq psi1 = my_psi(I,j) @@ -1663,7 +1666,7 @@ real function my_psi(ig,jg) r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only r = min(1.0, r) ! Flatten stream function in corners of box my_psi = 0.5*(1.0 - cos(dpi*r)) - my_psi = my_psi * (circular_max_u * G%US%m_to_L*G%len_lon*1e3 / dpi) ! len_lon is in km + my_psi = my_psi * (circular_max_u * G%len_lon * G%grid_unit_to_L / dpi) ! len_lon is in km end function my_psi end subroutine initialize_velocity_circular diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 858ca32f93..e608dbd1c2 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -49,10 +49,11 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) real :: min_depth ! The minimum ocean depth [Z ~> m] real :: shelf_depth ! The ocean depth on the shelf in the DOME configuration [Z ~> m] real :: slope ! The bottom slope in the DOME configuration [Z L-1 ~> nondim] - real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf [km] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] - real :: inflow_width ! The longitudinal width of the DOME inflow channel [km] - real :: km_to_L ! The conversion factor from the units of latitude to L [L km-1 ~> 1e3] + real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf in the same units as geolat, often [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] + real :: inflow_width ! The longitudinal width of the DOME inflow channel in the same units as geolat, often [km] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "DOME_initialize_topography" ! This subroutine's name. @@ -60,7 +61,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif call MOM_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5) @@ -75,15 +85,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) default=600.0, units="m", scale=US%m_to_Z) call get_param(param_file, mdl, "DOME_SHELF_EDGE_LAT", shelf_edge_lat, & "The latitude of the shelf edge in the DOME configuration.", & - default=600.0, units="km") + default=600.0, units="km", scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_WIDTH", inflow_width, & - "The longitudinal width of the DOME inflow channel.", units="km", default=100.0) + "The longitudinal width of the DOME inflow channel.", & + units="km", default=100.0, scale=km_to_grid_unit) do j=js,je ; do i=is,ie if (G%geoLatT(i,j) < shelf_edge_lat) then - D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*km_to_L, max_depth) + D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*G%grid_unit_to_L, max_depth) else if ((G%geoLonT(i,j) > inflow_lon) .AND. (G%geoLonT(i,j) < inflow_lon+inflow_width)) then D(i,j) = shelf_depth @@ -177,7 +188,6 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) real :: min_depth ! The minimum depth at which to apply damping [Z ~> m] real :: damp_W, damp_E ! Damping rates in the western and eastern sponges [T-1 ~> s-1] real :: peak_damping ! The maximum sponge damping rates as the edges [T-1 ~> s-1] - real :: km_to_L ! The conversion factor from the units of longitude to L [L km-1 ~> 1e3] real :: edge_dist ! The distance to an edge [L ~> m] real :: sponge_width ! The width of the sponges [L ~> m] character(len=40) :: mdl = "DOME_initialize_sponges" ! This subroutine's name. @@ -186,7 +196,8 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_sponges is only set to work with Cartesian axis units.") ! Set up sponges for the DOME configuration call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & @@ -196,7 +207,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) default=10.0, units="day-1", scale=1.0/(86400.0*US%s_to_T)) call get_param(PF, mdl, "DOME_SPONGE_WIDTH", sponge_width, & "The width of the the DOME sponges.", & - default=200.0, units="km", scale=km_to_L) + default=200.0, units="km", scale=1.0e3*US%m_to_L) ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0 wherever ! there is no sponge, and the subroutines that are called will automatically @@ -204,7 +215,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) Idamp(:,:) = 0.0 do j=js,je ; do i=is,ie ; if (depth_tot(i,j) > min_depth) then - edge_dist = (G%geoLonT(i,j) - G%west_lon) * km_to_L + edge_dist = (G%geoLonT(i,j) - G%west_lon) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_W = peak_damping elseif (edge_dist < sponge_width) then @@ -213,7 +224,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) damp_W = 0.0 endif - edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * km_to_L + edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_E = peak_damping elseif (edge_dist < sponge_width) then @@ -328,10 +339,12 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) ! properties [T-1 ~> s-1] real :: g_prime_tot ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2] real :: Def_Rad ! The deformation radius, based on fluid of thickness D_edge [L ~> m] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] real :: I_Def_Rad ! The inverse of the deformation radius in the same units as G%geoLon [km-1] real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile [nondim] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] character(len=32) :: name ! The name of a tracer field. character(len=40) :: mdl = "DOME_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntherm, ntr_id @@ -343,6 +356,17 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif + call get_param(PF, mdl, "DOME_INFLOW_THICKNESS", D_edge, & "The thickness of the dense DOME inflow at the inner edge.", & default=300.0, units="m", scale=US%m_to_Z) @@ -362,7 +386,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) "The value of the Coriolis parameter that is used to determine the DOME "//& "inflow properties.", units="s-1", default=f_0*US%s_to_T, scale=US%T_to_s) call get_param(PF, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) if (associated(tv%S) .or. associated(tv%T)) then call get_param(PF, mdl, "S_REF", S_ref, & units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) @@ -383,7 +407,9 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * (Rlay_Ref + 0.5*Rlay_range) * GV%RZ_to_H endif - I_Def_Rad = 1.0 / (1.0e-3*US%L_to_m*Def_Rad) + I_Def_Rad = 1.0 / ((1.0e-3*US%L_to_m*km_to_grid_unit) * Def_Rad) + ! This is mathematically equivalent to + ! I_Def_Rad = G%grid_unit_to_L / Def_Rad if (OBC%number_of_segments /= 1) then call MOM_error(WARNING, 'Error in DOME OBC segment setup', .true.) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d03a07e313..4bf7931856 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -59,7 +59,6 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) real :: ly ! domain width (across ice flow) [L ~> m] real :: bx, by ! The x- and y- contributions to the bathymetric profiles at a point [Z ~> m] real :: xtil ! x-positon normalized by the characteristic along-flow length scale [nondim] - real :: km_to_L ! The conversion factor from the axis units to L [L km-1 ~> 1e3] logical :: is_2D ! If true, use a 2D setup ! This include declares and sets the variable "version". # include "version_variable.h" @@ -93,7 +92,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) "Characteristic width of the side walls of the channel in the ISOMIP configuration.", & units="m", default=4.0e3, scale=US%m_to_L) - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "ISOMIP_initialization.F90: " //& + "ISOMIP_initialize_topography is only set to work with Cartesian axis units.") ! The following variables should be transformed into runtime parameters. b0 = -150.0*US%m_to_Z ; b2 = -728.8*US%m_to_Z ; b4 = 343.91*US%m_to_Z ; b6 = -50.57*US%m_to_Z @@ -101,8 +101,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) if (is_2D) then do j=js,je ; do i=is,ie ! For the 2D setup take a slice through the middle of the domain - xtil = G%geoLonT(i,j)*km_to_L / xbar - !xtil = 450.*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar + ! xtil = 450.0e3*US%m_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = 2.0 * dc / (1.0 + exp(2.0*wc / fc)) @@ -117,17 +117,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) ! 3D setup ! ===== TEST ===== !if (G%geoLonT(i,j)<500.) then - ! xtil = 500.*km_to_L / xbar + ! xtil = 500.0e3*US%m_to_L / xbar !else - ! xtil = G%geoLonT(i,j)*km_to_L / xbar + ! xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar !endif ! ===== TEST ===== - xtil = G%geoLonT(i,j)*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly - wc) / fc))) + & - (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly + wc) / fc))) + by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly - wc) / fc))) + & + (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly + wc) / fc))) D(i,j) = -max(bx+by, -bmax) if (D(i,j) > max_depth) D(i,j) = max_depth diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index aca19e86d0..e092f42508 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -349,7 +349,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) real :: fbench !< The benchmark 'f' value [T-1 ~> s-1] real :: fbench_fac !< A factor that is set to 0 to use the !! benchmark 'f' value [nondim] - real :: km_to_L !< The conversion factor from the units of latitude to L [L km-1 ~> 1e3] real :: rel_tau_fac !< A factor that is set to 0 to disable !! current relative stress calculation [nondim] @@ -359,7 +358,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "Idealized_Hurricane.F90: " //& + "idealized_hurricane_wind_forcing is only set to work with Cartesian axis units.") ! Allocate the forcing arrays, if necessary. call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., tau_mag=.true.) @@ -376,7 +376,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YC = CS%Hurr_cen_Y0 + (time_type_to_real(day)*US%s_to_T * CS%hurr_translation_spd * & sin(CS%hurr_translation_dir)) - if (CS%BR_Bench) then ! f reset to value used in generated wind for benchmark test fbench = CS%f_column @@ -403,8 +402,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCu(I,j)*km_to_L - YC - XX = G%geoLonCu(I,j)*km_to_L - XC + YY = G%geoLatCu(I,j) * G%grid_unit_to_L - YC + XX = G%geoLonCu(I,j) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%taux(I,j) = G%mask2dCu(I,j) * TX @@ -427,8 +426,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCv(i,J)*km_to_L - YC - XX = G%geoLonCv(i,J)*km_to_L - XC + YY = G%geoLatCv(i,J) * G%grid_unit_to_L - YC + XX = G%geoLonCv(i,J) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%tauy(i,J) = G%mask2dCv(i,J) * TY diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 1fc8a2f564..118d76ac93 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -215,10 +215,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) if (.not.associated(OBC)) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & 'Kelvin_set_OBC_data() was called but OBC type was not initialized!') + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & + "Kelvin_set_OBC_data() is only set to work with Cartesian axis units.") time_sec = US%s_to_T*time_type_to_real(Time) PI = 4.0*atan(1.0) - km_to_L_scale = 1000.0*US%m_to_L do j=jsd,jed ; do i=isd,ied depth_tot(i,j) = 0.0 @@ -237,6 +238,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0) ! Two wavelengths in domain omega = (4.0 * CS%H0 * N0) / (CS%mode * US%m_to_L*G%len_lon) + !### There is a bug here when len_lon is in km. This should be + ! omega = (4.0 * CS%H0 * N0) / (CS%mode * G%grid_unit_to_L*G%len_lon) endif sina = sin(CS%coast_angle) @@ -257,8 +260,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) jsd = segment%HI%jsd ; jed = segment%HI%jed JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do j=jsd,jed ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonCu(I,j) - y1 = km_to_L_scale * G%geoLatCu(I,j) + x1 = G%grid_unit_to_L * G%geoLonCu(I,j) + y1 = G%grid_unit_to_L * G%geoLatCu(I,j) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = -(x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -299,8 +302,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB+1,JedB-1 ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i+1,j) ) @@ -316,8 +319,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) isd = segment%HI%isd ; ied = segment%HI%ied JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do J=JsdB,JedB ; do i=isd,ied - x1 = km_to_L_scale * G%geoLonCv(i,J) - y1 = km_to_L_scale * G%geoLatCv(i,J) + x1 = G%grid_unit_to_L * G%geoLonCv(i,J) + y1 = G%grid_unit_to_L * G%geoLatCv(i,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -355,8 +358,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB,JedB ; do I=IsdB+1,IedB-1 - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i,j+1) ) diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 33c7641a00..09219eb845 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -250,6 +250,8 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just if (max_depth <= 0.0) call MOM_error(FATAL, & "Rossby_front_initialize_thickness, Rossby_front_initialize_velocity: "//& "This module requires a positive value of MAXIMUM_DEPTH.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Rossby_front_2d_initialization.F90: '// & + "dTdy() is only set to work with Cartesian axis units.") v(:,:,:) = 0.0 u(:,:,:) = 0.0 @@ -335,18 +337,16 @@ end function Hml real function dTdy( G, dT, lat, US ) type(ocean_grid_type), intent(in) :: G !< Grid structure real, intent(in) :: dT !< Top to bottom temperature difference [C ~> degC] - real, intent(in) :: lat !< Latitude in [km] + real, intent(in) :: lat !< Latitude in the same units as geoLat, often [km] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] real :: dHML ! The range of the mixed layer depths [Z ~> m] real :: dHdy ! The mixed layer depth gradient [Z L-1 ~> m m-1] - real :: km_to_L ! Horizontal axis unit conversion factor when AXIS_UNITS = 'k' (1000 m) [L km-1 ~> 1000] PI = 4.0 * atan(1.0) - km_to_L = 1.0e3*US%m_to_L dHML = 0.5 * ( HMLmax - HMLmin ) * G%max_depth - dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * km_to_L ) ) * cos( yPseudo(G, lat) ) + dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * G%grid_unit_to_L ) ) * cos( yPseudo(G, lat) ) dTdy = -( dT / G%max_depth ) * dHdy end function dTdy diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index 722a41b7e5..a734574995 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -79,10 +79,12 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US, param_file, jus if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_thickness: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_thickness() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq I_nz = 1.0 / real(nz) x0 = 2.0*G%len_lon/3.0 @@ -150,10 +152,12 @@ subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read) if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_velocity: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_velocity() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq x0 = 2.0*G%len_lon/3.0 y0 = 0.0