Skip to content

+Rossby_front_2d_init and BFB_init cleanup #248

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 2 commits into from
Nov 28, 2022
Merged

+Rossby_front_2d_init and BFB_init cleanup #248

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b8e85d9
+Rossby_front_2d_init and BFB_init cleanup
Hallberg-NOAA Nov 18, 2022
030aecf
Merge branch 'dev/gfdl' into Rossby_front_cleanup
marshallward Nov 25, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 24 additions & 44 deletions src/user/BFB_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
4 changes: 2 additions & 2 deletions src/user/BFB_surface_forcing.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
80 changes: 46 additions & 34 deletions src/user/Rossby_front_2d_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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', &
Expand Down Expand Up @@ -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
Expand All @@ -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.

Expand All @@ -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
Expand All @@ -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
Expand All @@ -226,31 +234,35 @@ 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
Hml = HMLmean + dHML * sin( yPseudo(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
Expand Down