21 implicit none ;
private
23 public bfb_buoyancy_forcing, bfb_surface_forcing_init
28 logical :: use_temperature
29 logical :: restorebuoy
50 subroutine bfb_buoyancy_forcing(state, fluxes, day, dt, G, US, CS)
53 type(
forcing),
intent(inout) :: fluxes
56 type(time_type),
intent(in) :: day
57 real,
intent(in) :: dt
67 real :: density_restore
70 real :: buoy_rest_const
72 integer :: i, j, is, ie, js, je
73 integer :: isd, ied, jsd, jed
75 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
76 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
80 if (cs%use_temperature)
then
96 if ( cs%use_temperature )
then
99 do j=js,je ;
do i=is,ie
102 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
103 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
106 fluxes%vprec(i,j) = 0.0
109 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
110 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
111 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
112 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
115 do j=js,je ;
do i=is,ie
118 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
122 if (cs%restorebuoy)
then
123 if (cs%use_temperature)
then
127 call mom_error(fatal,
"User_buoyancy_surface_forcing: " // &
128 "Temperature and salinity restoring used without modification." )
130 rhoxcp = cs%Rho0 * fluxes%C_p
131 do j=js,je ;
do i=is,ie
137 fluxes%heat_added(i,j) = (g%mask2dT(i,j) * (rhoxcp * cs%Flux_const)) * &
138 (temp_restore - state%SST(i,j))
139 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%Rho0*cs%Flux_const)) * &
140 ((salin_restore - state%SSS(i,j)) / &
141 (0.5 * (salin_restore + state%SSS(i,j))))
150 buoy_rest_const = -1.0 * (cs%G_Earth * us%m_to_Z*us%T_to_s*cs%Flux_const) / cs%Rho0
152 do j=js,je ;
do i=is,ie
155 if (g%geoLatT(i,j) < cs%lfrslat)
then
156 temp_restore = cs%SST_s
157 elseif (g%geoLatT(i,j) > cs%lfrnlat)
then
158 temp_restore = cs%SST_n
160 temp_restore = (cs%SST_s - cs%SST_n)/(cs%lfrslat - cs%lfrnlat) * &
161 (g%geoLatT(i,j) - cs%lfrslat) + cs%SST_s
164 density_restore = temp_restore*cs%drho_dt + cs%Rho0
166 fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
167 (density_restore - state%sfc_density(i,j))
172 end subroutine bfb_buoyancy_forcing
175 subroutine bfb_surface_forcing_init(Time, G, US, param_file, diag, CS)
176 type(time_type),
intent(in) :: time
180 type(
diag_ctrl),
target,
intent(in) :: diag
184 #include "version_variable.h"
185 character(len=40) :: mdl =
"BFB_surface_forcing"
187 if (
associated(cs))
then
188 call mom_error(warning,
"BFB_surface_forcing_init called with an associated "// &
189 "control structure.")
197 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
198 "If true, Temperature and salinity are used as state "//&
199 "variables.", default=.true.)
201 call get_param(param_file, mdl,
"G_EARTH", cs%G_Earth, &
202 "The gravitational acceleration of the Earth.", &
203 units=
"m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
204 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
205 "The mean ocean density used with BOUSSINESQ true to "//&
206 "calculate accelerations and the mass for conservation "//&
207 "properties, or with BOUSSINSEQ false to convert some "//&
208 "parameters from vertical units of m to kg m-2.", &
209 units=
"kg m-3", default=1035.0)
210 call get_param(param_file, mdl,
"LFR_SLAT", cs%lfrslat, &
211 "Southern latitude where the linear forcing ramp begins.", &
212 units=
"degrees", default = 20.0)
213 call get_param(param_file, mdl,
"LFR_NLAT", cs%lfrnlat, &
214 "Northern latitude where the linear forcing ramp ends.", &
215 units=
"degrees", default = 40.0)
216 call get_param(param_file, mdl,
"SST_S", cs%SST_s, &
217 "SST at the southern edge of the linear forcing ramp.", &
218 units=
"C", default = 20.0)
219 call get_param(param_file, mdl,
"SST_N", cs%SST_n, &
220 "SST at the northern edge of the linear forcing ramp.", &
221 units=
"C", default = 10.0)
222 call get_param(param_file, mdl,
"DRHO_DT", cs%drho_dt, &
223 "The rate of change of density with temperature.", &
224 units=
"kg m-3 K-1", default = -0.2)
225 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
226 "The background gustiness in the winds.", units=
"Pa", &
229 call get_param(param_file, mdl,
"RESTOREBUOY", cs%restorebuoy, &
230 "If true, the buoyancy fluxes drive the model back "//&
231 "toward some specified surface state with a rate "//&
232 "given by FLUXCONST.", default= .false.)
233 if (cs%restorebuoy)
then
234 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
235 "The constant that relates the restoring surface fluxes "//&
236 "to the relative surface anomalies (akin to a piston "//&
237 "velocity). Note the non-MKS units.", units=
"m day-1", &
238 fail_if_missing=.true.)
240 cs%Flux_const = cs%Flux_const / 86400.0
243 end subroutine bfb_surface_forcing_init