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