21 implicit none ;
private
23 public dumbbell_dynamic_forcing, dumbbell_buoyancy_forcing, dumbbell_surface_forcing_init
27 logical :: use_temperature
29 logical :: restorebuoy
38 real,
dimension(:,:),
allocatable :: &
40 real,
dimension(:,:),
allocatable :: &
50 subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, CS)
53 type(
forcing),
intent(inout) :: fluxes
56 type(time_type),
intent(in) :: day
57 real,
intent(in) :: dt
65 real :: density_restore
68 integer :: i, j, is, ie, js, je
69 integer :: isd, ied, jsd, jed
71 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
72 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
76 if (cs%use_temperature)
then
95 if ( cs%use_temperature )
then
98 do j=js,je ;
do i=is,ie
101 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
102 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
105 fluxes%vprec(i,j) = 0.0
108 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
109 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
110 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
111 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
114 do j=js,je ;
do i=is,ie
117 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
121 if (cs%use_temperature .and. cs%restorebuoy)
then
122 do j=js,je ;
do i=is,ie
125 if (cs%forcing_mask(i,j)>0.)
then
126 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%Rho0*cs%Flux_const)) * &
127 ((cs%S_restore(i,j) - state%SSS(i,j)) / &
128 (0.5 * (cs%S_restore(i,j) + state%SSS(i,j))))
134 end subroutine dumbbell_buoyancy_forcing
137 subroutine dumbbell_dynamic_forcing(state, fluxes, day, dt, G, CS)
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))
176 end subroutine dumbbell_dynamic_forcing
179 subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
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)
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.", units=
"m day-1", &
242 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
268 end subroutine dumbbell_surface_forcing_init