21 implicit none ;
private
27 logical :: use_temperature
29 logical :: restorebuoy
38 real,
dimension(:,:),
allocatable :: &
40 real,
dimension(:,:),
allocatable :: &
53 type(
forcing),
intent(inout) :: fluxes
56 type(time_type),
intent(in) :: day
57 real,
intent(in) :: dt
66 real :: density_restore
69 integer :: i, j, is, ie, js, je
70 integer :: isd, ied, jsd, jed
72 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
73 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
77 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%use_temperature .and. cs%restorebuoy)
then
123 do j=js,je ;
do i=is,ie
126 if (cs%forcing_mask(i,j)>0.)
then
127 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%Rho0*cs%Flux_const)) * &
128 ((cs%S_restore(i,j) - state%SSS(i,j)) / (0.5 * (cs%S_restore(i,j) + state%SSS(i,j))))
140 type(
forcing),
intent(inout) :: fluxes
143 type(time_type),
intent(in) :: day
144 real,
intent(in) :: dt
150 integer :: i, j, is, ie, js, je
151 integer :: isd, ied, jsd, jed
152 integer :: idays, isecs
153 real :: deg_rad, rdays
156 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
157 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
159 deg_rad = atan(1.0)*4.0/180.
161 call get_time(day,isecs,idays)
162 rdays = real(idays) + real(isecs)/8.64e4
169 do j=js,je ;
do i=is,ie
170 fluxes%p_surf(i,j) = cs%forcing_mask(i,j)* cs%slp_amplitude * &
171 g%mask2dT(i,j) * sin(deg_rad*(rdays/cs%slp_period))
172 fluxes%p_surf_full(i,j) = cs%forcing_mask(i,j) * cs%slp_amplitude * &
173 g%mask2dT(i,j) * sin(deg_rad*(rdays/cs%slp_period))
180 type(time_type),
intent(in) :: time
184 type(
diag_ctrl),
target,
intent(in) :: diag
189 real :: s_surf, s_range
192 #include "version_variable.h"
193 character(len=40) :: mdl =
"dumbbell_surface_forcing"
195 if (
associated(cs))
then
196 call mom_error(warning,
"dumbbell_surface_forcing_init called with an associated "// &
197 "control structure.")
205 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
206 "If true, Temperature and salinity are used as state "//&
207 "variables.", default=.true.)
209 call get_param(param_file, mdl,
"G_EARTH", cs%G_Earth, &
210 "The gravitational acceleration of the Earth.", &
211 units=
"m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
212 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
213 "The mean ocean density used with BOUSSINESQ true to "//&
214 "calculate accelerations and the mass for conservation "//&
215 "properties, or with BOUSSINSEQ false to convert some "//&
216 "parameters from vertical units of m to kg m-2.", &
217 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
218 call get_param(param_file, mdl,
"DUMBBELL_SLP_AMP", cs%slp_amplitude, &
219 "Amplitude of SLP forcing in reservoirs.", &
220 units=
"kg m2 s-1", default = 10000.0)
221 call get_param(param_file, mdl,
"DUMBBELL_SLP_PERIOD", cs%slp_period, &
222 "Periodicity of SLP forcing in reservoirs.", &
223 units=
"days", default = 1.0)
224 call get_param(param_file, mdl,
"DUMBBELL_SLP_PERIOD", cs%slp_period, &
225 "Periodicity of SLP forcing in reservoirs.", &
226 units=
"days", default = 1.0)
227 call get_param(param_file, mdl,
"INITIAL_SSS", s_surf, &
228 "Initial surface salinity", units=
"1e-3", default=34.0, do_not_log=.true.)
229 call get_param(param_file, mdl,
"INITIAL_S_RANGE", s_range, &
230 "Initial salinity range (bottom - surface)", units=
"1e-3", &
231 default=2., do_not_log=.true.)
233 call get_param(param_file, mdl,
"RESTOREBUOY", cs%restorebuoy, &
234 "If true, the buoyancy fluxes drive the model back "//&
235 "toward some specified surface state with a rate "//&
236 "given by FLUXCONST.", default= .false.)
237 if (cs%restorebuoy)
then
238 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
239 "The constant that relates the restoring surface fluxes "//&
240 "to the relative surface anomalies (akin to a piston "//&
241 "velocity). Note the non-MKS units.", &
242 units=
"m day-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
244 cs%Flux_const = cs%Flux_const / 86400.0
247 allocate(cs%forcing_mask(g%isd:g%ied, g%jsd:g%jed)); cs%forcing_mask(:,:)=0.0
248 allocate(cs%S_restore(g%isd:g%ied, g%jsd:g%jed))
253 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon - 0.5
254 y = ( g%geoLatT(i,j) - g%south_lat ) / g%len_lat - 0.5
255 cs%forcing_mask(i,j)=0
256 cs%S_restore(i,j) = s_surf
258 cs%forcing_mask(i,j) = 1
259 cs%S_restore(i,j) = s_surf + s_range
260 elseif ((x<-0.25))
then
261 cs%forcing_mask(i,j) = 1
262 cs%S_restore(i,j) = s_surf - s_range