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