MOM_surface_forcing_gfdl.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
6
7!#CTRL# use MOM_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts
8!#CTRL# use MOM_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end
9!#CTRL# use MOM_controlled_forcing, only : ctrl_forcing_CS
10use mom_coms, only : reproducing_sum, field_chksum
11use mom_constants, only : hlv, hlf
12use mom_coupler_types, only : coupler_2d_bc_type, coupler_type_write_chksums
13use mom_coupler_types, only : coupler_type_initialized, coupler_type_spawn
14use mom_coupler_types, only : coupler_type_copy_data
15use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
16use mom_cpu_clock, only : clock_subcomponent
17use mom_data_override, only : data_override_init, data_override
18use mom_diag_mediator, only : diag_ctrl, safe_alloc_ptr, time_type
19use mom_domains, only : pass_vector, pass_var, fill_symmetric_edges
20use mom_domains, only : agrid, bgrid_ne, cgrid_ne, to_all
21use mom_domains, only : to_north, to_east, omit_corners
22use mom_eos, only : gsw_sr_from_sp
23use mom_error_handler, only : mom_error, warning, fatal, is_root_pe, mom_mesg
24use mom_file_parser, only : get_param, log_param, log_version, param_file_type
27use mom_forcing_type, only : allocate_forcing_type, deallocate_forcing_type
28use mom_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing
30use mom_grid, only : ocean_grid_type
31use mom_interpolate, only : init_external_field, time_interp_external
32use mom_interpolate, only : time_interp_external_init
33use mom_interpolate, only : external_field
34use mom_io, only : slasher, write_version_number, mom_read_data
35use mom_io, only : read_netcdf_data
36use mom_io, only : stdout_if_root
37use mom_restart, only : register_restart_field, restart_init, mom_restart_cs
42use mom_variables, only : surface
45use iso_fortran_env, only : int64
46
47implicit none ; private
48
49#include <MOM_memory.h>
50
55
56! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
57! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
58! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
59! vary with the Boussinesq approximation, the Boussinesq variant is given first.
60
61!> surface_forcing_CS is a structure containing pointers to the forcing fields
62!! which may be used to drive MOM. All fluxes are positive downward.
63type, public :: surface_forcing_cs ; private
64 integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integer values
65 !! from MOM_domains) to indicate the staggering of
66 !! the winds that are being provided in calls to
67 !! update_ocean_model.
68 logical :: use_temperature !< If true, temp and saln used as state variables.
69 logical :: nonbous !< If true, this run is fully non-Boussinesq
70 real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim].
71
72 real :: rho0 !< Boussinesq reference density [R ~> kg m-3]
73 real :: area_surf = -1.0 !< Total ocean surface area [L2 ~> m2]
74 real :: latent_heat_fusion !< Latent heat of fusion [Q ~> J kg-1]
75 real :: latent_heat_vapor !< Latent heat of vaporization [Q ~> J kg-1]
76
77 real :: max_p_surf !< The maximum surface pressure that can be exerted by
78 !! the atmosphere and floating sea-ice [R L2 T-2 ~> Pa].
79 !! This is needed because the FMS coupling structure
80 !! does not limit the water that can be frozen out
81 !! of the ocean and the ice-ocean heat fluxes are
82 !! treated explicitly.
83 logical :: use_limited_p_ssh !< If true, return the sea surface height with
84 !! the correction for the atmospheric (and sea-ice)
85 !! pressure limited by max_p_surf instead of the
86 !! full atmospheric pressure. The default is true.
87 logical :: approx_net_mass_src !< If true, use the net mass sources from the ice-ocean boundary
88 !! type without any further adjustments to drive the ocean dynamics.
89 !! The actual net mass source may differ due to corrections.
90
91 real :: gust_const !< Constant unresolved background gustiness for ustar [R Z2 T-2 ~> Pa]
92 logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied from an input file.
93 real, pointer, dimension(:,:) :: &
94 bbl_tidal_dis => null() !< Tidal energy dissipation in the bottom boundary layer that can act as a
95 !! source of energy for bottom boundary layer mixing [R Z L2 T-3 ~> W m-2]
96 real, pointer, dimension(:,:) :: &
97 gust => null() !< A spatially varying unresolved background gustiness that
98 !! contributes to ustar [R Z2 T-2 ~> Pa]. gust is used when read_gust_2d is true.
99 real, pointer, dimension(:,:) :: &
100 ustar_tidal => null() !< Tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1]
101 real :: cd_tides !< Drag coefficient that applies to the tides [nondim]
102 real :: utide !< Constant tidal velocity to use if read_tideamp is false [Z T-1 ~> m s-1].
103 logical :: read_tideamp !< If true, spatially varying tidal amplitude read from a file.
104
105 logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts to damp surface
106 !! deflections (especially surface gravity waves). The default is false.
107 real :: g_earth !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
108 real :: kv_sea_ice !< Viscosity in sea-ice that resists sheared vertical motions [L4 Z-2 T-1 ~> m2 s-1]
109 real :: density_sea_ice !< Typical density of sea-ice [R ~> kg m-3]. The value is only used to convert
110 !! the ice pressure into appropriate units for use with Kv_sea_ice.
111 real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which sea-ice viscosity
112 !! becomes effective [R Z ~> kg m-2], typically of order 1000 kg m-2.
113 logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments
114 logical :: allow_carbon_flux_exchange !< If true, allows fluxes and diagnostics of carbon in runoff.
115
116 logical :: restore_salt !< If true, the coupled MOM driver adds a term to restore surface
117 !! salinity to a specified value.
118 logical :: restore_temp !< If true, the coupled MOM driver adds a term to restore sea
119 !! surface temperature to a specified value.
120 real :: flux_const_salt !< Piston velocity for surface salinity restoring [Z T-1 ~> m s-1]
121 real :: flux_const_temp !< Piston velocity for surface temperature restoring [Z T-1 ~> m s-1]
122 real :: rho_restore !< The density that is used to convert piston velocities into salt
123 !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
124 logical :: trestore_spear_ecda !< If true, modify restoring data wrt local SSS
125 real :: spear_dtf_ds !< The derivative of the freezing temperature with
126 !! salinity [C S-1 ~> degC ppt-1].
127 logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux
128 logical :: adjust_net_srestore_to_zero !< Adjust srestore to zero (for both salt_flux or vprec)
129 logical :: adjust_net_srestore_by_scaling !< Adjust srestore w/o moving zero contour
130 logical :: adjust_net_fresh_water_to_zero !< Adjust net surface fresh-water (with restoring) to zero
131 logical :: use_net_fw_adjustment_sign_bug !< Use the wrong sign when adjusting net FW
132 logical :: adjust_net_fresh_water_by_scaling !< Adjust net surface fresh-water w/o moving zero contour
133 logical :: mask_srestore_under_ice !< If true, use an ice mask defined by frazil criteria
134 !! for salinity restoring.
135 real :: ice_salt_concentration !< Salt concentration for sea ice [kg/kg]
136 logical :: mask_srestore_marginal_seas !< If true, then mask SSS restoring in marginal seas
137 logical :: max_delta_srestore_file !< If true, apply a 2-dimensional maximum delta salinity
138 !! when restoring. The file should be
139 !! in inputdir/max_delta_srestore.nc and the field
140 !! should be named 'max_delta_srestore'
141 real, pointer, dimension(:,:) :: max_delta_srestore_2d => null()
142 !< Maximum delta salinity used for restoring [S ~> ppt]
143 real :: max_delta_srestore !< Maximum delta salinity used for restoring [S ~> ppt]
144 real :: min_ratio_srestore !< Minimum fraction of restoring salinity to preserve [nondim]
145 real :: max_delta_trestore !< Maximum delta sst used for restoring [C ~> degC]
146 real, pointer, dimension(:,:) :: basin_mask => null() !< Mask for surface salinity restoring by basin [nondim]
147 integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
148 !! gustiness calculations. Values below 20190101 recover the answers
149 !! from the end of 2018, while higher values use a simpler expression
150 !! to calculate gustiness.
151 logical :: ustar_gustless_bug !< If true, include a bug in the time-averaging of the
152 !! gustless wind friction velocity.
153 logical :: check_no_land_fluxes !< Return warning if IOB flux over land is non-zero
154
155 type(diag_ctrl), pointer :: diag => null() !< Structure to regulate diagnostic output timing
156 character(len=200) :: inputdir !< Directory where NetCDF input files are
157 character(len=200) :: salt_restore_file !< Filename for salt restoring data
158 character(len=30) :: salt_restore_var_name !< Name of surface salinity in salt_restore_file
159 logical :: salt_restore_is_practical !< Specifies that the target salinity is practical and not absolute.
160 logical :: mask_srestore !< If true, apply a 2-dimensional mask to the surface
161 !! salinity restoring fluxes. The masking file should be
162 !! in inputdir/salt_restore_mask.nc and the field should
163 !! be named 'mask'
164 real, pointer, dimension(:,:) :: srestore_mask => null() !< mask for SSS restoring [nondim]
165 character(len=200) :: temp_restore_file !< Filename for sst restoring data
166 character(len=30) :: temp_restore_var_name !< Name of surface temperature in temp_restore_file
167 logical :: mask_trestore !< If true, apply a 2-dimensional mask to the surface
168 !! temperature restoring fluxes. The masking file should be
169 !! in inputdir/temp_restore_mask.nc and the field should
170 !! be named 'mask'
171 real, pointer, dimension(:,:) :: trestore_mask => null() !< Mask for SST restoring [nondim]
172 type(external_field) :: srestore_handle
173 !< Handle for time-interpolated salt restoration field
174 type(external_field) :: trestore_handle
175 !< Handle for time-interpolated temperature restoration field
176
177 type(forcing_diags), public :: handles !< Diagnostics handles
178
179!#CTRL# type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL()
180 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control structure
181 type(user_revise_forcing_cs), pointer :: urf_cs => null() !< A control structure for user forcing revisions
182end type surface_forcing_cs
183
184
185!> ice_ocean_boundary_type is a structure corresponding to forcing, but with the elements, units,
186!! and conventions that exactly conform to the use for MOM6-based coupled models.
187type, public :: ice_ocean_boundary_type
188 real, pointer, dimension(:,:) :: u_flux =>null() !< i-direction wind stress [Pa]
189 real, pointer, dimension(:,:) :: v_flux =>null() !< j-direction wind stress [Pa]
190 real, pointer, dimension(:,:) :: t_flux =>null() !< sensible heat flux [W m-2]
191 real, pointer, dimension(:,:) :: q_flux =>null() !< specific humidity flux [kg m-2 s-1]
192 real, pointer, dimension(:,:) :: salt_flux =>null() !< salt flux [kg m-2 s-1]
193 real, pointer, dimension(:,:) :: excess_salt =>null() !< salt left behind by brine rejection [kg m-2 s-1]
194 real, pointer, dimension(:,:) :: lw_flux =>null() !< long wave radiation [W m-2]
195 real, pointer, dimension(:,:) :: sw_flux_vis_dir =>null() !< direct visible sw radiation [W m-2]
196 real, pointer, dimension(:,:) :: sw_flux_vis_dif =>null() !< diffuse visible sw radiation [W m-2]
197 real, pointer, dimension(:,:) :: sw_flux_nir_dir =>null() !< direct Near InfraRed sw radiation [W m-2]
198 real, pointer, dimension(:,:) :: sw_flux_nir_dif =>null() !< diffuse Near InfraRed sw radiation [W m-2]
199 real, pointer, dimension(:,:) :: lprec =>null() !< mass flux of liquid precip [kg m-2 s-1]
200 real, pointer, dimension(:,:) :: fprec =>null() !< mass flux of frozen precip [kg m-2 s-1]
201 real, pointer, dimension(:,:) :: runoff =>null() !< mass flux of liquid runoff [kg m-2 s-1]
202 real, pointer, dimension(:,:) :: runoff_carbon =>null() !< mass flux of carbon in liquid runoff [kg m-2 s-1]
203 real, pointer, dimension(:,:) :: calving =>null() !< mass flux of frozen runoff [kg m-2 s-1]
204 real, pointer, dimension(:,:) :: stress_mag =>null() !< The time-mean magnitude of the stress on the ocean [Pa]
205 real, pointer, dimension(:,:) :: ustar_berg =>null() !< frictional velocity beneath icebergs [m s-1]
206 real, pointer, dimension(:,:) :: area_berg =>null() !< fractional area covered by icebergs [m2 m-2]
207 real, pointer, dimension(:,:) :: mass_berg =>null() !< mass of icebergs per unit ocean area [kg m-2]
208 real, pointer, dimension(:,:) :: runoff_hflx =>null() !< heat content of liquid runoff [W m-2]
209 real, pointer, dimension(:,:) :: calving_hflx =>null() !< heat content of frozen runoff [W m-2]
210 real, pointer, dimension(:,:) :: p =>null() !< pressure of overlying ice and atmosphere
211 !< on ocean surface [Pa]
212 real, pointer, dimension(:,:) :: mi =>null() !< mass of ice per unit ocean area [kg m-2]
213 real, pointer, dimension(:,:) :: ice_rigidity =>null() !< rigidity of the sea ice, sea-ice and
214 !! ice-shelves, expressed as a coefficient
215 !! for divergence damping, as determined
216 !! outside of the ocean model [m3 s-1]
217 real, pointer, dimension(:,:) :: shelf_sfc_mass_flux =>null() !< mass flux to surface of ice sheet [kg m-2 s-1]
218 integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT
219 type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of named fields
220 !! used for passive tracer fluxes.
221 integer :: wind_stagger = -999 !< A flag indicating the spatial discretization of wind stresses.
222 !! This flag may be set by the flux-exchange code, based on what
223 !! the sea-ice model is providing. Otherwise, the value from
224 !! the surface_forcing_CS is used.
226
227integer :: id_clock_forcing !< A CPU time clock
228
229contains
230
231!> This subroutine translates the Ice_ocean_boundary_type into a MOM
232!! thermodynamic forcing type, including changes of units, sign conventions,
233!! and putting the fields into arrays with MOM-standard halos.
234subroutine convert_iob_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, sfc_state)
236 target, intent(in) :: iob !< An ice-ocean boundary type with fluxes to drive
237 !! the ocean in a coupled model
238 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to all
239 !! possible mass, heat or salt flux forcing fields.
240 !! Unused fields have NULL ptrs.
241 integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
242 type(time_type), intent(in) :: time !< The time of the fluxes, used for interpolating the
243 !! salinity to the right time, when it is being restored.
244 real, intent(in) :: valid_time !< The amount of time over which these fluxes
245 !! should be applied [T ~> s].
246 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
247 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
248 type(surface_forcing_cs),pointer :: cs !< A pointer to the control structure returned by a
249 !! previous call to surface_forcing_init.
250 type(surface), intent(in) :: sfc_state !< A structure containing fields that describe the
251 !! surface state of the ocean.
252
253 real, dimension(SZI_(G),SZJ_(G)) :: &
254 data_restore, & ! The surface value toward which to restore [S ~> ppt] or [C ~> degC]
255 sst_anom, & ! Instantaneous sea surface temperature anomalies from a target value [C ~> degC]
256 sss_anom, & ! Instantaneous sea surface salinity anomalies from a target value [S ~> ppt]
257 sss_mean, & ! A (mean?) salinity about which to normalize local salinity
258 ! anomalies when calculating restorative precipitation anomalies [S ~> ppt]
259 net_fw, & ! The area integrated net freshwater flux into the ocean [R Z L2 T-1 ~> kg s-1]
260 net_fw2, & ! The area averaged net freshwater flux into the ocean [R Z T-1 ~> kg m-2 s-1]
261 work_sum, & ! A 2-d array that is used as the work space for global sums [L2 ~> m2] or [R Z L2 T-1 ~> kg s-1]
262 open_ocn_mask ! a binary field indicating where ice is present based on frazil criteria [nondim]
263
264 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
265 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
266 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
267
268 real :: delta_sss ! temporary storage for sss diff from restoring value [S ~> ppt]
269 real :: delta_sst ! temporary storage for sst diff from restoring value [C ~> degC]
270
271 real :: kg_m2_s_conversion ! A combination of unit conversion factors for rescaling
272 ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1]
273 real :: rhoxcp ! Reference density times heat capacity times unit scaling
274 ! factors [Q R C-1 ~> J m-3 degC-1]
275 real :: sign_for_net_fw_bug ! Should be +1. but an old bug can be recovered by using -1 [nondim]
276
277 call cpu_clock_begin(id_clock_forcing)
278
279 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
280 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
281 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
282 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
283 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
284 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
285 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
286
287 kg_m2_s_conversion = us%kg_m2s_to_RZ_T
288 if (cs%restore_temp) rhoxcp = cs%rho_restore * fluxes%C_p
289 open_ocn_mask(:,:) = 1.0
290 fluxes%vPrecGlobalAdj = 0.0
291 fluxes%vPrecGlobalScl = 0.0
292 fluxes%saltFluxGlobalAdj = 0.0
293 fluxes%saltFluxGlobalScl = 0.0
294 fluxes%netFWGlobalAdj = 0.0
295 fluxes%netFWGlobalScl = 0.0
296
297 ! allocation and initialization if this is the first time that this
298 ! flux type has been used.
299 if (fluxes%dt_buoy_accum < 0) then
300 call allocate_forcing_type(g, fluxes, water=.true., heat=.true., ustar=.not.cs%nonBous, press=.true., &
301 fix_accum_bug=.not.cs%ustar_gustless_bug, tau_mag=cs%nonBous,&
302 carbon=cs%allow_carbon_flux_exchange)
303
304 call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
305 call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
306 call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed)
307 call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed)
308
309 call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed)
310 call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed)
311 if (cs%use_limited_P_SSH) then
312 fluxes%p_surf_SSH => fluxes%p_surf
313 else
314 fluxes%p_surf_SSH => fluxes%p_surf_full
315 endif
316
317 call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed)
318 call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed)
319
320 call safe_alloc_ptr(fluxes%BBL_tidal_dis,isd,ied,jsd,jed)
321 call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed)
322
323 call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
324 call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
325
326 if (associated(iob%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed)
327
328 do j=js-2,je+2 ; do i=is-2,ie+2
329 fluxes%BBL_tidal_dis(i,j) = cs%BBL_tidal_dis(i,j)
330 fluxes%ustar_tidal(i,j) = cs%ustar_tidal(i,j)
331 enddo ; enddo
332
333
334 endif ! endif for allocation and initialization
335
336
337 if (((associated(iob%ustar_berg) .and. (.not.associated(fluxes%ustar_berg))) &
338 .or. (associated(iob%area_berg) .and. (.not.associated(fluxes%area_berg)))) &
339 .or. (associated(iob%mass_berg) .and. (.not.associated(fluxes%mass_berg)))) &
340 call allocate_forcing_type(g, fluxes, iceberg=.true.)
341
342 if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. &
343 coupler_type_initialized(iob%fluxes)) &
344 call coupler_type_spawn(iob%fluxes, fluxes%tr_fluxes, (/is,is,ie,ie/), (/js,js,je,je/))
345 ! It might prove valuable to use the same array extents as the rest of the
346 ! ocean model, rather than using haloless arrays, in which case the last line
347 ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/))
348
349 ! allocation and initialization on first call to this routine
350 if (cs%area_surf < 0.0) then
351 do j=js,je ; do i=is,ie
352 work_sum(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
353 enddo ; enddo
354 cs%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer, unscale=us%L_to_m**2)
355 endif ! endif for allocation and initialization
356
357
358 ! Indicate that there are new unused fluxes.
359 fluxes%fluxes_used = .false.
360 fluxes%dt_buoy_accum = valid_time
361
362 fluxes%heat_added(:,:) = 0.0
363 fluxes%salt_flux_added(:,:) = 0.0
364
365 do j=js,je ; do i=is,ie
366 fluxes%salt_flux(i,j) = 0.0
367 fluxes%vprec(i,j) = 0.0
368 enddo ; enddo
369
370 ! Salinity restoring logic
371 if (cs%restore_salt) then
372 call time_interp_external(cs%srestore_handle, time, data_restore, scale=us%ppt_to_S)
373 if (sfc_state%S_is_absS .and. cs%salt_restore_is_practical) then
374 !Adjust the salt restoring data to absolute
375 do j=js,je
376 do i=is,ie
377 data_restore(i,j) = gsw_sr_from_sp(data_restore(i,j))
378 enddo
379 enddo
380 endif
381 ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not)
382 open_ocn_mask(:,:) = 1.0
383 if (cs%mask_srestore_under_ice) then ! Do not restore under sea-ice
384 do j=js,je ; do i=is,ie
385 if (sfc_state%SST(i,j) <= cs%SPEAR_dTf_dS*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
386 enddo ; enddo
387 endif
388 if (cs%salt_restore_as_sflux) then
389 do j=js,je ; do i=is,ie
390 delta_sss = data_restore(i,j) - sfc_state%SSS(i,j)
391 if (sfc_state%SSS(i,j) >= data_restore(i,j)*cs%min_ratio_srestore) then
392 if (.not. cs%max_delta_srestore_file) then
393 delta_sss = sign(1.0,delta_sss) * min(abs(delta_sss), cs%max_delta_srestore)
394 else
395 if (abs(delta_sss) > abs(cs%max_delta_srestore_2d(i,j))) then
396 if (cs%max_delta_srestore_2d(i,j) < 0.0) then
397 delta_sss = 0.0 !turn off restoring
398 else !clip restoring
399 delta_sss = sign(1.0,delta_sss) * min(abs(delta_sss), cs%max_delta_srestore_2d(i,j))
400 endif
401 endif
402 endif !max_delta_srestore_file
403 endif !min_ratio_srestore
404 fluxes%salt_flux(i,j) = 1.e-3*us%S_to_ppt*g%mask2dT(i,j) * (cs%rho_restore*cs%Flux_const_salt)* &
405 (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j)) * delta_sss ! R Z T-1 ~> kg Salt m-2 s-1
406 enddo ; enddo
407 if (cs%adjust_net_srestore_to_zero) then
408 if (cs%adjust_net_srestore_by_scaling) then
409 call adjust_area_mean_to_zero(fluxes%salt_flux, g, fluxes%saltFluxGlobalScl, &
410 unit_scale=us%RZ_T_to_kg_m2s)
411 fluxes%saltFluxGlobalAdj = 0.
412 else
413 work_sum(is:ie,js:je) = g%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je) * g%mask2dT(is:ie,js:je)
414 fluxes%saltFluxGlobalAdj = reproducing_sum(work_sum(:,:), isr,ier, jsr,jer, unscale=us%RZL2_to_kg*us%s_to_T) &
415 / cs%area_surf
416 fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - &
417 fluxes%saltFluxGlobalAdj * g%mask2dT(is:ie,js:je)
418 endif
419 endif
420 fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic
421 else
422 do j=js,je ; do i=is,ie
423 if (g%mask2dT(i,j) > 0.0) then
424 delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
425 if (sfc_state%SSS(i,j) >= data_restore(i,j)*cs%min_ratio_srestore) then
426 if (.not. cs%max_delta_srestore_file) then
427 delta_sss = sign(1.0,delta_sss) * min(abs(delta_sss), cs%max_delta_srestore)
428 else
429 if (abs(delta_sss) > abs(cs%max_delta_srestore_2d(i,j))) then
430 if (cs%max_delta_srestore_2d(i,j) < 0.0) then
431 delta_sss = 0.0 !turn off restoring
432 else !clip restoring
433 delta_sss = sign(1.0,delta_sss) * min(abs(delta_sss), cs%max_delta_srestore_2d(i,j))
434 endif
435 endif
436 endif !max_delta_srestore_file
437 endif !min_ratio_srestore
438 fluxes%vprec(i,j) = (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j))* &
439 (cs%rho_restore*cs%Flux_const_salt) * &
440 delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j)))
441 endif
442 enddo ; enddo
443 if (cs%adjust_net_srestore_to_zero) then
444 if (cs%adjust_net_srestore_by_scaling) then
445 call adjust_area_mean_to_zero(fluxes%vprec, g, fluxes%vPrecGlobalScl, &
446 unit_scale=us%RZ_T_to_kg_m2s)
447 fluxes%vPrecGlobalAdj = 0.
448 else
449 work_sum(is:ie,js:je) = g%areaT(is:ie,js:je) * fluxes%vprec(is:ie,js:je)
450 fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer, unscale=us%RZL2_to_kg*us%s_to_T) &
451 / cs%area_surf
452 do j=js,je ; do i=is,ie
453 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%vPrecGlobalAdj ) * g%mask2dT(i,j)
454 enddo ; enddo
455 endif
456 endif
457 endif
458 endif
459
460 ! SST restoring logic
461 if (cs%restore_temp) then
462 call time_interp_external(cs%trestore_handle, time, data_restore, scale=us%degC_to_C)
463 if ( cs%trestore_SPEAR_ECDA ) then
464 do j=js,je ; do i=is,ie
465 if (abs(data_restore(i,j)+1.8*us%degC_to_C) < 0.0001*us%degC_to_C) then
466 data_restore(i,j) = cs%SPEAR_dTf_dS*sfc_state%SSS(i,j)
467 endif
468 enddo ; enddo
469 endif
470
471 do j=js,je ; do i=is,ie
472 delta_sst = data_restore(i,j) - sfc_state%SST(i,j)
473 delta_sst = sign(1.0,delta_sst) * min(abs(delta_sst), cs%max_delta_trestore)
474 fluxes%heat_added(i,j) = g%mask2dT(i,j) * cs%trestore_mask(i,j) * &
475 rhoxcp * delta_sst * cs%Flux_const_temp ! [Q R Z T-1 ~> W m-2]
476 enddo ; enddo
477 endif
478
479
480 ! obtain fluxes from IOB; note the staggering of indices
481 i0 = is - isc_bnd ; j0 = js - jsc_bnd
482 do j=js,je ; do i=is,ie
483
484 if (associated(iob%lprec)) then
485 fluxes%lprec(i,j) = kg_m2_s_conversion * iob%lprec(i-i0,j-j0) * g%mask2dT(i,j)
486 if (cs%check_no_land_fluxes) &
487 call check_mask_val_consistency(iob%lprec(i-i0,j-j0), g%mask2dT(i,j), i, j, 'lprec', g)
488 endif
489
490 if (associated(iob%fprec)) then
491 fluxes%fprec(i,j) = kg_m2_s_conversion * iob%fprec(i-i0,j-j0) * g%mask2dT(i,j)
492 if (cs%check_no_land_fluxes) &
493 call check_mask_val_consistency(iob%fprec(i-i0,j-j0), g%mask2dT(i,j), i, j, 'fprec', g)
494 endif
495
496 if (associated(iob%q_flux)) then
497 fluxes%evap(i,j) = - kg_m2_s_conversion * iob%q_flux(i-i0,j-j0) * g%mask2dT(i,j)
498 if (cs%check_no_land_fluxes) &
499 call check_mask_val_consistency(iob%q_flux(i-i0,j-j0), g%mask2dT(i,j), i, j, 'q_flux', g)
500 endif
501
502 if (associated(iob%runoff)) then
503 fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%runoff(i-i0,j-j0) * g%mask2dT(i,j)
504 if (cs%check_no_land_fluxes) &
505 call check_mask_val_consistency(iob%runoff(i-i0,j-j0), g%mask2dT(i,j), i, j, 'runoff', g)
506 endif
507
508 if (associated(iob%calving)) then
509 fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%calving(i-i0,j-j0) * g%mask2dT(i,j)
510 if (cs%check_no_land_fluxes) &
511 call check_mask_val_consistency(iob%calving(i-i0,j-j0), g%mask2dT(i,j), i, j, 'calving', g)
512 endif
513
514 if (associated(iob%shelf_sfc_mass_flux)) then
515 fluxes%shelf_sfc_mass_flux(i,j) = kg_m2_s_conversion * iob%shelf_sfc_mass_flux(i-i0,j-j0)
516 endif
517
518 if (associated(iob%ustar_berg)) then
519 fluxes%ustar_berg(i,j) = us%m_to_Z*us%T_to_s * iob%ustar_berg(i-i0,j-j0) * g%mask2dT(i,j)
520 if (cs%check_no_land_fluxes) &
521 call check_mask_val_consistency(iob%ustar_berg(i-i0,j-j0), g%mask2dT(i,j), i, j, 'ustar_berg', g)
522 endif
523
524 if (associated(iob%area_berg)) then
525 fluxes%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
526 if (cs%check_no_land_fluxes) &
527 call check_mask_val_consistency(iob%area_berg(i-i0,j-j0), g%mask2dT(i,j), i, j, 'area_berg', g)
528 endif
529
530 if (associated(iob%mass_berg)) then
531 fluxes%mass_berg(i,j) = us%m_to_Z*us%kg_m3_to_R * iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
532 if (cs%check_no_land_fluxes) &
533 call check_mask_val_consistency(iob%mass_berg(i-i0,j-j0), g%mask2dT(i,j), i, j, 'mass_berg', g)
534 endif
535
536 if (associated(iob%runoff_hflx)) then
537 fluxes%heat_content_lrunoff(i,j) = us%W_m2_to_QRZ_T * iob%runoff_hflx(i-i0,j-j0) * g%mask2dT(i,j)
538 if (cs%check_no_land_fluxes) &
539 call check_mask_val_consistency(iob%runoff_hflx(i-i0,j-j0), g%mask2dT(i,j), i, j, 'runoff_hflx', g)
540 endif
541
542 if (associated(iob%runoff_carbon) .and. cs%allow_carbon_flux_exchange) then
543 fluxes%carbon_content_lrunoff(i,j) = us%kg_m2s_to_RZ_T * iob%runoff_carbon(i-i0,j-j0) * g%mask2dT(i,j)
544 if (cs%check_no_land_fluxes) &
545 call check_mask_val_consistency(iob%runoff_carbon(i-i0,j-j0), g%mask2dT(i,j), i, j, 'runoff_carbon', g)
546 endif
547
548 if (associated(iob%calving_hflx)) then
549 fluxes%heat_content_frunoff(i,j) = us%W_m2_to_QRZ_T * iob%calving_hflx(i-i0,j-j0) * g%mask2dT(i,j)
550 if (cs%check_no_land_fluxes) &
551 call check_mask_val_consistency(iob%calving_hflx(i-i0,j-j0), g%mask2dT(i,j), i, j, 'calving_hflx', g)
552 endif
553
554 if (associated(iob%lw_flux)) then
555 fluxes%LW(i,j) = us%W_m2_to_QRZ_T * iob%lw_flux(i-i0,j-j0) * g%mask2dT(i,j)
556 if (cs%check_no_land_fluxes) &
557 call check_mask_val_consistency(iob%lw_flux(i-i0,j-j0), g%mask2dT(i,j), i, j, 'lw_flux', g)
558 endif
559
560 if (associated(iob%t_flux)) then
561 fluxes%sens(i,j) = -us%W_m2_to_QRZ_T* iob%t_flux(i-i0,j-j0) * g%mask2dT(i,j)
562 if (cs%check_no_land_fluxes) &
563 call check_mask_val_consistency(iob%t_flux(i-i0,j-j0), g%mask2dT(i,j), i, j, 't_flux', g)
564 endif
565
566 fluxes%latent(i,j) = 0.0
567 if (associated(iob%fprec)) then
568 fluxes%latent(i,j) = fluxes%latent(i,j) - iob%fprec(i-i0,j-j0)*kg_m2_s_conversion * cs%latent_heat_fusion
569 fluxes%latent_fprec_diag(i,j) = -g%mask2dT(i,j) * iob%fprec(i-i0,j-j0)*kg_m2_s_conversion * cs%latent_heat_fusion
570 endif
571 if (associated(iob%calving)) then
572 fluxes%latent(i,j) = fluxes%latent(i,j) - iob%calving(i-i0,j-j0)*kg_m2_s_conversion * cs%latent_heat_fusion
573 fluxes%latent_frunoff_diag(i,j) = -g%mask2dT(i,j) * iob%calving(i-i0,j-j0)*kg_m2_s_conversion * &
574 cs%latent_heat_fusion
575 endif
576 if (associated(iob%q_flux)) then
577 fluxes%latent(i,j) = fluxes%latent(i,j) - iob%q_flux(i-i0,j-j0)*kg_m2_s_conversion * cs%latent_heat_vapor
578 fluxes%latent_evap_diag(i,j) = -g%mask2dT(i,j) * iob%q_flux(i-i0,j-j0)*kg_m2_s_conversion * cs%latent_heat_vapor
579 endif
580
581 fluxes%latent(i,j) = g%mask2dT(i,j) * fluxes%latent(i,j)
582
583 if (associated(iob%sw_flux_vis_dir)) then
584 fluxes%sw_vis_dir(i,j) = g%mask2dT(i,j) * us%W_m2_to_QRZ_T * iob%sw_flux_vis_dir(i-i0,j-j0)
585 if (cs%check_no_land_fluxes) &
586 call check_mask_val_consistency(iob%sw_flux_vis_dir(i-i0,j-j0), g%mask2dT(i,j), i, j, 'sw_flux_vis_dir', g)
587 endif
588 if (associated(iob%sw_flux_vis_dif)) then
589 fluxes%sw_vis_dif(i,j) = g%mask2dT(i,j) * us%W_m2_to_QRZ_T * iob%sw_flux_vis_dif(i-i0,j-j0)
590 if (cs%check_no_land_fluxes) &
591 call check_mask_val_consistency(iob%sw_flux_vis_dif(i-i0,j-j0), g%mask2dT(i,j), i, j, 'sw_flux_vis_dif', g)
592 endif
593 if (associated(iob%sw_flux_nir_dir)) then
594 fluxes%sw_nir_dir(i,j) = g%mask2dT(i,j) * us%W_m2_to_QRZ_T * iob%sw_flux_nir_dir(i-i0,j-j0)
595 if (cs%check_no_land_fluxes) &
596 call check_mask_val_consistency(iob%sw_flux_nir_dir(i-i0,j-j0), g%mask2dT(i,j), i, j, 'sw_flux_nir_dir', g)
597 endif
598 if (associated(iob%sw_flux_nir_dif)) then
599 fluxes%sw_nir_dif(i,j) = g%mask2dT(i,j) * us%W_m2_to_QRZ_T * iob%sw_flux_nir_dif(i-i0,j-j0)
600 if (cs%check_no_land_fluxes) &
601 call check_mask_val_consistency(iob%sw_flux_nir_dif(i-i0,j-j0), g%mask2dT(i,j), i, j, 'sw_flux_nir_dif', g)
602 endif
603 if (cs%answer_date < 20190101) then
604 fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
605 fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)
606 else
607 fluxes%sw(i,j) = (fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j)) + &
608 (fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j))
609 endif
610
611 enddo ; enddo
612
613 ! applied surface pressure from atmosphere and cryosphere
614 if (associated(iob%p)) then
615 if (cs%max_p_surf >= 0.0) then
616 do j=js,je ; do i=is,ie
617 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * us%Pa_to_RL2_T2*iob%p(i-i0,j-j0)
618 fluxes%p_surf(i,j) = min(fluxes%p_surf_full(i,j),cs%max_p_surf)
619 if (cs%check_no_land_fluxes) &
620 call check_mask_val_consistency(iob%p(i-i0,j-j0), g%mask2dT(i,j), i, j, 'p', g)
621 enddo ; enddo
622 else
623 do j=js,je ; do i=is,ie
624 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * us%Pa_to_RL2_T2*iob%p(i-i0,j-j0)
625 fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
626 if (cs%check_no_land_fluxes) &
627 call check_mask_val_consistency(iob%p(i-i0,j-j0), g%mask2dT(i,j), i, j, 'p', g)
628 enddo ; enddo
629 endif
630 fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
631 endif
632
633 ! more salt restoring logic
634 if (associated(iob%salt_flux)) then
635 do j=js,je ; do i=is,ie
636 fluxes%salt_flux(i,j) = g%mask2dT(i,j)*(fluxes%salt_flux(i,j) - kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0))
637 fluxes%salt_flux_in(i,j) = g%mask2dT(i,j)*( -kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0) )
638 if (cs%check_no_land_fluxes) &
639 call check_mask_val_consistency(iob%salt_flux(i-i0,j-j0), g%mask2dT(i,j), i, j, 'salt_flux', g)
640 enddo ; enddo
641 endif
642 if (associated(iob%excess_salt)) then
643 do j=js,je ; do i=is,ie
644 fluxes%salt_left_behind(i,j) = g%mask2dT(i,j)*(kg_m2_s_conversion*iob%excess_salt(i-i0,j-j0))
645 enddo ; enddo
646 endif
647
648!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
649!#CTRL# do j=js,je ; do i=is,ie
650!#CTRL# SST_anom(i,j) = sfc_state%SST(i,j) - CS%T_Restore(i,j)
651!#CTRL# SSS_anom(i,j) = sfc_state%SSS(i,j) - CS%S_Restore(i,j)
652!#CTRL# SSS_mean(i,j) = 0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j))
653!#CTRL# enddo ; enddo
654!#CTRL# call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
655!#CTRL# fluxes%vprec, day, valid_time, G, US, CS%ctrl_forcing_CSp)
656!#CTRL# endif
657
658 ! adjust the NET fresh-water flux to zero, if flagged
659 if (cs%adjust_net_fresh_water_to_zero) then
660 sign_for_net_fw_bug = 1.
661 if (cs%use_net_FW_adjustment_sign_bug) sign_for_net_fw_bug = -1.
662 do j=js,je ; do i=is,ie
663 net_fw(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + &
664 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
665 (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * g%areaT(i,j)
666 ! The following contribution appears to be calculating the volume flux of sea-ice
667 ! melt. This calculation is clearly WRONG if either sea-ice has variable
668 ! salinity or the sea-ice is completely fresh.
669 ! Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system
670 ! is constant.
671 ! To do this correctly we will need a sea-ice melt field added to IOB. -AJA
672 if (associated(iob%salt_flux) .and. (cs%ice_salt_concentration>0.0)) &
673 net_fw(i,j) = net_fw(i,j) + sign_for_net_fw_bug * g%areaT(i,j) * &
674 (kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0) / cs%ice_salt_concentration)
675 net_fw2(i,j) = net_fw(i,j) / g%areaT(i,j)
676 enddo ; enddo
677
678 if (cs%adjust_net_fresh_water_by_scaling) then
679 call adjust_area_mean_to_zero(net_fw2, g, fluxes%netFWGlobalScl, unscale=us%RZ_T_to_kg_m2s)
680 do j=js,je ; do i=is,ie
681 fluxes%vprec(i,j) = fluxes%vprec(i,j) + &
682 (net_fw2(i,j) - net_fw(i,j)/g%areaT(i,j)) * g%mask2dT(i,j)
683 enddo ; enddo
684 else
685 fluxes%netFWGlobalAdj = reproducing_sum(net_fw(:,:), isr, ier, jsr, jer, unscale=us%RZL2_to_kg*us%s_to_T) / &
686 cs%area_surf
687 do j=js,je ; do i=is,ie
688 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%netFWGlobalAdj ) * g%mask2dT(i,j)
689 enddo ; enddo
690 endif
691
692 endif
693
694 ! Set the wind stresses and ustar.
695 if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) .and. associated(fluxes%tau_mag) &
696 .and. associated(fluxes%tau_mag_gustless) ) then
697 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, ustar=fluxes%ustar, &
698 mag_tau=fluxes%tau_mag, gustless_ustar=fluxes%ustar_gustless, &
699 gustless_mag_tau=fluxes%tau_mag_gustless)
700 else
701 if (associated(fluxes%ustar)) &
702 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, ustar=fluxes%ustar)
703 if (associated(fluxes%ustar_gustless)) &
704 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, gustless_ustar=fluxes%ustar_gustless)
705 if (associated(fluxes%tau_mag)) &
706 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, mag_tau=fluxes%tau_mag)
707 if (associated(fluxes%tau_mag_gustless)) &
708 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, gustless_mag_tau=fluxes%tau_mag_gustless)
709 endif
710
711 if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
712 coupler_type_initialized(iob%fluxes)) &
713 call coupler_type_copy_data(iob%fluxes, fluxes%tr_fluxes)
714
715 if (cs%allow_flux_adjustments) then
716 ! Apply adjustments to fluxes
717 call apply_flux_adjustments(g, us, cs, time, fluxes)
718 endif
719
720 ! Allow for user-written code to alter fluxes after all the above
721 call user_alter_forcing(sfc_state, fluxes, time, g, cs%urf_CS)
722
723 call cpu_clock_end(id_clock_forcing)
724
725end subroutine convert_iob_to_fluxes
726
727!> This subroutine translates the Ice_ocean_boundary_type into a MOM
728!! mechanical forcing type, including changes of units, sign conventions,
729!! and putting the fields into arrays with MOM-standard halos.
730subroutine convert_iob_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_forcing, reset_avg)
732 target, intent(in) :: iob !< An ice-ocean boundary type with fluxes to drive
733 !! the ocean in a coupled model
734 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
735 integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
736 type(time_type), intent(in) :: time !< The time of the fluxes, used for interpolating the
737 !! salinity to the right time, when it is being restored.
738 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
739 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
740 type(surface_forcing_cs),pointer :: cs !< A pointer to the control structure returned by a
741 !! previous call to surface_forcing_init.
742 real, optional, intent(in) :: dt_forcing !< A time interval over which to apply the
743 !! current value of ustar as a weighted running
744 !! average [T ~> s], or if 0 do not average ustar.
745 !! Missing is equivalent to 0.
746 logical, optional, intent(in) :: reset_avg !< If true, reset the time average.
747
748 ! Local variables
749 real, dimension(SZI_(G),SZJ_(G)) :: &
750 rigidity_at_h, & ! Ice rigidity at tracer points [L4 Z-1 T-1 ~> m3 s-1]
751 net_mass_src, & ! A temporary of net mass sources [R Z T-1 ~> kg m-2 s-1].
752 ustar_tmp, & ! A temporary array of ustar values [Z T-1 ~> m s-1].
753 tau_mag_tmp ! A temporary array of surface stress magnitudes [R Z2 T-2 ~> Pa]
754
755 real :: i_gearth ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
756 real :: kv_rho_ice ! (CS%Kv_sea_ice / CS%density_sea_ice) [L4 Z-2 T-1 R-1 ~> m5 s-1 kg-1]
757 real :: mass_ice ! mass of sea ice at a face [R Z ~> kg m-2]
758 real :: mass_eff ! effective mass of sea ice for rigidity [R Z ~> kg m-2]
759 real :: wt1, wt2 ! Relative weights of previous and current values of ustar [nondim].
760 real :: kg_m2_s_conversion ! A combination of unit conversion factors for rescaling
761 ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1]
762
763 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
764 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
765 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
766
767 call cpu_clock_begin(id_clock_forcing)
768
769 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
770 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
771 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
772 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
773 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
774 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
775 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
776 i0 = is - isc_bnd ; j0 = js - jsc_bnd
777
778 kg_m2_s_conversion = us%kg_m2s_to_RZ_T
779
780 ! allocation and initialization if this is the first time that this
781 ! mechanical forcing type has been used.
782 if (.not.forces%initialized) then
783 call allocate_mech_forcing(g, forces, stress=.true., ustar=.not.cs%nonBous, &
784 press=.true., tau_mag=cs%nonBous)
785
786 call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
787 call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
788 if (cs%use_limited_P_SSH) then
789 forces%p_surf_SSH => forces%p_surf
790 else
791 forces%p_surf_SSH => forces%p_surf_full
792 endif
793
794 if (cs%rigid_sea_ice) then
795 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
796 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
797 endif
798
799 forces%initialized = .true.
800 endif
801
802 if ( (associated(iob%area_berg) .and. (.not. associated(forces%area_berg))) .or. &
803 (associated(iob%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
804 call allocate_mech_forcing(g, forces, iceberg=.true.)
805
806 if (associated(iob%ice_rigidity)) then
807 rigidity_at_h(:,:) = 0.0
808 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
809 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
810 endif
811
812 forces%accumulate_rigidity = .true. ! Multiple components may contribute to rigidity.
813 if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
814 if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
815
816 ! Set the weights for forcing fields that use running time averages.
817 if (present(reset_avg)) then ; if (reset_avg) forces%dt_force_accum = 0.0 ; endif
818 wt1 = 0.0 ; wt2 = 1.0
819 if (present(dt_forcing)) then
820 if ((forces%dt_force_accum > 0.0) .and. (dt_forcing > 0.0)) then
821 wt1 = forces%dt_force_accum / (forces%dt_force_accum + dt_forcing)
822 wt2 = 1.0 - wt1
823 endif
824 if (dt_forcing > 0.0) then
825 forces%dt_force_accum = max(forces%dt_force_accum, 0.0) + dt_forcing
826 else
827 forces%dt_force_accum = 0.0 ! Reset the averaging time interval.
828 endif
829 else
830 forces%dt_force_accum = 0.0 ! Reset the averaging time interval.
831 endif
832
833 ! applied surface pressure from atmosphere and cryosphere
834 if (associated(iob%p)) then
835 if (cs%max_p_surf >= 0.0) then
836 do j=js,je ; do i=is,ie
837 forces%p_surf_full(i,j) = g%mask2dT(i,j) * us%Pa_to_RL2_T2*iob%p(i-i0,j-j0)
838 forces%p_surf(i,j) = min(forces%p_surf_full(i,j),cs%max_p_surf)
839 enddo ; enddo
840 else
841 do j=js,je ; do i=is,ie
842 forces%p_surf_full(i,j) = g%mask2dT(i,j) * us%Pa_to_RL2_T2*iob%p(i-i0,j-j0)
843 forces%p_surf(i,j) = forces%p_surf_full(i,j)
844 enddo ; enddo
845 endif
846 else
847 do j=js,je ; do i=is,ie
848 forces%p_surf_full(i,j) = 0.0
849 forces%p_surf(i,j) = 0.0
850 enddo ; enddo
851 endif
852 forces%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
853
854 ! Set the wind stresses and ustar.
855 if (wt1 <= 0.0) then
856 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, taux=forces%taux, tauy=forces%tauy, &
857 tau_halo=1)
858 if (associated(forces%ustar)) &
859 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, ustar=forces%ustar)
860 if (associated(forces%tau_mag)) &
861 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, mag_tau=forces%tau_mag)
862 else
863 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, taux=forces%taux, tauy=forces%tauy, &
864 tau_halo=1)
865 if (associated(forces%ustar)) then
866 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, ustar=ustar_tmp)
867 do j=js,je ; do i=is,ie
868 forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
869 enddo ; enddo
870 endif
871 if (associated(forces%tau_mag)) then
872 call extract_iob_stresses(iob, index_bounds, time, g, us, cs, mag_tau=tau_mag_tmp)
873 do j=js,je ; do i=is,ie
874 forces%tau_mag(i,j) = wt1*forces%tau_mag(i,j) + wt2*tau_mag_tmp(i,j)
875 enddo ; enddo
876 endif
877 endif
878
879 ! Find the net mass source in the input forcing without other adjustments.
880 if (cs%approx_net_mass_src .and. associated(forces%net_mass_src)) then
881 net_mass_src(:,:) = 0.0
882 i0 = is - isc_bnd ; j0 = js - jsc_bnd
883 do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
884 if (associated(iob%lprec)) &
885 net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * iob%lprec(i-i0,j-j0)
886 if (associated(iob%fprec)) &
887 net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * iob%fprec(i-i0,j-j0)
888 if (associated(iob%runoff)) &
889 net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * iob%runoff(i-i0,j-j0)
890 if (associated(iob%calving)) &
891 net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * iob%calving(i-i0,j-j0)
892 if (associated(iob%q_flux)) &
893 net_mass_src(i,j) = net_mass_src(i,j) - kg_m2_s_conversion * iob%q_flux(i-i0,j-j0)
894 endif ; enddo ; enddo
895 if (wt1 <= 0.0) then
896 do j=js,je ; do i=is,ie
897 forces%net_mass_src(i,j) = wt2*net_mass_src(i,j)
898 enddo ; enddo
899 else
900 do j=js,je ; do i=is,ie
901 forces%net_mass_src(i,j) = wt1*forces%net_mass_src(i,j) + wt2*net_mass_src(i,j)
902 enddo ; enddo
903 endif
904 forces%net_mass_src_set = .true.
905 else
906 forces%net_mass_src_set = .false.
907 endif
908
909 ! Obtain optional ice-berg related fluxes from the IOB type:
910 if (associated(iob%area_berg)) then ; do j=js,je ; do i=is,ie
911 forces%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
912 enddo ; enddo ; endif
913
914 if (associated(iob%mass_berg)) then ; do j=js,je ; do i=is,ie
915 forces%mass_berg(i,j) = us%m_to_Z*us%kg_m3_to_R * iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
916 enddo ; enddo ; endif
917
918 ! Obtain sea ice related dynamic fields
919 if (associated(iob%ice_rigidity)) then
920 do j=js,je ; do i=is,ie
921 rigidity_at_h(i,j) = us%m_to_L**3*us%Z_to_L*us%T_to_s * iob%ice_rigidity(i-i0,j-j0) * g%mask2dT(i,j)
922 enddo ; enddo
923 call pass_var(rigidity_at_h, g%Domain, halo=1)
924 do i=is-1,ie ; do j=js,je
925 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
926 min(rigidity_at_h(i,j), rigidity_at_h(i+1,j))
927 enddo ; enddo
928 do i=is,ie ; do j=js-1,je
929 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
930 min(rigidity_at_h(i,j), rigidity_at_h(i,j+1))
931 enddo ; enddo
932 endif
933
934 if (cs%rigid_sea_ice) then
935 call pass_var(forces%p_surf_full, g%Domain, halo=1)
936 i_gearth = 1.0 / cs%g_Earth
937 kv_rho_ice = (cs%Kv_sea_ice / cs%density_sea_ice)
938 do i=is-1,ie ; do j=js,je
939 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * i_gearth
940 mass_eff = 0.0
941 if (mass_ice > cs%rigid_sea_ice_mass) then
942 mass_eff = (mass_ice - cs%rigid_sea_ice_mass)**2 / (mass_ice + cs%rigid_sea_ice_mass)
943 endif
944 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * mass_eff
945 enddo ; enddo
946 do i=is,ie ; do j=js-1,je
947 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * i_gearth
948 mass_eff = 0.0
949 if (mass_ice > cs%rigid_sea_ice_mass) then
950 mass_eff = (mass_ice - cs%rigid_sea_ice_mass)**2 / (mass_ice + cs%rigid_sea_ice_mass)
951 endif
952 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * mass_eff
953 enddo ; enddo
954 endif
955
956 if (cs%allow_flux_adjustments) then
957 ! Apply adjustments to forces
958 call apply_force_adjustments(g, us, cs, time, forces)
959 endif
960
961!### ! Allow for user-written code to alter fluxes after all the above
962!### call user_alter_mech_forcing(forces, Time, G, CS%urf_CS)
963
964 call cpu_clock_end(id_clock_forcing)
965end subroutine convert_iob_to_forces
966
967
968!> This subroutine extracts the wind stresses and related fields like ustar from an
969!! Ice_ocean_boundary_type into optional argument arrays, including changes of units, sign
970!! conventions, and putting the fields into arrays with MOM-standard sized halos.
971subroutine extract_iob_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ustar, &
972 gustless_ustar, mag_tau, gustless_mag_tau, tau_halo)
974 target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive
975 !! the ocean in a coupled model
976 integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
977 type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the
978 !! salinity to the right time, when it is being restored.
979 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
980 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
981 type(surface_forcing_cs),pointer :: CS !< A pointer to the control structure returned by a
982 !! previous call to surface_forcing_init.
983 real, dimension(SZIB_(G),SZJ_(G)), &
984 optional, intent(inout) :: taux !< The zonal wind stresses on a C-grid [R Z L T-2 ~> Pa].
985 real, dimension(SZI_(G),SZJB_(G)), &
986 optional, intent(inout) :: tauy !< The meridional wind stresses on a C-grid [R Z L T-2 ~> Pa].
987 real, dimension(SZI_(G),SZJ_(G)), &
988 optional, intent(inout) :: ustar !< The surface friction velocity [Z T-1 ~> m s-1].
989 real, dimension(SZI_(G),SZJ_(G)), &
990 optional, intent(out) :: gustless_ustar !< The surface friction velocity without
991 !! any contributions from gustiness [Z T-1 ~> m s-1].
992 real, dimension(SZI_(G),SZJ_(G)), &
993 optional, intent(inout) :: mag_tau !< The magintude of the wind stress at tracer points
994 !! including subgridscale variability and gustiness [R Z2 T-2 ~> Pa]
995 real, dimension(SZI_(G),SZJ_(G)), &
996 optional, intent(out) :: gustless_mag_tau !< The magintude of the wind stress at tracer points
997 !! without any contributions from gustiness [R Z2 T-2 ~> Pa]
998 integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default.
999
1000 ! Local variables
1001 real, dimension(SZI_(G),SZJ_(G)) :: taux_in_A ! Zonal wind stresses [R Z L T-2 ~> Pa] at h points
1002 real, dimension(SZI_(G),SZJ_(G)) :: tauy_in_A ! Meridional wind stresses [R Z L T-2 ~> Pa] at h points
1003 real, dimension(SZIB_(G),SZJ_(G)) :: taux_in_C ! Zonal wind stresses [R Z L T-2 ~> Pa] at u points
1004 real, dimension(SZI_(G),SZJB_(G)) :: tauy_in_C ! Meridional wind stresses [R Z L T-2 ~> Pa] at v points
1005 real, dimension(SZIB_(G),SZJB_(G)) :: taux_in_B ! Zonal wind stresses [R Z L T-2 ~> Pa] at q points
1006 real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses [R Z L T-2 ~> Pa] at q points
1007
1008 real :: gustiness ! unresolved gustiness that contributes to ustar [R Z2 T-2 ~> Pa]
1009 real :: Irho0 ! Inverse of the Boussinesq mean density [R-1 ~> m3 kg-1]
1010 real :: taux2, tauy2 ! squared wind stresses [R2 Z2 L2 T-4 ~> Pa2]
1011 real :: tau_mag ! magnitude of the wind stress [R Z2 T-2 ~> Pa]
1012 real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
1013 real :: Pa_to_RZ2_T2 ! The combination of unit conversion factors used for mag_tau [R Z2 T-2 Pa-1 ~> 1]
1014
1015 logical :: do_ustar, do_gustless, do_tau_mag, do_gustless_tau_mag
1016 integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
1017 integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo
1018
1019 halo = 0 ; if (present(tau_halo)) halo = tau_halo
1020 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1021 ish = g%isc-halo ; ieh = g%iec+halo ; jsh = g%jsc-halo ; jeh = g%jec+halo
1022 isqh = g%IscB-halo ; ieqh = g%IecB+halo ; jsqh = g%JscB-halo ; jeqh = g%JecB+halo
1023 i0 = is - index_bounds(1) ; j0 = js - index_bounds(3)
1024
1025 irho0 = 1.0 / cs%Rho0
1026 stress_conversion = us%Pa_to_RLZ_T2 * cs%wind_stress_multiplier
1027
1028 do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)
1029 do_tau_mag = present(mag_tau) ; do_gustless_tau_mag = present(gustless_mag_tau)
1030
1031 wind_stagger = cs%wind_stagger
1032 if ((iob%wind_stagger == agrid) .or. (iob%wind_stagger == bgrid_ne) .or. &
1033 (iob%wind_stagger == cgrid_ne)) wind_stagger = iob%wind_stagger
1034
1035 if (associated(iob%u_flux).neqv.associated(iob%v_flux)) call mom_error(fatal,"extract_IOB_stresses: "//&
1036 "associated(IOB%u_flux) /= associated(IOB%v_flux !!!")
1037 if (present(taux).neqv.present(tauy)) call mom_error(fatal,"extract_IOB_stresses: "//&
1038 "present(taux) /= present(tauy) !!!")
1039
1040 ! Set surface momentum stress related fields as a function of staggering.
1041 if (present(taux) .or. present(tauy) .or. &
1042 ((do_ustar .or. do_tau_mag .or. do_gustless .or. do_gustless_tau_mag) &
1043 .and. .not.associated(iob%stress_mag)) ) then
1044
1045 if (wind_stagger == bgrid_ne) then
1046 taux_in_b(:,:) = 0.0 ; tauy_in_b(:,:) = 0.0
1047 if (associated(iob%u_flux).and.associated(iob%v_flux)) then
1048 do j=js,je ; do i=is,ie
1049 taux_in_b(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
1050 tauy_in_b(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
1051 enddo ; enddo
1052 endif
1053
1054 if (g%symmetric) call fill_symmetric_edges(taux_in_b, tauy_in_b, g%Domain, stagger=bgrid_ne)
1055 call pass_vector(taux_in_b, tauy_in_b, g%Domain, stagger=bgrid_ne, halo=max(1,halo))
1056
1057 if (present(taux).and.present(tauy)) then
1058 do j=jsh,jeh ; do i=isqh,ieqh
1059 taux(i,j) = 0.0
1060 if ((g%mask2dBu(i,j) + g%mask2dBu(i,j-1)) > 0.0) &
1061 taux(i,j) = (g%mask2dBu(i,j)*taux_in_b(i,j) + g%mask2dBu(i,j-1)*taux_in_b(i,j-1)) / &
1062 (g%mask2dBu(i,j) + g%mask2dBu(i,j-1))
1063 enddo ; enddo
1064 do j=jsqh,jeqh ; do i=ish,ieh
1065 tauy(i,j) = 0.0
1066 if ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j)) > 0.0) &
1067 tauy(i,j) = (g%mask2dBu(i,j)*tauy_in_b(i,j) + g%mask2dBu(i-1,j)*tauy_in_b(i-1,j)) / &
1068 (g%mask2dBu(i,j) + g%mask2dBu(i-1,j))
1069 enddo ; enddo
1070 endif
1071 elseif (wind_stagger == agrid) then
1072 taux_in_a(:,:) = 0.0 ; tauy_in_a(:,:) = 0.0
1073 if (associated(iob%u_flux).and.associated(iob%v_flux)) then
1074 do j=js,je ; do i=is,ie
1075 taux_in_a(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
1076 tauy_in_a(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
1077 enddo ; enddo
1078 endif
1079
1080 if (halo == 0) then
1081 call pass_vector(taux_in_a, tauy_in_a, g%Domain, to_all+omit_corners, stagger=agrid, halo=1)
1082 else
1083 call pass_vector(taux_in_a, tauy_in_a, g%Domain, stagger=agrid, halo=max(1,halo))
1084 endif
1085
1086 if (present(taux)) then ; do j=jsh,jeh ; do i=isqh,ieqh
1087 taux(i,j) = 0.0
1088 if ((g%mask2dT(i,j) + g%mask2dT(i+1,j)) > 0.0) &
1089 taux(i,j) = (g%mask2dT(i,j)*taux_in_a(i,j) + g%mask2dT(i+1,j)*taux_in_a(i+1,j)) / &
1090 (g%mask2dT(i,j) + g%mask2dT(i+1,j))
1091 enddo ; enddo ; endif
1092
1093 if (present(tauy)) then ; do j=jsqh,jeqh ; do i=ish,ieh
1094 tauy(i,j) = 0.0
1095 if ((g%mask2dT(i,j) + g%mask2dT(i,j+1)) > 0.0) &
1096 tauy(i,j) = (g%mask2dT(i,j)*tauy_in_a(i,j) + g%mask2dT(i,j+1)*tauy_in_a(i,j+1)) / &
1097 (g%mask2dT(i,j) + g%mask2dT(i,j+1))
1098 enddo ; enddo ; endif
1099
1100 else ! C-grid wind stresses.
1101 taux_in_c(:,:) = 0.0 ; tauy_in_c(:,:) = 0.0
1102 if (associated(iob%u_flux).and.associated(iob%v_flux)) then
1103 do j=js,je ; do i=is,ie
1104 taux_in_c(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
1105 tauy_in_c(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
1106 enddo ; enddo
1107 endif
1108
1109 if (g%symmetric) call fill_symmetric_edges(taux_in_c, tauy_in_c, g%Domain)
1110 call pass_vector(taux_in_c, tauy_in_c, g%Domain, halo=max(1,halo))
1111
1112 if (present(taux).and.present(tauy)) then
1113 do j=jsh,jeh ; do i=isqh,ieqh
1114 taux(i,j) = g%mask2dCu(i,j)*taux_in_c(i,j)
1115 enddo ; enddo
1116 do j=jsqh,jeqh ; do i=ish,ieh
1117 tauy(i,j) = g%mask2dCv(i,j)*tauy_in_c(i,j)
1118 enddo ; enddo
1119 endif
1120 endif ! endif for extracting wind stress fields with various staggerings
1121 endif
1122
1123 if (do_ustar .or. do_tau_mag .or. do_gustless .or. do_gustless_tau_mag) then
1124 ! Set surface friction velocity directly or as a function of staggering.
1125 ! ustar is required for the bulk mixed layer formulation and other turbulent mixing
1126 ! parametizations. The background gustiness (for example with a relatively small value
1127 ! of 0.02 Pa) is intended to give reasonable behavior in regions of very weak winds.
1128 if (associated(iob%stress_mag)) then
1129 pa_to_rz2_t2 = us%Pa_to_RLZ_T2 * us%L_to_Z
1130
1131 if (do_ustar .or. do_tau_mag) then ; do j=js,je ; do i=is,ie
1132 gustiness = cs%gust_const
1133 if (cs%read_gust_2d) then
1134 if ((wind_stagger == cgrid_ne) .or. &
1135 ((wind_stagger == agrid) .and. (g%mask2dT(i,j) > 0.0)) .or. &
1136 ((wind_stagger == bgrid_ne) .and. &
1137 (((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + &
1138 (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) > 0.0)) ) &
1139 gustiness = cs%gust(i,j)
1140 endif
1141 if (do_tau_mag) &
1142 mag_tau(i,j) = gustiness + pa_to_rz2_t2*iob%stress_mag(i-i0,j-j0)
1143 if (do_gustless_tau_mag) &
1144 gustless_mag_tau(i,j) = pa_to_rz2_t2*iob%stress_mag(i-i0,j-j0)
1145 if (do_ustar) &
1146 ustar(i,j) = sqrt(gustiness*irho0 + irho0*pa_to_rz2_t2*iob%stress_mag(i-i0,j-j0))
1147 enddo ; enddo ; endif
1148 if (cs%answer_date < 20190101) then
1149 if (do_gustless) then ; do j=js,je ; do i=is,ie
1150 gustless_ustar(i,j) = sqrt(pa_to_rz2_t2*iob%stress_mag(i-i0,j-j0) / cs%Rho0)
1151 enddo ; enddo ; endif
1152 else
1153 if (do_gustless) then ; do j=js,je ; do i=is,ie
1154 gustless_ustar(i,j) = sqrt(irho0 * pa_to_rz2_t2*iob%stress_mag(i-i0,j-j0))
1155 enddo ; enddo ; endif
1156 endif
1157 elseif (wind_stagger == bgrid_ne) then
1158 do j=js,je ; do i=is,ie
1159 tau_mag = 0.0 ; gustiness = cs%gust_const
1160 if (((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + &
1161 (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) > 0.0) then
1162 tau_mag = us%L_to_Z * sqrt(((g%mask2dBu(i,j)*((taux_in_b(i,j)**2) + (tauy_in_b(i,j)**2)) + &
1163 g%mask2dBu(i-1,j-1)*((taux_in_b(i-1,j-1)**2) + (tauy_in_b(i-1,j-1)**2))) + &
1164 (g%mask2dBu(i,j-1)*((taux_in_b(i,j-1)**2) + (tauy_in_b(i,j-1)**2)) + &
1165 g%mask2dBu(i-1,j)*((taux_in_b(i-1,j)**2) + (tauy_in_b(i-1,j)**2))) ) / &
1166 ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) )
1167 if (cs%read_gust_2d) gustiness = cs%gust(i,j)
1168 endif
1169 if (do_ustar) ustar(i,j) = sqrt(gustiness*irho0 + irho0 * tau_mag)
1170 if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
1171 if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
1172 if (cs%answer_date < 20190101) then
1173 if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / cs%Rho0)
1174 else
1175 if (do_gustless) gustless_ustar(i,j) = sqrt(irho0 * tau_mag)
1176 endif
1177 enddo ; enddo
1178 elseif (wind_stagger == agrid) then
1179 do j=js,je ; do i=is,ie
1180 tau_mag = g%mask2dT(i,j) * us%L_to_Z * sqrt((taux_in_a(i,j)**2) + (tauy_in_a(i,j)**2))
1181 gustiness = cs%gust_const
1182 if (cs%read_gust_2d .and. (g%mask2dT(i,j) > 0.0)) gustiness = cs%gust(i,j)
1183 if (do_ustar) ustar(i,j) = sqrt(gustiness*irho0 + irho0 * tau_mag)
1184 if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
1185 if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
1186 if (cs%answer_date < 20190101) then
1187 if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / cs%Rho0)
1188 else
1189 if (do_gustless) gustless_ustar(i,j) = sqrt(irho0 * tau_mag)
1190 endif
1191 enddo ; enddo
1192 else ! C-grid wind stresses.
1193 do j=js,je ; do i=is,ie
1194 taux2 = 0.0 ; tauy2 = 0.0
1195 if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0.0) &
1196 taux2 = (g%mask2dCu(i-1,j)*(taux_in_c(i-1,j)**2) + g%mask2dCu(i,j)*(taux_in_c(i,j)**2)) / &
1197 (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
1198 if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0.0) &
1199 tauy2 = (g%mask2dCv(i,j-1)*(tauy_in_c(i,j-1)**2) + g%mask2dCv(i,j)*(tauy_in_c(i,j)**2)) / &
1200 (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
1201 tau_mag = us%L_to_Z * sqrt(taux2 + tauy2)
1202
1203 gustiness = cs%gust_const
1204 if (cs%read_gust_2d) gustiness = cs%gust(i,j)
1205
1206 if (do_ustar) ustar(i,j) = sqrt(gustiness*irho0 + irho0 * tau_mag)
1207 if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
1208 if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
1209 if (cs%answer_date < 20190101) then
1210 if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / cs%Rho0)
1211 else
1212 if (do_gustless) gustless_ustar(i,j) = sqrt(irho0 * tau_mag)
1213 endif
1214 enddo ; enddo
1215 endif ! endif for wind friction velocity fields
1216 endif
1217
1218end subroutine extract_iob_stresses
1219
1220
1221!> Adds thermodynamic flux adjustments obtained via data_override
1222!! Component name is 'OCN'
1223!! Available adjustments are:
1224!! - hflx_adj (Heat flux into the ocean [W m-2])
1225!! - sflx_adj (Salt flux into the ocean [kg salt m-2 s-1])
1226!! - prcme_adj (Fresh water flux into the ocean [kg m-2 s-1])
1227subroutine apply_flux_adjustments(G, US, CS, Time, fluxes)
1228 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1229 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1230 type(surface_forcing_cs), pointer :: CS !< Surface forcing control structure
1231 type(time_type), intent(in) :: Time !< Model time structure
1232 type(forcing), intent(inout) :: fluxes !< Surface fluxes structure
1233
1234 ! Local variables
1235 real, dimension(G%isc:G%iec,G%jsc:G%jec) :: temp_at_h ! Various fluxes at h points
1236 ! [Q R Z T-1 ~> W m-2] or [R Z T-1 ~> kg m-2 s-1]
1237
1238 integer :: isc, iec, jsc, jec, i, j
1239 logical :: overrode_h
1240
1241 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1242
1243 call data_override(g%Domain, 'hflx_adj', temp_at_h, time, override=overrode_h, &
1244 scale=us%W_m2_to_QRZ_T)
1245
1246 if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
1247 fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j) * g%mask2dT(i,j)
1248 enddo ; enddo ; endif
1249 ! Not needed? ! if (overrode_h) call pass_var(fluxes%heat_added, G%Domain)
1250
1251 call data_override(g%Domain, 'sflx_adj', temp_at_h, time, override=overrode_h, &
1252 scale=us%kg_m2s_to_RZ_T)
1253
1254 if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
1255 fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + temp_at_h(i,j) * g%mask2dT(i,j)
1256 enddo ; enddo ; endif
1257 ! Not needed? ! if (overrode_h) call pass_var(fluxes%salt_flux_added, G%Domain)
1258
1259 call data_override(g%Domain, 'prcme_adj', temp_at_h, time, override=overrode_h, &
1260 scale=us%kg_m2s_to_RZ_T)
1261
1262 if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
1263 fluxes%vprec(i,j) = fluxes%vprec(i,j) + temp_at_h(i,j)* g%mask2dT(i,j)
1264 enddo ; enddo ; endif
1265 ! Not needed? ! if (overrode_h) call pass_var(fluxes%vprec, G%Domain)
1266end subroutine apply_flux_adjustments
1267
1268!> Adds mechanical forcing adjustments obtained via data_override
1269!! Component name is 'OCN'
1270!! Available adjustments are:
1271!! - taux_adj (Zonal wind stress delta, positive to the east [Pa])
1272!! - tauy_adj (Meridional wind stress delta, positive to the north [Pa])
1273subroutine apply_force_adjustments(G, US, CS, Time, forces)
1274 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1275 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1276 type(surface_forcing_cs), pointer :: CS !< Surface forcing control structure
1277 type(time_type), intent(in) :: Time !< Model time structure
1278 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
1279
1280 ! Local variables
1281 real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points [R Z L T-2 ~> Pa]
1282 real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points [R Z L T-2 ~> Pa]
1283
1284 integer :: isc, iec, jsc, jec, i, j
1285 real :: dLonDx, dLonDy ! The change in longitude across the cell in the x- and y-directions [degrees_E]
1286 real :: rDlon ! The magnitude of the change in longitude [degrees_E] and then its inverse [degrees_E-1]
1287 real :: cosA, sinA ! The cosine and sine of the angle between the grid and true north [nondim]
1288 real :: zonal_tau, merid_tau ! True zonal and meridional wind stresses [R Z L T-2 ~> Pa]
1289 logical :: overrode_x, overrode_y
1290
1291 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1292
1293 tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0
1294 ! Either reads data or leaves contents unchanged
1295 overrode_x = .false. ; overrode_y = .false.
1296 call data_override(g%Domain, 'taux_adj', tempx_at_h(isc:iec,jsc:jec), time, &
1297 override=overrode_x, scale=us%Pa_to_RLZ_T2)
1298 call data_override(g%Domain, 'tauy_adj', tempy_at_h(isc:iec,jsc:jec), time, &
1299 override=overrode_y, scale=us%Pa_to_RLZ_T2)
1300
1301 if (overrode_x .or. overrode_y) then
1302 if (.not. (overrode_x .and. overrode_y)) call mom_error(fatal,"apply_flux_adjustments: "//&
1303 "Both taux_adj and tauy_adj must be specified, or neither, in data_table")
1304
1305 ! Rotate winds
1306 call pass_vector(tempx_at_h, tempy_at_h, g%Domain, to_all, agrid, halo=1)
1307 do j=jsc-1,jec+1 ; do i=isc-1,iec+1
1308 dlondx = g%geoLonCu(i,j) - g%geoLonCu(i-1,j)
1309 dlondy = g%geoLonCv(i,j) - g%geoLonCv(i,j-1)
1310 rdlon = sqrt( dlondx * dlondx + dlondy * dlondy )
1311 if (rdlon > 0.) rdlon = 1. / rdlon
1312 cosa = dlondx * rdlon
1313 sina = dlondy * rdlon
1314 zonal_tau = tempx_at_h(i,j)
1315 merid_tau = tempy_at_h(i,j)
1316 tempx_at_h(i,j) = cosa * zonal_tau - sina * merid_tau
1317 tempy_at_h(i,j) = sina * zonal_tau + cosa * merid_tau
1318 enddo ; enddo
1319
1320 ! Average to C-grid locations
1321 do j=jsc,jec ; do i=isc-1,iec
1322 forces%taux(i,j) = forces%taux(i,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) )
1323 enddo ; enddo
1324
1325 do j=jsc-1,jec ; do i=isc,iec
1326 forces%tauy(i,j) = forces%tauy(i,j) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) )
1327 enddo ; enddo
1328 endif ! overrode_x .or. overrode_y
1329
1330end subroutine apply_force_adjustments
1331
1332!> Save any restart files associated with the surface forcing.
1333subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
1334 filename_suffix)
1335 type(surface_forcing_cs), pointer :: cs !< A pointer to the control structure returned
1336 !! by a previous call to surface_forcing_init
1337 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
1338 type(time_type), intent(in) :: time !< The current model time
1339 character(len=*), intent(in) :: directory !< The directory into which to write the
1340 !! restart files
1341 logical, optional, intent(in) :: time_stamped !< If true, the restart file names include
1342 !! a unique time stamp. The default is false.
1343 character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-
1344 !! stamp) to append to the restart file names.
1345
1346 if (.not.associated(cs)) return
1347 if (.not.associated(cs%restart_CSp)) return
1348 call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1349
1350end subroutine forcing_save_restart
1351
1352!> Initialize the surface forcing, including setting parameters and allocating permanent memory.
1353subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
1354 type(time_type), intent(in) :: time !< The current model time
1355 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1356 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1357 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1358 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
1359 !! diagnostic output
1360 type(surface_forcing_cs), pointer :: cs !< A pointer that is set to point to the control
1361 !! structure for this module
1362 integer, optional, intent(in) :: wind_stagger !< If present, the staggering of the winds
1363 !! that are being provided in calls to update_ocean_model
1364
1365 ! Local variables
1366 real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1].
1367 real, dimension(SZI_(G),SZJ_(G)) :: &
1368 utide_2d ! A 2d array of RMS tidal velocities [Z T-1 ~> m s-1].
1369 real :: flux_const_dflt ! A default piston velocity for restoring surface properties [m day-1]
1370 logical :: boussinesq ! If true, this run is fully Boussinesq
1371 logical :: semi_boussinesq ! If true, this run is partially non-Boussinesq
1372 real :: rho_tke_tidal ! The constant bottom density used to translate tidal amplitudes into
1373 ! the tidal bottom TKE input used with INT_TIDE_DISSIPATION, times the
1374 ! factor rescaling from the units of TKE to those of mean kinetic
1375 ! energy [R L2 Z-2 ~> kg m-3]
1376 logical :: new_sim ! False if this simulation was started from a restart file
1377 ! or other equivalent files.
1378 logical :: iceberg_flux_diags ! If true, diagnostics of fluxes from icebergs are available.
1379 logical :: fix_ustar_gustless_bug ! If false, include a bug using an older run-time parameter.
1380 logical :: test_value ! This is used to determine whether a logical parameter is being set explicitly.
1381 logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly.
1382 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1383 type(time_type) :: time_frc
1384 type(directories) :: dirs ! A structure containing relevant directory paths and input filenames.
1385 character(len=200) :: tideamp_file, gust_file, salt_file, temp_file ! Input file names.
1386 ! This include declares and sets the variable "version".
1387# include "version_variable.h"
1388 character(len=40) :: mdl = "MOM_surface_forcing" ! This module's name.
1389 character(len=48) :: stagger
1390 character(len=80) :: varnam
1391 character(len=240) :: flnam
1392 character(len=240) :: basin_file
1393 integer :: i, j, isd, ied, jsd, jed
1394
1395 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1396
1397 if (associated(cs)) then
1398 call mom_error(warning, "surface_forcing_init called with an associated "// &
1399 "control structure.")
1400 return
1401 endif
1402 allocate(cs)
1403
1404 id_clock_forcing=cpu_clock_id('Ocean surface forcing', grain=clock_subcomponent)
1405 call cpu_clock_begin(id_clock_forcing)
1406
1407 cs%diag => diag
1408
1409 call write_version_number(version)
1410 ! Read all relevant parameters and write them to the model log.
1411 call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.)
1412
1413 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, &
1414 "The directory in which all input files are found.", &
1415 default=".")
1416 cs%inputdir = slasher(cs%inputdir)
1417 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
1418 "If true, Temperature and salinity are used as state "//&
1419 "variables.", default=.true.)
1420 call get_param(param_file, mdl, "BOUSSINESQ", boussinesq, &
1421 "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
1422 call get_param(param_file, mdl, "SEMI_BOUSSINESQ", semi_boussinesq, &
1423 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
1424 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
1425 "height changes. This only applies if BOUSSINESQ is false.", &
1426 default=.true., do_not_log=.true.)
1427 cs%nonBous = .not.(boussinesq .or. semi_boussinesq)
1428 call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
1429 "The mean ocean density used with BOUSSINESQ true to "//&
1430 "calculate accelerations and the mass for conservation "//&
1431 "properties, or with BOUSSINESQ false to convert some "//&
1432 "parameters from vertical units of m to kg m-2.", &
1433 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R) ! (, do_not_log=CS%nonBous)
1434 call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1435 "The latent heat of fusion.", units="J/kg", default=hlf, scale=us%J_kg_to_Q)
1436 call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1437 "The latent heat of fusion.", units="J/kg", default=hlv, scale=us%J_kg_to_Q)
1438 call get_param(param_file, mdl, "MAX_P_SURF", cs%max_p_surf, &
1439 "The maximum surface pressure that can be exerted by the "//&
1440 "atmosphere and floating sea-ice or ice shelves. This is "//&
1441 "needed because the FMS coupling structure does not "//&
1442 "limit the water that can be frozen out of the ocean and "//&
1443 "the ice-ocean heat fluxes are treated explicitly. No "//&
1444 "limit is applied if a negative value is used.", &
1445 units="Pa", default=-1.0, scale=us%Pa_to_RL2_T2)
1446 call get_param(param_file, mdl, "RESTORE_SALINITY", cs%restore_salt, &
1447 "If true, the coupled driver will add a globally-balanced "//&
1448 "fresh-water flux that drives sea-surface salinity "//&
1449 "toward specified values.", default=.false.)
1450 call get_param(param_file, mdl, "RESTORE_TEMPERATURE", cs%restore_temp, &
1451 "If true, the coupled driver will add a "//&
1452 "heat flux that drives sea-surface temperature "//&
1453 "toward specified values.", default=.false.)
1454 call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", &
1455 cs%adjust_net_srestore_to_zero, &
1456 "If true, adjusts the salinity restoring seen to zero "//&
1457 "whether restoring is via a salt flux or virtual precip.",&
1458 default=cs%restore_salt)
1459 call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_BY_SCALING", &
1460 cs%adjust_net_srestore_by_scaling, &
1461 "If true, adjustments to salt restoring to achieve zero net are "//&
1462 "made by scaling values without moving the zero contour.",&
1463 default=.false.)
1464 call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_TO_ZERO", &
1465 cs%adjust_net_fresh_water_to_zero, &
1466 "If true, adjusts the net fresh-water forcing seen "//&
1467 "by the ocean (including restoring) to zero.", default=.false.)
1468 if (cs%adjust_net_fresh_water_to_zero) &
1469 call get_param(param_file, mdl, "USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
1470 cs%use_net_FW_adjustment_sign_bug, &
1471 "If true, use the wrong sign for the adjustment to "//&
1472 "the net fresh-water.", default=.false.)
1473 call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", &
1474 cs%adjust_net_fresh_water_by_scaling, &
1475 "If true, adjustments to net fresh water to achieve zero net are "//&
1476 "made by scaling values without moving the zero contour.",&
1477 default=.false.)
1478 call get_param(param_file, mdl, "ICE_SALT_CONCENTRATION", &
1479 cs%ice_salt_concentration, &
1480 "The assumed sea-ice salinity needed to reverse engineer the "//&
1481 "melt flux (or ice-ocean fresh-water flux).", &
1482 units="kg/kg", default=0.005)
1483 call get_param(param_file, mdl, "USE_LIMITED_PATM_SSH", cs%use_limited_P_SSH, &
1484 "If true, return the sea surface height with the "//&
1485 "correction for the atmospheric (and sea-ice) pressure "//&
1486 "limited by max_p_surf instead of the full atmospheric "//&
1487 "pressure.", default=.true.)
1488 call get_param(param_file, mdl, "APPROX_NET_MASS_SRC", cs%approx_net_mass_src, &
1489 "If true, use the net mass sources from the ice-ocean "//&
1490 "boundary type without any further adjustments to drive "//&
1491 "the ocean dynamics. The actual net mass source may differ "//&
1492 "due to internal corrections.", default=.false.)
1493
1494 if (present(wind_stagger)) then
1495 if (wind_stagger == agrid) then ; stagger = 'AGRID'
1496 elseif (wind_stagger == bgrid_ne) then ; stagger = 'BGRID_NE'
1497 elseif (wind_stagger == cgrid_ne) then ; stagger = 'CGRID_NE'
1498 else ; stagger = 'UNKNOWN' ; call mom_error(fatal,"surface_forcing_init: WIND_STAGGER = "// &
1499 trim(stagger)// "is invalid.") ; endif
1500 call log_param(param_file, mdl, "WIND_STAGGER", stagger, &
1501 "The staggering of the input wind stress field "//&
1502 "from the coupler that is actually used.")
1503 cs%wind_stagger = wind_stagger
1504 else
1505 call get_param(param_file, mdl, "WIND_STAGGER", stagger, &
1506 "A case-insensitive character string to indicate the "//&
1507 "staggering of the input wind stress field. Valid "//&
1508 "values are 'A', 'B', or 'C'.", default="C")
1509 if (uppercase(stagger(1:1)) == 'A') then ; cs%wind_stagger = agrid
1510 elseif (uppercase(stagger(1:1)) == 'B') then ; cs%wind_stagger = bgrid_ne
1511 elseif (uppercase(stagger(1:1)) == 'C') then ; cs%wind_stagger = cgrid_ne
1512 else ; call mom_error(fatal,"surface_forcing_init: WIND_STAGGER = "// &
1513 trim(stagger)//" is invalid.") ; endif
1514 endif
1515
1516 call get_param(param_file, mdl, "WIND_STRESS_MULTIPLIER", cs%wind_stress_multiplier, &
1517 "A factor multiplying the wind-stress given to the ocean by the "//&
1518 "coupler. This is used for testing and should be =1.0 for any "//&
1519 "production runs.", units="nondim", default=1.0)
1520
1521 if (cs%restore_salt) then
1522 call get_param(param_file, mdl, "FLUXCONST", flux_const_dflt, &
1523 "The constant that relates the restoring surface fluxes to the relative "//&
1524 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
1525 units="m day-1", default=0.0)
1526 call get_param(param_file, mdl, "FLUXCONST_SALT", cs%Flux_const_salt, &
1527 "The constant that relates the restoring surface salt fluxes to the relative "//&
1528 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
1529 units="m day-1", default=flux_const_dflt, scale=us%m_to_Z*us%T_to_s)
1530 ! Finish converting CS%Flux_const_salt from m day-1 to [Z T-1 ~> m s-1]. Ideally this would be
1531 ! included in the scale factors above, but doing so would change answers because a/b /= a*(1/b).
1532 cs%Flux_const_salt = cs%Flux_const_salt / 86400.0
1533 call get_param(param_file, mdl, "SALT_RESTORE_FILE", cs%salt_restore_file, &
1534 "A file in which to find the surface salinity to use for restoring.", &
1535 default="salt_restore.nc")
1536 call get_param(param_file, mdl, "SALT_RESTORE_VARIABLE", cs%salt_restore_var_name, &
1537 "The name of the surface salinity variable to read from "//&
1538 "SALT_RESTORE_FILE for restoring salinity.", &
1539 default="salt")
1540 call get_param(param_file, mdl, "SALT_RESTORE_PRACTICAL_SALINITY", cs%salt_restore_is_practical, &
1541 "Specifies if the restoring surface salinity variable is practical salinity. If this "//&
1542 "flag is set to false it is assumed that the salinity is absolute salinity.", default=.false.)
1543 call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", cs%salt_restore_as_sflux, &
1544 "If true, the restoring of salinity is applied as a salt "//&
1545 "flux instead of as a freshwater flux.", default=.false.)
1546 call get_param(param_file, mdl, "MAX_DELTA_SRESTORE_FROM_FILE", cs%max_delta_srestore_file, &
1547 "If true, read a file MAX_DELTA_SRESTORE_FILE containing the field "//&
1548 "MAX_DELTA_SRESTORE_VARNAME for the maximum salinity difference used in "//&
1549 "restoring terms. Where the field's value is negative turn off restoring when "//&
1550 "the salinity difference magnitude exceeds abs(value).", default=.false.)
1551 if (.not. cs%max_delta_srestore_file) then
1552 call get_param(param_file, mdl, "MAX_DELTA_SRESTORE", cs%max_delta_srestore, &
1553 "The maximum salinity difference used in restoring terms.", &
1554 units="PSU or g kg-1", default=999.0, scale=us%ppt_to_S)
1555 else
1556 call get_param(param_file, mdl, "MAX_DELTA_SRESTORE_FILE", flnam, &
1557 "The path to the file containing the maximum salinity difference field.", &
1558 default="max_delta_srestore.nc")
1559 flnam = trim(cs%inputdir) // trim(flnam)
1560 call get_param(param_file, mdl, "MAX_DELTA_SRESTORE_VARNAME", varnam, &
1561 "The name of the maximum salinity difference variable in the input file.", &
1562 default="max_delta_srestore")
1563 cs%max_delta_srestore = 999.0
1564 call safe_alloc_ptr(cs%max_delta_srestore_2d,isd,ied,jsd,jed)
1565 call mom_read_data(flnam,varnam, cs%max_delta_srestore_2d, g%domain, timelevel=1)
1566 endif
1567 call get_param(param_file, mdl, "MIN_RATIO_SRESTORE", cs%min_ratio_srestore, &
1568 "Turn off MAX_DELTA_SRESTORE where the ratio of SSS to restoring salinity "//&
1569 "is less than this value.", units="nondim", default=0.0)
1570 call get_param(param_file, mdl, "MASK_SRESTORE_UNDER_ICE", cs%mask_srestore_under_ice, &
1571 "If true, disables SSS restoring under sea-ice based on a frazil "//&
1572 "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", &
1573 default=.false.)
1574 call get_param(param_file, mdl, "MASK_SRESTORE_MARGINAL_SEAS", &
1575 cs%mask_srestore_marginal_seas, &
1576 "If true, disable SSS restoring in marginal seas. Only used when "//&
1577 "RESTORE_SALINITY is True.", default=.false.)
1578 call get_param(param_file, mdl, "BASIN_FILE", basin_file, &
1579 "A file in which to find the basin masks, in variable 'basin'.", &
1580 default="basin.nc")
1581 basin_file = trim(cs%inputdir) // trim(basin_file)
1582 call safe_alloc_ptr(cs%basin_mask,isd,ied,jsd,jed) ; cs%basin_mask(:,:) = 1.0
1583 if (cs%mask_srestore_marginal_seas) then
1584 call mom_read_data(basin_file,'basin',cs%basin_mask,g%domain, timelevel=1)
1585 do j=jsd,jed ; do i=isd,ied
1586 if (cs%basin_mask(i,j) >= 6.0) then ; cs%basin_mask(i,j) = 0.0
1587 else ; cs%basin_mask(i,j) = 1.0 ; endif
1588 enddo ; enddo
1589 endif
1590 call get_param(param_file, mdl, "MASK_SRESTORE", cs%mask_srestore, &
1591 "If true, read a file (salt_restore_mask) containing "//&
1592 "a mask for SSS restoring.", default=.false.)
1593 endif
1594
1595 if (cs%restore_temp) then
1596 call get_param(param_file, mdl, "FLUXCONST", flux_const_dflt, &
1597 "The constant that relates the restoring surface fluxes to the relative "//&
1598 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
1599 units="m day-1", default=0.0)
1600 call get_param(param_file, mdl, "FLUXCONST_TEMP", cs%Flux_const_temp, &
1601 "The constant that relates the restoring surface temperature fluxes to the relative "//&
1602 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
1603 units="m day-1", default=flux_const_dflt, scale=us%m_to_Z*us%T_to_s)
1604 ! Finish converting CS%Flux_const_temp from [m day-1] to [Z T-1 ~> m s-1]. Ideally this would be
1605 ! included in the scale factors above, but doing so would change answers because a/b /= a*(1/b).
1606 cs%Flux_const_temp = cs%Flux_const_temp / 86400.0
1607 call get_param(param_file, mdl, "SST_RESTORE_FILE", cs%temp_restore_file, &
1608 "A file in which to find the surface temperature to use for restoring.", &
1609 default="temp_restore.nc")
1610 call get_param(param_file, mdl, "SST_RESTORE_VARIABLE", cs%temp_restore_var_name, &
1611 "The name of the surface temperature variable to read from "//&
1612 "SST_RESTORE_FILE for restoring sst.", &
1613 default="temp")
1614
1615 call get_param(param_file, mdl, "MAX_DELTA_TRESTORE", cs%max_delta_trestore, &
1616 "The maximum sst difference used in restoring terms.", &
1617 units="degC ", default=999.0, scale=us%degC_to_C)
1618 call get_param(param_file, mdl, "MASK_TRESTORE", cs%mask_trestore, &
1619 "If true, read a file (temp_restore_mask) containing "//&
1620 "a mask for SST restoring.", default=.false.)
1621
1622 call get_param(param_file, mdl, "SPEAR_ECDA_SST_RESTORE_TFREEZE", cs%trestore_SPEAR_ECDA, &
1623 "If true, modify SST restoring field using SSS state. This only modifies the "//&
1624 "restoring data that is within 0.0001degC of -1.8degC.", default=.false.)
1625 else
1626 cs%trestore_SPEAR_ECDA = .false. ! Needed to toggle logging of SPEAR_DTFREEZE_DS
1627 endif
1628 call get_param(param_file, mdl, "SPEAR_DTFREEZE_DS", cs%SPEAR_dTf_dS, &
1629 "The derivative of the freezing temperature with salinity.", &
1630 units="degC ppt-1", default=-0.054, scale=us%degC_to_C*us%S_to_ppt, &
1631 do_not_log=.not.cs%trestore_SPEAR_ECDA)
1632 call get_param(param_file, mdl, "RESTORE_FLUX_RHO", cs%rho_restore, &
1633 "The density that is used to convert piston velocities into salt or heat "//&
1634 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
1635 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &
1636 do_not_log=.not.(cs%restore_temp.or.cs%restore_salt))
1637
1638 ! Optionally read tidal amplitude from input file [Z T-1 ~> m s-1] on model grid.
1639 ! Otherwise use default tidal amplitude for bottom frictionally-generated
1640 ! dissipation. Default cd_tides is chosen to yield approx 1 TWatt of
1641 ! work done against tides globally using OSU tidal amplitude.
1642 ! Note that the slightly unusual length scaling is deliberate, because the tidal
1643 ! amplitudes are used to set the friction velocity.
1644 call get_param(param_file, mdl, "CD_TIDES", cs%cd_tides, &
1645 "The drag coefficient that applies to the tides.", &
1646 units="nondim", default=1.0e-4)
1647 call get_param(param_file, mdl, "READ_TIDEAMP", cs%read_TIDEAMP, &
1648 "If true, read a file (given by TIDEAMP_FILE) containing "//&
1649 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1650 if (cs%read_TIDEAMP) then
1651 call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
1652 "The path to the file containing the spatially varying "//&
1653 "tidal amplitudes with INT_TIDE_DISSIPATION.", &
1654 default="tideamp.nc")
1655 cs%utide=0.0
1656 else
1657 call get_param(param_file, mdl, "UTIDE", cs%utide, &
1658 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1659 units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1660 endif
1661 call get_param(param_file, mdl, "TKE_TIDAL_RHO", rho_tke_tidal, &
1662 "The constant bottom density used to translate tidal amplitudes into the tidal "//&
1663 "bottom TKE input used with INT_TIDE_DISSIPATION.", &
1664 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R*us%Z_to_L**2, &
1665 do_not_log=.not.(cs%read_TIDEAMP.or.(cs%utide>0.0)))
1666
1667 call safe_alloc_ptr(cs%BBL_tidal_dis,isd,ied,jsd,jed)
1668 call safe_alloc_ptr(cs%ustar_tidal,isd,ied,jsd,jed)
1669
1670 if (cs%read_TIDEAMP) then
1671 tideamp_file = trim(cs%inputdir) // trim(tideamp_file)
1672 ! NOTE: There are certain cases where FMS is unable to read this file, so
1673 ! we use read_netCDF_data in place of MOM_read_data.
1674 utide_2d(:,:) = 0.0
1675 call read_netcdf_data(tideamp_file, 'tideamp', utide_2d, g%Domain, &
1676 rescale=us%m_to_Z*us%T_to_s)
1677 do j=jsd,jed ; do i=isd,ied
1678 utide = utide_2d(i,j)
1679 cs%BBL_tidal_dis(i,j) = g%mask2dT(i,j)*rho_tke_tidal*cs%cd_tides*(utide*utide*utide)
1680 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1681 enddo ; enddo
1682 else
1683 do j=jsd,jed ; do i=isd,ied
1684 utide = cs%utide
1685 cs%BBL_tidal_dis(i,j) = rho_tke_tidal*cs%cd_tides*(utide*utide*utide)
1686 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1687 enddo ; enddo
1688 endif
1689
1690 call time_interp_external_init()
1691
1692 ! Optionally read a x-y gustiness field in place of a global constant.
1693 call get_param(param_file, mdl, "READ_GUST_2D", cs%read_gust_2d, &
1694 "If true, use a 2-dimensional gustiness supplied from "//&
1695 "an input file", default=.false.)
1696 call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
1697 "The background gustiness in the winds.", &
1698 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2*us%L_to_Z)
1699 if (cs%read_gust_2d) then
1700 call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
1701 "The file in which the wind gustiness is found in "//&
1702 "variable gustiness.", fail_if_missing=.true.)
1703
1704 call safe_alloc_ptr(cs%gust,isd,ied,jsd,jed)
1705 gust_file = trim(cs%inputdir) // trim(gust_file)
1706 ! NOTE: There are certain cases where FMS is unable to read this file, so
1707 ! we use read_netCDF_data in place of MOM_read_data.
1708 call read_netcdf_data(gust_file, 'gustiness', cs%gust, g%Domain, &
1709 rescale=us%Pa_to_RLZ_T2*us%L_to_Z) ! units in file should be [Pa]
1710 endif
1711 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1712 "This sets the default value for the various _ANSWER_DATE parameters.", &
1713 default=99991231)
1714 call get_param(param_file, mdl, "SURFACE_FORCING_ANSWER_DATE", cs%answer_date, &
1715 "The vintage of the order of arithmetic and expressions in the gustiness "//&
1716 "calculations. Values below 20190101 recover the answers from the end "//&
1717 "of 2018, while higher values use a simpler expression to calculate gustiness.", &
1718 default=default_answer_date)
1719
1720 call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", cs%ustar_gustless_bug, &
1721 "If true include a bug in the time-averaging of the gustless wind friction velocity", &
1722 default=.false., do_not_log=.true.)
1723 ! This is used to test whether USTAR_GUSTLESS_BUG is being actively set.
1724 call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", test_value, default=.true., do_not_log=.true.)
1725 explicit_bug = cs%ustar_gustless_bug .eqv. test_value
1726 call get_param(param_file, mdl, "FIX_USTAR_GUSTLESS_BUG", fix_ustar_gustless_bug, &
1727 "If true correct a bug in the time-averaging of the gustless wind friction velocity", &
1728 default=.true., do_not_log=.true.)
1729 call get_param(param_file, mdl, "FIX_USTAR_GUSTLESS_BUG", test_value, default=.false., do_not_log=.true.)
1730 explicit_fix = fix_ustar_gustless_bug .eqv. test_value
1731
1732 if (explicit_bug .and. explicit_fix .and. (fix_ustar_gustless_bug .eqv. cs%ustar_gustless_bug)) then
1733 ! USTAR_GUSTLESS_BUG is being explicitly set, and should not be changed.
1734 call mom_error(fatal, "USTAR_GUSTLESS_BUG and FIX_USTAR_GUSTLESS_BUG are both being set "//&
1735 "with inconsistent values. FIX_USTAR_GUSTLESS_BUG is an obsolete "//&
1736 "parameter and should be removed.")
1737 elseif (explicit_fix) then
1738 call mom_error(warning, "FIX_USTAR_GUSTLESS_BUG is an obsolete parameter. "//&
1739 "Use USTAR_GUSTLESS_BUG instead (noting that it has the opposite sense).")
1740 cs%ustar_gustless_bug = .not.fix_ustar_gustless_bug
1741 endif
1742 call log_param(param_file, mdl, "USTAR_GUSTLESS_BUG", cs%ustar_gustless_bug, &
1743 "If true include a bug in the time-averaging of the gustless wind friction velocity", &
1744 default=.false.)
1745
1746
1747! See whether sufficiently thick sea ice should be treated as rigid.
1748 call get_param(param_file, mdl, "USE_RIGID_SEA_ICE", cs%rigid_sea_ice, &
1749 "If true, sea-ice is rigid enough to exert a "//&
1750 "nonhydrostatic pressure that resist vertical motion.", &
1751 default=.false.)
1752 if (cs%rigid_sea_ice) then
1753 call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
1754 "The gravitational acceleration of the Earth.", &
1755 units="m s-2", default=9.80, scale=us%Z_to_m*us%m_s_to_L_T**2)
1756 call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", cs%density_sea_ice, &
1757 "A typical density of sea ice, used with the kinematic "//&
1758 "viscosity, when USE_RIGID_SEA_ICE is true.", &
1759 units="kg m-3", default=900.0, scale=us%kg_m3_to_R)
1760 call get_param(param_file, mdl, "SEA_ICE_VISCOSITY", cs%Kv_sea_ice, &
1761 "The kinematic viscosity of sufficiently thick sea ice "//&
1762 "for use in calculating the rigidity of sea ice.", &
1763 units="m2 s-1", default=1.0e9, scale=us%Z_to_L**2*us%m_to_L**2*us%T_to_s)
1764 call get_param(param_file, mdl, "SEA_ICE_RIGID_MASS", cs%rigid_sea_ice_mass, &
1765 "The mass of sea-ice per unit area at which the sea-ice "//&
1766 "starts to exhibit rigidity", &
1767 units="kg m-2", default=1000.0, scale=us%kg_m3_to_R*us%m_to_Z)
1768 endif
1769
1770 call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
1771 "If true, makes available diagnostics of fluxes from icebergs "//&
1772 "as seen by MOM6.", default=.false.)
1773 call get_param(param_file, mdl, "ALLOW_CARBON_FLUX_EXCHANGE", cs%allow_carbon_flux_exchange, &
1774 "If true, makes available fluxes and diagnostics of carbon in runoff "//&
1775 "within MOM6.", default=.false.)
1776 call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles, &
1777 use_berg_fluxes=iceberg_flux_diags, &
1778 use_carbon_runoff=cs%allow_carbon_flux_exchange)
1779
1780 call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", cs%allow_flux_adjustments, &
1781 "If true, allows flux adjustments to specified via the "//&
1782 "data_table using the component name 'OCN'.", default=.false.)
1783
1784 call get_param(param_file, mdl, "CHECK_NO_LAND_FLUXES", cs%check_no_land_fluxes, &
1785 "If true, checks that values from IOB fluxes are zero "//&
1786 "above land points (i.e. G%mask2dT = 0).", default=.false., &
1787 debuggingparam=.true.)
1788
1789 call data_override_init(g%Domain)
1790
1791 if (cs%restore_salt) then
1792 salt_file = trim(cs%inputdir) // trim(cs%salt_restore_file)
1793 cs%srestore_handle = init_external_field(salt_file, cs%salt_restore_var_name, mom_domain=g%Domain)
1794 call safe_alloc_ptr(cs%srestore_mask,isd,ied,jsd,jed) ; cs%srestore_mask(:,:) = 1.0
1795 if (cs%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes
1796 flnam = trim(cs%inputdir) // 'salt_restore_mask.nc'
1797 call mom_read_data(flnam,'mask', cs%srestore_mask, g%domain, timelevel=1)
1798 endif
1799 endif
1800
1801 if (cs%restore_temp) then
1802 temp_file = trim(cs%inputdir) // trim(cs%temp_restore_file)
1803 cs%trestore_handle = init_external_field(temp_file, cs%temp_restore_var_name, mom_domain=g%Domain)
1804 call safe_alloc_ptr(cs%trestore_mask,isd,ied,jsd,jed) ; cs%trestore_mask(:,:) = 1.0
1805 if (cs%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes
1806 flnam = trim(cs%inputdir) // 'temp_restore_mask.nc'
1807 call mom_read_data(flnam, 'mask', cs%trestore_mask, g%domain, timelevel=1)
1808 endif
1809 endif
1810
1811 ! Set up any restart fields associated with the forcing.
1812 call restart_init(param_file, cs%restart_CSp, "MOM_forcing.res")
1813!#CTRL# call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, &
1814!#CTRL# CS%restart_CSp)
1815 call restart_init_end(cs%restart_CSp)
1816
1817 if (associated(cs%restart_CSp)) then
1818 call get_mom_input(dirs=dirs)
1819
1820 new_sim = .false.
1821 if ((dirs%input_filename(1:1) == 'n') .and. &
1822 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1823 if (.not.new_sim) then
1824 call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1825 g, cs%restart_CSp)
1826 endif
1827 endif
1828
1829!#CTRL# call controlled_forcing_init(Time, G, US, param_file, diag, CS%ctrl_forcing_CSp)
1830
1831 call user_revise_forcing_init(param_file, cs%urf_CS)
1832
1833 call cpu_clock_end(id_clock_forcing)
1834end subroutine surface_forcing_init
1835
1836!> Clean up and deallocate any memory associated with this module and its children.
1837subroutine surface_forcing_end(CS, fluxes)
1838 type(surface_forcing_cs), pointer :: CS !< A pointer to the control structure returned by
1839 !! a previous call to surface_forcing_init, it will
1840 !! be deallocated here.
1841 type(forcing), optional, intent(inout) :: fluxes !< A structure containing pointers to all
1842 !! possible mass, heat or salt flux forcing fields.
1843 !! If present, it will be deallocated here.
1844
1845 if (present(fluxes)) call deallocate_forcing_type(fluxes)
1846
1847!#CTRL# call controlled_forcing_end(CS%ctrl_forcing_CSp)
1848
1849 if (associated(cs)) deallocate(cs)
1850 cs => null()
1851
1852end subroutine surface_forcing_end
1853
1854!> Write out a set of messages with checksums of the fields in an ice_ocean_boundary type
1855subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
1856
1857 character(len=*), intent(in) :: id !< An identifying string for this call
1858 integer, intent(in) :: timestep !< The number of elapsed timesteps
1860 intent(in) :: iobt !< An ice-ocean boundary type with fluxes to drive the
1861 !! ocean in a coupled model whose checksums are reported
1862 ! Local variables
1863 integer(kind=int64) :: chks ! A checksum for the field
1864 logical :: root ! True only on the root PE
1865 integer :: outunit ! The output unit to write to
1866
1867 root = is_root_pe()
1868 outunit = stdout_if_root()
1869
1870 if (root) write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
1871 chks = field_chksum( iobt%u_flux ) ; if (root) write(outunit,100) 'iobt%u_flux ', chks
1872 chks = field_chksum( iobt%v_flux ) ; if (root) write(outunit,100) 'iobt%v_flux ', chks
1873 chks = field_chksum( iobt%t_flux ) ; if (root) write(outunit,100) 'iobt%t_flux ', chks
1874 chks = field_chksum( iobt%q_flux ) ; if (root) write(outunit,100) 'iobt%q_flux ', chks
1875 chks = field_chksum( iobt%salt_flux ) ; if (root) write(outunit,100) 'iobt%salt_flux ', chks
1876 chks = field_chksum( iobt%lw_flux ) ; if (root) write(outunit,100) 'iobt%lw_flux ', chks
1877 chks = field_chksum( iobt%sw_flux_vis_dir) ; if (root) write(outunit,100) 'iobt%sw_flux_vis_dir', chks
1878 chks = field_chksum( iobt%sw_flux_vis_dif) ; if (root) write(outunit,100) 'iobt%sw_flux_vis_dif', chks
1879 chks = field_chksum( iobt%sw_flux_nir_dir) ; if (root) write(outunit,100) 'iobt%sw_flux_nir_dir', chks
1880 chks = field_chksum( iobt%sw_flux_nir_dif) ; if (root) write(outunit,100) 'iobt%sw_flux_nir_dif', chks
1881 chks = field_chksum( iobt%lprec ) ; if (root) write(outunit,100) 'iobt%lprec ', chks
1882 chks = field_chksum( iobt%fprec ) ; if (root) write(outunit,100) 'iobt%fprec ', chks
1883 chks = field_chksum( iobt%runoff ) ; if (root) write(outunit,100) 'iobt%runoff ', chks
1884 chks = field_chksum( iobt%calving ) ; if (root) write(outunit,100) 'iobt%calving ', chks
1885 chks = field_chksum( iobt%p ) ; if (root) write(outunit,100) 'iobt%p ', chks
1886 if (associated(iobt%shelf_sfc_mass_flux)) then
1887 chks = field_chksum( iobt%shelf_sfc_mass_flux ) ; if (root) write(outunit,100) 'iobt%shelf_sfc_mass_flux ',&
1888 chks
1889 endif
1890 if (associated(iobt%ustar_berg)) then
1891 chks = field_chksum( iobt%ustar_berg ) ; if (root) write(outunit,100) 'iobt%ustar_berg ', chks
1892 endif
1893 if (associated(iobt%area_berg)) then
1894 chks = field_chksum( iobt%area_berg ) ; if (root) write(outunit,100) 'iobt%area_berg ', chks
1895 endif
1896 if (associated(iobt%mass_berg)) then
1897 chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks
1898 endif
1899 if (associated(iobt%excess_salt)) then
1900 chks = field_chksum( iobt%excess_salt ) ; if (root) write(outunit,100) 'iobt%excess_salt ', chks
1901 endif
1902100 FORMAT(" CHECKSUM::",a20," = ",z20)
1903
1904 call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')
1905
1906end subroutine ice_ocn_bnd_type_chksum
1907
1908!> Check the values passed by IOB over land are zero
1909subroutine check_mask_val_consistency(val, mask, i, j, varname, G)
1910
1911 real, intent(in) :: val !< value of flux/variable passed by IOB [various]
1912 real, intent(in) :: mask !< value of ocean mask [nondim]
1913 integer, intent(in) :: i !< model grid cell indices
1914 integer, intent(in) :: j !< model grid cell indices
1915 character(len=*), intent(in) :: varname !< variable name
1916 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1917 ! Local variables
1918 character(len=48) :: ci, cj !< model local grid cell indices as strings
1919 character(len=48) :: ciglo, cjglo !< model global grid cell indices as strings
1920 character(len=48) :: cval !< value to be displayed
1921 character(len=256) :: error_message !< error message to be displayed
1922
1923 if ((mask == 0.) .and. (val /= 0.)) then
1924 write(ci, '(I8)') i
1925 write(cj, '(I8)') j
1926 write(ciglo, '(I8)') i + g%HI%idg_offset
1927 write(cjglo, '(I8)') j + g%HI%jdg_offset
1928 write(cval, '(E22.16)') val
1929 error_message = "MOM_surface_forcing: found non-zero value (="//trim(cval)//") over land "//&
1930 "for variable "//trim(varname)//" at local point (i, j) = ("//trim(ci)//", "//trim(cj)//&
1931 ", global point (iglo, jglo) = ("//trim(ciglo)//", "//trim(cjglo)//")"
1932 call mom_error(warning, error_message)
1933 endif
1934
1935end subroutine
1936
1937end module mom_surface_forcing_gfdl