MOM_surface_forcing.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Functions that calculate the surface wind stresses and fluxes of buoyancy
6!! or temperature/salinity and fresh water, in ocean-only (solo) mode.
7!!
8!! These functions are called every time step, even if the wind stresses
9!! or buoyancy fluxes are constant in time - in that case these routines
10!! return quickly without doing anything. In addition, any I/O of forcing
11!! fields is controlled by surface_forcing_init, located in this file.
13
14use mom_constants, only : hlv, hlf
15use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
16use mom_cpu_clock, only : clock_module
17use mom_data_override, only : data_override_init, data_override
18use mom_diag_mediator, only : post_data, query_averaging_enabled
19use mom_diag_mediator, only : diag_ctrl, safe_alloc_ptr
20use mom_domains, only : pass_var, pass_vector, agrid, to_south, to_west, to_all
21use mom_domains, only : fill_symmetric_edges, cgrid_ne
22use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
24use mom_file_parser, only : get_param, log_param, log_version, param_file_type
30use mom_forcing_type, only : allocate_forcing_type, deallocate_forcing_type
31use mom_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing
32use mom_grid, only : ocean_grid_type
34use mom_io, only : file_exists, mom_read_data, mom_read_vector, slasher
35use mom_io, only : read_netcdf_data, east_face, north_face, num_timelevels
36use mom_restart, only : register_restart_field, restart_init, mom_restart_cs
38use mom_time_manager, only : time_type, operator(+), operator(/), operator(*)
39use mom_time_manager, only : set_time, get_time, get_date, time_to_real
42use mom_variables, only : surface
61
62implicit none ; private
63
64#include <MOM_memory.h>
65
66public set_forcing
69
70!> Structure containing pointers to the forcing fields that may be used to drive MOM.
71!! All fluxes are positive into the ocean.
72type, public :: surface_forcing_cs ; private
73
74 logical :: use_temperature !< if true, temp & salinity used as state variables
75 logical :: restorebuoy !< if true, use restoring surface buoyancy forcing
76 logical :: adiabatic !< if true, no diapycnal mass fluxes or surface buoyancy forcing
77 logical :: nonbous !< If true, this run is fully non-Boussinesq
78 logical :: variable_winds !< if true, wind stresses vary with time
79 logical :: variable_buoyforce !< if true, buoyancy forcing varies with time.
80 real :: south_lat !< southern latitude of the domain [degrees_N] or [km] or [m]
81 real :: len_lat !< domain length in latitude [degrees_N] or [km] or [m]
82
83 real :: rho0 !< Boussinesq reference density [R ~> kg m-3]
84 real :: g_earth !< gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
85 real :: flux_const = 0.0 !< piston velocity for surface restoring [Z T-1 ~> m s-1]
86 real :: flux_const_t = 0.0 !< piston velocity for surface temperature restoring [Z T-1 ~> m s-1]
87 real :: flux_const_s = 0.0 !< piston velocity for surface salinity restoring [Z T-1 ~> m s-1]
88 real :: rho_restore !< The density that is used to convert piston velocities into salt
89 !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
90 real :: latent_heat_fusion !< latent heat of fusion times [Q ~> J kg-1]
91 real :: latent_heat_vapor !< latent heat of vaporization [Q ~> J kg-1]
92 real :: tau_x0 !< Constant zonal wind stress used in the WIND_CONFIG="const"
93 !! forcing [R L Z T-2 ~> Pa]
94 real :: tau_y0 !< Constant meridional wind stress used in the WIND_CONFIG="const"
95 !! forcing [R L Z T-2 ~> Pa]
96 real :: taux_mag !< Peak magnitude of the zonal wind stress for several analytic
97 !! profiles [R L Z T-2 ~> Pa]
98
99 real :: gust_const !< constant unresolved background gustiness for ustar [R Z2 T-2 ~> Pa]
100 logical :: read_gust_2d !< if true, use 2-dimensional gustiness supplied from a file
101 real, pointer :: gust(:,:) => null() !< spatially varying unresolved background gustiness [R L Z T-2 ~> Pa]
102 !! gust is used when read_gust_2d is true.
103
104 real, pointer :: t_restore(:,:) => null() !< temperature to damp (restore) the SST to [C ~> degC]
105 real, pointer :: s_restore(:,:) => null() !< salinity to damp (restore) the SSS [S ~> ppt]
106 real, pointer :: dens_restore(:,:) => null() !< density to damp (restore) surface density [R ~> kg m-3]
107
108 ! if WIND_CONFIG=='gyres' then use the following as = A, B, C and n respectively for
109 ! taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L)
110 real :: gyres_taux_const !< A constant wind stress [R L Z T-2 ~> Pa].
111 real :: gyres_taux_sin_amp !< The amplitude of cosine wind stress gyres [R L Z T-2 ~> Pa], if WIND_CONFIG=='gyres'
112 real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [R L Z T-2 ~> Pa], if WIND_CONFIG=='gyres'
113 real :: gyres_taux_n_pis !< The number of sine lobes in the basin if WIND_CONFIG=='gyres' [nondim]
114 integer :: answer_date !< This 8-digit integer gives the approximate date with which the order
115 !! of arithmetic and expressions were added to the code.
116 !! Dates before 20190101 use original answers.
117 !! Dates after 20190101 use a form of the gyre wind stresses that are
118 !! rotationally invariant and more likely to be the same between compilers.
119 logical :: ustar_gustless_bug !< If true, include a bug in the time-averaging of the
120 !! gustless wind friction velocity.
121 logical :: use_marbl_tracers !< If true, allocate memory for forcing needed by MARBL
122 ! if WIND_CONFIG=='scurves' then use the following to define a piecewise scurve profile
123 real :: scurves_ydata(20) = 90. !< Latitudes of scurve nodes [degreesN]
124 real :: scurves_taux(20) = 0. !< Zonal wind stress values at scurve nodes [R L Z T-2 ~> Pa]
125
126 real :: t_north !< Target temperatures at north used in buoyancy_forcing_linear [C ~> degC]
127 real :: t_south !< Target temperatures at south used in buoyancy_forcing_linear [C ~> degC]
128 real :: s_north !< Target salinity at north used in buoyancy_forcing_linear [S ~> ppt]
129 real :: s_south !< Target salinity at south used in buoyancy_forcing_linear [S ~> ppt]
130
131 logical :: first_call_set_forcing = .true. !< True until after the first call to set_forcing
132 logical :: archaic_omip_file = .true. !< If true use the variable names and data fields from
133 !! a very old version of the OMIP forcing
134 logical :: dataoverrideisinitialized = .false. !< If true, data override has been initialized
135
136 real :: wind_scale !< value by which wind-stresses are scaled [nondim]
137 real :: constantheatforcing !< value used for sensible heat flux when buoy_config="const" [Q R Z T-1 ~> W m-2]
138
139 character(len=8) :: wind_stagger !< A character indicating how the wind stress components
140 !! are staggered in WIND_FILE. Valid values are A or C for now.
141 type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null() !< A pointer to the structure
142 !! that is used to orchestrate the calling of tracer packages
143!#CTRL# type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL()
144 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure
145
146 type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
147
148 character(len=200) :: inputdir !< directory where NetCDF input files are.
149 character(len=200) :: wind_config !< indicator for wind forcing type (2gyre, USER, FILE..)
150 character(len=200) :: wind_file !< if wind_config is "file", file to use
151 character(len=200) :: buoy_config !< indicator for buoyancy forcing type
152
153 character(len=200) :: longwave_file = '' !< The file from which the longwave heat flux is read
154 character(len=200) :: shortwave_file = '' !< The file from which the shortwave heat flux is read
155 character(len=200) :: evaporation_file = '' !< The file from which the evaporation is read
156 character(len=200) :: sensibleheat_file = '' !< The file from which the sensible heat flux is read
157 character(len=200) :: latentheat_file = '' !< The file from which the latent heat flux is read
158
159 character(len=200) :: rain_file = '' !< The file from which the rainfall is read
160 character(len=200) :: snow_file = '' !< The file from which the snowfall is read
161 character(len=200) :: runoff_file = '' !< The file from which the runoff is read
162
163 character(len=200) :: longwaveup_file = '' !< The file from which the upward longwave heat flux is read
164 character(len=200) :: shortwaveup_file = '' !< The file from which the upward shortwave heat flux is read
165
166 character(len=200) :: sstrestore_file = '' !< The file from which to read the sea surface
167 !! temperature to restore toward
168 character(len=200) :: salinityrestore_file = '' !< The file from which to read the sea surface
169 !! salinity to restore toward
170
171 character(len=80) :: stress_x_var = '' !< X-wind stress variable name in the input file
172 character(len=80) :: stress_y_var = '' !< Y-wind stress variable name in the input file
173 character(len=80) :: ustar_var = '' !< ustar variable name in the input file
174 character(len=80) :: lw_var = '' !< longwave heat flux variable name in the input file
175 character(len=80) :: sw_var = '' !< shortwave heat flux variable name in the input file
176 character(len=80) :: latent_var = '' !< latent heat flux variable name in the input file
177 character(len=80) :: sens_var = '' !< sensible heat flux variable name in the input file
178 character(len=80) :: evap_var = '' !< evaporation variable name in the input file
179 character(len=80) :: rain_var = '' !< rainfall variable name in the input file
180 character(len=80) :: snow_var = '' !< snowfall variable name in the input file
181 character(len=80) :: lrunoff_var = '' !< liquid runoff variable name in the input file
182 character(len=80) :: frunoff_var = '' !< frozen runoff variable name in the input file
183 character(len=80) :: sst_restore_var = '' !< target sea surface temperature variable name in the input file
184 character(len=80) :: sss_restore_var = '' !< target sea surface salinity variable name in the input file
185
186 ! These variables relate model times to time levels in the various forcing files.
187 integer :: wind_days_per_rec = 0 !< If positive the number of days of wind stress per forcing file time level,
188 !! or if negative the number of time levels per day. If 31 change forcing
189 !! monthly, or if 0 the model will guess the right value based on the file size.
190 integer :: sw_days_per_rec = 0 !< If positive the number of days shortwave heat flux per forcing file time level,
191 !! or if negative the number of time levels per day. If 31 change forcing
192 !! monthly, or if 0 the model will guess the right value based on the file size.
193 integer :: lw_days_per_rec = 0 !< If positive the number of days longwave heat flux per forcing file time level,
194 !! or if negative the number of time levels per day. If 31 change forcing
195 !! monthly, or if 0 the model will guess the right value based on the file size.
196 integer :: latent_days_per_rec = 0 !< If positive the number of days latent heat flux per forcing file time level,
197 !! or if negative the number of time levels per day. If 31 change forcing
198 !! monthly, or if 0 the model will guess the right value based on the file size.
199 integer :: sens_days_per_rec = 0 !< If positive the number of days sensible heat flux per forcing file time level,
200 !! or if negative the number of time levels per day. If 31 change forcing
201 !! monthly, or if 0 the model will guess the right value based on the file size.
202 integer :: evap_days_per_rec = 0 !< If positive the number of days evaporation per forcing file time level,
203 !! or if negative the number of time levels per day. If 31 change forcing
204 !! monthly, or if 0 the model will guess the right value based on the file size.
205 integer :: precip_days_per_rec = 0 !< If positive the number of days precipitation per forcing file time level,
206 !! or if negative the number of time levels per day. If 31 change forcing
207 !! monthly, or if 0 the model will guess the right value based on the file size.
208 integer :: runoff_days_per_rec = 0 !< If positive the number of days runoff per forcing file time level,
209 !! or if negative the number of time levels per day. If 31 change forcing
210 !! monthly, or if 0 the model will guess the right value based on the file size.
211 integer :: sst_days_per_rec = 0 !< If positive the number of days target SST per forcing file time level,
212 !! or if negative the number of time levels per day. If 31 change forcing
213 !! monthly, or if 0 the model will guess the right value based on the file size.
214 integer :: sss_days_per_rec = 0 !< If positive the number of days target SSS per forcing file time level,
215 !! or if negative the number of time levels per day. If 31 change forcing
216 !! monthly, or if 0 the model will guess the right value based on the file size.
217
218 ! These variables give the number of time levels in the various forcing files.
219 integer :: wind_nlev = -1 !< The number of time levels in the file of wind stress
220 integer :: sw_nlev = -1 !< The number of time levels in the file of shortwave heat flux
221 integer :: lw_nlev = -1 !< The number of time levels in the file of longwave heat flux
222 integer :: latent_nlev = -1 !< The number of time levels in the file of latent heat flux
223 integer :: sens_nlev = -1 !< The number of time levels in the file of sensible heat flux
224 integer :: evap_nlev = -1 !< The number of time levels in the file of evaporation
225 integer :: precip_nlev = -1 !< The number of time levels in the file of precipitation
226 integer :: runoff_nlev = -1 !< The number of time levels in the file of runoff
227 integer :: sst_nlev = -1 !< The number of time levels in the file of target SST
228 integer :: sss_nlev = -1 !< The number of time levels in the file of target SSS
229
230 ! These variables give the last time level read for the various forcing files.
231 integer :: wind_last_lev = -1 !< The last time level read of wind stress
232 integer :: sw_last_lev = -1 !< The last time level read of shortwave heat flux
233 integer :: lw_last_lev = -1 !< The last time level read of longwave heat flux
234 integer :: latent_last_lev = -1 !< The last time level read of latent heat flux
235 integer :: sens_last_lev = -1 !< The last time level read of sensible heat flux
236 integer :: evap_last_lev = -1 !< The last time level read of evaporation
237 integer :: precip_last_lev = -1 !< The last time level read of precipitation
238 integer :: runoff_last_lev = -1 !< The last time level read of runoff
239 integer :: sst_last_lev = -1 !< The last time level read of target SST
240 integer :: sss_last_lev = -1 !< The last time level read of target SSS
241
242 type(forcing_diags), public :: handles !< A structure with diagnostics handles
243
244 !>@{ Control structures for named forcing packages
245 type(user_revise_forcing_cs), pointer :: urf_cs => null()
246 type(user_surface_forcing_cs), pointer :: user_forcing_csp => null()
247 type(bfb_surface_forcing_cs), pointer :: bfb_forcing_csp => null()
248 type(dumbbell_surface_forcing_cs), pointer :: dumbbell_forcing_csp => null()
249 type(meso_surface_forcing_cs), pointer :: meso_forcing_csp => null()
250 type(idealized_hurricane_cs), pointer :: idealized_hurricane_csp => null()
251 type(scm_cvmix_tests_cs), pointer :: scm_cvmix_tests_csp => null()
252 type(marbl_forcing_cs), pointer :: marbl_forcing_csp => null()
253 !>@}
254
255end type surface_forcing_cs
256
257integer :: id_clock_forcing !< A CPU time clock
258
259contains
260
261!> Calls subroutines in this file to get surface forcing fields.
262!!
263!! It also allocates and initializes the fields in the forcing and mech_forcing types
264!! the first time it is called.
265subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US, CS)
266 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
267 !! describe the surface state of the ocean.
268 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
269 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
270 type(time_type), intent(in) :: day_start !< The start time of the fluxes
271 type(time_type), intent(in) :: day_interval !< Length of time over which these fluxes applied
272 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
273 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
274 type(surface_forcing_cs), pointer :: cs !< pointer to control structure returned by
275 !! a previous surface_forcing_init call
276 ! Local variables
277 real :: dt ! length of time over which fluxes applied [T ~> s]
278 type(time_type) :: day_center ! central time of the fluxes.
279 integer :: isd, ied, jsd, jed
280 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
281
282 call cpu_clock_begin(id_clock_forcing)
283 call calltree_enter("set_forcing, MOM_surface_forcing.F90")
284
285 day_center = day_start + day_interval/2
286 dt = time_to_real(day_interval, scale=us%s_to_T)
287
288 if (cs%first_call_set_forcing) then
289 ! Allocate memory for the mechanical and thermodynamic forcing fields.
290 call allocate_mech_forcing(g, forces, stress=.true., ustar=.not.cs%nonBous, press=.true., tau_mag=cs%nonBous)
291
292 call allocate_forcing_type(g, fluxes, ustar=.not.cs%nonBous, marbl=cs%use_marbl_tracers, tau_mag=cs%nonBous, &
293 fix_accum_bug=.not.cs%ustar_gustless_bug)
294 if (trim(cs%buoy_config) /= "NONE") then
295 if ( cs%use_temperature ) then
296 call allocate_forcing_type(g, fluxes, water=.true., heat=.true., press=.true.)
297 if (cs%restorebuoy) then
298 call safe_alloc_ptr(cs%T_Restore,isd, ied, jsd, jed)
299 call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
300 call safe_alloc_ptr(cs%S_Restore, isd, ied, jsd, jed)
301 endif
302 else ! CS%use_temperature false.
303 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
304
305 if (cs%restorebuoy) call safe_alloc_ptr(cs%Dens_Restore, isd, ied, jsd, jed)
306 endif ! endif for CS%use_temperature
307 endif
308 endif
309
310 ! calls to various wind options
311 if (cs%variable_winds .or. cs%first_call_set_forcing) then
312
313 if (trim(cs%wind_config) == "file") then
314 call wind_forcing_from_file(sfc_state, forces, day_center, g, us, cs)
315 elseif (trim(cs%wind_config) == "data_override") then
316 call wind_forcing_by_data_override(sfc_state, forces, day_center, g, us, cs)
317 elseif (trim(cs%wind_config) == "2gyre") then
318 call wind_forcing_2gyre(sfc_state, forces, day_center, g, us, cs)
319 elseif (trim(cs%wind_config) == "1gyre") then
320 call wind_forcing_1gyre(sfc_state, forces, day_center, g, us, cs)
321 elseif (trim(cs%wind_config) == "gyres") then
322 call wind_forcing_gyres(sfc_state, forces, day_center, g, us, cs)
323 elseif (trim(cs%wind_config) == "zero") then
324 call wind_forcing_const(sfc_state, forces, 0., 0., day_center, g, us, cs)
325 elseif (trim(cs%wind_config) == "const") then
326 call wind_forcing_const(sfc_state, forces, cs%tau_x0, cs%tau_y0, day_center, g, us, cs)
327 elseif (trim(cs%wind_config) == "Neverworld" .or. trim(cs%wind_config) == "Neverland") then
328 call neverworld_wind_forcing(sfc_state, forces, day_center, g, us, cs)
329 elseif (trim(cs%wind_config) == "scurves") then
330 call scurve_wind_forcing(sfc_state, forces, day_center, g, us, cs)
331 elseif (trim(cs%wind_config) == "ideal_hurr") then
332 call idealized_hurricane_wind_forcing(sfc_state, forces, day_center, g, us, cs%idealized_hurricane_CSp)
333 elseif (trim(cs%wind_config) == "SCM_ideal_hurr") then
334 call mom_error(fatal, "MOM_surface_forcing (set_forcing): "//&
335 'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option.')
336 elseif (trim(cs%wind_config) == "SCM_CVmix_tests") then
337 call scm_cvmix_tests_wind_forcing(sfc_state, forces, day_center, g, us, cs%SCM_CVmix_tests_CSp)
338 elseif (trim(cs%wind_config) == "USER") then
339 call user_wind_forcing(sfc_state, forces, day_center, g, us, cs%user_forcing_CSp)
340 elseif (cs%variable_winds .and. .not.cs%first_call_set_forcing) then
341 call mom_error(fatal, &
342 "MOM_surface_forcing: Variable winds defined with no wind config")
343 else
344 call mom_error(fatal, &
345 "MOM_surface_forcing:Unrecognized wind config "//trim(cs%wind_config))
346 endif
347 endif
348
349 ! calls to various buoyancy forcing options
350 if (cs%restorebuoy .and. .not.cs%variable_buoyforce) then
351 call mom_error(fatal, "With RESTOREBUOY = True, VARIABLE_BUOYFORCE = True should be used. "//&
352 "Otherwise, this can lead to diverging solutions when a simulation "//&
353 "is continued using a restart file.")
354 endif
355
356 if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
357 (.not.cs%adiabatic)) then
358 if (trim(cs%buoy_config) == "file") then
359 call buoyancy_forcing_from_files(sfc_state, fluxes, day_center, dt, g, us, cs)
360 elseif (trim(cs%buoy_config) == "data_override") then
361 call buoyancy_forcing_from_data_override(sfc_state, fluxes, day_center, dt, g, us, cs)
362 elseif (trim(cs%buoy_config) == "zero") then
363 call buoyancy_forcing_zero(sfc_state, fluxes, day_center, dt, g, cs)
364 elseif (trim(cs%buoy_config) == "const") then
365 call buoyancy_forcing_const(sfc_state, fluxes, day_center, dt, g, us, cs)
366 elseif (trim(cs%buoy_config) == "linear") then
367 call buoyancy_forcing_linear(sfc_state, fluxes, day_center, dt, g, us, cs)
368 elseif (trim(cs%buoy_config) == "MESO") then
369 call meso_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%MESO_forcing_CSp)
370 elseif (trim(cs%buoy_config) == "SCM_CVmix_tests") then
371 call scm_cvmix_tests_buoyancy_forcing(sfc_state, fluxes, day_center, g, us, cs%SCM_CVmix_tests_CSp)
372 elseif (trim(cs%buoy_config) == "USER") then
373 call user_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%user_forcing_CSp)
374 elseif (trim(cs%buoy_config) == "BFB") then
375 call bfb_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%BFB_forcing_CSp)
376 elseif (trim(cs%buoy_config) == "dumbbell") then
377 call dumbbell_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%dumbbell_forcing_CSp)
378 elseif (trim(cs%buoy_config) == "NONE") then
379 call mom_mesg("MOM_surface_forcing: buoyancy forcing has been set to omitted.")
380 elseif (cs%variable_buoyforce .and. .not.cs%first_call_set_forcing) then
381 call mom_error(fatal, &
382 "MOM_surface_forcing: Variable buoy defined with no buoy config.")
383 else
384 call mom_error(fatal, &
385 "MOM_surface_forcing: Unrecognized buoy config "//trim(cs%buoy_config))
386 endif
387 endif
388
389 if (cs%use_marbl_tracers) then
390 call marbl_forcing_from_data_override(fluxes, day_center, g, us, cs)
391 endif
392
393 if (associated(cs%tracer_flow_CSp)) then
394 call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, g, us, cs%Rho0, &
395 cs%tracer_flow_CSp)
396 endif
397
398 ! Allow for user-written code to alter the fluxes after all the above
399 call user_alter_forcing(sfc_state, fluxes, day_center, g, cs%urf_CS)
400
401 ! Fields that exist in both the forcing and mech_forcing types must be copied.
402 if (cs%variable_winds .or. cs%first_call_set_forcing) then
403 call copy_common_forcing_fields(forces, fluxes, g)
404 call set_derived_forcing_fields(forces, fluxes, g, us, cs%Rho0)
405 endif
406
407 if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
408 (.not.cs%adiabatic)) then
409 call set_net_mass_forcing(fluxes, forces, g, us)
410 endif
411
412 cs%first_call_set_forcing = .false.
413
414 call cpu_clock_end(id_clock_forcing)
415 call calltree_leave("set_forcing")
416
417end subroutine set_forcing
418
419!> Sets the surface wind stresses to constant values
420subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS)
421 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
422 !! describe the surface state of the ocean.
423 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
424 real, intent(in) :: tau_x0 !< The zonal wind stress [R Z L T-2 ~> Pa]
425 real, intent(in) :: tau_y0 !< The meridional wind stress [R Z L T-2 ~> Pa]
426 type(time_type), intent(in) :: day !< The time of the fluxes
427 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
428 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
429 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
430 !! a previous surface_forcing_init call
431 ! Local variables
432 real :: mag_tau ! Magnitude of the wind stress [R Z2 T-2 ~> Pa]
433 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
434
435 call calltree_enter("wind_forcing_const, MOM_surface_forcing.F90")
436 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
437 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
438
439 mag_tau = us%L_to_Z * sqrt( tau_x0**2 + tau_y0**2)
440
441 ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa].
442 do j=js,je ; do i=is-1,ieq
443 forces%taux(i,j) = tau_x0
444 enddo ; enddo
445
446 do j=js-1,jeq ; do i=is,ie
447 forces%tauy(i,j) = tau_y0
448 enddo ; enddo
449
450 if (cs%read_gust_2d) then
451 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
452 forces%ustar(i,j) = sqrt( ( mag_tau + cs%gust(i,j) ) / cs%Rho0 )
453 enddo ; enddo ; endif
454 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
455 forces%tau_mag(i,j) = mag_tau + cs%gust(i,j)
456 enddo ; enddo ; endif
457 else
458 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
459 forces%ustar(i,j) = sqrt( ( mag_tau + cs%gust_const ) / cs%Rho0 )
460 enddo ; enddo ; endif
461 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
462 forces%tau_mag(i,j) = mag_tau + cs%gust_const
463 enddo ; enddo ; endif
464 endif
465
466 call calltree_leave("wind_forcing_const")
467end subroutine wind_forcing_const
468
469
470!> Sets the surface wind stresses to set up two idealized gyres.
471subroutine wind_forcing_2gyre(sfc_state, forces, day, G, US, CS)
472 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
473 !! describe the surface state of the ocean.
474 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
475 type(time_type), intent(in) :: day !< The time of the fluxes
476 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
477 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
478 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
479 !! a previous surface_forcing_init call
480 ! Local variables
481 real :: PI ! A common irrational number, 3.1415926535... [nondim]
482 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
483
484 call calltree_enter("wind_forcing_2gyre, MOM_surface_forcing.F90")
485 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
486 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
487
488 pi = 4.0*atan(1.0)
489
490 ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa].
491 do j=js,je ; do i=is-1,ieq
492 forces%taux(i,j) = cs%taux_mag * (1.0 - cos(2.0*pi*(g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat))
493 enddo ; enddo
494
495 do j=js-1,jeq ; do i=is,ie
496 forces%tauy(i,j) = 0.0
497 enddo ; enddo
498
499 if (associated(forces%ustar)) call stresses_to_ustar(forces, g, us, cs)
500
501 call calltree_leave("wind_forcing_2gyre")
502end subroutine wind_forcing_2gyre
503
504
505!> Sets the surface wind stresses to set up a single idealized gyre.
506subroutine wind_forcing_1gyre(sfc_state, forces, day, G, US, CS)
507 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
508 !! describe the surface state of the ocean.
509 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
510 type(time_type), intent(in) :: day !< The time of the fluxes
511 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
512 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
513 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
514 !! a previous surface_forcing_init call
515 ! Local variables
516 real :: PI ! A common irrational number, 3.1415926535... [nondim]
517 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
518
519 call calltree_enter("wind_forcing_1gyre, MOM_surface_forcing.F90")
520 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
521 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
522
523 pi = 4.0*atan(1.0)
524
525 ! Set the steady surface wind stresses, in units of [R Z L T-2 ~> Pa].
526 do j=js,je ; do i=is-1,ieq
527 forces%taux(i,j) = cs%taux_mag * cos(pi*(g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat)
528 enddo ; enddo
529
530 do j=js-1,jeq ; do i=is,ie
531 forces%tauy(i,j) = 0.0
532 enddo ; enddo
533
534 if (associated(forces%ustar)) call stresses_to_ustar(forces, g, us, cs)
535
536 call calltree_leave("wind_forcing_1gyre")
537end subroutine wind_forcing_1gyre
538
539!> Sets the surface wind stresses to set up idealized gyres.
540subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
541 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
542 !! describe the surface state of the ocean.
543 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
544 type(time_type), intent(in) :: day !< The time of the fluxes
545 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
546 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
547 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
548 !! a previous surface_forcing_init call
549 ! Local variables
550 real :: PI ! A common irrational number, 3.1415926535... [nondim]
551 real :: y ! The latitude relative to the south normalized by the domain extent [nondim]
552 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
553
554 call calltree_enter("wind_forcing_gyres, MOM_surface_forcing.F90")
555 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
556 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
557
558 pi = 4.0*atan(1.0)
559
560 ! steady surface wind stresses [R L Z T-2 ~> Pa]
561 do j=js-1,je+1 ; do i=is-1,ieq
562 y = (g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat
563 forces%taux(i,j) = cs%gyres_taux_const + &
564 ( cs%gyres_taux_sin_amp*sin(cs%gyres_taux_n_pis*pi*y) &
565 + cs%gyres_taux_cos_amp*cos(cs%gyres_taux_n_pis*pi*y) )
566 enddo ; enddo
567
568 do j=js-1,jeq ; do i=is-1,ie+1
569 forces%tauy(i,j) = 0.0
570 enddo ; enddo
571
572 ! set the friction velocity
573 if (cs%answer_date < 20190101) then
574 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
575 forces%tau_mag(i,j) = cs%gust_const + us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
576 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
577 enddo ; enddo ; endif
578 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
579 forces%ustar(i,j) = sqrt( (cs%gust_const/cs%Rho0) + &
580 us%L_to_Z * sqrt(0.5*((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2) + &
581 (forces%taux(i-1,j)**2) + (forces%taux(i,j)**2)))/cs%Rho0 )
582 enddo ; enddo ; endif
583 else
584 call stresses_to_ustar(forces, g, us, cs)
585 endif
586
587 call calltree_leave("wind_forcing_gyres")
588end subroutine wind_forcing_gyres
589
590!> Sets the surface wind stresses, forces%taux and forces%tauy for the
591!! Neverworld forcing configuration.
592subroutine neverworld_wind_forcing(sfc_state, forces, day, G, US, CS)
593 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
594 !! describe the surface state of the ocean.
595 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
596 type(time_type), intent(in) :: day !< Time used for determining the fluxes.
597 type(ocean_grid_type), intent(inout) :: G !< Grid structure.
598 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
599 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
600 !! a previous surface_forcing_init call
601 ! Local variables
602 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
603 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
604 real :: PI ! A common irrational number, 3.1415926535... [nondim]
605 real :: y ! The latitude relative to the south normalized by the domain extent [nondim]
606 real :: tau_max ! The magnitude of the wind stress [R Z L T-2 ~> Pa]
607 real :: off ! An offset in the relative latitude [nondim]
608
609 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
610 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
611 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
612 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
613
614 ! Allocate the forcing arrays, if necessary.
615 call allocate_mech_forcing(g, forces, stress=.true.)
616
617 ! Set the surface wind stresses, in units of [R Z L T-2 ~> Pa]. A positive taux
618 ! accelerates the ocean to the (pseudo-)east.
619
620 ! The i-loop extends to is-1 so that taux can be used later in the
621 ! calculation of ustar - otherwise the lower bound would be Isq.
622 pi = 4.0*atan(1.0)
623
624 forces%taux(:,:) = 0.0
625 tau_max = cs%taux_mag
626 off = 0.02
627 do j=js,je ; do i=is-1,ieq
628 y = (g%geoLatT(i,j)-g%south_lat)/g%len_lat
629
630 if (y <= 0.29) then
631 forces%taux(i,j) = forces%taux(i,j) + tau_max * ( (1/0.29)*y - ( 1/(2*pi) )*sin( (2*pi*y) / 0.29 ) )
632 endif
633 if ((y > 0.29) .and. (y <= (0.8-off))) then
634 forces%taux(i,j) = forces%taux(i,j) + tau_max *(0.35+0.65*cos(pi*(y-0.29)/(0.51-off)) )
635 endif
636 if ((y > (0.8-off)) .and. (y <= (1-off))) then
637 forces%taux(i,j) = forces%taux(i,j) + tau_max *( 1.5*( (y-1+off) - (0.1/pi)*sin(10.0*pi*(y-0.8+off)) ) )
638 endif
639 forces%taux(i,j) = g%mask2dCu(i,j) * forces%taux(i,j)
640 enddo ; enddo
641
642 do j=js-1,jeq ; do i=is,ie
643 forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0
644 enddo ; enddo
645
646 ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1]. ustar is always positive.
647 if (associated(forces%ustar)) call stresses_to_ustar(forces, g, us, cs)
648
649end subroutine neverworld_wind_forcing
650
651!> Sets the zonal wind stresses to a piecewise series of s-curves.
652subroutine scurve_wind_forcing(sfc_state, forces, day, G, US, CS)
653 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
654 !! describe the surface state of the ocean.
655 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
656 type(time_type), intent(in) :: day !< Time used for determining the fluxes.
657 type(ocean_grid_type), intent(inout) :: G !< Grid structure.
658 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
659 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
660 !! a previous surface_forcing_init call
661 ! Local variables
662 integer :: i, j, kseg
663 real :: y_curve ! The latitude relative to the southern end of a curve segment [degreesN]
664 real :: L_curve ! The latitudinal extent of a curve segment [degreesN]
665! real :: ydata(7) = (/ -70., -45., -15., 0., 15., 45., 70. /)
666! real :: taudt(7) = (/ 0., 0.2, -0.1, -0.02, -0.1, 0.1, 0. /)
667
668 ! Allocate the forcing arrays, if necessary.
669 call allocate_mech_forcing(g, forces, stress=.true.)
670
671 kseg = 1
672 do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
673 ! Find segment k s.t. ydata(k)<= G%geoLatCu(I,j) < ydata(k+1)
674 do while (g%geoLatCu(i,j) >= cs%scurves_ydata(kseg+1) .and. kseg<6) ! Should this be kseg<19?
675 kseg = kseg+1
676 enddo
677 do while (g%geoLatCu(i,j) < cs%scurves_ydata(kseg) .and. kseg>1)
678 kseg = kseg-1
679 enddo
680
681 y_curve = g%geoLatCu(i,j) - cs%scurves_ydata(kseg)
682 l_curve = cs%scurves_ydata(kseg+1) - cs%scurves_ydata(kseg)
683 forces%taux(i,j) = cs%scurves_taux(kseg) + &
684 (cs%scurves_taux(kseg+1) - cs%scurves_taux(kseg)) * scurve(y_curve, l_curve)
685 forces%taux(i,j) = g%mask2dCu(i,j) * forces%taux(i,j)
686 enddo ; enddo
687
688 do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
689 forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0
690 enddo ; enddo
691
692 ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1]. ustar is always positive.
693 if (associated(forces%ustar)) call stresses_to_ustar(forces, g, us, cs)
694
695end subroutine scurve_wind_forcing
696
697!> Returns the value of a cosine-bell function evaluated at x/L
698real function scurve(x,L)
699 real , intent(in) :: x !< non-dimensional position [nondim]
700 real , intent(in) :: l !< non-dimensional width [nondim]
701 real :: s ! The evaluated function value [nondim]
702
703 s = x/l
704 scurve = (3. - 2.*s) * (s*s)
705end function scurve
706
707! Sets the surface wind stresses from input files.
708subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
709 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
710 !! describe the surface state of the ocean.
711 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
712 type(time_type), intent(in) :: day !< The time of the fluxes
713 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
714 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
715 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
716 !! a previous surface_forcing_init call
717 ! Local variables
718 character(len=200) :: filename ! The name of the input file.
719 real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R L Z T-2 ~> Pa]
720 real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R L Z T-2 ~> Pa]
721 real :: ustar_loc(SZI_(G),SZJ_(G)) ! The local value of ustar [Z T-1 ~> m s-1]
722 real :: tau_mag ! The magnitude of the wind stress including any contributions from
723 ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
724 integer :: time_lev ! The time level that is used for a field.
725 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
726 logical :: read_Ustar
727
728 call calltree_enter("wind_forcing_from_file, MOM_surface_forcing.F90")
729 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
730 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
731
732 time_lev = get_file_time_level(day, cs%wind_nlev, cs%wind_days_per_rec)
733
734 if (time_lev /= cs%wind_last_lev) then
735 filename = trim(cs%wind_file)
736 read_ustar = (len_trim(cs%ustar_var) > 0)
737! if (is_root_pe()) &
738! write(*,'("Wind_forcing Reading time level ",I," last was ",I,".")')&
739! time_lev-1,CS%wind_last_lev-1
740 select case ( uppercase(cs%wind_stagger(1:1)) )
741 case ("A")
742 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
743 call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
744 temp_x(:,:), temp_y(:,:), g%Domain, stagger=agrid, &
745 timelevel=time_lev, scale=us%Pa_to_RLZ_T2)
746
747 call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
748 do j=js,je ; do i=is-1,ieq
749 forces%taux(i,j) = 0.5 * cs%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
750 enddo ; enddo
751 do j=js-1,jeq ; do i=is,ie
752 forces%tauy(i,j) = 0.5 * cs%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
753 enddo ; enddo
754
755 if (.not.read_ustar) then
756 if (cs%read_gust_2d) then
757 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
758 forces%tau_mag(i,j) = cs%gust(i,j) + us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
759 enddo ; enddo ; endif
760 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
761 tau_mag = cs%gust(i,j) + us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
762 forces%ustar(i,j) = sqrt(tau_mag / cs%Rho0)
763 enddo ; enddo ; endif
764 else
765 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
766 forces%tau_mag(i,j) = cs%gust_const + us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))
767 enddo ; enddo ; endif
768 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
769 forces%ustar(i,j) = sqrt( cs%gust_const/cs%Rho0 + &
770 us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) / cs%Rho0 )
771 enddo ; enddo ; endif
772 endif
773 endif
774 case ("C")
775 if (g%symmetric) then
776 if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
777 " wind_forcing_from_file with C-grid input and symmetric memory "//&
778 " called with a non-associated auxiliary domain in the grid type.")
779 ! Read the data as though symmetric memory were not being used, and
780 ! then translate it appropriately.
781 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
782 call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
783 temp_x(:,:), temp_y(:,:), &
784 g%Domain_aux, stagger=cgrid_ne, timelevel=time_lev, &
785 scale=us%Pa_to_RLZ_T2)
786 do j=js,je ; do i=is,ie
787 forces%taux(i,j) = cs%wind_scale * temp_x(i,j)
788 forces%tauy(i,j) = cs%wind_scale * temp_y(i,j)
789 enddo ; enddo
790 call fill_symmetric_edges(forces%taux, forces%tauy, g%Domain, stagger=cgrid_ne)
791 else
792 call mom_read_vector(filename, cs%stress_x_var, cs%stress_y_var, &
793 forces%taux(:,:), forces%tauy(:,:), &
794 g%Domain, stagger=cgrid_ne, timelevel=time_lev, &
795 scale=us%Pa_to_RLZ_T2)
796
797 if (cs%wind_scale /= 1.0) then
798 do j=js,je ; do i=isq,ieq
799 forces%taux(i,j) = cs%wind_scale * forces%taux(i,j)
800 enddo ; enddo
801 do j=jsq,jeq ; do i=is,ie
802 forces%tauy(i,j) = cs%wind_scale * forces%tauy(i,j)
803 enddo ; enddo
804 endif
805 endif
806
807 call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
808 if (.not.read_ustar) then
809 if (cs%read_gust_2d) then
810 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
811 forces%tau_mag(i,j) = cs%gust(i,j) + &
812 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
813 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
814 enddo ; enddo ; endif
815 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
816 tau_mag = cs%gust(i,j) + &
817 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
818 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
819 forces%ustar(i,j) = sqrt( tau_mag / cs%Rho0 )
820 enddo ; enddo ; endif
821 else
822 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
823 forces%tau_mag(i,j) = cs%gust_const + &
824 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
825 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
826 enddo ; enddo ; endif
827 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
828 forces%ustar(i,j) = sqrt( cs%gust_const/cs%Rho0 + &
829 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
830 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))/cs%Rho0 )
831 enddo ; enddo ; endif
832 endif
833 endif
834 case default
835 call mom_error(fatal, "wind_forcing_from_file: Unrecognized stagger "//&
836 trim(cs%wind_stagger)//" is not 'A' or 'C'.")
837 end select
838
839 if (read_ustar) then
840 call mom_read_data(filename, cs%Ustar_var, ustar_loc(:,:), &
841 g%Domain, timelevel=time_lev, scale=us%m_to_Z*us%T_to_s)
842 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
843 forces%tau_mag(i,j) = cs%Rho0 * ustar_loc(i,j)**2
844 enddo ; enddo ; endif
845 if (associated(forces%ustar)) then ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
846 forces%ustar(i,j) = ustar_loc(i,j)
847 enddo ; enddo ; endif
848 endif
849
850 cs%wind_last_lev = time_lev
851
852 endif ! time_lev /= CS%wind_last_lev
853
854 call calltree_leave("wind_forcing_from_file")
855end subroutine wind_forcing_from_file
856
857
858! Sets the surface wind stresses via the data override facility.
859subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
860 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
861 !! describe the surface state of the ocean.
862 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
863 type(time_type), intent(in) :: day !< The time of the fluxes
864 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
865 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
866 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
867 !! a previous surface_forcing_init call
868 ! Local variables
869 real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R Z L T-2 ~> Pa].
870 real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R Z L T-2 ~> Pa].
871 real :: ustar_prev(SZI_(G),SZJ_(G)) ! The pre-override value of ustar [Z T-1 ~> m s-1]
872 real :: ustar_loc(SZI_(G),SZJ_(G)) ! The value of ustar, perhaps altered by data override [Z T-1 ~> m s-1]
873 real :: tau_mag ! The magnitude of the wind stress including any contributions from
874 ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
875 integer :: i, j
876
877 call calltree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90")
878
879 if (.not.cs%dataOverrideIsInitialized) then
880 call allocate_mech_forcing(g, forces, stress=.true., ustar=.not.cs%nonBous, press=.true., tau_mag=cs%nonBous)
881 call data_override_init(g%Domain)
882 cs%dataOverrideIsInitialized = .true.
883 endif
884
885 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
886 ! CS%wind_scale is ignored here because it is not set in this mode.
887 call data_override(g%Domain, 'taux', temp_x, day, scale=us%Pa_to_RLZ_T2)
888 call data_override(g%Domain, 'tauy', temp_y, day, scale=us%Pa_to_RLZ_T2)
889 call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
890 do j=g%jsc,g%jec ; do i=g%isc-1,g%IecB
891 forces%taux(i,j) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
892 enddo ; enddo
893 do j=g%jsc-1,g%JecB ; do i=g%isc,g%iec
894 forces%tauy(i,j) = 0.5 * (temp_y(i,j) + temp_y(i,j+1))
895 enddo ; enddo
896
897 if (cs%read_gust_2d) then
898 call data_override(g%Domain, 'gust', cs%gust, day, scale=us%Pa_to_RLZ_T2*us%L_to_Z)
899 if (associated(forces%tau_mag)) then ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
900 forces%tau_mag(i,j) = us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + cs%gust(i,j)
901 enddo ; enddo ; endif
902 do j=g%jsc,g%jec ; do i=g%isc,g%iec
903 tau_mag = us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + cs%gust(i,j)
904 ustar_loc(i,j) = sqrt( tau_mag / cs%Rho0 )
905 enddo ; enddo
906 else
907 if (associated(forces%tau_mag)) then
908 do j=g%jsc,g%jec ; do i=g%isc,g%iec
909 forces%tau_mag(i,j) = us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + cs%gust_const
910 ! ustar_loc(i,j) = sqrt( forces%tau_mag(i,j) / CS%Rho0 )
911 enddo ; enddo
912 endif
913 do j=g%jsc,g%jec ; do i=g%isc,g%iec
914 ustar_loc(i,j) = sqrt(us%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))/cs%Rho0 + &
915 cs%gust_const/cs%Rho0)
916 enddo ; enddo
917 endif
918
919 ! Give the data override the option to modify the newly calculated forces%ustar.
920 ustar_prev(:,:) = ustar_loc(:,:)
921 call data_override(g%Domain, 'ustar', ustar_loc, day, scale=us%m_to_Z*us%T_to_s)
922
923 ! Only reset values where data override of ustar has occurred
924 if (associated(forces%tau_mag)) then
925 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (ustar_prev(i,j) /= ustar_loc(i,j)) then
926 forces%tau_mag(i,j) = cs%Rho0 * ustar_loc(i,j)**2
927 endif ; enddo ; enddo
928 endif
929
930 if (associated(forces%ustar)) then ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
931 forces%ustar(i,j) = ustar_loc(i,j)
932 enddo ; enddo ; endif
933
934 call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
935
936 call calltree_leave("wind_forcing_by_data_override")
938
939!> Translate the wind stresses into the friction velocity, including effects of background gustiness.
940subroutine stresses_to_ustar(forces, G, US, CS)
941 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
942 type(ocean_grid_type), intent(in) :: G !< Grid structure.
943 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
944 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
945 !! a previous surface_forcing_init call
946 ! Local variables
947 real :: I_rho ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]
948 real :: tau_mag ! The magnitude of the wind stress including any contributions from
949 ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa]
950 integer :: i, j, is, ie, js, je
951
952 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
953
954 i_rho = 1.0 / cs%Rho0
955
956 if (cs%read_gust_2d) then
957 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
958 forces%tau_mag(i,j) = cs%gust(i,j) + &
959 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
960 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
961 enddo ; enddo ; endif
962 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
963 tau_mag = cs%gust(i,j) + &
964 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
965 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
966 forces%ustar(i,j) = sqrt( tau_mag * i_rho )
967 enddo ; enddo ; endif
968 else
969 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
970 forces%tau_mag(i,j) = cs%gust_const + &
971 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
972 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
973 enddo ; enddo ; endif
974 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
975 tau_mag = cs%gust_const + &
976 us%L_to_Z * sqrt(0.5*(((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)) + &
977 ((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2))))
978 forces%ustar(i,j) = sqrt( tau_mag * i_rho )
979 enddo ; enddo ; endif
980 endif
981
982end subroutine stresses_to_ustar
983
984!> Specifies zero surface buoyancy fluxes from input files.
985subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS)
986 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
987 !! describe the surface state of the ocean.
988 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
989 type(time_type), intent(in) :: day !< The time of the fluxes
990 real, intent(in) :: dt !< The amount of time over which
991 !! the fluxes apply [T ~> s]
992 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
993 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
994 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
995 !! a previous surface_forcing_init call
996 ! Local variables
997 real, dimension(SZI_(G),SZJ_(G)) :: &
998 temp ! A 2-d temporary work array in various units of [Q R Z T-1 ~> W m-2] or
999 ! [R Z T-1 ~> kg m-2 s-1]
1000!#CTRL# real, dimension(SZI_(G),SZJ_(G)) :: &
1001!#CTRL# SST_anom, & ! Instantaneous sea surface temperature anomalies from a
1002!#CTRL# ! target (observed) value [C ~> degC].
1003!#CTRL# SSS_anom, & ! Instantaneous sea surface salinity anomalies from a target
1004!#CTRL# ! (observed) value [S ~> ppt].
1005!#CTRL# SSS_mean ! A (mean?) salinity about which to normalize local salinity
1006!#CTRL# ! anomalies when calculating restorative precipitation anomalies [S ~> ppt].
1007
1008 real :: rhoXcp ! reference density times heat capacity [Q R C-1 ~> J m-3 degC-1]
1009
1010 logical :: fluxes_changed ! True if any of the fluxes might have been altered
1011 integer :: time_lev ! time level that for a field
1012
1013 integer :: i, j, is, ie, js, je
1014
1015 call calltree_enter("buoyancy_forcing_from_files, MOM_surface_forcing.F90")
1016
1017 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1018
1019 if (cs%use_temperature) rhoxcp = cs%rho_restore * fluxes%C_p
1020
1021 ! Read the buoyancy forcing file
1022 fluxes_changed = .false.
1023
1024 ! longwave
1025 time_lev = get_file_time_level(day, cs%LW_nlev, cs%LW_days_per_rec)
1026 if (time_lev /= cs%LW_last_lev) then
1027 call mom_read_data(cs%longwave_file, cs%LW_var, fluxes%lw(:,:), &
1028 g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1029 if (cs%archaic_OMIP_file) then
1030 call mom_read_data(cs%longwaveup_file, "lwup_sfc", temp(:,:), g%Domain, &
1031 timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1032 do j=js,je ; do i=is,ie ; fluxes%LW(i,j) = fluxes%LW(i,j) - temp(i,j) ; enddo ; enddo
1033 endif
1034 cs%LW_last_lev = time_lev ; fluxes_changed = .true.
1035 endif
1036
1037 ! evaporation
1038 if ( (cs%evap_nlev /= cs%LW_nlev) .or. (cs%evap_days_per_rec /= cs%LW_days_per_rec) ) &
1039 time_lev = get_file_time_level(day, cs%evap_nlev, cs%evap_days_per_rec)
1040 if (time_lev /= cs%evap_last_lev) then
1041 if (cs%archaic_OMIP_file) then
1042 call mom_read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
1043 g%Domain, timelevel=time_lev, scale=-us%kg_m2s_to_RZ_T)
1044 do j=js,je ; do i=is,ie
1045 fluxes%latent(i,j) = cs%latent_heat_vapor*fluxes%evap(i,j)
1046 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
1047 enddo ; enddo
1048 else
1049 call mom_read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
1050 g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T)
1051 endif
1052 cs%evap_last_lev = time_lev ; fluxes_changed = .true.
1053 endif
1054
1055 if ( (cs%latent_nlev /= cs%evap_nlev) .or. (cs%latent_days_per_rec /= cs%evap_days_per_rec) ) &
1056 time_lev = get_file_time_level(day, cs%latent_nlev, cs%latent_days_per_rec)
1057 if (time_lev /= cs%latent_last_lev) then
1058 if (.not.cs%archaic_OMIP_file) then
1059 call mom_read_data(cs%latentheat_file, cs%latent_var, fluxes%latent(:,:), &
1060 g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1061 do j=js,je ; do i=is,ie
1062 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
1063 enddo ; enddo
1064 endif
1065 cs%latent_last_lev = time_lev ; fluxes_changed = .true.
1066 endif
1067
1068 if ( (cs%sens_nlev /= cs%latent_nlev) .or. (cs%sens_days_per_rec /= cs%latent_days_per_rec) ) &
1069 time_lev = get_file_time_level(day, cs%sens_nlev, cs%sens_days_per_rec)
1070 if (time_lev /= cs%sens_last_lev) then
1071 if (cs%archaic_OMIP_file) then
1072 call mom_read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
1073 g%Domain, timelevel=time_lev, scale=-us%W_m2_to_QRZ_T)
1074 else
1075 call mom_read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
1076 g%Domain, timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1077 endif
1078 cs%sens_last_lev = time_lev ; fluxes_changed = .true.
1079 endif
1080
1081 if ( (cs%SW_nlev /= cs%sens_nlev) .or. (cs%SW_days_per_rec /= cs%sens_days_per_rec) ) &
1082 time_lev = get_file_time_level(day, cs%SW_nlev, cs%SW_days_per_rec)
1083 if (time_lev /= cs%SW_last_lev) then
1084 call mom_read_data(cs%shortwave_file, cs%SW_var, fluxes%sw(:,:), g%Domain, &
1085 timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1086 if (cs%archaic_OMIP_file) then
1087 call mom_read_data(cs%shortwaveup_file, "swup_sfc", temp(:,:), g%Domain, &
1088 timelevel=time_lev, scale=us%W_m2_to_QRZ_T)
1089 do j=js,je ; do i=is,ie
1090 fluxes%sw(i,j) = fluxes%sw(i,j) - temp(i,j)
1091 enddo ; enddo
1092 endif
1093 cs%SW_last_lev = time_lev ; fluxes_changed = .true.
1094 endif
1095
1096 if ( (cs%precip_nlev /= cs%SW_nlev) .or. (cs%precip_days_per_rec /= cs%SW_days_per_rec) ) &
1097 time_lev = get_file_time_level(day, cs%precip_nlev, cs%precip_days_per_rec)
1098 if (time_lev /= cs%precip_last_lev) then
1099 call mom_read_data(cs%snow_file, cs%snow_var, &
1100 fluxes%fprec(:,:), g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T)
1101 call mom_read_data(cs%rain_file, cs%rain_var, &
1102 fluxes%lprec(:,:), g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T)
1103 if (cs%archaic_OMIP_file) then
1104 do j=js,je ; do i=is,ie
1105 fluxes%lprec(i,j) = fluxes%lprec(i,j) - fluxes%fprec(i,j)
1106 enddo ; enddo
1107 endif
1108 cs%precip_last_lev = time_lev ; fluxes_changed = .true.
1109 endif
1110
1111 if ( (cs%runoff_nlev /= cs%precip_nlev) .or. (cs%runoff_days_per_rec /= cs%precip_days_per_rec) ) &
1112 time_lev = get_file_time_level(day, cs%runoff_nlev, cs%runoff_days_per_rec)
1113 if (time_lev /= cs%runoff_last_lev) then
1114 if (cs%archaic_OMIP_file) then
1115 call mom_read_data(cs%runoff_file, cs%lrunoff_var, temp(:,:), &
1116 g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T*us%m_to_L**2)
1117 do j=js,je ; do i=is,ie
1118 fluxes%lrunoff(i,j) = temp(i,j)*g%IareaT(i,j)
1119 enddo ; enddo
1120 call mom_read_data(cs%runoff_file, cs%frunoff_var, temp(:,:), &
1121 g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T*us%m_to_L**2)
1122 do j=js,je ; do i=is,ie
1123 fluxes%frunoff(i,j) = temp(i,j)*g%IareaT(i,j)
1124 enddo ; enddo
1125 else
1126 call mom_read_data(cs%runoff_file, cs%lrunoff_var, fluxes%lrunoff(:,:), &
1127 g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T)
1128 call mom_read_data(cs%runoff_file, cs%frunoff_var, fluxes%frunoff(:,:), &
1129 g%Domain, timelevel=time_lev, scale=us%kg_m2s_to_RZ_T)
1130 endif
1131 cs%runoff_last_lev = time_lev ; fluxes_changed = .true.
1132 endif
1133
1134 ! Read the SST and SSS fields for damping.
1135 if (cs%restorebuoy) then !#CTRL# .or. associated(CS%ctrl_forcing_CSp)) then
1136 time_lev = get_file_time_level(day, cs%SST_nlev, cs%SST_days_per_rec)
1137 if (time_lev /= cs%SST_last_lev) then
1138 call mom_read_data(cs%SSTrestore_file, cs%SST_restore_var, &
1139 cs%T_Restore(:,:), g%Domain, timelevel=time_lev, scale=us%degC_to_C)
1140 cs%SST_last_lev = time_lev
1141 endif
1142
1143 if ( (cs%SSS_nlev /= cs%SST_nlev) .or. (cs%SSS_days_per_rec /= cs%SST_days_per_rec) ) &
1144 time_lev = get_file_time_level(day, cs%SSS_nlev, cs%SSS_days_per_rec)
1145 if (time_lev /= cs%SSS_last_lev) then
1146 call mom_read_data(cs%salinityrestore_file, cs%SSS_restore_var, &
1147 cs%S_Restore(:,:), g%Domain, timelevel=time_lev, scale=us%ppt_to_S)
1148 cs%SSS_last_lev = time_lev
1149 endif
1150 endif
1151
1152 if (fluxes_changed) then
1153 ! mask out land points and compute heat content of water fluxes
1154 ! assume liquid precipitation enters ocean at SST
1155 ! assume frozen precipitation enters ocean at 0degC
1156 ! assume liquid runoff enters ocean at SST
1157 ! assume solid runoff (calving) enters ocean at 0degC
1158 ! mass leaving the ocean has heat_content determined in MOM_diabatic_driver.F90
1159 do j=js,je ; do i=is,ie
1160 fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
1161 fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
1162 fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
1163 fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
1164 fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
1165 fluxes%lw(i,j) = fluxes%lw(i,j) * g%mask2dT(i,j)
1166 fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
1167 fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
1168 fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
1169
1170 fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
1171 fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
1172 fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
1173 enddo ; enddo
1174 endif ! fluxes have changed and need to be masked
1175
1176
1177 ! restoring surface boundary fluxes
1178 if (cs%restorebuoy) then
1179
1180 if (cs%use_temperature) then
1181 do j=js,je ; do i=is,ie
1182 if (g%mask2dT(i,j) > 0.0) then
1183 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1184 ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
1185 fluxes%vprec(i,j) = - (cs%rho_restore*cs%Flux_const_S) * &
1186 (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
1187 (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
1188 else
1189 fluxes%heat_added(i,j) = 0.0
1190 fluxes%vprec(i,j) = 0.0
1191 endif
1192 enddo ; enddo
1193 else
1194 do j=js,je ; do i=is,ie
1195 if (g%mask2dT(i,j) > 0.0) then
1196 fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
1197 (cs%G_Earth * cs%Flux_const / cs%rho_restore)
1198 else
1199 fluxes%buoy(i,j) = 0.0
1200 endif
1201 enddo ; enddo
1202 endif
1203
1204 else ! not RESTOREBUOY
1205 if (.not.cs%use_temperature) then
1206 call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
1207 "The fluxes need to be defined without RESTOREBUOY.")
1208 endif
1209
1210 endif ! end RESTOREBUOY
1211
1212!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
1213!#CTRL# do j=js,je ; do i=is,ie
1214!#CTRL# SST_anom(i,j) = sfc_state%SST(i,j) - CS%T_Restore(i,j)
1215!#CTRL# SSS_anom(i,j) = sfc_state%SSS(i,j) - CS%S_Restore(i,j)
1216!#CTRL# SSS_mean(i,j) = 0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j))
1217!#CTRL# enddo ; enddo
1218!#CTRL# call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
1219!#CTRL# fluxes%vprec, day, dt, G, US, CS%ctrl_forcing_CSp)
1220!#CTRL# endif
1221
1222 call calltree_leave("buoyancy_forcing_from_files")
1223end subroutine buoyancy_forcing_from_files
1224
1225!> Specifies zero surface buoyancy fluxes from data over-ride.
1226subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, US, CS)
1227 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
1228 !! describe the surface state of the ocean.
1229 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
1230 type(time_type), intent(in) :: day !< The time of the fluxes
1231 real, intent(in) :: dt !< The amount of time over which
1232 !! the fluxes apply [T ~> s]
1233 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1234 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1235 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
1236 !! a previous surface_forcing_init call
1237 ! Local variables
1238!#CTRL# real, dimension(SZI_(G),SZJ_(G)) :: &
1239!#CTRL# SST_anom, & ! Instantaneous sea surface temperature anomalies from a
1240!#CTRL# ! target (observed) value [C ~> degC].
1241!#CTRL# SSS_anom, & ! Instantaneous sea surface salinity anomalies from a target
1242!#CTRL# ! (observed) value [S ~> ppt].
1243!#CTRL# SSS_mean ! A (mean?) salinity about which to normalize local salinity
1244!#CTRL# ! anomalies when calculating restorative precipitation anomalies [S ~> ppt].
1245 real :: rhoXcp ! The mean density times the heat capacity [Q R C-1 ~> J m-3 degC-1].
1246 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1247
1248 call calltree_enter("buoyancy_forcing_from_data_override, MOM_surface_forcing.F90")
1249
1250 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1251 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1252
1253 if (cs%use_temperature) rhoxcp = cs%rho_restore * fluxes%C_p
1254
1255 if (.not.cs%dataOverrideIsInitialized) then
1256 call data_override_init(g%Domain)
1257 cs%dataOverrideIsInitialized = .true.
1258 endif
1259
1260 call data_override(g%Domain, 'lw', fluxes%lw, day, scale=us%W_m2_to_QRZ_T)
1261 call data_override(g%Domain, 'sw', fluxes%sw, day, scale=us%W_m2_to_QRZ_T)
1262
1263 ! The normal MOM6 sign conventions are that fluxes%evap and fluxes%sens are positive into the
1264 ! ocean but evap and sens are normally positive quantities in the files.
1265 call data_override(g%Domain, 'evap', fluxes%evap, day, scale=-us%kg_m2s_to_RZ_T)
1266 call data_override(g%Domain, 'sens', fluxes%sens, day, scale=-us%W_m2_to_QRZ_T)
1267
1268 do j=js,je ; do i=is,ie
1269 fluxes%latent(i,j) = cs%latent_heat_vapor*fluxes%evap(i,j)
1270 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
1271 enddo ; enddo
1272
1273 call data_override(g%Domain, 'snow', fluxes%fprec, day, scale=us%kg_m2s_to_RZ_T)
1274 call data_override(g%Domain, 'rain', fluxes%lprec, day, scale=us%kg_m2s_to_RZ_T)
1275 call data_override(g%Domain, 'runoff', fluxes%lrunoff, day, scale=us%kg_m2s_to_RZ_T)
1276 call data_override(g%Domain, 'calving', fluxes%frunoff, day, scale=us%kg_m2s_to_RZ_T)
1277
1278! Read the SST and SSS fields for damping.
1279 if (cs%restorebuoy) then !#CTRL# .or. associated(CS%ctrl_forcing_CSp)) then
1280 call data_override(g%Domain, 'SST_restore', cs%T_restore, day, scale=us%degC_to_C)
1281 call data_override(g%Domain, 'SSS_restore', cs%S_restore, day, scale=us%ppt_to_S)
1282 endif
1283
1284 ! restoring boundary fluxes
1285 if (cs%restorebuoy) then
1286 if (cs%use_temperature) then
1287 do j=js,je ; do i=is,ie
1288 if (g%mask2dT(i,j) > 0.0) then
1289 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1290 ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
1291 fluxes%vprec(i,j) = - (cs%rho_restore*cs%Flux_const_S) * &
1292 (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
1293 (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
1294 else
1295 fluxes%heat_added(i,j) = 0.0
1296 fluxes%vprec(i,j) = 0.0
1297 endif
1298 enddo ; enddo
1299 else
1300 do j=js,je ; do i=is,ie
1301 if (g%mask2dT(i,j) > 0.0) then
1302 fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
1303 (cs%G_Earth * cs%Flux_const / cs%rho_restore)
1304 else
1305 fluxes%buoy(i,j) = 0.0
1306 endif
1307 enddo ; enddo
1308 endif
1309 else ! not RESTOREBUOY
1310 if (.not.cs%use_temperature) then
1311 call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
1312 "The fluxes need to be defined without RESTOREBUOY.")
1313 endif
1314 endif ! end RESTOREBUOY
1315
1316
1317 ! mask out land points and compute heat content of water fluxes
1318 ! assume liquid precip enters ocean at SST
1319 ! assume frozen precip enters ocean at 0degC
1320 ! assume liquid runoff enters ocean at SST
1321 ! assume solid runoff (calving) enters ocean at 0degC
1322 ! mass leaving ocean has heat_content determined in MOM_diabatic_driver.F90
1323 do j=js,je ; do i=is,ie
1324 fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
1325 fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
1326 fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
1327 fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
1328 fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
1329 fluxes%lw(i,j) = fluxes%lw(i,j) * g%mask2dT(i,j)
1330 fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
1331 fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
1332 fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
1333
1334 fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
1335 fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
1336 fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
1337 enddo ; enddo
1338
1339!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
1340!#CTRL# do j=js,je ; do i=is,ie
1341!#CTRL# SST_anom(i,j) = sfc_state%SST(i,j) - CS%T_Restore(i,j)
1342!#CTRL# SSS_anom(i,j) = sfc_state%SSS(i,j) - CS%S_Restore(i,j)
1343!#CTRL# SSS_mean(i,j) = 0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j))
1344!#CTRL# enddo ; enddo
1345!#CTRL# call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
1346!#CTRL# fluxes%vprec, day, dt, G, US, CS%ctrl_forcing_CSp)
1347!#CTRL# endif
1348
1349 call calltree_leave("buoyancy_forcing_from_data_override")
1351
1352!> This subroutine specifies zero surface buoyancy fluxes
1353subroutine buoyancy_forcing_zero(sfc_state, fluxes, day, dt, G, CS)
1354 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
1355 !! describe the surface state of the ocean.
1356 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
1357 type(time_type), intent(in) :: day !< The time of the fluxes
1358 real, intent(in) :: dt !< The amount of time over which
1359 !! the fluxes apply [T ~> s]
1360 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1361 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
1362 !! a previous surface_forcing_init call
1363 ! Local variables
1364 integer :: i, j, is, ie, js, je
1365
1366 call calltree_enter("buoyancy_forcing_zero, MOM_surface_forcing.F90")
1367 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1368
1369 if (cs%use_temperature) then
1370 do j=js,je ; do i=is,ie
1371 fluxes%evap(i,j) = 0.0
1372 fluxes%lprec(i,j) = 0.0
1373 fluxes%fprec(i,j) = 0.0
1374 fluxes%vprec(i,j) = 0.0
1375 fluxes%lrunoff(i,j) = 0.0
1376 fluxes%frunoff(i,j) = 0.0
1377 fluxes%lw(i,j) = 0.0
1378 fluxes%latent(i,j) = 0.0
1379 fluxes%sens(i,j) = 0.0
1380 fluxes%sw(i,j) = 0.0
1381 fluxes%latent_evap_diag(i,j) = 0.0
1382 fluxes%latent_fprec_diag(i,j) = 0.0
1383 fluxes%latent_frunoff_diag(i,j) = 0.0
1384 enddo ; enddo
1385 else
1386 do j=js,je ; do i=is,ie
1387 fluxes%buoy(i,j) = 0.0
1388 enddo ; enddo
1389 endif
1390
1391 call calltree_leave("buoyancy_forcing_zero")
1392end subroutine buoyancy_forcing_zero
1393
1394
1395!> Sets up spatially and temporally constant surface heat fluxes.
1396subroutine buoyancy_forcing_const(sfc_state, fluxes, day, dt, G, US, CS)
1397 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
1398 !! describe the surface state of the ocean.
1399 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
1400 type(time_type), intent(in) :: day !< The time of the fluxes
1401 real, intent(in) :: dt !< The amount of time over which
1402 !! the fluxes apply [T ~> s]
1403 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1404 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1405 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
1406 !! a previous surface_forcing_init call
1407 ! Local variables
1408 integer :: i, j, is, ie, js, je
1409 call calltree_enter("buoyancy_forcing_const, MOM_surface_forcing.F90")
1410 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1411
1412 if (cs%use_temperature) then
1413 do j=js,je ; do i=is,ie
1414 fluxes%evap(i,j) = 0.0
1415 fluxes%lprec(i,j) = 0.0
1416 fluxes%fprec(i,j) = 0.0
1417 fluxes%vprec(i,j) = 0.0
1418 fluxes%lrunoff(i,j) = 0.0
1419 fluxes%frunoff(i,j) = 0.0
1420 fluxes%lw(i,j) = 0.0
1421 fluxes%latent(i,j) = 0.0
1422 fluxes%sens(i,j) = cs%constantHeatForcing * g%mask2dT(i,j)
1423 fluxes%sw(i,j) = 0.0
1424 fluxes%latent_evap_diag(i,j) = 0.0
1425 fluxes%latent_fprec_diag(i,j) = 0.0
1426 fluxes%latent_frunoff_diag(i,j) = 0.0
1427 enddo ; enddo
1428 else
1429 do j=js,je ; do i=is,ie
1430 fluxes%buoy(i,j) = 0.0
1431 enddo ; enddo
1432 endif
1433
1434 call calltree_leave("buoyancy_forcing_const")
1435end subroutine buoyancy_forcing_const
1436
1437!> Sets surface fluxes of heat and salinity by restoring to temperature and
1438!! salinity profiles that vary linearly with latitude.
1439subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS)
1440 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
1441 !! describe the surface state of the ocean.
1442 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
1443 type(time_type), intent(in) :: day !< The time of the fluxes
1444 real, intent(in) :: dt !< The amount of time over which
1445 !! the fluxes apply [T ~> s]
1446 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1447 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1448 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
1449 !! a previous surface_forcing_init call
1450 ! Local variables
1451 real :: y ! The latitude relative to the south normalized by the domain extent [nondim]
1452 real :: T_restore ! The temperature towards which to restore [C ~> degC]
1453 real :: S_restore ! The salinity towards which to restore [S ~> ppt]
1454 integer :: i, j, is, ie, js, je
1455
1456 call calltree_enter("buoyancy_forcing_linear, MOM_surface_forcing.F90")
1457 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1458
1459 ! This case has no surface buoyancy forcing.
1460 if (cs%use_temperature) then
1461 do j=js,je ; do i=is,ie
1462 fluxes%evap(i,j) = 0.0
1463 fluxes%lprec(i,j) = 0.0
1464 fluxes%fprec(i,j) = 0.0
1465 fluxes%vprec(i,j) = 0.0
1466 fluxes%lrunoff(i,j) = 0.0
1467 fluxes%frunoff(i,j) = 0.0
1468 fluxes%lw(i,j) = 0.0
1469 fluxes%latent(i,j) = 0.0
1470 fluxes%sens(i,j) = 0.0
1471 fluxes%sw(i,j) = 0.0
1472 fluxes%latent_evap_diag(i,j) = 0.0
1473 fluxes%latent_fprec_diag(i,j) = 0.0
1474 fluxes%latent_frunoff_diag(i,j) = 0.0
1475 enddo ; enddo
1476 else
1477 do j=js,je ; do i=is,ie
1478 fluxes%buoy(i,j) = 0.0
1479 enddo ; enddo
1480 endif
1481
1482 if (cs%restorebuoy) then
1483 if (cs%use_temperature) then
1484 do j=js,je ; do i=is,ie
1485 y = (g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat
1486 t_restore = cs%T_south + (cs%T_north-cs%T_south)*y
1487 s_restore = cs%S_south + (cs%S_north-cs%S_south)*y
1488 if (g%mask2dT(i,j) > 0.0) then
1489 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1490 ((t_restore - sfc_state%SST(i,j)) * ((cs%rho_restore * fluxes%C_p) * cs%Flux_const))
1491 fluxes%vprec(i,j) = - (cs%rho_restore*cs%Flux_const) * &
1492 (s_restore - sfc_state%SSS(i,j)) / &
1493 (0.5*(sfc_state%SSS(i,j) + s_restore))
1494 else
1495 fluxes%heat_added(i,j) = 0.0
1496 fluxes%vprec(i,j) = 0.0
1497 endif
1498 enddo ; enddo
1499 else
1500 call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
1501 "RESTOREBUOY to linear not written yet.")
1502 !do j=js,je ; do i=is,ie
1503 ! if (G%mask2dT(i,j) > 0.0) then
1504 ! fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
1505 ! (CS%G_Earth * CS%Flux_const / CS%rho_restore)
1506 ! else
1507 ! fluxes%buoy(i,j) = 0.0
1508 ! endif
1509 !enddo ; enddo
1510 endif
1511 else ! not RESTOREBUOY
1512 if (.not.cs%use_temperature) then
1513 call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
1514 "The fluxes need to be defined without RESTOREBUOY.")
1515 endif
1516 endif ! end RESTOREBUOY
1517
1518 call calltree_leave("buoyancy_forcing_linear")
1519end subroutine buoyancy_forcing_linear
1520
1521!> Return a time record to read from a file based on the model time, the number of time records in
1522!! that file and the number of time records per model day.
1523function get_file_time_level(Time, nlev_file, days_per_rec) result (time_lev)
1524 type(time_type), intent(in) :: time !< The time of the fluxes
1525 integer , intent(in) :: nlev_file !< The number of time records in a forcing file
1526 integer , intent(in) :: days_per_rec !< If positive, the number of days spanned by each
1527 !! time record in a file, if negative the number
1528 !! time records per day, or if 0 determine this
1529 !! by guessing based on the number of records in
1530 !! the file. If this is 31, the time levels will
1531 !! be based on the months of the calendar.
1532 !! Setting this larger than 1000000 will always
1533 !! cause the time level to be set to 1.
1534 integer :: time_lev !< The time level in a file that will be read.
1535
1536 ! Local variables
1537 integer :: days, seconds ! The number of days and seconds since the start of the calendar
1538 integer :: year, month, day, hour, minute, second ! The components of the model time
1539 integer :: recs_per_day ! The number of file time records per day
1540 integer :: recs ! The number of time levels into the file to read without
1541 ! taking the periodicity of the file into account.
1542
1543 if ( (days_per_rec >= 1000000) .or. &
1544 ( (days_per_rec == 0) .and. .not.((nlev_file == 12) .or. (nlev_file == 365)) ) ) then
1545 ! The second condition above is to recreate the existing behavior, but it should perhaps be
1546 ! phased out.
1547 time_lev = 1
1548 elseif ( (days_per_rec == 31) .or. ((days_per_rec == 0) .and. (nlev_file == 12)) ) then
1549 call get_date(time, year, month, day, hour, minute, second)
1550 time_lev = month
1551 else
1552 call get_time(time, seconds, days)
1553 if ( (days_per_rec == 0) .or. (abs(days_per_rec) == 1) ) then
1554 recs = days
1555 elseif (days_per_rec < 0) then
1556 recs_per_day = -days_per_rec
1557 recs = days * recs_per_day + ( (recs_per_day*set_time(seconds, 0)) / set_time(0, 1) )
1558 ! When integer rounding in the time-type arithmetic is considered, the line above is equivalent to:
1559 ! seconds_per_day = set_time(0, 1) / set_time(1, 0)
1560 ! recs = days * recs_per_day + floor(real(recs_per_day*seconds) / real(seconds_per_day))
1561 else
1562 recs = days / days_per_rec
1563 endif
1564 time_lev = recs - nlev_file*floor(real(recs) / real(nlev_file)) + 1
1565 endif
1566
1567end function get_file_time_level
1568
1569!> Sets the necessary MARBL forcings via the data override facility.
1570subroutine marbl_forcing_from_data_override(fluxes, day, G, US, CS)
1571 type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
1572 type(time_type), intent(in) :: day !< The time of the fluxes
1573 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1574 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1575 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
1576 !! a previous surface_forcing_init call
1577 ! Local variables
1578 real, pointer, dimension(:,:) :: atm_co2_prog =>null() !< Prognostic atmospheric CO2 concentration [ppm]
1579 real, pointer, dimension(:,:) :: atm_co2_diag =>null() !< Diagnostic atmospheric CO2 concentration [ppm]
1580 real, pointer, dimension(:,:) :: atm_fine_dust_flux =>null() !< Fine dust flux from atmosphere
1581 !! [R Z T-1 ~> kg m-2 s-1]
1582 real, pointer, dimension(:,:) :: atm_coarse_dust_flux =>null() !< Coarse dust flux from atmosphere
1583 !! [R Z T-1 ~> kg m-2 s-1]
1584 real, pointer, dimension(:,:) :: seaice_dust_flux =>null() !< Dust flux from seaice
1585 !! [R Z T-1 ~> kg m-2 s-1]
1586 real, pointer, dimension(:,:) :: atm_bc_flux =>null() !< Black carbon flux from atmosphere
1587 !! [R Z T-1 ~> kg m-2 s-1]
1588 real, pointer, dimension(:,:) :: seaice_bc_flux =>null() !< Black carbon flux from seaice
1589 !! [R Z T-1 ~> kg m-2 s-1]
1590 real, pointer, dimension(:,:) :: nhx_dep =>null() !< Nitrogen deposition
1591 !! [R Z T-1 ~> kg m-2 s-1]
1592 real, pointer, dimension(:,:) :: noy_dep =>null() !< Nitrogen deposition
1593 !! [R Z T-1 ~> kg m-2 s-1]
1594 integer :: isc, iec, jsc, jec
1595
1596 ! Necessary null pointers for arguments to convert_driver_fields_to_forcings()
1597 ! Since they are null, MARBL will not use multiple ice categories
1598 real, pointer, dimension(:,:) :: afracr =>null()
1599 real, pointer, dimension(:,:) :: swnet_afracr =>null()
1600 real, pointer, dimension(:,:,:) :: swpen_ifrac_n =>null()
1601 real, pointer, dimension(:,:,:) :: ifrac_n =>null()
1602
1603 call calltree_enter("MARBL_forcing_from_data_override, MOM_surface_forcing.F90")
1604
1605 if (.not.cs%dataOverrideIsInitialized) then
1606 call data_override_init(g%Domain)
1607 cs%dataOverrideIsInitialized = .true.
1608 endif
1609
1610 ! Allocate memory for pointers
1611 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1612 allocate ( atm_co2_prog(isc:iec,jsc:jec), &
1613 atm_co2_diag(isc:iec,jsc:jec), &
1614 atm_fine_dust_flux(isc:iec,jsc:jec), &
1615 atm_coarse_dust_flux(isc:iec,jsc:jec), &
1616 seaice_dust_flux(isc:iec,jsc:jec), &
1617 atm_bc_flux(isc:iec,jsc:jec), &
1618 seaice_bc_flux(isc:iec,jsc:jec), &
1619 nhx_dep(isc:iec,jsc:jec), &
1620 noy_dep(isc:iec,jsc:jec), &
1621 source=0.0)
1622
1623
1624 ! fluxes used directly as MARBL inputs
1625 ! (should be scaled)
1626 call data_override(g%Domain, 'ice_fraction', fluxes%ice_fraction, day)
1627 call data_override(g%Domain, 'u10_sqr', fluxes%u10_sqr, day, scale=us%m_s_to_L_T**2)
1628
1629 ! fluxes used to compute MARBL inputs
1630 ! These are kept in physical units, and will be scaled appropriately in
1631 ! convert_driver_fields_to_forcings()
1632 call data_override(g%Domain, 'atm_co2_prog', atm_co2_prog, day)
1633 call data_override(g%Domain, 'atm_co2_diag', atm_co2_diag, day)
1634 call data_override(g%Domain, 'atm_fine_dust_flux', atm_fine_dust_flux, day)
1635 call data_override(g%Domain, 'atm_coarse_dust_flux', atm_coarse_dust_flux, day)
1636 call data_override(g%Domain, 'atm_bc_flux', atm_bc_flux, day)
1637 call data_override(g%Domain, 'seaice_dust_flux', seaice_dust_flux, day)
1638 call data_override(g%Domain, 'seaice_bc_flux', seaice_bc_flux, day)
1639 call data_override(g%Domain, 'nhx_dep', nhx_dep, day)
1640 call data_override(g%Domain, 'noy_dep', noy_dep, day)
1641
1642 call convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, &
1643 seaice_dust_flux, atm_bc_flux, seaice_bc_flux, &
1644 nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, &
1645 afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, &
1646 day, g, us, 0, 0, fluxes, cs%marbl_forcing_CSp)
1647
1648 deallocate ( atm_co2_prog, &
1649 atm_co2_diag, &
1650 atm_fine_dust_flux, &
1651 atm_coarse_dust_flux, &
1652 seaice_dust_flux, &
1653 atm_bc_flux, &
1654 seaice_bc_flux, &
1655 nhx_dep, &
1656 noy_dep)
1657
1658 call calltree_leave("MARBL_forcing_from_data_override")
1659
1661
1662!> Save a restart file for the forcing fields
1663subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
1664 filename_suffix)
1665 type(surface_forcing_cs), pointer :: cs !< pointer to control structure returned by
1666 !! a previous surface_forcing_init call
1667 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
1668 type(time_type), intent(in) :: time !< model time at this call; needed for mpp_write calls
1669 character(len=*), intent(in) :: directory !< directory into which to write these restart files
1670 logical, optional, intent(in) :: time_stamped !< If true, the restart file names
1671 !! include a unique time stamp; the default is false.
1672 character(len=*), optional, intent(in) :: filename_suffix !< optional suffix (e.g., a time-stamp)
1673 !! to append to the restart file name
1674
1675 if (.not.associated(cs)) return
1676 if (.not.associated(cs%restart_CSp)) return
1677
1678 call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1679
1680end subroutine forcing_save_restart
1681
1682!> Initialize the surface forcing module
1683subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_CSp)
1684 type(time_type), intent(in) :: time !< The current model time
1685 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1686 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1687 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1688 type(diag_ctrl), target, intent(inout) :: diag !< structure used to regulate diagnostic output
1689 type(surface_forcing_cs), pointer :: cs !< pointer to control structure returned by
1690 !! a previous surface_forcing_init call
1691 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< Forcing for tracers?
1692
1693 ! Local variables
1694 type(directories) :: dirs
1695 logical :: new_sim
1696 type(time_type) :: time_frc
1697 ! This include declares and sets the variable "version".
1698# include "version_variable.h"
1699 real :: flux_const_default ! The unscaled value of FLUXCONST [m day-1]
1700 logical :: boussinesq ! If true, this run is fully Boussinesq
1701 logical :: semi_boussinesq ! If true, this run is partially non-Boussinesq
1702 logical :: fix_ustar_gustless_bug ! If false, include a bug using an older run-time parameter.
1703 logical :: test_value ! This is used to determine whether a logical parameter is being set explicitly.
1704 logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly.
1705 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1706 character(len=40) :: mdl = "MOM_surface_forcing" ! This module's name.
1707 character(len=200) :: filename, gust_file ! The name of the gustiness input file.
1708
1709 if (associated(cs)) then
1710 call mom_error(warning, "surface_forcing_init called with an associated "// &
1711 "control structure.")
1712 return
1713 endif
1714 allocate(cs)
1715
1716 id_clock_forcing=cpu_clock_id('(Ocean surface forcing)', grain=clock_module)
1717 call cpu_clock_begin(id_clock_forcing)
1718
1719 cs%diag => diag
1720 if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
1721
1722 ! Read all relevant parameters and write them to the model log.
1723 call log_version(param_file, mdl, version, '')
1724 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
1725 "If true, Temperature and salinity are used as state "//&
1726 "variables.", default=.true.)
1727 call get_param(param_file, "MOM", "BOUSSINESQ", boussinesq, &
1728 "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
1729 call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_boussinesq, &
1730 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
1731 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
1732 "height changes. This only applies if BOUSSINESQ is false.", &
1733 default=.true., do_not_log=.true.)
1734 cs%nonBous = .not.(boussinesq .or. semi_boussinesq)
1735 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, &
1736 "The directory in which all input files are found.", &
1737 default=".")
1738 cs%inputdir = slasher(cs%inputdir)
1739
1740 call get_param(param_file, mdl, "ADIABATIC", cs%adiabatic, &
1741 "There are no diapycnal mass fluxes if ADIABATIC is "//&
1742 "true. This assumes that KD = KDML = 0.0 and that "//&
1743 "there is no buoyancy forcing, but makes the model "//&
1744 "faster by eliminating subroutine calls.", default=.false.)
1745 call get_param(param_file, mdl, "VARIABLE_WINDS", cs%variable_winds, &
1746 "If true, the winds vary in time after the initialization.", &
1747 default=.true.)
1748 call get_param(param_file, mdl, "VARIABLE_BUOYFORCE", cs%variable_buoyforce, &
1749 "If true, the buoyancy forcing varies in time after the "//&
1750 "initialization of the model.", default=.true.)
1751
1752 ! Determine parameters related to the buoyancy forcing.
1753 call get_param(param_file, mdl, "BUOY_CONFIG", cs%buoy_config, &
1754 "The character string that indicates how buoyancy forcing is specified. Valid "//&
1755 "options include (file), (data_override), (zero), (const), (linear), (MESO), "//&
1756 "(SCM_CVmix_tests), (BFB), (dumbbell), (USER) and (NONE).", default="zero")
1757 if (trim(cs%buoy_config) == "file") then
1758 call get_param(param_file, mdl, "ARCHAIC_OMIP_FORCING_FILE", cs%archaic_OMIP_file, &
1759 "If true, use the forcing variable decomposition from "//&
1760 "the old German OMIP prescription that predated CORE. If "//&
1761 "false, use the variable groupings available from MOM "//&
1762 "output diagnostics of forcing variables.", default=.true.)
1763 if (cs%archaic_OMIP_file) then
1764 call get_param(param_file, mdl, "LONGWAVEDOWN_FILE", cs%longwave_file, &
1765 "The file with the downward longwave heat flux, in "//&
1766 "variable lwdn_sfc.", fail_if_missing=.true.)
1767 call get_param(param_file, mdl, "LONGWAVEUP_FILE", cs%longwaveup_file, &
1768 "The file with the upward longwave heat flux, in "//&
1769 "variable lwup_sfc.", fail_if_missing=.true.)
1770 call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
1771 "The file with the evaporative moisture flux, in "//&
1772 "variable evap.", fail_if_missing=.true.)
1773 call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1774 "The file with the sensible heat flux, in "//&
1775 "variable shflx.", fail_if_missing=.true.)
1776 call get_param(param_file, mdl, "SHORTWAVEUP_FILE", cs%shortwaveup_file, &
1777 "The file with the upward shortwave heat flux.", &
1778 fail_if_missing=.true.)
1779 call get_param(param_file, mdl, "SHORTWAVEDOWN_FILE", cs%shortwave_file, &
1780 "The file with the downward shortwave heat flux.", &
1781 fail_if_missing=.true.)
1782 call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
1783 "The file with the downward frozen precip flux, in "//&
1784 "variable snow.", fail_if_missing=.true.)
1785 call get_param(param_file, mdl, "PRECIP_FILE", cs%rain_file, &
1786 "The file with the downward total precip flux, in "//&
1787 "variable precip.", fail_if_missing=.true.)
1788 call get_param(param_file, mdl, "FRESHDISCHARGE_FILE", cs%runoff_file, &
1789 "The file with the fresh and frozen runoff/calving fluxes, "//&
1790 "invariables disch_w and disch_s.", fail_if_missing=.true.)
1791
1792 ! These variable names are hard-coded, per the archaic OMIP conventions.
1793 cs%latentheat_file = cs%evaporation_file ; cs%latent_var = "evap"
1794 cs%LW_var = "lwdn_sfc" ; cs%SW_var = "swdn_sfc" ; cs%sens_var = "shflx"
1795 cs%evap_var = "evap" ; cs%rain_var = "precip" ; cs%snow_var = "snow"
1796 cs%lrunoff_var = "disch_w" ; cs%frunoff_var = "disch_s"
1797
1798 else
1799 call get_param(param_file, mdl, "LONGWAVE_FILE", cs%longwave_file, &
1800 "The file with the longwave heat flux, in the variable "//&
1801 "given by LONGWAVE_FORCING_VAR.", fail_if_missing=.true.)
1802 call get_param(param_file, mdl, "LONGWAVE_FORCING_VAR", cs%LW_var, &
1803 "The variable with the longwave forcing field.", default="LW")
1804 call get_param(param_file, mdl, "LONGWAVE_FILE_DAYS_PER_RECORD", cs%LW_days_per_rec, &
1805 "If positive the number of days of longwave fluxes per time level in LONGWAVE_FILE, "//&
1806 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1807 "or if 0 the model will guess the right value based on the file size.", &
1808 default=0)
1809
1810 call get_param(param_file, mdl, "SHORTWAVE_FILE", cs%shortwave_file, &
1811 "The file with the shortwave heat flux, in the variable "//&
1812 "given by SHORTWAVE_FORCING_VAR.", fail_if_missing=.true.)
1813 call get_param(param_file, mdl, "SHORTWAVE_FORCING_VAR", cs%SW_var, &
1814 "The variable with the shortwave forcing field.", default="SW")
1815 call get_param(param_file, mdl, "SHORTWAVE_FILE_DAYS_PER_RECORD", cs%SW_days_per_rec, &
1816 "If positive the number of days of shortwave fluxes per time level in SHORTWAVE_FILE, "//&
1817 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1818 "or if 0 the model will guess the right value based on the file size.", &
1819 default=cs%LW_days_per_rec)
1820
1821 call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
1822 "The file with the evaporative moisture flux, in the "//&
1823 "variable given by EVAP_FORCING_VAR.", fail_if_missing=.true.)
1824 call get_param(param_file, mdl, "EVAP_FORCING_VAR", cs%evap_var, &
1825 "The variable with the evaporative moisture flux.", &
1826 default="evap")
1827 call get_param(param_file, mdl, "EVAPORATION_FILE_DAYS_PER_RECORD", cs%evap_days_per_rec, &
1828 "If positive the number of days of evaporation per time level in EVAPORATION_FILE, "//&
1829 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1830 "or if 0 the model will guess the right value based on the file size.", &
1831 default=cs%LW_days_per_rec)
1832
1833 call get_param(param_file, mdl, "LATENTHEAT_FILE", cs%latentheat_file, &
1834 "The file with the latent heat flux, in the variable "//&
1835 "given by LATENT_FORCING_VAR.", fail_if_missing=.true.)
1836 call get_param(param_file, mdl, "LATENT_FORCING_VAR", cs%latent_var, &
1837 "The variable with the latent heat flux.", default="latent")
1838 call get_param(param_file, mdl, "LATENTHEAT_FILE_DAYS_PER_RECORD", cs%latent_days_per_rec, &
1839 "If positive the number of days of latent heat fluxes per time level in LATENTHEAT_FILE, "//&
1840 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1841 "or if 0 the model will guess the right value based on the file size.", &
1842 default=cs%LW_days_per_rec)
1843
1844 call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1845 "The file with the sensible heat flux, in the variable "//&
1846 "given by SENSIBLE_FORCING_VAR.", fail_if_missing=.true.)
1847 call get_param(param_file, mdl, "SENSIBLE_FORCING_VAR", cs%sens_var, &
1848 "The variable with the sensible heat flux.", default="sensible")
1849 call get_param(param_file, mdl, "SENSIBLEHEAT_FILE_DAYS_PER_RECORD", cs%sens_days_per_rec, &
1850 "If positive the number of days of sensible heat fluxes per time level in SENSIBLEHEAT_FILE, "//&
1851 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1852 "or if 0 the model will guess the right value based on the file size.", &
1853 default=cs%LW_days_per_rec)
1854
1855 call get_param(param_file, mdl, "RAIN_FILE", cs%rain_file, &
1856 "The file with the liquid precipitation flux, in the "//&
1857 "variable given by RAIN_FORCING_VAR.", fail_if_missing=.true.)
1858 call get_param(param_file, mdl, "RAIN_FORCING_VAR", cs%rain_var, &
1859 "The variable with the liquid precipitation flux.", &
1860 default="liq_precip")
1861 call get_param(param_file, mdl, "RAIN_FILE_DAYS_PER_RECORD", cs%precip_days_per_rec, &
1862 "If positive the number of days of rain fluxes per time level in RAIN_FILE, "//&
1863 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1864 "or if 0 the model will guess the right value based on the file size.", &
1865 default=cs%LW_days_per_rec)
1866 call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
1867 "The file with the frozen precipitation flux, in the "//&
1868 "variable given by SNOW_FORCING_VAR.", fail_if_missing=.true.)
1869 call get_param(param_file, mdl, "SNOW_FORCING_VAR", cs%snow_var, &
1870 "The variable with the frozen precipitation flux.", &
1871 default="froz_precip")
1872 call get_param(param_file, mdl, "SHORTWAVE_FILE_DAYS_PER_RECORD", cs%SW_days_per_rec, &
1873 "If positive the number of days of shortwave fluxes per time level in SHORTWAVE_FILE, "//&
1874 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1875 "or if 0 the model will guess the right value based on the file size.", &
1876 default=cs%LW_days_per_rec)
1877
1878 call get_param(param_file, mdl, "RUNOFF_FILE", cs%runoff_file, &
1879 "The file with the fresh and frozen runoff/calving "//&
1880 "fluxes, in variables given by LIQ_RUNOFF_FORCING_VAR "//&
1881 "and FROZ_RUNOFF_FORCING_VAR.", fail_if_missing=.true.)
1882 call get_param(param_file, mdl, "LIQ_RUNOFF_FORCING_VAR", cs%lrunoff_var, &
1883 "The variable with the liquid runoff flux.", &
1884 default="liq_runoff")
1885 call get_param(param_file, mdl, "FROZ_RUNOFF_FORCING_VAR", cs%frunoff_var, &
1886 "The variable with the frozen runoff flux.", &
1887 default="froz_runoff")
1888 call get_param(param_file, mdl, "RUNOFF_FILE_DAYS_PER_RECORD", cs%SW_days_per_rec, &
1889 "If positive the number of days of runoff per time level in RUNOFF_FILE, "//&
1890 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1891 "or if 0 the model will guess the right value based on the file size.", &
1892 default=0)
1893 endif
1894
1895 call get_param(param_file, mdl, "SSTRESTORE_FILE", cs%SSTrestore_file, &
1896 "The file with the SST toward which to restore in the "//&
1897 "variable given by SST_RESTORE_VAR.", fail_if_missing=.true.)
1898 call get_param(param_file, mdl, "SALINITYRESTORE_FILE", cs%salinityrestore_file, &
1899 "The file with the surface salinity toward which to "//&
1900 "restore in the variable given by SSS_RESTORE_VAR.", &
1901 fail_if_missing=.true.)
1902 if (cs%archaic_OMIP_file) then
1903 cs%SST_restore_var = "TEMP" ; cs%SSS_restore_var = "SALT"
1904 else
1905 call get_param(param_file, mdl, "SST_RESTORE_VAR", cs%SST_restore_var, &
1906 "The variable with the SST toward which to restore.", &
1907 default="SST")
1908 call get_param(param_file, mdl, "SSTRESTORE_FILE_DAYS_PER_RECORD", cs%SST_days_per_rec, &
1909 "If positive the number of days of SST per time level in SSTRESTORE_FILE, "//&
1910 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1911 "or if 0 the model will guess the right value based on the file size.", &
1912 default=0)
1913 call get_param(param_file, mdl, "SSS_RESTORE_VAR", cs%SSS_restore_var, &
1914 "The variable with the SSS toward which to restore.", &
1915 default="SSS")
1916 call get_param(param_file, mdl, "SALINITYRESTORE_FILE_DAYS_PER_RECORD", cs%SSS_days_per_rec, &
1917 "If positive the number of days of salinity per time level in SALINITYRESTORE_FILE, "//&
1918 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1919 "or if 0 the model will guess the right value based on the file size.", &
1920 default=cs%SST_days_per_rec)
1921 endif
1922
1923 ! Add inputdir to the file names.
1924 cs%shortwave_file = trim(cs%inputdir)//trim(cs%shortwave_file)
1925 cs%longwave_file = trim(cs%inputdir)//trim(cs%longwave_file)
1926 cs%sensibleheat_file = trim(cs%inputdir)//trim(cs%sensibleheat_file)
1927 cs%latentheat_file = trim(cs%inputdir)//trim(cs%latentheat_file)
1928 cs%evaporation_file = trim(cs%inputdir)//trim(cs%evaporation_file)
1929 cs%snow_file = trim(cs%inputdir)//trim(cs%snow_file)
1930 cs%rain_file = trim(cs%inputdir)//trim(cs%rain_file)
1931 cs%runoff_file = trim(cs%inputdir)//trim(cs%runoff_file)
1932
1933 cs%shortwaveup_file = trim(cs%inputdir)//trim(cs%shortwaveup_file)
1934 cs%longwaveup_file = trim(cs%inputdir)//trim(cs%longwaveup_file)
1935
1936 cs%SSTrestore_file = trim(cs%inputdir)//trim(cs%SSTrestore_file)
1937 cs%salinityrestore_file = trim(cs%inputdir)//trim(cs%salinityrestore_file)
1938 elseif (trim(cs%buoy_config) == "const") then
1939 call get_param(param_file, mdl, "SENSIBLE_HEAT_FLUX", cs%constantHeatForcing, &
1940 "A constant heat forcing (positive into ocean) applied "//&
1941 "through the sensible heat flux field. ", &
1942 units='W/m2', scale=us%W_m2_to_QRZ_T, fail_if_missing=.true.)
1943 endif
1944
1945 ! Determine parameters related to the wind forcing.
1946 call get_param(param_file, mdl, "WIND_CONFIG", cs%wind_config, &
1947 "The character string that indicates how wind forcing is specified. Valid "//&
1948 "options include (file), (data_override), (2gyre), (1gyre), (gyres), (zero), "//&
1949 "(const), (Neverworld), (scurves), (ideal_hurr), (SCM_CVmix_tests) and (USER).", &
1950 default="zero")
1951 if (trim(cs%wind_config) == "file") then
1952 call get_param(param_file, mdl, "WIND_FILE", cs%wind_file, &
1953 "The file in which the wind stresses are found in "//&
1954 "variables STRESS_X and STRESS_Y.", fail_if_missing=.true.)
1955 call get_param(param_file, mdl, "WINDSTRESS_X_VAR",cs%stress_x_var, &
1956 "The name of the x-wind stress variable in WIND_FILE.", &
1957 default="STRESS_X")
1958 call get_param(param_file, mdl, "WINDSTRESS_Y_VAR", cs%stress_y_var, &
1959 "The name of the y-wind stress variable in WIND_FILE.", &
1960 default="STRESS_Y")
1961 call get_param(param_file, mdl, "WIND_STAGGER",cs%wind_stagger, &
1962 "A character indicating how the wind stress components "//&
1963 "are staggered in WIND_FILE. This may be A or C for now.", &
1964 default="C")
1965 call get_param(param_file, mdl, "WINDSTRESS_SCALE", cs%wind_scale, &
1966 "A value by which the wind stresses in WIND_FILE are rescaled.", &
1967 default=1.0, units="nondim")
1968 call get_param(param_file, mdl, "USTAR_FORCING_VAR", cs%ustar_var, &
1969 "The name of the friction velocity variable in WIND_FILE "//&
1970 "or blank to get ustar from the wind stresses plus the "//&
1971 "gustiness.", default=" ")
1972 cs%wind_file = trim(cs%inputdir) // trim(cs%wind_file)
1973 call get_param(param_file, mdl, "WIND_FILE_DAYS_PER_RECORD", cs%wind_days_per_rec, &
1974 "If positive the number of days of wind stress per time level in WIND_FILE, "//&
1975 "or if negative the number of time levels per day. If 31 change forcing monthly, "//&
1976 "or if 0 the model will guess the right value based on the file size.", &
1977 default=0)
1978 endif
1979 if (trim(cs%wind_config) == "gyres") then
1980 call get_param(param_file, mdl, "TAUX_CONST", cs%gyres_taux_const, &
1981 "With the gyres wind_config, the constant offset in the "//&
1982 "zonal wind stress profile: "//&
1983 " A in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1984 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2)
1985 call get_param(param_file, mdl, "TAUX_SIN_AMP", cs%gyres_taux_sin_amp, &
1986 "With the gyres wind_config, the sine amplitude in the "//&
1987 "zonal wind stress profile: "//&
1988 " B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1989 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2)
1990 call get_param(param_file, mdl, "TAUX_COS_AMP", cs%gyres_taux_cos_amp, &
1991 "With the gyres wind_config, the cosine amplitude in "//&
1992 "the zonal wind stress profile: "//&
1993 " C in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1994 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2)
1995 call get_param(param_file, mdl, "TAUX_N_PIS",cs%gyres_taux_n_pis, &
1996 "With the gyres wind_config, the number of gyres in "//&
1997 "the zonal wind stress profile: "//&
1998 " n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1999 units="nondim", default=0.0)
2000 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2001 "This sets the default value for the various _ANSWER_DATE parameters.", &
2002 default=99991231)
2003 call get_param(param_file, mdl, "WIND_GYRES_ANSWER_DATE", cs%answer_date, &
2004 "The vintage of the expressions used to set gyre wind stresses. "//&
2005 "Values below 20190101 recover the answers from the end of 2018, "//&
2006 "while higher values use a form of the gyre wind stresses that are "//&
2007 "rotationally invariant and more likely to be the same between compilers.", &
2008 default=default_answer_date)
2009 else
2010 cs%answer_date = 20190101
2011 endif
2012 if (trim(cs%wind_config) == "scurves") then
2013 call get_param(param_file, mdl, "WIND_SCURVES_LATS", cs%scurves_ydata, &
2014 "A list of latitudes defining a piecewise scurve profile "//&
2015 "for zonal wind stress.", &
2016 units="degrees N", fail_if_missing=.true.)
2017 call get_param(param_file, mdl, "WIND_SCURVES_TAUX", cs%scurves_taux, &
2018 "A list of zonal wind stress values at latitudes "//&
2019 "WIND_SCURVES_LATS defining a piecewise scurve profile.", &
2020 units="Pa", scale=us%Pa_to_RLZ_T2, fail_if_missing=.true.)
2021 endif
2022 if (trim(cs%wind_config) == "2gyre") then
2023 call get_param(param_file, mdl, "TAUX_MAGNITUDE", cs%taux_mag, &
2024 "The peak zonal wind stress when WIND_CONFIG = 2gyre.", &
2025 units="Pa", default=0.1, scale=us%Pa_to_RLZ_T2)
2026 endif
2027 if (trim(cs%wind_config) == "1gyre") then
2028 call get_param(param_file, mdl, "TAUX_MAGNITUDE", cs%taux_mag, &
2029 "The peak zonal wind stress when WIND_CONFIG = 1gyre.", &
2030 units="Pa", default=-0.2, scale=us%Pa_to_RLZ_T2)
2031 endif
2032 if (trim(cs%wind_config) == "Neverworld" .or. trim(cs%wind_config) == "Neverland") then
2033 call get_param(param_file, mdl, "TAUX_MAGNITUDE", cs%taux_mag, &
2034 "The peak zonal wind stress when WIND_CONFIG = Neverworld.", &
2035 units="Pa", default=0.2, scale=us%Pa_to_RLZ_T2)
2036 endif
2037
2038 if ((trim(cs%wind_config) == "2gyre") .or. &
2039 (trim(cs%wind_config) == "1gyre") .or. &
2040 (trim(cs%wind_config) == "gyres") .or. &
2041 (trim(cs%buoy_config) == "linear")) then
2042 cs%south_lat = g%south_lat
2043 cs%len_lat = g%len_lat
2044 endif
2045
2046 call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
2047 "The mean ocean density used with BOUSSINESQ true to "//&
2048 "calculate accelerations and the mass for conservation "//&
2049 "properties, or with BOUSSINESQ false to convert some "//&
2050 "parameters from vertical units of m to kg m-2.", &
2051 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R) ! (, do_not_log=CS%nonBous)
2052 call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
2053 "If true, the buoyancy fluxes drive the model back toward some "//&
2054 "specified surface state with a rate given by FLUXCONST.", default=.false.)
2055 call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
2056 "The latent heat of fusion.", default=hlf, &
2057 units="J/kg", scale=us%J_kg_to_Q)
2058 call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
2059 "The latent heat of fusion.", default=hlv, units="J/kg", scale=us%J_kg_to_Q)
2060 if (cs%restorebuoy) then
2061 ! These three variables use non-standard time units, but are rescaled as they are read.
2062 call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
2063 "The constant that relates the restoring surface fluxes to the relative "//&
2064 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
2065 default=0.0, units="m day-1", scale=us%m_to_Z*us%T_to_s/86400.0)
2066
2067 if (cs%use_temperature) then
2068 call get_param(param_file, mdl, "FLUXCONST", flux_const_default, &
2069 default=0.0, units="m day-1", do_not_log=.true.)
2070 call get_param(param_file, mdl, "FLUXCONST_T", cs%Flux_const_T, &
2071 "The constant that relates the restoring surface temperature flux to the "//&
2072 "relative surface anomaly (akin to a piston velocity). Note the non-MKS units.", &
2073 units="m day-1", scale=us%m_to_Z*us%T_to_s/86400.0, default=flux_const_default)
2074 call get_param(param_file, mdl, "FLUXCONST_S", cs%Flux_const_S, &
2075 "The constant that relates the restoring surface salinity flux to the "//&
2076 "relative surface anomaly (akin to a piston velocity). Note the non-MKS units.", &
2077 units="m day-1", scale=us%m_to_Z*us%T_to_s/86400.0, default=flux_const_default)
2078 endif
2079
2080 if (trim(cs%buoy_config) == "linear") then
2081 call get_param(param_file, mdl, "SST_NORTH", cs%T_north, &
2082 "With buoy_config linear, the sea surface temperature "//&
2083 "at the northern end of the domain toward which to "//&
2084 "to restore.", units="degC", default=0.0, scale=us%degC_to_C)
2085 call get_param(param_file, mdl, "SST_SOUTH", cs%T_south, &
2086 "With buoy_config linear, the sea surface temperature "//&
2087 "at the southern end of the domain toward which to "//&
2088 "to restore.", units="degC", default=0.0, scale=us%degC_to_C)
2089 call get_param(param_file, mdl, "SSS_NORTH", cs%S_north, &
2090 "With buoy_config linear, the sea surface salinity "//&
2091 "at the northern end of the domain toward which to "//&
2092 "to restore.", units="ppt", default=35.0, scale=us%ppt_to_S)
2093 call get_param(param_file, mdl, "SSS_SOUTH", cs%S_south, &
2094 "With buoy_config linear, the sea surface salinity "//&
2095 "at the southern end of the domain toward which to "//&
2096 "to restore.", units="ppt", default=35.0, scale=us%ppt_to_S)
2097 endif
2098 call get_param(param_file, mdl, "RESTORE_FLUX_RHO", cs%rho_restore, &
2099 "The density that is used to convert piston velocities into salt or heat "//&
2100 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
2101 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &
2102 do_not_log=(((cs%Flux_const==0.0).and.(cs%Flux_const_T==0.0).and.(cs%Flux_const_S==0.0))&
2103 .or.(.not.cs%restorebuoy)))
2104 endif
2105 call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
2106 "The gravitational acceleration of the Earth.", &
2107 units="m s-2", default=9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
2108
2109 call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
2110 "The background gustiness in the winds.", &
2111 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2*us%L_to_Z)
2112
2113 call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", cs%ustar_gustless_bug, &
2114 "If true include a bug in the time-averaging of the gustless wind friction velocity", &
2115 default=.false., do_not_log=.true.)
2116 ! This is used to test whether USTAR_GUSTLESS_BUG is being actively set.
2117 call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", test_value, default=.true., do_not_log=.true.)
2118 explicit_bug = cs%ustar_gustless_bug .eqv. test_value
2119 call get_param(param_file, mdl, "FIX_USTAR_GUSTLESS_BUG", fix_ustar_gustless_bug, &
2120 "If true correct a bug in the time-averaging of the gustless wind friction velocity", &
2121 default=.true., do_not_log=.true.)
2122 call get_param(param_file, mdl, "FIX_USTAR_GUSTLESS_BUG", test_value, default=.false., do_not_log=.true.)
2123 explicit_fix = fix_ustar_gustless_bug .eqv. test_value
2124
2125 if (explicit_bug .and. explicit_fix .and. (fix_ustar_gustless_bug .eqv. cs%ustar_gustless_bug)) then
2126 ! USTAR_GUSTLESS_BUG is being explicitly set, and should not be changed.
2127 call mom_error(fatal, "USTAR_GUSTLESS_BUG and FIX_USTAR_GUSTLESS_BUG are both being set "//&
2128 "with inconsistent values. FIX_USTAR_GUSTLESS_BUG is an obsolete "//&
2129 "parameter and should be removed.")
2130 elseif (explicit_fix) then
2131 call mom_error(warning, "FIX_USTAR_GUSTLESS_BUG is an obsolete parameter. "//&
2132 "Use USTAR_GUSTLESS_BUG instead (noting that it has the opposite sense).")
2133 cs%ustar_gustless_bug = .not.fix_ustar_gustless_bug
2134 endif
2135 call log_param(param_file, mdl, "USTAR_GUSTLESS_BUG", cs%ustar_gustless_bug, &
2136 "If true include a bug in the time-averaging of the gustless wind friction velocity", &
2137 default=.false.)
2138
2139 call get_param(param_file, mdl, "READ_GUST_2D", cs%read_gust_2d, &
2140 "If true, use a 2-dimensional gustiness supplied from "//&
2141 "an input file", default=.false.)
2142 if (cs%read_gust_2d) then
2143 call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
2144 "The file in which the wind gustiness is found in "//&
2145 "variable gustiness.", fail_if_missing=.true.)
2146 call safe_alloc_ptr(cs%gust,g%isd,g%ied,g%jsd,g%jed)
2147 filename = trim(cs%inputdir) // trim(gust_file)
2148 ! NOTE: There are certain cases where FMS is unable to read this file, so
2149 ! we use read_netCDF_data in place of MOM_read_data.
2150 call read_netcdf_data(filename, 'gustiness', cs%gust, g%Domain, &
2151 rescale=us%Pa_to_RLZ_T2*us%L_to_Z) ! units in file should be [Pa]
2152 endif
2153 call get_param(param_file, mdl, "USE_MARBL_TRACERS", cs%use_marbl_tracers, &
2154 default=.false., do_not_log=.true.)
2155
2156! All parameter settings are now known.
2157
2158 if (trim(cs%wind_config) == "USER" .or. trim(cs%buoy_config) == "USER" ) then
2159 call user_surface_forcing_init(time, g, us, param_file, diag, cs%user_forcing_CSp)
2160 elseif (trim(cs%buoy_config) == "BFB" ) then
2161 call bfb_surface_forcing_init(time, g, us, param_file, diag, cs%BFB_forcing_CSp)
2162 elseif (trim(cs%buoy_config) == "dumbbell" ) then
2163 call dumbbell_surface_forcing_init(time, g, us, param_file, diag, cs%dumbbell_forcing_CSp)
2164 elseif (trim(cs%wind_config) == "MESO" .or. trim(cs%buoy_config) == "MESO" ) then
2165 call meso_surface_forcing_init(time, g, us, param_file, diag, cs%MESO_forcing_CSp)
2166 elseif (trim(cs%wind_config) == "ideal_hurr") then
2167 call idealized_hurricane_wind_init(time, g, us, param_file, cs%idealized_hurricane_CSp)
2168 elseif (trim(cs%wind_config) == "SCM_ideal_hurr") then
2169 call mom_error(fatal, "MOM_surface_forcing (surface_forcing_init): "//&
2170 'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option. '//&
2171 'To obtain mathematically equivalent results set '//&
2172 'WIND_CONFIG = "ideal_hurr", IDL_HURR_SCM = True and IDL_HURR_X0 = 6.48e+05.')
2173 elseif (trim(cs%wind_config) == "const") then
2174 call get_param(param_file, mdl, "CONST_WIND_TAUX", cs%tau_x0, &
2175 "With wind_config const, this is the constant zonal wind-stress", &
2176 units="Pa", scale=us%Pa_to_RLZ_T2, fail_if_missing=.true.)
2177 call get_param(param_file, mdl, "CONST_WIND_TAUY", cs%tau_y0, &
2178 "With wind_config const, this is the constant meridional wind-stress", &
2179 units="Pa", scale=us%Pa_to_RLZ_T2, fail_if_missing=.true.)
2180 elseif (trim(cs%wind_config) == "SCM_CVmix_tests" .or. &
2181 trim(cs%buoy_config) == "SCM_CVmix_tests") then
2182 call scm_cvmix_tests_surface_forcing_init(time, g, param_file, cs%SCM_CVmix_tests_CSp)
2183 endif
2184
2185 ! Set up MARBL forcing control structure
2186 call marbl_forcing_init(g, us, param_file, diag, time, cs%inputdir, cs%use_marbl_tracers, &
2187 cs%marbl_forcing_CSp)
2188
2189 call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles)
2190
2191 ! Set up any restart fields associated with the forcing.
2192 call restart_init(param_file, cs%restart_CSp, "MOM_forcing.res")
2193!#CTRL# call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, &
2194!#CTRL# CS%restart_CSp)
2195 call restart_init_end(cs%restart_CSp)
2196
2197 if (associated(cs%restart_CSp)) then
2198 call get_mom_input(dirs=dirs)
2199
2200 new_sim = .false.
2201 if ((dirs%input_filename(1:1) == 'n') .and. &
2202 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
2203 if (.not.new_sim) then
2204 call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
2205 g, cs%restart_CSp)
2206 endif
2207 endif
2208
2209 ! Determine how many time levels are in each forcing variable.
2210 if (trim(cs%buoy_config) == "file") then
2211 cs%SW_nlev = num_timelevels(cs%shortwave_file, cs%SW_var, min_dims=3)
2212 cs%LW_nlev = num_timelevels(cs%longwave_file, cs%LW_var, min_dims=3)
2213 cs%latent_nlev = num_timelevels(cs%latentheat_file, cs%latent_var, 3)
2214 cs%sens_nlev = num_timelevels(cs%sensibleheat_file, cs%sens_var, min_dims=3)
2215
2216 cs%evap_nlev = num_timelevels(cs%evaporation_file, cs%evap_var, min_dims=3)
2217 cs%precip_nlev = num_timelevels(cs%rain_file, cs%rain_var, min_dims=3)
2218 cs%runoff_nlev = num_timelevels(cs%runoff_file, cs%lrunoff_var, 3)
2219
2220 cs%SST_nlev = num_timelevels(cs%SSTrestore_file, cs%SST_restore_var, 3)
2221 cs%SSS_nlev = num_timelevels(cs%salinityrestore_file, cs%SSS_restore_var, 3)
2222 endif
2223
2224 if (trim(cs%wind_config) == "file") &
2225 cs%wind_nlev = num_timelevels(cs%wind_file, cs%stress_x_var, min_dims=3)
2226
2227!#CTRL# call controlled_forcing_init(Time, G, US, param_file, diag, CS%ctrl_forcing_CSp)
2228
2229 call user_revise_forcing_init(param_file, cs%urf_CS)
2230
2231 call cpu_clock_end(id_clock_forcing)
2232end subroutine surface_forcing_init
2233
2234
2235!> Deallocate memory associated with the surface forcing module
2236subroutine surface_forcing_end(CS, fluxes)
2237 type(surface_forcing_cs), pointer :: CS !< pointer to control structure returned by
2238 !! a previous surface_forcing_init call
2239 type(forcing), optional, intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
2240
2241 if (present(fluxes)) call deallocate_forcing_type(fluxes)
2242
2243!#CTRL# call controlled_forcing_end(CS%ctrl_forcing_CSp)
2244
2245 if (associated(cs)) deallocate(cs)
2246 cs => null()
2247
2248 call calltree_leave("MARBL_forcing_from_data_override, MOM_surface_forcing.F90")
2249end subroutine surface_forcing_end
2250
2251end module mom_surface_forcing