From b8e85d9da0b5b40588e50324ef0089f74e895622 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 18 Nov 2022 07:22:10 -0500 Subject: [PATCH] +Rossby_front_2d_init and BFB_init cleanup Completed the dimensional rescaling of internal variables in the Rossby_front_initialization and BFB_initialization modules, including the addition of comments describing numerous internal variables and their units. The BFB module was more extensively modified to place logging of the module name and version before the get_param calls for this module, following the MOM6 pattern, to correct spelling errors in get_param descriptions, and to use the grid mask in the grid file rather than a reexamination of the minimum depth to determine the land-mask. The internal subroutine write_BFB_log is no longer needed and has been folded into BFB_set_coord. All answers are bitwise identical, but there are minor changes in the MOM_parameter_doc files for the buoy_forced_basin test case. --- src/user/BFB_initialization.F90 | 68 +++++++----------- src/user/BFB_surface_forcing.F90 | 4 +- src/user/Rossby_front_2d_initialization.F90 | 80 ++++++++++++--------- 3 files changed, 72 insertions(+), 80 deletions(-) diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index 68a6b6530b..3efc908ffb 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -24,10 +24,6 @@ module BFB_initialization ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units ! vary with the Boussinesq approximation, the Boussinesq variant is given first. -!> Unsafe model variable -!! \todo Remove this module variable -logical :: first_call = .true. - contains !> This subroutine specifies the vertical coordinate in terms of temperature at the surface and at the bottom. @@ -42,17 +38,22 @@ subroutine BFB_set_coord(Rlay, g_prime, GV, US, param_file) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - real :: drho_dt, SST_s, T_bot, rho_top, rho_bot - integer :: k, nz - character(len=40) :: mdl = "BFB_set_coord" ! This subroutine's name. + real :: dRho_dT ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] + real :: SST_s, T_bot ! Temperatures at the surface and seafloor [C ~> degC] + real :: rho_top, rho_bot ! Densities at the surface and seafloor [R ~> kg m-3] + integer :: k, nz + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=40) :: mdl = "BFB_initialization" ! This module's name. + call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "DRHO_DT", drho_dt, & "Rate of change of density with temperature.", & - units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R) + units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC) call get_param(param_file, mdl, "SST_S", SST_s, & - "SST at the suothern edge of the domain.", units="C", default=20.0) + "SST at the southern edge of the domain.", units="degC", default=20.0, scale=US%degC_to_C) call get_param(param_file, mdl, "T_BOT", T_bot, & - "Bottom Temp", units="C", default=5.0) + "Bottom temperature", units="degC", default=5.0, scale=US%degC_to_C) rho_top = GV%Rho0 + drho_dt*SST_s rho_bot = GV%Rho0 + drho_dt*T_bot nz = GV%ke @@ -64,15 +65,11 @@ subroutine BFB_set_coord(Rlay, g_prime, GV, US, param_file) else g_prime(k) = GV%g_Earth endif - !Rlay(:) = 0.0 - !g_prime(:) = 0.0 enddo - if (first_call) call write_BFB_log(param_file) - end subroutine BFB_set_coord -!> This subroutine sets up the sponges for the southern bouundary of the domain. Maximum damping occurs +!> This subroutine sets up the sponges for the southern boundary of the domain. Maximum damping occurs !! within 2 degrees lat of the boundary. The damping linearly decreases northward over the next 2 degrees. subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, param_file, CSp, h) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure @@ -92,29 +89,27 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for eta, in depth units [Z ~> m]. real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: H0(SZK_(GV)) ! Resting layer thicknesses in depth units [Z ~> m]. - real :: min_depth ! The minimum ocean depth in depth units [Z ~> m]. real :: slat ! The southern latitude of the domain [degrees_N] real :: wlon ! The western longitude of the domain [degrees_E] real :: lenlat ! The latitudinal length of the domain [degrees_N] real :: lenlon ! The longitudinal length of the domain [degrees_E] real :: nlat ! The northern latitude of the domain [degrees_N] real :: max_damping ! The maximum damping rate [T-1 ~> s-1] + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "BFB_initialize_sponges_southonly" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz 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 - eta(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 - ! 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 set up the sponges only where Idamp is positive ! and mask2dT is 1. -! Set up sponges for DOME configuration - call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) + ! Set up sponges for this configuration + ! call log_version(param_file, mdl, version) slat = G%south_lat lenlat = G%len_lat @@ -124,12 +119,14 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept do k=1,nz ; H0(k) = -G%max_depth * real(k-1) / real(nz) ; enddo ! Use for meridional thickness profile initialization -! do k=1,nz ; H0(k) = -G%max_depth * real(k-1) / real(nz-1) ; enddo + ! do k=1,nz ; H0(k) = -G%max_depth * real(k-1) / real(nz-1) ; enddo max_damping = 1.0 / (86400.0*US%s_to_T) + eta(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 + do j=js,je ; do i=is,ie - if (depth_tot(i,j) <= min_depth) then ; Idamp(i,j) = 0.0 + if (G%mask2dT(i,j) <= 0.0) then ; Idamp(i,j) = 0.0 elseif (G%geoLatT(i,j) < slat+2.0) then ; Idamp(i,j) = max_damping elseif (G%geoLatT(i,j) < slat+4.0) then Idamp(i,j) = max_damping * (slat+4.0-G%geoLatT(i,j))/2.0 @@ -140,16 +137,16 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept ! depth space for Boussinesq or non-Boussinesq models. ! This section is used for uniform thickness initialization - do k = 1,nz; eta(i,j,k) = H0(k); enddo + do k=1,nz ; eta(i,j,k) = H0(k) ; enddo - ! The below section is used for meridional temperature profile thickness initiation - ! do k = 1,nz; eta(i,j,k) = H0(k); enddo + ! The below section is used for meridional temperature profile thickness initialization + ! do k=1,nz ; eta(i,j,k) = H0(k) ; enddo ! if (G%geoLatT(i,j) > 40.0) then ! do k = 1,nz ! eta(i,j,k) = -G%Angstrom_Z*(k-1) ! enddo ! elseif (G%geoLatT(i,j) > 20.0) then - ! do k = 1,nz + ! do k=1,nz ! eta(i,j,k) = min(H0(k) + (G%geoLatT(i,j) - 20.0)*(G%max_depth - nz*G%Angstrom_Z)/20.0, & ! -(k-1)*G%Angstrom_Z) ! enddo @@ -166,23 +163,6 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept ! By default, momentum is advected vertically within the sponge, but ! ! momentum is typically not damped within the sponge. ! - if (first_call) call write_BFB_log(param_file) - end subroutine BFB_initialize_sponges_southonly -!> Write output about the parameter values being used. -subroutine write_BFB_log(param_file) - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - -! This include declares and sets the variable "version". -#include "version_variable.h" - character(len=40) :: mdl = "BFB_initialization" ! This module's name. - - call log_version(param_file, mdl, version) - first_call = .false. - -end subroutine write_BFB_log - end module BFB_initialization diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 6f16bdd6f0..38361ab070 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -212,10 +212,10 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS) units="degrees", default=40.0) call get_param(param_file, mdl, "SST_S", CS%SST_s, & "SST at the southern edge of the linear forcing ramp.", & - units="C", default=20.0, scale=US%degC_to_C) + units="degC", default=20.0, scale=US%degC_to_C) call get_param(param_file, mdl, "SST_N", CS%SST_n, & "SST at the northern edge of the linear forcing ramp.", & - units="C", default=10.0, scale=US%degC_to_C) + units="degC", default=10.0, scale=US%degC_to_C) call get_param(param_file, mdl, "DRHO_DT", CS%drho_dt, & "The rate of change of density with temperature.", & units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC) diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 2d0dcb85e5..d3d1ad2368 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -47,9 +47,13 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read !! parameters without changing h. integer :: i, j, k, is, ie, js, je, nz - real :: Tz, Dml, eta, stretch, h0 - real :: min_thickness, T_range - real :: dRho_dT ! The partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] + real :: Tz ! Vertical temperature gradient [C Z-1 ~> degC m-1] + real :: Dml ! Mixed layer depth [Z ~> m] + real :: eta ! An interface height depth [Z ~> m] + real :: stretch ! A nondimensional stretching factor [nondim] + real :: h0 ! The stretched thickness per layer [Z ~> m] + real :: T_range ! Range of temperatures over the vertical [C ~> degC] + real :: dRho_dT ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -59,13 +63,12 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read if (.not.just_read) call log_version(param_file, mdl, version, "") ! Read parameters needed to set thickness - call get_param(param_file, mdl, "MIN_THICKNESS", min_thickness, & - 'Minimum layer thickness',units='m',default=1.e-3, do_not_log=just_read) call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file, mdl, "T_RANGE", T_range, 'Initial temperature range', & - units='C', default=0.0, do_not_log=just_read) - call get_param(param_file, mdl, "DRHO_DT", dRho_dT, default=-0.2, scale=US%kg_m3_to_R, do_not_log=.true.) + units='C', default=0.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(param_file, mdl, "DRHO_DT", dRho_dT, & + units="kg m-3 degC-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC, do_not_log=.true.) if (just_read) return ! All run-time parameters have been read, so return. @@ -76,7 +79,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read case (REGRIDDING_LAYER, REGRIDDING_RHO) do j = G%jsc,G%jec ; do i = G%isc,G%iec Dml = Hml( G, G%geoLatT(i,j) ) - eta = -( -dRho_DT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) + eta = -( -dRho_dT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz @@ -87,7 +90,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read case (REGRIDDING_ZSTAR, REGRIDDING_SIGMA) do j = G%jsc,G%jec ; do i = G%isc,G%iec Dml = Hml( G, G%geoLatT(i,j) ) - eta = -( -dRho_DT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) + eta = -( -dRho_dT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz @@ -118,15 +121,18 @@ subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, & !! only read parameters without changing T & S. integer :: i, j, k, is, ie, js, je, nz - real :: T_ref, S_ref ! Reference salinity and temerature within surface layer - real :: T_range ! Range of salinities and temperatures over the vertical - real :: zc, zi, dTdz + real :: T_ref ! Reference temperature within the surface layer [C ~> degC] + real :: S_ref ! Reference salinity within the surface layer [S ~> [ppt] + real :: T_range ! Range of temperatures over the vertical [C ~> degC] + real :: zc ! Position of the middle of the cell [Z ~> m] + real :: zi ! Bottom interface position relative to the sea surface [H ~> m or kg m-2] + real :: dTdz ! Vertical temperature gradient [C Z-1 ~> degC m-1] character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file, mdl, "S_REF", S_ref, 'Reference salinity', & default=35.0, units='1e-3', scale=US%ppt_to_S, do_not_log=just_read) call get_param(param_file, mdl,"T_REF",T_ref,'Reference temperature', & @@ -169,12 +175,13 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just logical, intent(in) :: just_read !< If present and true, this call will only !! read parameters without setting u & v. - real :: T_range ! Range of salinities and temperatures over the vertical - real :: dUdT ! Factor to convert dT/dy into dU/dz, g*alpha/f [L2 Z-1 T-1 degC-1 ~> m s-1 degC-1] - real :: dRho_dT ! The partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] - real :: Dml, zi, zc, zm ! Depths [Z ~> m]. + real :: T_range ! Range of temperatures over the vertical [C ~> degC] + real :: dUdT ! Factor to convert dT/dy into dU/dz, g*alpha/f [L2 Z-1 T-1 C-1 ~> m s-1 degC-1] + real :: dRho_dT ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] + real :: Dml ! Mixed layer depth [Z ~> m] + real :: zi, zc, zm ! Depths [Z ~> m]. real :: f ! The local Coriolis parameter [T-1 ~> s-1] - real :: Ty ! The meridional temperature gradient [degC L-1 ~> degC m-1] + real :: Ty ! The meridional temperature gradient [C L-1 ~> degC m-1] real :: hAtU ! Interpolated layer thickness [Z ~> m]. integer :: i, j, k, is, ie, js, je, nz character(len=40) :: verticalCoordinate @@ -184,8 +191,9 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file, mdl, "T_RANGE", T_range, 'Initial temperature range', & - units='C', default=0.0, do_not_log=just_read) - call get_param(param_file, mdl, "DRHO_DT", dRho_dT, default=-0.2, scale=US%kg_m3_to_R, do_not_log=.true.) + units='C', default=0.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(param_file, mdl, "DRHO_DT", dRho_dT, & + units='kg m-3 degC-1', default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC, do_not_log=.true.) if (just_read) return ! All run-time parameters have been read, so return. @@ -197,7 +205,7 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just dUdT = 0.0 ; if (abs(f) > 0.0) & dUdT = ( GV%g_Earth*dRho_dT ) / ( f * GV%Rho0 ) Dml = Hml( G, G%geoLatT(i,j) ) - Ty = US%L_to_m*dTdy( G, T_range, G%geoLatT(i,j) ) + Ty = dTdy( G, T_range, G%geoLatT(i,j), US ) zi = 0. do k = 1, nz hAtU = 0.5*(h(i,j,k)+h(i+1,j,k)) * GV%H_to_Z @@ -212,12 +220,12 @@ end subroutine Rossby_front_initialize_velocity !> Pseudo coordinate across domain used by Hml() and dTdy() !! returns a coordinate from -PI/2 .. PI/2 squashed towards the -!! center of the domain. +!! center of the domain [radians]. real function yPseudo( G, lat ) type(ocean_grid_type), intent(in) :: G !< Grid structure - real, intent(in) :: lat !< Latitude + real, intent(in) :: lat !< Latitude in arbitrary units, often [km] ! Local - real :: PI + real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] PI = 4.0 * atan(1.0) yPseudo = ( ( lat - G%south_lat ) / G%len_lat ) - 0.5 ! -1/2 .. 1/.2 @@ -226,12 +234,12 @@ end function yPseudo !> Analytic prescription of mixed layer depth in 2d Rossby front test, -!! in the same units as G%max_depth +!! in the same units as G%max_depth (usually [Z ~> m]) real function Hml( G, lat ) type(ocean_grid_type), intent(in) :: G !< Grid structure - real, intent(in) :: lat !< Latitude + real, intent(in) :: lat !< Latitude in arbitrary units, often [km] ! Local - real :: dHML, HMLmean + real :: dHML, HMLmean ! The range and mean of the mixed layer depths [Z ~> m] dHML = 0.5 * ( HMLmax - HMLmin ) * G%max_depth HMLmean = 0.5 * ( HMLmin + HMLmax ) * G%max_depth @@ -239,18 +247,22 @@ real function Hml( G, lat ) end function Hml -!> Analytic prescription of mixed layer temperature gradient in 2d Rossby front test -real function dTdy( G, dT, lat ) +!> Analytic prescription of mixed layer temperature gradient in [C L-1 ~> degC m-1] in 2d Rossby front test +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 - real, intent(in) :: lat !< Latitude + real, intent(in) :: dT !< Top to bottom temperature difference [C ~> degC] + real, intent(in) :: lat !< Latitude in [km] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local - real :: PI, dHML, dHdy - real :: km = 1.e3 ! AXIS_UNITS = 'k' (1000 m) + 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] 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 ) ) * cos( yPseudo(G, lat) ) + dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * km_to_L ) ) * cos( yPseudo(G, lat) ) dTdy = -( dT / G%max_depth ) * dHdy end function dTdy