21 implicit none ;
private
28 logical :: use_temperature
29 logical :: restorebuoy
53 type(
forcing),
intent(inout) :: fluxes
56 type(time_type),
intent(in) :: day
57 real,
intent(in) :: dt
67 real :: density_restore
71 real :: buoy_rest_const
73 integer :: i, j, is, ie, js, je
74 integer :: isd, ied, jsd, jed
76 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
77 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
81 if (cs%use_temperature)
then
97 if ( cs%use_temperature )
then
100 do j=js,je ;
do i=is,ie
103 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
104 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
107 fluxes%vprec(i,j) = 0.0
110 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
111 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
112 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
113 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
116 do j=js,je ;
do i=is,ie
119 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
123 if (cs%restorebuoy)
then
124 if (cs%use_temperature)
then
128 call mom_error(fatal,
"User_buoyancy_surface_forcing: " // &
129 "Temperature and salinity restoring used without modification." )
131 rhoxcp = us%R_to_kg_m3*us%Z_to_m*us%s_to_T * cs%Rho0 * fluxes%C_p
132 do j=js,je ;
do i=is,ie
138 fluxes%heat_added(i,j) = (g%mask2dT(i,j) * (rhoxcp * cs%Flux_const)) * &
139 (temp_restore - state%SST(i,j))
140 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%Rho0*cs%Flux_const)) * &
141 ((salin_restore - state%SSS(i,j)) / (0.5 * (salin_restore + state%SSS(i,j))))
150 buoy_rest_const = -1.0 * (cs%G_Earth * 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 - us%kg_m3_to_R*state%sfc_density(i,j))
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, scale=us%kg_m3_to_R)
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, scale=us%kg_m3_to_R)
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.", &
238 units=
"m day-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
240 cs%Flux_const = cs%Flux_const / 86400.0