diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index a1c300c94f..be00de8779 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -97,7 +97,7 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric) integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0). logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully !! symmetric computational domain. - real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> nondim] or [nondim] + real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> 1] or [1] integer :: hs logical :: sym diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 7592dc8477..9c9ebf4960 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -171,10 +171,11 @@ module MOM_grid ! initialization routines (but not all) real :: south_lat !< The latitude (or y-coordinate) of the first v-line real :: west_lon !< The longitude (or x-coordinate) of the first u-line - real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain - real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain - real :: Rad_Earth = 6.378e6 !< The radius of the planet [m]. - real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m]. + real :: len_lat !< The latitudinal (or y-coord) extent of physical domain + real :: len_lon !< The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth !< The radius of the planet [m] + real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] + real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m] end type ocean_grid_type contains diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index d3447f6590..58063f0669 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -126,7 +126,8 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) 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 - oG%Rad_Earth = dG%Rad_Earth ; oG%max_depth = dG%max_depth + oG%Rad_Earth = dG%Rad_Earth ; oG%Rad_Earth_L = dG%Rad_Earth_L + oG%max_depth = dG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(oG%areaT, oG%Domain) @@ -272,7 +273,8 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) 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 - dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth + dG%Rad_Earth = oG%Rad_Earth ; dG%Rad_Earth_L = oG%Rad_Earth_L + dG%max_depth = oG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(dG%areaT, dG%Domain) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index f6f7205dcf..3b6fb0c510 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -392,9 +392,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci Temp_int, Salt_int ! Layer and cell integrated heat and salt [J] and [g Salt]. real :: HL2_to_kg ! A conversion factor from a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1] real :: KE_scale_factor ! The combination of unit rescaling factors in the kinetic energy - ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or nondim] + ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or 1] real :: PE_scale_factor ! The combination of unit rescaling factors in the potential energy - ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> nondim] + ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> 1] integer :: num_nc_fields ! The number of fields that will actually go into ! the NetCDF file. integer :: i, j, k, is, ie, js, je, ns, nz, m, Isq, Ieq, Jsq, Jeq, isr, ier, jsr, jer diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index c7db67ee17..b67056cd80 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -163,10 +163,11 @@ module MOM_dyn_horgrid ! initialization routines (but not all) real :: south_lat !< The latitude (or y-coordinate) of the first v-line real :: west_lon !< The longitude (or x-coordinate) of the first u-line - real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain - real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain - real :: Rad_Earth = 6.378e6 !< The radius of the planet [m]. - real :: max_depth !< The maximum depth of the ocean [Z ~> m]. + real :: len_lat !< The latitudinal (or y-coord) extent of physical domain + real :: len_lon !< The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth !< The radius of the planet [m] + real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] + real :: max_depth !< The maximum depth of the ocean [Z ~> m] end type dyn_horgrid_type contains @@ -361,6 +362,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global G%Rad_Earth = G_in%Rad_Earth + G%Rad_Earth_L = G_in%Rad_Earth_L G%max_depth = G_in%max_depth call set_derived_dyn_horgrid(G, US) @@ -406,12 +408,8 @@ subroutine set_derived_dyn_horgrid(G, US) type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Various inverse grid spacings and derived areas are calculated within this ! subroutine. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [L m-1 ~> nondim] integer :: i, j, isd, ied, jsd, jed integer :: IsdB, IedB, JsdB, JedB - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 81e7b66d7a..f67a977d27 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -26,14 +26,14 @@ module MOM_grid_initialize ! vary with the Boussinesq approximation, the Boussinesq variant is given first. !> Global positioning system (aka container for information to describe the grid) -type, public :: GPS ; private +type, private :: GPS ; private real :: len_lon !< The longitudinal or x-direction length of the domain. real :: len_lat !< The latitudinal or y-direction length of the domain. real :: west_lon !< The western longitude of the domain or the equivalent !! starting value for the x-axis. real :: south_lat !< The southern latitude of the domain or the equivalent !! starting value for the y-axis. - real :: Rad_Earth !< The radius of the Earth [m]. + real :: Rad_Earth_L !< The radius of the Earth in rescaled units [L ~> m] real :: Lat_enhance_factor !< The amount by which the meridional resolution !! is enhanced within LAT_EQ_ENHANCE of the equator. real :: Lat_eq_enhance !< The latitude range to the north and south of the equator @@ -59,12 +59,18 @@ 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" + real :: m_to_L ! A length unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A length unit conversion factor [m L-1 ~> 1] + ! This include declares and sets the variable "version". +# include "version_variable.h" logical :: debug character(len=256) :: config + m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L + L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m + call callTree_enter("set_grid_metrics(), MOM_grid_initialize.F90") call log_version(param_file, "MOM_grid_init", version, "") call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, & @@ -82,6 +88,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%Rad_Earth_L = -1.0*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) case ("cartesian"); call set_grid_metrics_cartesian(G, param_file, US) @@ -93,6 +100,15 @@ subroutine set_grid_metrics(G, param_file, US) case default ; call MOM_error(FATAL, "MOM_grid_init: set_grid_metrics "//& "Unrecognized grid configuration "//trim(config)) end select + if (G%Rad_Earth_L <= 0.0) then + ! The grid metrics were set with an option that does not explicitly initialize Rad_Earth. + ! ### Rad_Earth should be read as in: + ! call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & + ! "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) + ! but for now it is being set via a hard-coded value to reproduce current behavior. + G%Rad_Earth_L = 6.378e6*m_to_L + endif + G%Rad_Earth = L_to_m*G%Rad_Earth_L ! Calculate derived metrics (i.e. reciprocals and products) call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") @@ -113,8 +129,8 @@ subroutine grid_metrics_chksum(parent, G, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] integer :: halo m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m @@ -181,7 +197,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ real, dimension(:,:), allocatable :: tmpGlbl - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=200) :: filename, grid_file, inputdir character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic" type(MOM_domain_type), pointer :: SGdom => NULL() ! Supergrid domain @@ -386,11 +402,10 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) integer :: niglobal, njglobal real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) - real :: dx_everywhere, dy_everywhere ! Grid spacings [m]. - real :: I_dx, I_dy ! Inverse grid spacings [m-1]. + real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m]. + real :: I_dx, I_dy ! Inverse grid spacings [L-1 ~> m-1]. real :: PI - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=80) :: units_temp character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_cartesian" @@ -402,7 +417,6 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) call callTree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.F90") m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m PI = 4.0*atan(1.0) call get_param(param_file, mdl, "AXIS_UNITS", units_temp, & @@ -423,8 +437,8 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) call get_param(param_file, mdl, "LENLON", G%len_lon, & "The longitudinal or x-direction length of the domain.", & units=units_temp, fail_if_missing=.true.) - call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & + "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" @@ -462,14 +476,14 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) enddo if (units_temp(1:1) == 'k') then ! Axes are measured in km. - dx_everywhere = 1000.0 * G%len_lon / (REAL(niglobal)) - dy_everywhere = 1000.0 * G%len_lat / (REAL(njglobal)) + dx_everywhere = 1000.0*m_to_L * G%len_lon / (REAL(niglobal)) + dy_everywhere = 1000.0*m_to_L * G%len_lat / (REAL(njglobal)) elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. - dx_everywhere = G%len_lon / (REAL(niglobal)) - dy_everywhere = G%len_lat / (REAL(njglobal)) + dx_everywhere = m_to_L*G%len_lon / (REAL(niglobal)) + dy_everywhere = m_to_L*G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. - dx_everywhere = G%Rad_Earth * G%len_lon * PI / (180.0 * niglobal) - dy_everywhere = G%Rad_Earth * G%len_lat * PI / (180.0 * njglobal) + dx_everywhere = G%Rad_Earth_L * G%len_lon * PI / (180.0 * niglobal) + dy_everywhere = G%Rad_Earth_L * G%len_lat * PI / (180.0 * njglobal) endif I_dx = 1.0 / dx_everywhere ; I_dy = 1.0 / dy_everywhere @@ -477,30 +491,30 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) do J=JsdB,JedB ; do I=IsdB,IedB G%geoLonBu(I,J) = grid_lonB(I) ; G%geoLatBu(I,J) = grid_latB(J) - G%dxBu(I,J) = m_to_L*dx_everywhere ; G%IdxBu(I,J) = L_to_m*I_dx - G%dyBu(I,J) = m_to_L*dy_everywhere ; G%IdyBu(I,J) = L_to_m*I_dy - G%areaBu(I,J) = m_to_L**2*dx_everywhere * dy_everywhere ; G%IareaBu(I,J) = L_to_m**2*I_dx * I_dy + G%dxBu(I,J) = dx_everywhere ; G%IdxBu(I,J) = I_dx + G%dyBu(I,J) = dy_everywhere ; G%IdyBu(I,J) = I_dy + G%areaBu(I,J) = dx_everywhere * dy_everywhere ; G%IareaBu(I,J) = I_dx * I_dy enddo ; enddo do j=jsd,jed ; do i=isd,ied G%geoLonT(i,j) = grid_lonT(i) ; G%geoLatT(i,j) = grid_LatT(j) - G%dxT(i,j) = m_to_L*dx_everywhere ; G%IdxT(i,j) = L_to_m*I_dx - G%dyT(i,j) = m_to_L*dy_everywhere ; G%IdyT(i,j) = L_to_m*I_dy - G%areaT(i,j) = m_to_L**2*dx_everywhere * dy_everywhere ; G%IareaT(i,j) = L_to_m**2*I_dx * I_dy + G%dxT(i,j) = dx_everywhere ; G%IdxT(i,j) = I_dx + G%dyT(i,j) = dy_everywhere ; G%IdyT(i,j) = I_dy + G%areaT(i,j) = dx_everywhere * dy_everywhere ; G%IareaT(i,j) = I_dx * I_dy enddo ; enddo do j=jsd,jed ; do I=IsdB,IedB G%geoLonCu(I,j) = grid_lonB(I) ; G%geoLatCu(I,j) = grid_LatT(j) - G%dxCu(I,j) = m_to_L*dx_everywhere ; G%IdxCu(I,j) = L_to_m*I_dx - G%dyCu(I,j) = m_to_L*dy_everywhere ; G%IdyCu(I,j) = L_to_m*I_dy + G%dxCu(I,j) = dx_everywhere ; G%IdxCu(I,j) = I_dx + G%dyCu(I,j) = dy_everywhere ; G%IdyCu(I,j) = I_dy enddo ; enddo do J=JsdB,JedB ; do i=isd,ied G%geoLonCv(i,J) = grid_lonT(i) ; G%geoLatCv(i,J) = grid_latB(J) - G%dxCv(i,J) = m_to_L*dx_everywhere ; G%IdxCv(i,J) = L_to_m*I_dx - G%dyCv(i,J) = m_to_L*dy_everywhere ; G%IdyCv(i,J) = L_to_m*I_dy + G%dxCv(i,J) = dx_everywhere ; G%IdxCv(i,J) = I_dx + G%dyCv(i,J) = dy_everywhere ; G%IdyCv(i,J) = I_dy enddo ; enddo call callTree_leave("set_grid_metrics_cartesian()") @@ -527,7 +541,7 @@ subroutine set_grid_metrics_spherical(G, param_file, US) real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) real :: dLon,dLat,latitude,longitude,dL_di - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical" is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -555,8 +569,8 @@ subroutine set_grid_metrics_spherical(G, param_file, US) call get_param(param_file, mdl, "LENLON", G%len_lon, & "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) - call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & + "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) dLon = G%len_lon/G%Domain%niglobal dLat = G%len_lat/G%Domain%njglobal @@ -600,9 +614,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxBu(I,J) = m_to_L*G%Rad_Earth * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di -! G%dxBu(I,J) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) - G%dyBu(I,J) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxBu(I,J) = G%Rad_Earth_L * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di +! G%dxBu(I,J) = G%Rad_Earth_L * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) + G%dyBu(I,J) = G%Rad_Earth_L * dLat*PI_180 G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) enddo ; enddo @@ -612,9 +626,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxCv(i,J) = m_to_L*G%Rad_Earth * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di -! G%dxCv(i,J) = m_to_L*G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) - G%dyCv(i,J) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxCv(i,J) = G%Rad_Earth_L * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di +! G%dxCv(i,J) = G%Rad_Earth_L * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) + G%dyCv(i,J) = G%Rad_Earth_L * dLat*PI_180 enddo ; enddo do j=jsd,jed ; do I=IsdB,IedB @@ -623,9 +637,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxCu(I,j) = m_to_L*G%Rad_Earth * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di -! G%dxCu(I,j) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( latitude ) - G%dyCu(I,j) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxCu(I,j) = G%Rad_Earth_L * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di +! G%dxCu(I,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) + G%dyCu(I,j) = G%Rad_Earth_L * dLat*PI_180 enddo ; enddo do j=jsd,jed ; do i=isd,ied @@ -634,13 +648,13 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxT(i,j) = m_to_L*G%Rad_Earth * COS( G%geoLatT(i,j)*PI_180 ) * dL_di -! G%dxT(i,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude ) - G%dyT(i,j) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxT(i,j) = G%Rad_Earth_L * COS( G%geoLatT(i,j)*PI_180 ) * dL_di +! G%dxT(i,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) + G%dyT(i,j) = G%Rad_Earth_L * dLat*PI_180 ! latitude = G%geoLatCv(i,J)*PI_180 ! In radians ! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians -! G%areaT(i,j) = m_to_L**2 * Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) +! G%areaT(i,j) = Rad_Earth_L**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) G%areaT(i,j) = G%dxT(i,j) * G%dyT(i,j) enddo ; enddo @@ -677,7 +691,7 @@ subroutine set_grid_metrics_mercator(G, param_file, US) real :: fnRef ! fnRef is the value of Int_dj_dy or ! Int_dj_dy at a latitude or longitude that is real :: jRef, iRef ! being set to be at grid index jRef or iRef. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] integer :: itt1, itt2 logical :: debug = .FALSE., simple_area = .true. integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB @@ -713,11 +727,12 @@ subroutine set_grid_metrics_mercator(G, param_file, US) call get_param(param_file, mdl, "LENLON", GP%len_lon, & "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) - call get_param(param_file, mdl, "RAD_EARTH", GP%Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + call get_param(param_file, mdl, "RAD_EARTH", GP%Rad_Earth_L, & + "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) G%south_lat = GP%south_lat ; G%len_lat = GP%len_lat G%west_lon = GP%west_lon ; G%len_lon = GP%len_lon - G%Rad_Earth = GP%Rad_Earth + G%Rad_Earth_L = GP%Rad_Earth_L + call get_param(param_file, mdl, "ISOTROPIC", GP%isotropic, & "If true, an isotropic grid on a sphere (also known as "//& "a Mercator grid) is used. With an isotropic grid, the "//& @@ -826,8 +841,8 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do J=JsdB,JedB ; do I=IsdB,IedB G%geoLonBu(I,J) = xq(I,J)*180.0/PI G%geoLatBu(I,J) = yq(I,J)*180.0/PI - G%dxBu(I,J) = m_to_L*ds_di(xq(I,J), yq(I,J), GP) - G%dyBu(I,J) = m_to_L*ds_dj(xq(I,J), yq(I,J), GP) + G%dxBu(I,J) = ds_di(xq(I,J), yq(I,J), GP) + G%dyBu(I,J) = ds_dj(xq(I,J), yq(I,J), GP) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) G%IareaBu(I,J) = 1.0 / (G%areaBu(I,J)) @@ -836,8 +851,8 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do j=jsd,jed ; do i=isd,ied G%geoLonT(i,j) = xh(i,j)*180.0/PI G%geoLatT(i,j) = yh(i,j)*180.0/PI - G%dxT(i,j) = m_to_L*ds_di(xh(i,j), yh(i,j), GP) - G%dyT(i,j) = m_to_L*ds_dj(xh(i,j), yh(i,j), GP) + G%dxT(i,j) = ds_di(xh(i,j), yh(i,j), GP) + G%dyT(i,j) = ds_dj(xh(i,j), yh(i,j), GP) G%areaT(i,j) = G%dxT(i,j)*G%dyT(i,j) G%IareaT(i,j) = 1.0 / (G%areaT(i,j)) @@ -846,20 +861,20 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do j=jsd,jed ; do I=IsdB,IedB G%geoLonCu(I,j) = xu(I,j)*180.0/PI G%geoLatCu(I,j) = yu(I,j)*180.0/PI - G%dxCu(I,j) = m_to_L*ds_di(xu(I,j), yu(I,j), GP) - G%dyCu(I,j) = m_to_L*ds_dj(xu(I,j), yu(I,j), GP) + G%dxCu(I,j) = ds_di(xu(I,j), yu(I,j), GP) + G%dyCu(I,j) = ds_dj(xu(I,j), yu(I,j), GP) enddo ; enddo do J=JsdB,JedB ; do i=isd,ied G%geoLonCv(i,J) = xv(i,J)*180.0/PI G%geoLatCv(i,J) = yv(i,J)*180.0/PI - G%dxCv(i,J) = m_to_L*ds_di(xv(i,J), yv(i,J), GP) - G%dyCv(i,J) = m_to_L*ds_dj(xv(i,J), yv(i,J), GP) + G%dxCv(i,J) = ds_di(xv(i,J), yv(i,J), GP) + G%dyCv(i,J) = ds_dj(xv(i,J), yv(i,J), GP) enddo ; enddo if (.not.simple_area) then do j=JsdB+1,jed ; do i=IsdB+1,ied - G%areaT(I,J) = m_to_L**2*GP%Rad_Earth**2 * & + G%areaT(I,J) = GP%Rad_Earth_L**2 * & (dL(xq(I-1,J-1),xq(I-1,J),yq(I-1,J-1),yq(I-1,J)) + & (dL(xq(I-1,J),xq(I,J),yq(I-1,J),yq(I,J)) + & (dL(xq(I,J),xq(I,J-1),yq(I,J),yq(I,J-1)) + & @@ -884,31 +899,31 @@ subroutine set_grid_metrics_mercator(G, param_file, US) end subroutine set_grid_metrics_mercator -!> This function returns the grid spacing in the logical x direction. +!> This function returns the grid spacing in the logical x direction in [L ~> m]. function ds_di(x, y, GP) real, intent(in) :: x !< The longitude in question real, intent(in) :: y !< The latitude in question type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: ds_di - ! Local variables - ds_di = GP%Rad_Earth * cos(y) * dx_di(x,GP) + real :: ds_di ! The returned grid spacing [L ~> m] + + ds_di = GP%Rad_Earth_L * cos(y) * dx_di(x,GP) ! In general, this might be... - ! ds_di = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + & + ! ds_di = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + & ! dy_di(x,y,GP)*dy_di(x,y,GP)) end function ds_di -!> This function returns the grid spacing in the logical y direction. +!> This function returns the grid spacing in the logical y direction in [L ~> m]. function ds_dj(x, y, GP) real, intent(in) :: x !< The longitude in question real, intent(in) :: y !< The latitude in question type(GPS), intent(in) :: GP !< A structure of grid parameters - ! Local variables - real :: ds_dj - ds_dj = GP%Rad_Earth * dy_dj(y,GP) + real :: ds_dj ! The returned grid spacing [L ~> m] + + ds_dj = GP%Rad_Earth_L * dy_dj(y,GP) ! In general, this might be... - ! ds_dj = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + & + ! ds_dj = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + & ! dy_dj(x,y,GP)*dy_dj(x,y,GP)) end function ds_dj @@ -1199,8 +1214,7 @@ subroutine initialize_masks(G, PF, US) 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 :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + real :: m_to_Z_scale ! A unit conversion factor from m to Z [Z m-1 ~> 1] real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. real :: min_depth ! The minimum ocean depth in the same units as G%bathyT [Z ~> m]. real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m]. @@ -1209,7 +1223,6 @@ subroutine initialize_masks(G, PF, US) call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") m_to_Z_scale = 1.0 ; if (present(US)) m_to_Z_scale = US%m_to_Z - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & "If MASKING_DEPTH is unspecified, then anything shallower than "//& diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 838e33b68d..49a9bfa6f2 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -95,7 +95,7 @@ subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G, US) type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables integer :: i,j - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] real :: f1, f2 m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L @@ -299,13 +299,13 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config. ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. + real :: m_to_Z ! A dimensional rescaling factor [Z m-1 ~> 1] + real :: m_to_L ! A dimensional rescaling factor [L m-1 ~> 1] real :: min_depth ! The minimum depth [Z ~> m]. 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 [m]. - real :: Dedge ! The depth [Z ~> m], at the basin edge -! real :: south_lat, west_lon, len_lon, len_lat, Rad_earth + real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m] + real :: expdecay ! A decay scale of associated with the sloping boundaries [L ~> m] + real :: Dedge ! The depth at the basin edge [Z ~> m] integer :: i, j, is, ie, js, je, isd, ied, jsd, jed character(len=40) :: mdl = "initialize_topography_named" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -316,6 +316,7 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth "TOPO_CONFIG = "//trim(topog_config), 5) m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z + m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) @@ -326,23 +327,9 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth 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, 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.) -! call get_param(param_file, mdl, "LENLAT", len_lat, & -! "The latitudinal length of the domain.", units="degrees", & -! fail_if_missing=.true.) -! call get_param(param_file, mdl, "WESTLON", west_lon, & -! "The western longitude of the domain.", units="degrees", & -! default=0.0) -! call get_param(param_file, mdl, "LENLON", len_lon, & -! "The longitudinal length of the domain.", units="degrees", & -! fail_if_missing=.true.) -! call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth, & -! "The radius of the Earth.", units="m", default=6.378e6) call get_param(param_file, mdl, "TOPOG_SLOPE_SCALE", expdecay, & "The exponential decay scale used in defining some of "//& - "the named topographies.", units="m", default=400000.0) + "the named topographies.", units="m", default=400000.0, scale=m_to_L) endif @@ -352,30 +339,30 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth do i=is,ie ; do j=js,je ; D(i,j) = max_depth ; enddo ; enddo elseif (trim(topog_config) == "spoon") then D0 = (max_depth - Dedge) / & - ((1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay))) * & - (1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay)))) + ((1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay))) * & + (1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay)))) do i=is,ie ; do j=js,je ! This sets a bowl shaped (sort of) bottom topography, with a ! ! maximum depth of max_depth. ! D(i,j) = Dedge + D0 * & (sin(PI * (G%geoLonT(i,j) - (G%west_lon)) / G%len_lon) * & - (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))*G%Rad_earth*PI / & + (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))*G%Rad_Earth_L*PI / & (180.0*expdecay)) )) enddo ; enddo elseif (trim(topog_config) == "bowl") then D0 = (max_depth - Dedge) / & - ((1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay))) * & - (1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay)))) + ((1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay))) * & + (1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay)))) ! This sets a bowl shaped (sort of) bottom topography, with a ! maximum depth of max_depth. do i=is,ie ; do j=js,je D(i,j) = Dedge + D0 * & (sin(PI * (G%geoLonT(i,j) - G%west_lon) / G%len_lon) * & - ((1.0 - exp(-(G%geoLatT(i,j) - G%south_lat)*G%Rad_Earth*PI/ & + ((1.0 - exp(-(G%geoLatT(i,j) - G%south_lat)*G%Rad_Earth_L*PI/ & (180.0*expdecay))) * & (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))* & - G%Rad_Earth*PI/(180.0*expdecay))))) + G%Rad_Earth_L*PI/(180.0*expdecay))))) enddo ; enddo elseif (trim(topog_config) == "halfpipe") then D0 = max_depth - Dedge @@ -511,10 +498,13 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) ! This subroutine sets up the Coriolis parameter for a beta-plane integer :: I, J 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 m-1 ~> s-1 m-1] - real :: beta_lat_ref ! The reference latitude for the beta plane [degrees/km/m/cm] - real :: y_scl, Rad_Earth - real :: T_to_s ! A time unit conversion factor + 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/km/m/cm] + 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 :: T_to_s ! A time unit conversion factor [s T-1 ~> 1] + real :: m_to_L ! A length unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A length unit conversion factor [m L-1 ~> 1] real :: PI character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name. character(len=200) :: axis_units @@ -523,28 +513,30 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") T_to_s = 1.0 ; if (present(US)) T_to_s = US%T_to_s + m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L + L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m call get_param(param_file, mdl, "F_0", f_0, & "The reference value of the Coriolis parameter with the "//& "betaplane option.", units="s-1", default=0.0, scale=T_to_s) call get_param(param_file, mdl, "BETA", beta, & "The northward gradient of the Coriolis parameter with "//& - "the betaplane option.", units="m-1 s-1", default=0.0, scale=T_to_s) + "the betaplane option.", units="m-1 s-1", default=0.0, scale=T_to_s*L_to_m) call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees") PI = 4.0*atan(1.0) select case (axis_units(1:1)) case ("d") - call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth_L, & + "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) beta_lat_ref_units = "degrees" - y_scl = PI * Rad_Earth/ 180. + y_scl = PI * Rad_Earth_L / 180. case ("k") beta_lat_ref_units = "kilometers" - y_scl = 1.E3 + y_scl = 1.E3 * m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1. + y_scl = 1. * m_to_L case default ; call MOM_error(FATAL, & " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units)) end select @@ -645,8 +637,8 @@ subroutine reset_face_lengths_named(G, param_file, name, US) ! Local variables character(len=256) :: mesg ! Message for error messages. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] real :: dx_2 = -1.0, dy_2 = -1.0 real :: pi_180 integer :: option = -1 @@ -773,8 +765,8 @@ subroutine reset_face_lengths_file(G, param_file, US) 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 - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -851,8 +843,8 @@ subroutine reset_face_lengths_list(G, param_file, US) integer, allocatable, dimension(:) :: & u_line_no, v_line_no, & ! The line numbers in lines of u- and v-face lines u_line_used, v_line_used ! The number of times each u- and v-line is used. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] + real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] real :: lat, lon ! The latitude and longitude of a point. real :: len_lon ! The periodic range of longitudes, usually 360 degrees. real :: len_lat ! The range of latitudes, usually 180 degrees. @@ -1023,7 +1015,7 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. & ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then - G%dy_Cu(I,j) = G%mask2dCu(I,j) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0)) + G%dy_Cu(I,j) = G%mask2dCu(I,j) * min(G%dyCu(I,j), max(m_to_L*u_width(npt), 0.0)) if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain if ( G%mask2dCu(I,j) == 0.0 ) then write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",& @@ -1053,7 +1045,7 @@ subroutine reset_face_lengths_list(G, param_file, US) (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. & ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. & ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then - G%dx_Cv(i,J) = G%mask2dCv(i,J) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0)) + G%dx_Cv(i,J) = G%mask2dCv(i,J) * min(G%dxCv(i,J), max(m_to_L*v_width(npt), 0.0)) if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain if ( G%mask2dCv(i,J) == 0.0 ) then write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index b9ceb85cc5..5342c621dd 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -641,7 +641,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & real :: surfFricVel, surfBuoyFlux real :: sigma, sigmaRatio - real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim] + real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: dh ! The local thickness used for calculating interface positions [m] real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m] @@ -948,7 +948,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: zBottomMinusOffset ! Height of bottom plus a little bit [m] real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. real :: hTot ! Running sum of thickness used in the surface layer average [m] - real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim] + real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: delH ! Thickness of a layer [m] real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index ab430fa3b9..379275196e 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -1196,7 +1196,7 @@ subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & real :: b1Kd ! Temporary array [nondim] real :: ColHt_chg ! The change in column thickness [Z ~> m]. real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m]. - real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> 1 or m3 kg-1]. + real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> nondim or m3 kg-1] real :: dT_k, dT_km1 ! Temperature changes in layers k and k-1 [degC] real :: dS_k, dS_km1 ! Salinity changes in layers k and k-1 [ppt] real :: I_Kr_denom ! Temporary array [H-2 ~> m-2 or m4 kg-2] diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 4bed5dc12e..25b06d51fa 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -639,7 +639,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs real :: LA ! The value of the Langmuir number [nondim] real :: LAmod ! The modified Langmuir number by convection [nondim] real :: hbs_here ! The local minimum of hb_hs and MixLen_shape, times a - ! conversion factor from H to Z [Z H-1 ~> 1 or m3 kg-1]. + ! conversion factor from H to Z [Z H-1 ~> nondim or m3 kg-1]. real :: nstar_FC ! The fraction of conv_PErel that can be converted to mixing [nondim]. real :: TKE_reduc ! The fraction by which TKE and other energy fields are ! reduced to support mixing [nondim]. between 0 and 1. @@ -1637,7 +1637,7 @@ subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & real :: b1Kd ! Temporary array [nondim] real :: ColHt_chg ! The change in column thickness [Z ~> m]. real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m]. - real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> 1 or m3 kg-1]. + real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> nondim or m3 kg-1] real :: dT_k, dT_km1 ! Temperature changes in layers k and k-1 [degC] real :: dS_k, dS_km1 ! Salinity changes in layers k and k-1 [ppt] real :: I_Kr_denom ! Temporary array [H-2 ~> m-2 or m4 kg-2] diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index f3efb1f661..9884eec7a2 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -507,7 +507,7 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith ! Local variables character(len=200) :: inputdir, geo_file, filename, geotherm_var real :: geo_scale ! A constant heat flux or dimensionally rescaled geothermal flux scaling factor - ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> 1] + ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> nondim] integer :: i, j, isd, ied, jsd, jed, id isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index f60986fb74..83558b9a08 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -178,7 +178,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! Rho0 divided by G_Earth and the conversion ! from m to thickness units [H R ~> kg m-2 or kg2 m-5]. real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion - ! factor from lateral lengths to vertical depths [Z L-1 ~> 1]. + ! factor from lateral lengths to vertical depths [Z L-1 ~> nondim]. real :: cdrag_sqrt ! Square root of the drag coefficient [nondim]. real :: oldfn ! The integrated energy required to ! entrain up to the bottom of the layer, @@ -257,7 +257,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) real :: Cell_width ! The transverse width of the velocity cell [L ~> m]. real :: Rayleigh ! A nondimensional value that is multiplied by the layer's ! velocity magnitude to give the Rayleigh drag velocity, times - ! a lateral to vertical distance conversion factor [Z L-1 ~> 1]. + ! a lateral to vertical distance conversion factor [Z L-1 ~> nondim]. real :: gam ! The ratio of the change in the open interface width ! to the open interface width atop a cell [nondim]. real :: BBL_frac ! The fraction of a layer's drag that goes into the @@ -951,10 +951,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! set kv_bbl to the bound and recompute bbl_thick to be consistent ! but with a ridiculously large upper bound on thickness (for Cd u*=0) kv_bbl = CS%Kv_BBL_min - if (cdrag_sqrt*ustar(i)*BBL_visc_frac*G%Rad_Earth*US%m_to_Z > kv_bbl) then + if (cdrag_sqrt*ustar(i)*BBL_visc_frac*G%Rad_Earth_L*US%L_to_Z > kv_bbl) then bbl_thick_Z = kv_bbl / ( cdrag_sqrt*ustar(i)*BBL_visc_frac ) else - bbl_thick_Z = G%Rad_Earth * US%m_to_Z + bbl_thick_Z = G%Rad_Earth_L * US%L_to_Z endif else kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_Z*BBL_visc_frac @@ -983,10 +983,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! set kv_bbl to the bound and recompute bbl_thick to be consistent ! but with a ridiculously large upper bound on thickness (for Cd u*=0) kv_bbl = CS%Kv_BBL_min - if (cdrag_sqrt*ustar(i)*G%Rad_Earth*US%m_to_Z > kv_bbl) then + if (cdrag_sqrt*ustar(i)*G%Rad_Earth_L*US%L_to_Z > kv_bbl) then bbl_thick_Z = kv_bbl / ( cdrag_sqrt*ustar(i) ) else - bbl_thick_Z = G%Rad_Earth * US%m_to_Z + bbl_thick_Z = G%Rad_Earth_L * US%L_to_Z endif else kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_Z @@ -1218,7 +1218,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) ! Rho0 divided by G_Earth and the conversion ! from m to thickness units [H R ~> kg m-2 or kg2 m-5]. real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion - ! factor from lateral lengths to vertical depths [Z L-1 ~> 1]. + ! factor from lateral lengths to vertical depths [Z L-1 ~> nondim] real :: cdrag_sqrt ! Square root of the drag coefficient [nondim]. real :: oldfn ! The integrated energy required to ! entrain up to the bottom of the layer, diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index a44750c1fc..4278594913 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -310,7 +310,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag ! Add other user-provided calls here. if (CS%use_USER_tracer_example) & - call USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS%USER_tracer_example_CSp, & + call USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS%USER_tracer_example_CSp, & sponge_CSp) if (CS%use_DOME_tracer) & call initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS%DOME_tracer_CSp, & diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index b58e45b366..10551ea247 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -134,13 +134,14 @@ end function USER_register_tracer_example !> This subroutine initializes the NTR tracer fields in tr(:,:,:,:) !! and it sets up the tracer output. -subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & +subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & sponge_CSp) logical, intent(in) :: restart !< .true. if the fields have already !! been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate @@ -163,7 +164,7 @@ subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & real, pointer :: tr_ptr(:,:,:) => NULL() real :: PI ! 3.1415926... calculated as 4*atan(1) real :: tr_y ! Initial zonally uniform tracer concentrations. - real :: dist2 ! The distance squared from a line [m2]. + real :: dist2 ! The distance squared from a line [L2 ~> m2]. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB, lntr @@ -196,9 +197,9 @@ subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & ! This sets a stripe of tracer across the basin. PI = 4.0*atan(1.0) do j=js,je - dist2 = (G%Rad_Earth * PI / 180.0)**2 * & + dist2 = (G%Rad_Earth_L * PI / 180.0)**2 * & (G%geoLatT(i,j) - 40.0) * (G%geoLatT(i,j) - 40.0) - tr_y = 0.5*exp(-dist2/(1.0e5*1.0e5)) + tr_y = 0.5 * exp( -dist2 / (1.0e5*US%m_to_L)**2 ) do k=1,nz ; do i=is,ie ! This adds the stripes of tracer to every layer. @@ -282,10 +283,10 @@ subroutine tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, ! Local variables real :: hold0(SZI_(G)) ! The original topmost layer thickness, - ! with surface mass fluxes added back, m. - real :: b1(SZI_(G)) ! b1 and c1 are variables used by the - real :: c1(SZI_(G),SZK_(GV)) ! tridiagonal solver. - real :: d1(SZI_(G)) ! d1=1-c1 is used by the tridiagonal solver. + ! with surface mass fluxes added back [H ~> m or kg m-2]. + real :: b1(SZI_(G)) ! b1 is a variable used by the tridiagonal solver [H ~> m or kg m-2]. + real :: c1(SZI_(G),SZK_(GV)) ! c1 is a variable used by the tridiagonal solver [nondim]. + real :: d1(SZI_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. @@ -374,7 +375,7 @@ function USER_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index) !! stocks calculated here. ! Local variables - real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or nondim] + real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1] integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke