21 implicit none ;
private
34 logical :: use_temperature
35 logical :: restorebuoy
55 type(time_type),
intent(in) :: day
62 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
66 call mom_error(fatal,
"User_wind_surface_forcing: " // &
67 "User forcing routine called without modification." )
69 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
70 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
80 do j=js,je ;
do i=is-1,ieq
82 forces%taux(i,j) = g%mask2dCu(i,j) * 0.0*us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
84 do j=js-1,jeq ;
do i=is,ie
85 forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0
90 if (
associated(forces%ustar))
then ;
do j=js,je ;
do i=is,ie
92 forces%ustar(i,j) = g%mask2dT(i,j) * sqrt((cs%gust_const + &
93 sqrt(0.5*(forces%taux(i-1,j)**2 + forces%taux(i,j)**2) + &
94 0.5*(forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2))) * (us%L_to_Z/cs%Rho0))
105 type(
forcing),
intent(inout) :: fluxes
106 type(time_type),
intent(in) :: day
107 real,
intent(in) :: dt
130 real :: salin_restore
131 real :: density_restore
135 real :: buoy_rest_const
138 integer :: i, j, is, ie, js, je
139 integer :: isd, ied, jsd, jed
141 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
142 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
143 rho0_mks = cs%Rho0 * us%R_to_kg_m3
147 call mom_error(fatal,
"User_buoyancy_surface_forcing: " // &
148 "User forcing routine called without modification." )
152 if (cs%use_temperature)
then
153 call safe_alloc_ptr(fluxes%evap, isd, ied, jsd, jed)
154 call safe_alloc_ptr(fluxes%lprec, isd, ied, jsd, jed)
155 call safe_alloc_ptr(fluxes%fprec, isd, ied, jsd, jed)
156 call safe_alloc_ptr(fluxes%lrunoff, isd, ied, jsd, jed)
157 call safe_alloc_ptr(fluxes%frunoff, isd, ied, jsd, jed)
158 call safe_alloc_ptr(fluxes%vprec, isd, ied, jsd, jed)
160 call safe_alloc_ptr(fluxes%sw, isd, ied, jsd, jed)
161 call safe_alloc_ptr(fluxes%lw, isd, ied, jsd, jed)
162 call safe_alloc_ptr(fluxes%latent, isd, ied, jsd, jed)
163 call safe_alloc_ptr(fluxes%sens, isd, ied, jsd, jed)
165 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
171 if ( cs%use_temperature )
then
174 do j=js,je ;
do i=is,ie
177 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
178 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
181 fluxes%vprec(i,j) = 0.0
184 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
185 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
186 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
187 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
190 do j=js,je ;
do i=is,ie
193 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
197 if (cs%restorebuoy)
then
198 if (cs%use_temperature)
then
199 call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
202 call mom_error(fatal,
"User_buoyancy_surface_forcing: " // &
203 "Temperature and salinity restoring used without modification." )
205 rhoxcp = rho0_mks * fluxes%C_p
206 do j=js,je ;
do i=is,ie
212 fluxes%heat_added(i,j) = (g%mask2dT(i,j) * (rhoxcp * us%Z_to_m*us%s_to_T*cs%Flux_const)) * &
213 (temp_restore - sfc_state%SST(i,j))
214 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%Rho0*cs%Flux_const)) * &
215 ((salin_restore - sfc_state%SSS(i,j)) / (0.5 * (salin_restore + sfc_state%SSS(i,j))))
220 call mom_error(fatal,
"User_buoyancy_surface_forcing: " // &
221 "Buoyancy restoring used without modification." )
224 buoy_rest_const = -1.0 * (cs%G_Earth * cs%Flux_const) / rho0_mks
225 do j=js,je ;
do i=is,ie
228 density_restore = 1030.0
230 fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
231 (density_restore - sfc_state%sfc_density(i,j))
240 type(time_type),
intent(in) :: time
244 type(
diag_ctrl),
target,
intent(in) :: diag
249 #include "version_variable.h"
250 character(len=40) :: mdl =
"user_surface_forcing"
252 if (
associated(cs))
then
253 call mom_error(warning,
"USER_surface_forcing_init called with an associated "// &
254 "control structure.")
262 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
263 "If true, Temperature and salinity are used as state "//&
264 "variables.", default=.true.)
266 call get_param(param_file, mdl,
"G_EARTH", cs%G_Earth, &
267 "The gravitational acceleration of the Earth.", &
268 units=
"m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
269 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
270 "The mean ocean density used with BOUSSINESQ true to "//&
271 "calculate accelerations and the mass for conservation "//&
272 "properties, or with BOUSSINSEQ false to convert some "//&
273 "parameters from vertical units of m to kg m-2.", &
274 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
275 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
276 "The background gustiness in the winds.", &
277 units=
"Pa", default=0.02, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
279 call get_param(param_file, mdl,
"RESTOREBUOY", cs%restorebuoy, &
280 "If true, the buoyancy fluxes drive the model back "//&
281 "toward some specified surface state with a rate "//&
282 "given by FLUXCONST.", default= .false.)
283 if (cs%restorebuoy)
then
284 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
285 "The constant that relates the restoring surface fluxes "//&
286 "to the relative surface anomalies (akin to a piston "//&
287 "velocity). Note the non-MKS units.", &
288 units=
"m day-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
290 cs%Flux_const = cs%Flux_const / 86400.0