BFB_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initialization of the boundary-forced-basing configuration
7
8use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
9use mom_file_parser, only : get_param, log_version, param_file_type
10use mom_get_input, only : directories
11use mom_grid, only : ocean_grid_type
13use mom_tracer_registry, only : tracer_registry_type
17implicit none ; private
18
19#include <MOM_memory.h>
20
21public bfb_set_coord
23
24! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
25! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
26! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
27! vary with the Boussinesq approximation, the Boussinesq variant is given first.
28
29contains
30
31!> This subroutine specifies the vertical coordinate in terms of temperature at the surface and at the bottom.
32!! This case is set up in such a way that the temperature of the topmost layer is equal to the SST at the
33!! southern edge of the domain. The temperatures are then converted to densities of the top and bottom layers
34!! and linearly interpolated for the intermediate layers.
35subroutine bfb_set_coord(Rlay, g_prime, GV, US, param_file)
36 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
37 real, dimension(GV%ke), intent(out) :: rlay !< Layer potential density [R ~> kg m-3].
38 real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each
39 !! interface [L2 Z-1 T-2 ~> m s-2].
40 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
41 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
42
43 real :: rho_t0_s0 ! The density at T=0, S=0 [R ~> kg m-3]
44 real :: drho_dt ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
45 real :: drho_ds ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
46 real :: sst_s, t_bot ! Temperatures at the surface and seafloor [C ~> degC]
47 real :: s_ref ! Reference salinity [S ~> ppt]
48 real :: rho_top, rho_bot ! Densities at the surface and seafloor [R ~> kg m-3]
49 integer :: k, nz
50 ! This include declares and sets the variable "version".
51# include "version_variable.h"
52 character(len=40) :: mdl = "BFB_initialization" ! This module's name.
53
54 call log_version(param_file, mdl, version, "")
55 call get_param(param_file, mdl, "DRHO_DT", drho_dt, &
56 "The partial derivative of density with temperature.", &
57 units="kg m-3 K-1", default=-0.2, scale=us%kg_m3_to_R*us%C_to_degC)
58 call get_param(param_file, mdl, "DRHO_DS", drho_ds, &
59 "The partial derivative of density with salinity.", &
60 units="kg m-3 PSU-1", default=0.8, scale=us%kg_m3_to_R*us%S_to_ppt)
61 call get_param(param_file, mdl, "RHO_T0_S0", rho_t0_s0, &
62 "The density at T=0, S=0.", units="kg m-3", default=1000.0, scale=us%kg_m3_to_R)
63 call get_param(param_file, mdl, "SST_S", sst_s, &
64 "SST at the southern edge of the domain.", &
65 units="degC", default=20.0, scale=us%degC_to_C)
66 call get_param(param_file, mdl, "T_BOT", t_bot, &
67 "Bottom temperature", units="degC", default=5.0, scale=us%degC_to_C)
68 call get_param(param_file, mdl, "S_REF", s_ref, &
69 "The initial salinities.", units="PSU", default=35.0, scale=us%ppt_to_S)
70 rho_top = (rho_t0_s0 + drho_ds*s_ref) + drho_dt*sst_s
71 rho_bot = (rho_t0_s0 + drho_ds*s_ref) + drho_dt*t_bot
72 nz = gv%ke
73
74 do k = 1,nz
75 rlay(k) = (rho_bot - rho_top)/(nz-1)*real(k-1) + rho_top
76 if (k==1) then
77 g_prime(k) = gv%g_Earth
78 elseif (gv%Boussinesq) then
79 g_prime(k) = (rlay(k) - rlay(k-1)) * gv%g_Earth / gv%Rho0
80 else
81 g_prime(k) = (rlay(k) - rlay(k-1)) * gv%g_Earth / (0.5*(rlay(k) + rlay(k-1)))
82 endif
83 enddo
84
85end subroutine bfb_set_coord
86
87!> This subroutine sets up the sponges for the southern boundary of the domain. Maximum damping occurs
88!! within 2 degrees lat of the boundary. The damping linearly decreases northward over the next 2 degrees.
89subroutine bfb_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, param_file, CSp, h)
90 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
91 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
92 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
93 logical, intent(in) :: use_temperature !< If true, temperature and salinity are used as
94 !! state variables.
95 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
96 real, dimension(SZI_(G),SZJ_(G)), &
97 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
98 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
99 type(sponge_cs), pointer :: csp !< A pointer to the sponge control structure
100 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
101 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
102
103 ! Local variables
104 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for eta, in depth units [Z ~> m].
105 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate [T-1 ~> s-1]
106 real :: h0(szk_(gv)) ! Resting layer thicknesses in depth units [Z ~> m].
107 real :: slat ! The southern latitude of the domain [degrees_N]
108 real :: wlon ! The western longitude of the domain [degrees_E]
109 real :: lenlat ! The latitudinal length of the domain [degrees_N]
110 real :: lenlon ! The longitudinal length of the domain [degrees_E]
111 real :: nlat ! The northern latitude of the domain [degrees_N]
112 real :: max_damping ! The maximum damping rate [T-1 ~> s-1]
113 ! This include declares and sets the variable "version".
114# include "version_variable.h"
115 character(len=40) :: mdl = "BFB_initialize_sponges_southonly" ! This subroutine's name.
116 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
117
118 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
119 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
120
121! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0
122! wherever there is no sponge, and the subroutines that are called
123! will automatically set up the sponges only where Idamp is positive
124! and mask2dT is 1.
125
126 ! Set up sponges for this configuration
127 ! call log_version(param_file, mdl, version)
128
129 slat = g%south_lat
130 lenlat = g%len_lat
131 wlon = g%west_lon
132 lenlon = g%len_lon
133 nlat = slat + lenlat
134 do k=1,nz ; h0(k) = -g%max_depth * real(k-1) / real(nz) ; enddo
135
136 ! Use for meridional thickness profile initialization
137 ! do k=1,nz ; H0(k) = -G%max_depth * real(k-1) / real(nz-1) ; enddo
138
139 max_damping = 1.0 / (86400.0*us%s_to_T)
140
141 eta(:,:,:) = 0.0 ; idamp(:,:) = 0.0
142
143 do j=js,je ; do i=is,ie
144 if (g%mask2dT(i,j) <= 0.0) then ; idamp(i,j) = 0.0
145 elseif (g%geoLatT(i,j) < slat+2.0) then ; idamp(i,j) = max_damping
146 elseif (g%geoLatT(i,j) < slat+4.0) then
147 idamp(i,j) = max_damping * (slat+4.0-g%geoLatT(i,j))/2.0
148 else ; idamp(i,j) = 0.0
149 endif
150
151 ! These will be streched inside of apply_sponge, so they can be in
152 ! depth space for Boussinesq or non-Boussinesq models.
153
154 ! This section is used for uniform thickness initialization
155 do k=1,nz ; eta(i,j,k) = h0(k) ; enddo
156
157 ! The below section is used for meridional temperature profile thickness initialization
158 ! do k=1,nz ; eta(i,j,k) = H0(k) ; enddo
159 ! if (G%geoLatT(i,j) > 40.0) then
160 ! do k = 1,nz
161 ! eta(i,j,k) = -G%Angstrom_Z*(k-1)
162 ! enddo
163 ! elseif (G%geoLatT(i,j) > 20.0) then
164 ! do k=1,nz
165 ! eta(i,j,k) = min(H0(k) + (G%geoLatT(i,j) - 20.0)*(G%max_depth - nz*G%Angstrom_Z)/20.0, &
166 ! -(k-1)*G%Angstrom_Z)
167 ! enddo
168 ! endif
169 eta(i,j,nz+1) = -g%max_depth
170
171 enddo ; enddo
172
173! This call sets up the damping rates and interface heights.
174! This sets the inverse damping timescale fields in the sponges. !
175 call initialize_sponge(idamp, eta, g, param_file, csp, gv)
176
177! Now register all of the fields which are damped in the sponge. !
178! By default, momentum is advected vertically within the sponge, but !
179! momentum is typically not damped within the sponge. !
180
182
183end module bfb_initialization