MOM6
mom_surface_forcing_mct.F90
Go to the documentation of this file.
2 
3 ! This file is part of MOM6. See LICENSE.md for the license.
4 
5 use mom_coms, only : reproducing_sum
6 use mom_constants, only : hlv, hlf
7 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
8 use mom_cpu_clock, only : clock_subcomponent
10 use mom_diag_mediator, only : safe_alloc_ptr, time_type
12 use mom_domains, only : global_field_sum, bitwise_exact_sum
13 use mom_domains, only : agrid, bgrid_ne, cgrid_ne, to_all
14 use mom_domains, only : to_north, to_east, omit_corners
15 use mom_error_handler, only : mom_error, warning, fatal, is_root_pe, mom_mesg
22 use mom_grid, only : ocean_grid_type
23 use mom_io, only : slasher, write_version_number, mom_read_data
29 use mom_variables, only : surface
32 
33 use coupler_types_mod, only : coupler_2d_bc_type, coupler_type_write_chksums
34 use coupler_types_mod, only : coupler_type_initialized, coupler_type_spawn
35 use coupler_types_mod, only : coupler_type_copy_data
36 use data_override_mod, only : data_override_init, data_override
37 use fms_mod, only : stdout
38 use mpp_mod, only : mpp_chksum
39 use time_interp_external_mod, only : init_external_field, time_interp_external
40 use time_interp_external_mod, only : time_interp_external_init
41 
42 implicit none ; private
43 
44 #include <MOM_memory.h>
45 
51 
54 private surface_forcing_end
55 
56 !> Contains pointers to the forcing fields which may be used to drive MOM.
57 !! All fluxes are positive downward.
58 type, public :: surface_forcing_cs ; private
59  integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integer values
60  !! from MOM_domains) to indicate the staggering of
61  !! the winds that are being provided in calls to
62  !! update_ocean_model.
63  logical :: use_temperature !! If true, temp and saln used as state variables
64  real :: wind_stress_multiplier!< A multiplier applied to incoming wind stress (nondim).
65 
66  real :: rho0 !< Boussinesq reference density [R ~> kg m-3]
67  real :: area_surf = -1.0 !< total ocean surface area [m2]
68  real :: latent_heat_fusion !< latent heat of fusion [J kg-1]
69  real :: latent_heat_vapor !< latent heat of vaporization [J kg-1]
70 
71  real :: max_p_surf !< maximum surface pressure that can be
72  !! exerted by the atmosphere and floating sea-ice,
73  !! [Pa]. This is needed because the FMS coupling
74  !! structure does not limit the water that can be
75  !! frozen out of the ocean and the ice-ocean heat
76  !! fluxes are treated explicitly.
77  logical :: use_limited_p_ssh !< If true, return the sea surface height with
78  !! the correction for the atmospheric (and sea-ice)
79  !! pressure limited by max_p_surf instead of the
80  !! full atmospheric pressure. The default is true.
81  real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
82  logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied
83  !! from an input file.
84  real, pointer, dimension(:,:) :: &
85  tke_tidal => null(), & !< turbulent kinetic energy introduced to the
86  !! bottom boundary layer by drag on the tidal flows [R Z3 T-3 ~> W m-2]
87  gust => null(), & !< spatially varying unresolved background
88  !! gustiness that contributes to ustar [R L Z T-1 ~> Pa].
89  !! gust is used when read_gust_2d is true.
90  ustar_tidal => null() !< tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1]
91  real :: cd_tides !< drag coefficient that applies to the tides (nondimensional)
92  real :: utide !< constant tidal velocity to use if read_tideamp
93  !! is false [Z T-1 ~> m s-1]
94  logical :: read_tideamp !< If true, spatially varying tidal amplitude read from a file.
95  logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts
96  !! to damp surface deflections (especially surface
97  !! gravity waves). The default is false.
98  real :: g_earth !< Gravitational acceleration [m s-2]
99  real :: kv_sea_ice !! viscosity in sea-ice that resists sheared vertical motions [m2 s-1]
100  real :: density_sea_ice !< typical density of sea-ice [kg m-3]. The value is
101  !! only used to convert the ice pressure into
102  !! appropriate units for use with Kv_sea_ice.
103  real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which
104  !! sea-ice viscosity becomes effective, in kg m-2,
105  !! typically of order 1000 [kg m-2].
106  logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments
107  real :: flux_const !< piston velocity for surface restoring [m/s]
108  logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux
109  logical :: adjust_net_srestore_to_zero !< adjust srestore to zero (for both salt_flux or vprec)
110  logical :: adjust_net_srestore_by_scaling !< adjust srestore w/o moving zero contour
111  logical :: adjust_net_fresh_water_to_zero !< adjust net surface fresh-water (w/ restoring) to zero
112  logical :: use_net_fw_adjustment_sign_bug !< use the wrong sign when adjusting net FW
113  logical :: adjust_net_fresh_water_by_scaling !< adjust net surface fresh-water w/o moving zero contour
114  logical :: mask_srestore_under_ice !< If true, use an ice mask defined by frazil
115  !< criteria for salinity restoring.
116  real :: ice_salt_concentration !< salt concentration for sea ice [kg/kg]
117  logical :: mask_srestore_marginal_seas !< if true, then mask SSS restoring in marginal seas
118  real :: max_delta_srestore !< maximum delta salinity used for restoring
119  real :: max_delta_trestore !< maximum delta sst used for restoring
120  real, pointer, dimension(:,:) :: basin_mask => null() !< mask for SSS restoring by basin
121  type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing
122  character(len=200) :: inputdir !< directory where NetCDF input files are
123  character(len=200) :: salt_restore_file !< filename for salt restoring data
124  character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file
125  logical :: mask_srestore !< if true, apply a 2-dimensional mask to the surface
126  !< salinity restoring fluxes. The masking file should be
127  !< in inputdir/salt_restore_mask.nc and the field should
128  !! be named 'mask'
129  real, pointer, dimension(:,:) :: srestore_mask => null() !< mask for SSS restoring
130  character(len=200) :: temp_restore_file !< filename for sst restoring data
131  character(len=30) :: temp_restore_var_name !< name of surface temperature in temp_restore_file
132  logical :: mask_trestore !< if true, apply a 2-dimensional mask to the surface
133  !! temperature restoring fluxes. The masking file should be
134  !! in inputdir/temp_restore_mask.nc and the field should
135  !! be named 'mask'
136  real, pointer, dimension(:,:) :: trestore_mask => null() !< mask for SST restoring
137  integer :: id_srestore = -1 !< id number for time_interp_external.
138  integer :: id_trestore = -1 !< id number for time_interp_external.
139 
140  type(forcing_diags), public :: handles !< diagnostics handles
141  type(mom_restart_cs), pointer :: restart_csp => null() !< restart pointer
142  type(user_revise_forcing_cs), pointer :: urf_cs => null() !< user revise pointer
143 end type surface_forcing_cs
144 
145 !> Structure corresponding to forcing, but with the elements, units, and conventions
146 !! that exactly conform to the use for MOM-based coupled models.
147 type, public :: ice_ocean_boundary_type
148  real, pointer, dimension(:,:) :: rofl_flux =>null() !< liquid runoff [kg/m2/s]
149  real, pointer, dimension(:,:) :: rofi_flux =>null() !< ice runoff [kg/m2/s]
150  real, pointer, dimension(:,:) :: u_flux =>null() !< i-direction wind stress [Pa]
151  real, pointer, dimension(:,:) :: v_flux =>null() !< j-direction wind stress [Pa]
152  real, pointer, dimension(:,:) :: t_flux =>null() !< sensible heat flux [W/m2]
153  real, pointer, dimension(:,:) :: q_flux =>null() !< specific humidity flux [kg/m2/s]
154  real, pointer, dimension(:,:) :: salt_flux =>null() !< salt flux [kg/m2/s]
155  real, pointer, dimension(:,:) :: seaice_melt_heat =>null() !< sea ice and snow melt heat flux [W/m2]
156  real, pointer, dimension(:,:) :: seaice_melt =>null() !< water flux due to sea ice and snow melting [kg/m2/s]
157  real, pointer, dimension(:,:) :: lw_flux =>null() !< long wave radiation [W/m2]
158  real, pointer, dimension(:,:) :: sw_flux_vis_dir =>null() !< direct visible sw radiation [W/m2]
159  real, pointer, dimension(:,:) :: sw_flux_vis_dif =>null() !< diffuse visible sw radiation [W/m2]
160  real, pointer, dimension(:,:) :: sw_flux_nir_dir =>null() !< direct Near InfraRed sw radiation [W/m2]
161  real, pointer, dimension(:,:) :: sw_flux_nir_dif =>null() !< diffuse Near InfraRed sw radiation [W/m2]
162  real, pointer, dimension(:,:) :: lprec =>null() !< mass flux of liquid precip [kg/m2/s]
163  real, pointer, dimension(:,:) :: fprec =>null() !< mass flux of frozen precip [kg/m2/s]
164  real, pointer, dimension(:,:) :: runoff =>null() !< mass flux of liquid runoff [kg/m2/s]
165  real, pointer, dimension(:,:) :: calving =>null() !< mass flux of frozen runoff [kg/m2/s]
166  real, pointer, dimension(:,:) :: ustar_berg =>null() !< frictional velocity beneath icebergs [m/s]
167  real, pointer, dimension(:,:) :: area_berg =>null() !< area covered by icebergs[m2/m2]
168  real, pointer, dimension(:,:) :: mass_berg =>null() !< mass of icebergs(kg/m2)
169  real, pointer, dimension(:,:) :: runoff_hflx =>null() !< heat content of liquid runoff [W/m2]
170  real, pointer, dimension(:,:) :: calving_hflx =>null() !< heat content of frozen runoff [W/m2]
171  real, pointer, dimension(:,:) :: p =>null() !< pressure of overlying ice and atmosphere
172  !< on ocean surface [Pa]
173  real, pointer, dimension(:,:) :: mi =>null() !< mass of ice [kg/m2]
174  real, pointer, dimension(:,:) :: ice_rigidity =>null() !< rigidity of the sea ice, sea-ice and
175  !! ice-shelves, expressed as a coefficient
176  !! for divergence damping, as determined
177  !! outside of the ocean model in [m3/s]
178  integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT
179  type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of
180  !! named fields used for passive tracer fluxes.
181  integer :: wind_stagger = -999 !< A flag indicating the spatial discretization of
182  !! wind stresses. This flag may be set by the
183  !! flux-exchange code, based on what the sea-ice
184  !! model is providing. Otherwise, the value from
185  !! the surface_forcing_CS is used.
187 
189 
190 contains
191 
192 !> This subroutine translates the Ice_ocean_boundary_type into a MOM
193 !! thermodynamic forcing type, including changes of units, sign conventions,
194 !! and putting the fields into arrays with MOM-standard halos.
195 subroutine convert_iob_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, &
196  sfc_state, restore_salt, restore_temp)
197 
198  type(ice_ocean_boundary_type), &
199  target, intent(in) :: iob !< An ice-ocean boundary type with fluxes to drive
200  !! the ocean in a coupled model
201 
202  type(forcing), intent(inout) :: fluxes !< A structure containing pointers to all
203  !! possible mass, heat or salt flux forcing fields.
204  !! Unused fields have NULL ptrs.
205  integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
206  type(time_type), intent(in) :: time !< The time of the fluxes, used for interpolating the
207  !! salinity to the right time, when it is being restored.
208  real, intent(in) :: valid_time !< The amount of time over which these fluxes
209  !! should be applied [s].
210  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
211  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
212  type(surface_forcing_cs),pointer :: cs !< A pointer to the control structure returned by a
213  !! previous call to surface_forcing_init.
214  type(surface), intent(in) :: sfc_state !< A structure containing fields that describe the
215  !! surface state of the ocean.
216  logical, optional, intent(in) :: restore_salt !< If true, salinity is restored to a target value.
217  logical, optional, intent(in) :: restore_temp !< If true, temperature is restored to a target value.
218 
219  ! local variables
220  real, dimension(SZI_(G),SZJ_(G)) :: &
221  data_restore, & !< The surface value toward which to restore [g/kg or degC]
222  sst_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C]
223  sss_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg]
224  sss_mean, & !< A (mean?) salinity about which to normalize local salinity
225  !! anomalies when calculating restorative precipitation anomalies [g/kg]
226  pme_adj, & !< The adjustment to PminusE that will cause the salinity
227  !! to be restored toward its target value [kg/(m^2 * s)]
228  net_fw, & !< The area integrated net freshwater flux into the ocean [kg/s]
229  net_fw2, & !< The area integrated net freshwater flux into the ocean [kg/s]
230  work_sum, & !< A 2-d array that is used as the work space for a global
231  !! sum, used with units of m2 or [kg/s]
232  open_ocn_mask !< a binary field indicating where ice is present based on frazil criteria
233 
234  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
235  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
236  integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
237 
238  logical :: restore_salinity !< local copy of the argument restore_salt, if it
239  !! is present, or false (no restoring) otherwise.
240  logical :: restore_sst !< local copy of the argument restore_temp, if it
241  !! is present, or false (no restoring) otherwise.
242  real :: delta_sss !< temporary storage for sss diff from restoring value
243  real :: delta_sst !< temporary storage for sst diff from restoring value
244 
245  real :: kg_m2_s_conversion !< A combination of unit conversion factors for rescaling
246  !! mass fluxes [R Z s m2 kg-1 T-1 ~> 1].
247  real :: c_p !< heat capacity of seawater ( J/(K kg) )
248  real :: sign_for_net_fw_bug !< Should be +1. but an old bug can be recovered by using -1.
249 
250  call cpu_clock_begin(id_clock_forcing)
251 
252  isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
253  jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
254  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
255  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
256  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
257  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
258  isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
259 
260  kg_m2_s_conversion = us%kg_m3_to_R*us%m_to_Z*us%T_to_s
261  c_p = fluxes%C_p
262  open_ocn_mask(:,:) = 1.0
263  pme_adj(:,:) = 0.0
264  fluxes%vPrecGlobalAdj = 0.0
265  fluxes%vPrecGlobalScl = 0.0
266  fluxes%saltFluxGlobalAdj = 0.0
267  fluxes%saltFluxGlobalScl = 0.0
268  fluxes%netFWGlobalAdj = 0.0
269  fluxes%netFWGlobalScl = 0.0
270 
271  restore_salinity = .false.
272  if (present(restore_salt)) restore_salinity = restore_salt
273  restore_sst = .false.
274  if (present(restore_temp)) restore_sst = restore_temp
275 
276  ! allocation and initialization if this is the first time that this
277  ! flux type has been used.
278  if (fluxes%dt_buoy_accum < 0) then
279  call allocate_forcing_type(g, fluxes, water=.true., heat=.true., &
280  ustar=.true., press=.true.)
281 
282  call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
283  call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
284  call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed)
285  call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed)
286 
287  call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed)
288  call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed)
289  if (cs%use_limited_P_SSH) then
290  fluxes%p_surf_SSH => fluxes%p_surf
291  else
292  fluxes%p_surf_SSH => fluxes%p_surf_full
293  endif
294 
295  call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed)
296  call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed)
297  call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
298 
299  call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed)
300  call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed)
301 
302  if (cs%allow_flux_adjustments) then
303  call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
304  call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
305  endif
306 
307  do j=js-2,je+2 ; do i=is-2,ie+2
308  fluxes%TKE_tidal(i,j) = cs%TKE_tidal(i,j)
309  fluxes%ustar_tidal(i,j) = cs%ustar_tidal(i,j)
310  enddo; enddo
311 
312  if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
313 
314  endif ! endif for allocation and initialization
315 
316 
317  if (((associated(iob%ustar_berg) .and. (.not.associated(fluxes%ustar_berg))) &
318  .or. (associated(iob%area_berg) .and. (.not.associated(fluxes%area_berg)))) &
319  .or. (associated(iob%mass_berg) .and. (.not.associated(fluxes%mass_berg)))) &
320  call allocate_forcing_type(g, fluxes, iceberg=.true.)
321 
322  if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. &
323  coupler_type_initialized(iob%fluxes)) &
324  call coupler_type_spawn(iob%fluxes, fluxes%tr_fluxes, &
325  (/is,is,ie,ie/), (/js,js,je,je/))
326  ! It might prove valuable to use the same array extents as the rest of the
327  ! ocean model, rather than using haloless arrays, in which case the last line
328  ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/))
329 
330  ! allocation and initialization on first call to this routine
331  if (cs%area_surf < 0.0) then
332  do j=js,je ; do i=is,ie
333  work_sum(i,j) = us%L_to_m**2*g%areaT(i,j) * g%mask2dT(i,j)
334  enddo; enddo
335  cs%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer)
336  endif ! endif for allocation and initialization
337 
338 
339  ! Indicate that there are new unused fluxes.
340  fluxes%fluxes_used = .false.
341  fluxes%dt_buoy_accum = us%s_to_T*valid_time
342 
343  if (cs%allow_flux_adjustments) then
344  fluxes%heat_added(:,:)=0.0
345  fluxes%salt_flux_added(:,:)=0.0
346  endif
347 
348  do j=js,je ; do i=is,ie
349  fluxes%salt_flux(i,j) = 0.0
350  fluxes%vprec(i,j) = 0.0
351  enddo; enddo
352 
353  ! Salinity restoring logic
354  if (restore_salinity) then
355  call time_interp_external(cs%id_srestore,time,data_restore)
356  ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not)
357  open_ocn_mask(:,:) = 1.0
358  if (cs%mask_srestore_under_ice) then ! Do not restore under sea-ice
359  do j=js,je ; do i=is,ie
360  if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
361  enddo; enddo
362  endif
363  if (cs%salt_restore_as_sflux) then
364  do j=js,je ; do i=is,ie
365  delta_sss = data_restore(i,j)- sfc_state%SSS(i,j)
366  delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
367  fluxes%salt_flux(i,j) = 1.e-3*g%mask2dT(i,j) * (cs%Rho0*us%m_to_Z*us%T_to_s*cs%Flux_const)* &
368  (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
369  enddo; enddo
370  if (cs%adjust_net_srestore_to_zero) then
371  if (cs%adjust_net_srestore_by_scaling) then
372  call adjust_area_mean_to_zero(fluxes%salt_flux, g, fluxes%saltFluxGlobalScl, &
373  unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
374  fluxes%saltFluxGlobalAdj = 0.
375  else
376  work_sum(is:ie,js:je) = us%L_to_m**2*us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
377  g%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je)
378  fluxes%saltFluxGlobalAdj = reproducing_sum(work_sum(:,:), isr,ier, jsr,jer)/cs%area_surf
379  fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - kg_m2_s_conversion * fluxes%saltFluxGlobalAdj
380  endif
381  endif
382  fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic
383  else
384  do j=js,je ; do i=is,ie
385  if (g%mask2dT(i,j) > 0.5) then
386  delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
387  delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
388  fluxes%vprec(i,j) = (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j))* &
389  (us%m_to_Z*us%T_to_s * cs%Rho0*cs%Flux_const) * &
390  delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j)))
391  endif
392  enddo; enddo
393  if (cs%adjust_net_srestore_to_zero) then
394  if (cs%adjust_net_srestore_by_scaling) then
395  call adjust_area_mean_to_zero(fluxes%vprec, g, fluxes%vPrecGlobalScl, &
396  unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
397  fluxes%vPrecGlobalAdj = 0.
398  else
399  work_sum(is:ie,js:je) = us%L_to_m**2*g%areaT(is:ie,js:je) * &
400  us%R_to_kg_m3*us%Z_to_m*us%s_to_T*fluxes%vprec(is:ie,js:je)
401  fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / cs%area_surf
402  do j=js,je ; do i=is,ie
403  fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion*fluxes%vPrecGlobalAdj ) * g%mask2dT(i,j)
404  enddo; enddo
405  endif
406  endif
407  endif
408  endif
409 
410  ! SST restoring logic
411  if (restore_sst) then
412  call time_interp_external(cs%id_trestore,time,data_restore)
413  do j=js,je ; do i=is,ie
414  delta_sst = data_restore(i,j)- sfc_state%SST(i,j)
415  delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),cs%max_delta_trestore)
416  fluxes%heat_added(i,j) = g%mask2dT(i,j) * cs%trestore_mask(i,j) * &
417  (us%R_to_kg_m3*cs%Rho0*fluxes%C_p) * delta_sst * cs%Flux_const ! W m-2
418  enddo; enddo
419  endif
420 
421  ! obtain fluxes from IOB; note the staggering of indices
422  i0 = 0; j0 = 0
423  do j=js,je ; do i=is,ie
424  ! liquid precipitation (rain)
425  if (associated(iob%lprec)) &
426  fluxes%lprec(i,j) = kg_m2_s_conversion * iob%lprec(i-i0,j-j0) * g%mask2dT(i,j)
427 
428  ! frozen precipitation (snow)
429  if (associated(iob%fprec)) &
430  fluxes%fprec(i,j) = kg_m2_s_conversion * iob%fprec(i-i0,j-j0) * g%mask2dT(i,j)
431 
432  ! evaporation
433  if (associated(iob%q_flux)) &
434  fluxes%evap(i,j) = kg_m2_s_conversion * iob%q_flux(i-i0,j-j0) * g%mask2dT(i,j)
435 
436  ! liquid runoff flux
437  if (associated(iob%rofl_flux)) then
438  fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%rofl_flux(i-i0,j-j0) * g%mask2dT(i,j)
439  else if (associated(iob%runoff)) then
440  fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%runoff(i-i0,j-j0) * g%mask2dT(i,j)
441  end if
442 
443  ! ice runoff flux
444  if (associated(iob%rofi_flux)) then
445  fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%rofi_flux(i-i0,j-j0) * g%mask2dT(i,j)
446  else if (associated(iob%calving)) then
447  fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%calving(i-i0,j-j0) * g%mask2dT(i,j)
448  end if
449 
450  if (associated(iob%ustar_berg)) &
451  fluxes%ustar_berg(i,j) = us%m_to_Z * iob%ustar_berg(i-i0,j-j0) * g%mask2dT(i,j)
452 
453  if (associated(iob%area_berg)) &
454  fluxes%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
455 
456  if (associated(iob%mass_berg)) &
457  fluxes%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
458 
459  ! GMM, cime does not not have an equivalent for heat_content_lrunoff and
460  ! heat_content_frunoff. I am seeting these to zero for now.
461  if (associated(fluxes%heat_content_lrunoff)) &
462  fluxes%heat_content_lrunoff(i,j) = 0.0 * g%mask2dT(i,j)
463  if (associated(fluxes%heat_content_frunoff)) &
464  fluxes%heat_content_frunoff(i,j) = 0.0 * g%mask2dT(i,j)
465 
466  if (associated(iob%calving_hflx)) &
467  fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * iob%calving_hflx(i-i0,j-j0) * g%mask2dT(i,j)
468 
469  ! longwave radiation, sum up and down (W/m2)
470  if (associated(iob%lw_flux)) &
471  fluxes%LW(i,j) = iob%lw_flux(i-i0,j-j0) * g%mask2dT(i,j)
472 
473  ! sensible heat flux (W/m2)
474  if (associated(iob%t_flux)) &
475  fluxes%sens(i,j) = iob%t_flux(i-i0,j-j0) * g%mask2dT(i,j)
476 
477  ! sea ice and snow melt heat flux [W/m2]
478  if (associated(iob%seaice_melt_heat)) &
479  fluxes%seaice_melt_heat(i,j) = g%mask2dT(i,j) * iob%seaice_melt_heat(i-i0,j-j0)
480 
481  ! water flux due to sea ice and snow melt [kg/m2/s]
482  if (associated(iob%seaice_melt)) &
483  fluxes%seaice_melt(i,j) = g%mask2dT(i,j) * kg_m2_s_conversion * iob%seaice_melt(i-i0,j-j0)
484 
485  ! latent heat flux (W/m^2)
486  fluxes%latent(i,j) = 0.0
487  ! contribution from frozen ppt
488  if (associated(iob%fprec)) then
489  fluxes%latent(i,j) = fluxes%latent(i,j) + iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
490  fluxes%latent_fprec_diag(i,j) = g%mask2dT(i,j) * iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
491  endif
492  ! contribution from frozen runoff
493  if (associated(fluxes%frunoff)) then
494  fluxes%latent(i,j) = fluxes%latent(i,j) + iob%rofi_flux(i-i0,j-j0)*cs%latent_heat_fusion
495  fluxes%latent_frunoff_diag(i,j) = g%mask2dT(i,j) * iob%rofi_flux(i-i0,j-j0)*cs%latent_heat_fusion
496  endif
497  ! contribution from evaporation
498  if (associated(iob%q_flux)) then
499  fluxes%latent(i,j) = fluxes%latent(i,j) + iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
500  fluxes%latent_evap_diag(i,j) = g%mask2dT(i,j) * iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
501  endif
502  fluxes%latent(i,j) = g%mask2dT(i,j) * fluxes%latent(i,j)
503 
504  if (associated(iob%sw_flux_vis_dir)) &
505  fluxes%sw_vis_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dir(i-i0,j-j0)
506 
507  if (associated(iob%sw_flux_vis_dif)) &
508  fluxes%sw_vis_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dif(i-i0,j-j0)
509 
510  if (associated(iob%sw_flux_nir_dir)) &
511  fluxes%sw_nir_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dir(i-i0,j-j0)
512 
513  if (associated(iob%sw_flux_nir_dif)) &
514  fluxes%sw_nir_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dif(i-i0,j-j0)
515 
516  fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
517  fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)
518 
519  enddo; enddo
520 
521  ! applied surface pressure from atmosphere and cryosphere
522  if (associated(iob%p)) then
523  if (cs%max_p_surf >= 0.0) then
524  do j=js,je ; do i=is,ie
525  fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
526  fluxes%p_surf(i,j) = min(fluxes%p_surf_full(i,j),cs%max_p_surf)
527  enddo; enddo
528  else
529  do j=js,je ; do i=is,ie
530  fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
531  fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
532  enddo; enddo
533  endif
534  fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
535  endif
536 
537  if (associated(iob%salt_flux)) then
538  do j=js,je ; do i=is,ie
539  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))
540  fluxes%salt_flux_in(i,j) = g%mask2dT(i,j)*( kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0) )
541  enddo ; enddo
542  endif
543 
544  ! adjust the NET fresh-water flux to zero, if flagged
545  if (cs%adjust_net_fresh_water_to_zero) then
546  sign_for_net_fw_bug = 1.
547  if (cs%use_net_FW_adjustment_sign_bug) sign_for_net_fw_bug = -1.
548  do j=js,je ; do i=is,ie
549  net_fw(i,j) = us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
550  (((fluxes%lprec(i,j) + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
551  (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
552  (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * us%L_to_m**2*g%areaT(i,j)
553 
554  net_fw2(i,j) = net_fw(i,j) / (us%L_to_m**2*g%areaT(i,j))
555  enddo; enddo
556 
557  if (cs%adjust_net_fresh_water_by_scaling) then
558  call adjust_area_mean_to_zero(net_fw2, g, fluxes%netFWGlobalScl)
559  do j=js,je ; do i=is,ie
560  fluxes%vprec(i,j) = fluxes%vprec(i,j) + kg_m2_s_conversion * &
561  (net_fw2(i,j) - net_fw(i,j)/(us%L_to_m**2*g%areaT(i,j))) * g%mask2dT(i,j)
562  enddo; enddo
563  else
564  fluxes%netFWGlobalAdj = reproducing_sum(net_fw(:,:), isr, ier, jsr, jer) / cs%area_surf
565  do j=js,je ; do i=is,ie
566  fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion * fluxes%netFWGlobalAdj ) * g%mask2dT(i,j)
567  enddo; enddo
568  endif
569  endif
570 
571  if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
572  coupler_type_initialized(iob%fluxes)) &
573  call coupler_type_copy_data(iob%fluxes, fluxes%tr_fluxes)
574 
575  if (cs%allow_flux_adjustments) then
576  ! Apply adjustments to fluxes
577  call apply_flux_adjustments(g, us, cs, time, fluxes)
578  endif
579 
580  ! Allow for user-written code to alter fluxes after all the above
581  call user_alter_forcing(sfc_state, fluxes, time, g, cs%urf_CS)
582 
583  call cpu_clock_end(id_clock_forcing)
584 
585 end subroutine convert_iob_to_fluxes
586 
587 !> This subroutine translates the Ice_ocean_boundary_type into a MOM
588 !! mechanical forcing type, including changes of units, sign conventions,
589 !! and putting the fields into arrays with MOM-standard halos.
590 subroutine convert_iob_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
592  target, intent(in) :: iob !< An ice-ocean boundary type with fluxes to drive
593  !! the ocean in a coupled model
594  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
595  integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB.
596  type(time_type), intent(in) :: time !< The time of the fluxes, used for interpolating the
597  !! salinity to the right time, when it is being restored.
598  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
599  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
600  type(surface_forcing_cs),pointer :: cs !< A pointer to the control structure returned by a
601  !! previous call to surface_forcing_init.
602 
603  ! local variables
604  real, dimension(SZIB_(G),SZJB_(G)) :: &
605  taux_at_q, & !< Zonal wind stresses at q points [R Z L T-2 ~> Pa]
606  tauy_at_q !< Meridional wind stresses at q points [R Z L T-2 ~> Pa]
607 
608  real, dimension(SZI_(G),SZJ_(G)) :: &
609  rigidity_at_h, & !< Ice rigidity at tracer points [m3 s-1]
610  taux_at_h, & !< Zonal wind stresses at h points [R Z L T-2 ~> Pa]
611  tauy_at_h !< Meridional wind stresses at h points [R Z L T-2 ~> Pa]
612 
613  real :: gustiness !< unresolved gustiness that contributes to ustar [R Z L T-2 ~> Pa]
614  real :: irho0 !< inverse of the mean density in [Z L-1 R-1 ~> m3 kg-1]
615  real :: taux2, tauy2 !< squared wind stresses [R2 Z2 L2 T-4 ~> Pa2]
616  real :: tau_mag !< magnitude of the wind stress [R Z L T-2 ~> Pa]
617  real :: pa_conversion ! A unit conversion factor from Pa to the internal wind stress units [R Z L T-2 Pa-1 ~> 1]
618  real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
619  real :: i_gearth !< 1.0 / G%G_Earth [s2 m-1]
620  real :: kv_rho_ice !< (CS%kv_sea_ice / CS%density_sea_ice) [m5 s-1 kg-1]
621  real :: mass_ice !< mass of sea ice at a face [kg m-2]
622  real :: mass_eff !< effective mass of sea ice for rigidity [kg m-2]
623 
624  integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
625  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
626  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
627  integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
628 
629  call cpu_clock_begin(id_clock_forcing)
630 
631  isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
632  jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
633  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
634  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
635  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
636  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
637  isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
638  !i0 = is - isc_bnd ; j0 = js - jsc_bnd
639  i0 = 0; j0 = 0
640 
641  irho0 = us%L_to_Z / cs%Rho0
642  pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
643  stress_conversion = pa_conversion * cs%wind_stress_multiplier
644 
645  ! allocation and initialization if this is the first time that this
646  ! mechanical forcing type has been used.
647  if (.not.forces%initialized) then
648 
649  call allocate_mech_forcing(g, forces, stress=.true., ustar=.true., press=.true.)
650 
651  call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
652  call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
653 
654  if (cs%rigid_sea_ice) then
655  call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
656  call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
657  endif
658 
659  forces%initialized = .true.
660  endif
661 
662  if ( (associated(iob%area_berg) .and. (.not. associated(forces%area_berg))) .or. &
663  (associated(iob%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
664  call allocate_mech_forcing(g, forces, iceberg=.true.)
665  if (associated(iob%ice_rigidity)) then
666  rigidity_at_h(:,:) = 0.0
667  call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
668  call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
669  endif
670 
671  forces%accumulate_rigidity = .true. ! Multiple components may contribute to rigidity.
672  if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
673  if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
674 
675  !applied surface pressure from atmosphere and cryosphere
676  if (cs%use_limited_P_SSH) then
677  forces%p_surf_SSH => forces%p_surf
678  else
679  forces%p_surf_SSH => forces%p_surf_full
680  endif
681  if (associated(iob%p)) then
682  if (cs%max_p_surf >= 0.0) then
683  do j=js,je ; do i=is,ie
684  forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
685  forces%p_surf(i,j) = min(forces%p_surf_full(i,j),cs%max_p_surf)
686  enddo ; enddo
687  else
688  do j=js,je ; do i=is,ie
689  forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
690  forces%p_surf(i,j) = forces%p_surf_full(i,j)
691  enddo ; enddo
692  endif
693  else
694  do j=js,je ; do i=is,ie
695  forces%p_surf_full(i,j) = 0.0
696  forces%p_surf(i,j) = 0.0
697  enddo ; enddo
698  endif
699  forces%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
700 
701  ! GMM, CIME uses AGRID. All the BGRID_NE code can be cleaned later
702  wind_stagger = agrid
703 
704  if (wind_stagger == bgrid_ne) then
705  ! This is necessary to fill in the halo points.
706  taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0
707  endif
708  if (wind_stagger == agrid) then
709  ! This is necessary to fill in the halo points.
710  taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0
711  endif
712 
713  ! obtain fluxes from IOB; note the staggering of indices
714  do j=js,je ; do i=is,ie
715  if (associated(iob%area_berg)) &
716  forces%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
717 
718  if (associated(iob%mass_berg)) &
719  forces%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
720 
721  if (associated(iob%ice_rigidity)) &
722  rigidity_at_h(i,j) = iob%ice_rigidity(i-i0,j-j0) * g%mask2dT(i,j)
723 
724  if (wind_stagger == bgrid_ne) then
725  if (associated(iob%u_flux)) taux_at_q(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
726  if (associated(iob%v_flux)) tauy_at_q(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
727  elseif (wind_stagger == agrid) then
728  if (associated(iob%u_flux)) taux_at_h(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
729  if (associated(iob%v_flux)) tauy_at_h(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
730  else ! C-grid wind stresses.
731  if (associated(iob%u_flux)) forces%taux(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
732  if (associated(iob%v_flux)) forces%tauy(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
733  endif
734 
735  enddo ; enddo
736 
737  ! surface momentum stress related fields as function of staggering
738  if (wind_stagger == bgrid_ne) then
739  if (g%symmetric) &
740  call fill_symmetric_edges(taux_at_q, tauy_at_q, g%Domain, stagger=bgrid_ne)
741  call pass_vector(taux_at_q, tauy_at_q, g%Domain, stagger=bgrid_ne, halo=1)
742 
743  do j=js,je ; do i=isq,ieq
744  forces%taux(i,j) = 0.0
745  if ((g%mask2dBu(i,j) + g%mask2dBu(i,j-1)) > 0) &
746  forces%taux(i,j) = (g%mask2dBu(i,j)*taux_at_q(i,j) + &
747  g%mask2dBu(i,j-1)*taux_at_q(i,j-1)) / &
748  (g%mask2dBu(i,j) + g%mask2dBu(i,j-1))
749  enddo; enddo
750 
751  do j=jsq,jeq ; do i=is,ie
752  forces%tauy(i,j) = 0.0
753  if ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j)) > 0) &
754  forces%tauy(i,j) = (g%mask2dBu(i,j)*tauy_at_q(i,j) + &
755  g%mask2dBu(i-1,j)*tauy_at_q(i-1,j)) / &
756  (g%mask2dBu(i,j) + g%mask2dBu(i-1,j))
757  enddo; enddo
758 
759  ! ustar is required for the bulk mixed layer formulation. The background value
760  ! of 0.02 Pa is a relatively small value intended to give reasonable behavior
761  ! in regions of very weak winds.
762 
763  do j=js,je ; do i=is,ie
764  tau_mag = 0.0 ; gustiness = cs%gust_const
765  if (((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + &
766  (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) > 0) then
767  tau_mag = sqrt(((g%mask2dBu(i,j)*(taux_at_q(i,j)**2 + tauy_at_q(i,j)**2) + &
768  g%mask2dBu(i-1,j-1)*(taux_at_q(i-1,j-1)**2 + tauy_at_q(i-1,j-1)**2)) + &
769  (g%mask2dBu(i,j-1)*(taux_at_q(i,j-1)**2 + tauy_at_q(i,j-1)**2) + &
770  g%mask2dBu(i-1,j)*(taux_at_q(i-1,j)**2 + tauy_at_q(i-1,j)**2)) ) / &
771  ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) )
772  if (cs%read_gust_2d) gustiness = cs%gust(i,j)
773  endif
774  forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0*tau_mag)
775  enddo; enddo
776 
777  elseif (wind_stagger == agrid) then
778  call pass_vector(taux_at_h, tauy_at_h, g%Domain, to_all+omit_corners, stagger=agrid, halo=1)
779 
780  do j=js,je ; do i=isq,ieq
781  forces%taux(i,j) = 0.0
782  if ((g%mask2dT(i,j) + g%mask2dT(i+1,j)) > 0) &
783  forces%taux(i,j) = (g%mask2dT(i,j)*taux_at_h(i,j) + &
784  g%mask2dT(i+1,j)*taux_at_h(i+1,j)) / &
785  (g%mask2dT(i,j) + g%mask2dT(i+1,j))
786  enddo; enddo
787 
788  do j=jsq,jeq ; do i=is,ie
789  forces%tauy(i,j) = 0.0
790  if ((g%mask2dT(i,j) + g%mask2dT(i,j+1)) > 0) &
791  forces%tauy(i,j) = (g%mask2dT(i,j)*tauy_at_h(i,j) + &
792  g%mask2dT(i,j+1)*tauy_at_h(i,j+1)) / &
793  (g%mask2dT(i,j) + g%mask2dT(i,j+1))
794  enddo; enddo
795 
796  do j=js,je ; do i=is,ie
797  gustiness = cs%gust_const
798  if (cs%read_gust_2d .and. (g%mask2dT(i,j) > 0)) gustiness = cs%gust(i,j)
799  forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0 * g%mask2dT(i,j) * &
800  sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
801  enddo; enddo
802 
803  else ! C-grid wind stresses.
804  if (g%symmetric) &
805  call fill_symmetric_edges(forces%taux, forces%tauy, g%Domain)
806  call pass_vector(forces%taux, forces%tauy, g%Domain, halo=1)
807 
808  do j=js,je ; do i=is,ie
809  taux2 = 0.0
810  if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0) &
811  taux2 = (g%mask2dCu(i-1,j)*forces%taux(i-1,j)**2 + &
812  g%mask2dCu(i,j)*forces%taux(i,j)**2) / (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
813 
814  tauy2 = 0.0
815  if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0) &
816  tauy2 = (g%mask2dCv(i,j-1)*forces%tauy(i,j-1)**2 + &
817  g%mask2dCv(i,j)*forces%tauy(i,j)**2) / (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
818 
819  if (cs%read_gust_2d) then
820  forces%ustar(i,j) = sqrt(cs%gust(i,j)*irho0 + irho0*sqrt(taux2 + tauy2))
821  else
822  forces%ustar(i,j) = sqrt(cs%gust_const*irho0 + irho0*sqrt(taux2 + tauy2))
823  endif
824  enddo; enddo
825 
826  endif ! endif for wind related fields
827 
828  ! sea ice related dynamic fields
829  if (associated(iob%ice_rigidity)) then
830  call pass_var(rigidity_at_h, g%Domain, halo=1)
831  do i=is-1,ie ; do j=js,je
832  forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
833  min(rigidity_at_h(i,j), rigidity_at_h(i+1,j))
834  enddo ; enddo
835  do i=is,ie ; do j=js-1,je
836  forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
837  min(rigidity_at_h(i,j), rigidity_at_h(i,j+1))
838  enddo ; enddo
839  endif
840 
841  if (cs%rigid_sea_ice) then
842  call pass_var(forces%p_surf_full, g%Domain, halo=1)
843  i_gearth = 1.0 / cs%G_Earth
844  kv_rho_ice = (cs%kv_sea_ice / cs%density_sea_ice)
845  do i=is-1,ie ; do j=js,je
846  mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * i_gearth
847  mass_eff = 0.0
848  if (mass_ice > cs%rigid_sea_ice_mass) then
849  mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
850  (mass_ice + cs%rigid_sea_ice_mass)
851  endif
852  forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * mass_eff
853  enddo ; enddo
854  do i=is,ie ; do j=js-1,je
855  mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * i_gearth
856  mass_eff = 0.0
857  if (mass_ice > cs%rigid_sea_ice_mass) then
858  mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
859  (mass_ice + cs%rigid_sea_ice_mass)
860  endif
861  forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * mass_eff
862  enddo ; enddo
863  endif
864 
865  if (cs%allow_flux_adjustments) then
866  ! Apply adjustments to forces
867  call apply_force_adjustments(g, us, cs, time, forces)
868  endif
869 
870 !### ! Allow for user-written code to alter fluxes after all the above
871 !### call user_alter_mech_forcing(forces, Time, G, CS%urf_CS)
872 
873  call cpu_clock_end(id_clock_forcing)
874 end subroutine convert_iob_to_forces
875 
876 !> Adds thermodynamic flux adjustments obtained via data_override
877 !! Component name is 'OCN'
878 !! Available adjustments are:
879 !! - hflx_adj (Heat flux into the ocean, in W m-2)
880 !! - sflx_adj (Salt flux into the ocean, in kg salt m-2 s-1)
881 !! - prcme_adj (Fresh water flux into the ocean, in kg m-2 s-1)
882 subroutine apply_flux_adjustments(G, US, CS, Time, fluxes)
883  type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
884  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
885  type(surface_forcing_cs), pointer :: cs !< Surface forcing control structure
886  type(time_type), intent(in) :: time !< Model time structure
887  type(forcing), intent(inout) :: fluxes !< Surface fluxes structure
888 
889  ! Local variables
890  real, dimension(SZI_(G),SZJ_(G)) :: temp_at_h !< Fluxes at h points [W m-2 or kg m-2 s-1]
891 
892  integer :: isc, iec, jsc, jec, i, j
893  logical :: overrode_h
894 
895  isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
896 
897  overrode_h = .false.
898  call data_override('OCN', 'hflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
899 
900  if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
901  fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j)* g%mask2dT(i,j)
902  enddo ; enddo ; endif
903  ! Not needed? ! if (overrode_h) call pass_var(fluxes%heat_added, G%Domain)
904 
905  overrode_h = .false.
906  call data_override('OCN', 'sflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
907 
908  if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
909  fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + &
910  us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
911  enddo ; enddo ; endif
912  ! Not needed? ! if (overrode_h) call pass_var(fluxes%salt_flux_added, G%Domain)
913 
914  overrode_h = .false.
915  call data_override('OCN', 'prcme_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
916 
917  if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
918  fluxes%vprec(i,j) = fluxes%vprec(i,j) + us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
919  enddo ; enddo ; endif
920  ! Not needed? ! if (overrode_h) call pass_var(fluxes%vprec, G%Domain)
921 
922 end subroutine apply_flux_adjustments
923 
924 !> Adds mechanical forcing adjustments obtained via data_override
925 !! Component name is 'OCN'
926 !! Available adjustments are:
927 !! - taux_adj (Zonal wind stress delta, positive to the east, in Pa)
928 !! - tauy_adj (Meridional wind stress delta, positive to the north, in Pa)
929 subroutine apply_force_adjustments(G, US, CS, Time, forces)
930  type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
931  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
932  type(surface_forcing_cs), pointer :: cs !< Surface forcing control structure
933  type(time_type), intent(in) :: time !< Model time structure
934  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
935 
936  ! Local variables
937  real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h !< Delta to zonal wind stress at h points [R Z L T-2 ~> Pa]
938  real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h !< Delta to meridional wind stress at h points [R Z L T-2 ~> Pa]
939 
940  integer :: isc, iec, jsc, jec, i, j
941  real :: dlondx, dlondy, rdlon, cosa, sina, zonal_tau, merid_tau
942  real :: pa_conversion ! A unit conversion factor from Pa to the internal units [R Z L T-2 Pa-1 ~> 1]
943  logical :: overrode_x, overrode_y
944 
945  isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
946  pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
947 
948  tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0
949  ! Either reads data or leaves contents unchanged
950  overrode_x = .false. ; overrode_y = .false.
951  call data_override('OCN', 'taux_adj', tempx_at_h(isc:iec,jsc:jec), time, override=overrode_x)
952  call data_override('OCN', 'tauy_adj', tempy_at_h(isc:iec,jsc:jec), time, override=overrode_y)
953 
954  if (overrode_x .or. overrode_y) then
955  if (.not. (overrode_x .and. overrode_y)) call mom_error(fatal,"apply_flux_adjustments: "//&
956  "Both taux_adj and tauy_adj must be specified, or neither, in data_table")
957 
958  ! Rotate winds
959  call pass_vector(tempx_at_h, tempy_at_h, g%Domain, to_all, agrid, halo=1)
960  do j=jsc-1,jec+1 ; do i=isc-1,iec+1
961  dlondx = g%geoLonCu(i,j) - g%geoLonCu(i-1,j)
962  dlondy = g%geoLonCv(i,j) - g%geoLonCv(i,j-1)
963  rdlon = sqrt( dlondx * dlondx + dlondy * dlondy )
964  if (rdlon > 0.) rdlon = 1. / rdlon
965  cosa = dlondx * rdlon
966  sina = dlondy * rdlon
967  zonal_tau = pa_conversion * tempx_at_h(i,j)
968  merid_tau = pa_conversion * tempy_at_h(i,j)
969  tempx_at_h(i,j) = cosa * zonal_tau - sina * merid_tau
970  tempy_at_h(i,j) = sina * zonal_tau + cosa * merid_tau
971  enddo ; enddo
972 
973  ! Average to C-grid locations
974  do j=jsc,jec ; do i=isc-1,iec
975  forces%taux(i,j) = forces%taux(i,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) )
976  enddo ; enddo
977 
978  do j=jsc-1,jec ; do i=isc,iec
979  forces%tauy(i,j) = forces%tauy(i,j) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) )
980  enddo ; enddo
981  endif ! overrode_x .or. overrode_y
982 
983 end subroutine apply_force_adjustments
984 
985 !> Save any restart files associated with the surface forcing.
986 subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
987  filename_suffix)
988  type(surface_forcing_cs), pointer :: cs !< A pointer to the control structure returned
989  !! by a previous call to surface_forcing_init
990  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
991  type(time_type), intent(in) :: time !< The current model time
992  character(len=*), intent(in) :: directory !< The directory into which to write the
993  !! restart files
994  logical, optional, intent(in) :: time_stamped !< If true, the restart file names include
995  !! a unique time stamp. The default is false.
996  character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-
997  !! stamp) to append to the restart file names.
998 
999  if (.not.associated(cs)) return
1000  if (.not.associated(cs%restart_CSp)) return
1001  call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1002 
1003 end subroutine forcing_save_restart
1004 
1005 !> Initialize the surface forcing, including setting parameters and allocating permanent memory.
1006 subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp)
1007  type(time_type), intent(in) :: time !< The current model time
1008  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1009  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1010  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1011  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
1012  !! diagnostic output
1013  type(surface_forcing_cs), pointer :: cs !< A pointer that is set to point to the control
1014  !! structure for this module
1015  logical, optional, intent(in) :: restore_salt !< If present and true surface salinity
1016  !! restoring will be applied in this model.
1017  logical, optional, intent(in) :: restore_temp !< If present and true surface temperature
1018  !! restoring will be applied in this model.
1019 
1020  ! Local variables
1021  real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1].
1022  type(directories) :: dirs
1023  logical :: new_sim, iceberg_flux_diags
1024  type(time_type) :: time_frc
1025  character(len=200) :: tideamp_file, gust_file, salt_file, temp_file ! Input file names.
1026 ! This include declares and sets the variable "version".
1027 #include "version_variable.h"
1028  character(len=40) :: mdl = "MOM_surface_forcing_mct" ! This module's name.
1029  character(len=48) :: stagger
1030  character(len=48) :: flnam
1031  character(len=240) :: basin_file
1032  integer :: i, j, isd, ied, jsd, jed
1033 
1034  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1035 
1036  if (associated(cs)) then
1037  call mom_error(warning, "surface_forcing_init called with an associated "// &
1038  "control structure.")
1039  return
1040  endif
1041  allocate(cs)
1042 
1043  id_clock_forcing=cpu_clock_id('Ocean surface forcing', grain=clock_subcomponent)
1044  call cpu_clock_begin(id_clock_forcing)
1045 
1046  cs%diag => diag
1047 
1048  call write_version_number(version)
1049  ! Read all relevant parameters and write them to the model log.
1050  call log_version(param_file, mdl, version, "")
1051 
1052  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, &
1053  "The directory in which all input files are found.", &
1054  default=".")
1055  cs%inputdir = slasher(cs%inputdir)
1056  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
1057  "If true, Temperature and salinity are used as state "//&
1058  "variables.", default=.true.)
1059  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
1060  "The mean ocean density used with BOUSSINESQ true to "//&
1061  "calculate accelerations and the mass for conservation "//&
1062  "properties, or with BOUSSINSEQ false to convert some "//&
1063  "parameters from vertical units of m to kg m-2.", &
1064  units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
1065  call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1066  "The latent heat of fusion.", units="J/kg", default=hlf)
1067  call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1068  "The latent heat of fusion.", units="J/kg", default=hlv)
1069  call get_param(param_file, mdl, "MAX_P_SURF", cs%max_p_surf, &
1070  "The maximum surface pressure that can be exerted by the "//&
1071  "atmosphere and floating sea-ice or ice shelves. This is "//&
1072  "needed because the FMS coupling structure does not "//&
1073  "limit the water that can be frozen out of the ocean and "//&
1074  "the ice-ocean heat fluxes are treated explicitly. No "//&
1075  "limit is applied if a negative value is used.", units="Pa", &
1076  default=-1.0)
1077  call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", &
1078  cs%adjust_net_srestore_to_zero, &
1079  "If true, adjusts the salinity restoring seen to zero "//&
1080  "whether restoring is via a salt flux or virtual precip.",&
1081  default=restore_salt)
1082  call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_BY_SCALING", &
1083  cs%adjust_net_srestore_by_scaling, &
1084  "If true, adjustments to salt restoring to achieve zero net are "//&
1085  "made by scaling values without moving the zero contour.",&
1086  default=.false.)
1087  call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_TO_ZERO", &
1088  cs%adjust_net_fresh_water_to_zero, &
1089  "If true, adjusts the net fresh-water forcing seen "//&
1090  "by the ocean (including restoring) to zero.", default=.false.)
1091  if (cs%adjust_net_fresh_water_to_zero) &
1092  call get_param(param_file, mdl, "USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
1093  cs%use_net_FW_adjustment_sign_bug, &
1094  "If true, use the wrong sign for the adjustment to "//&
1095  "the net fresh-water.", default=.false.)
1096  call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", &
1097  cs%adjust_net_fresh_water_by_scaling, &
1098  "If true, adjustments to net fresh water to achieve zero net are "//&
1099  "made by scaling values without moving the zero contour.",&
1100  default=.false.)
1101  call get_param(param_file, mdl, "ICE_SALT_CONCENTRATION", &
1102  cs%ice_salt_concentration, &
1103  "The assumed sea-ice salinity needed to reverse engineer the "//&
1104  "melt flux (or ice-ocean fresh-water flux).", &
1105  units="kg/kg", default=0.005)
1106  call get_param(param_file, mdl, "USE_LIMITED_PATM_SSH", cs%use_limited_P_SSH, &
1107  "If true, return the sea surface height with the "//&
1108  "correction for the atmospheric (and sea-ice) pressure "//&
1109  "limited by max_p_surf instead of the full atmospheric "//&
1110  "pressure.", default=.true.)
1111 
1112  call get_param(param_file, mdl, "WIND_STAGGER", stagger, &
1113  "A case-insensitive character string to indicate the "//&
1114  "staggering of the input wind stress field. Valid "//&
1115  "values are 'A', 'B', or 'C'.", default="C")
1116  if (uppercase(stagger(1:1)) == 'A') then ; cs%wind_stagger = agrid
1117  elseif (uppercase(stagger(1:1)) == 'B') then ; cs%wind_stagger = bgrid_ne
1118  elseif (uppercase(stagger(1:1)) == 'C') then ; cs%wind_stagger = cgrid_ne
1119  else ; call mom_error(fatal,"surface_forcing_init: WIND_STAGGER = "// &
1120  trim(stagger)//" is invalid.") ; endif
1121  call get_param(param_file, mdl, "WIND_STRESS_MULTIPLIER", cs%wind_stress_multiplier, &
1122  "A factor multiplying the wind-stress given to the ocean by the "//&
1123  "coupler. This is used for testing and should be =1.0 for any "//&
1124  "production runs.", default=1.0)
1125 
1126  if (restore_salt) then
1127  call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
1128  "The constant that relates the restoring surface fluxes "//&
1129  "to the relative surface anomalies (akin to a piston "//&
1130  "velocity). Note the non-MKS units.", units="m day-1", &
1131  fail_if_missing=.true.)
1132  call get_param(param_file, mdl, "SALT_RESTORE_FILE", cs%salt_restore_file, &
1133  "A file in which to find the surface salinity to use for restoring.", &
1134  default="salt_restore.nc")
1135  call get_param(param_file, mdl, "SALT_RESTORE_VARIABLE", cs%salt_restore_var_name, &
1136  "The name of the surface salinity variable to read from "//&
1137  "SALT_RESTORE_FILE for restoring salinity.", &
1138  default="salt")
1139 ! Convert CS%Flux_const from m day-1 to m s-1.
1140  cs%Flux_const = cs%Flux_const / 86400.0
1141 
1142  call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", cs%salt_restore_as_sflux, &
1143  "If true, the restoring of salinity is applied as a salt "//&
1144  "flux instead of as a freshwater flux.", default=.false.)
1145  call get_param(param_file, mdl, "MAX_DELTA_SRESTORE", cs%max_delta_srestore, &
1146  "The maximum salinity difference used in restoring terms.", &
1147  units="PSU or g kg-1", default=999.0)
1148  call get_param(param_file, mdl, "MASK_SRESTORE_UNDER_ICE", &
1149  cs%mask_srestore_under_ice, &
1150  "If true, disables SSS restoring under sea-ice based on a frazil "//&
1151  "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", &
1152  default=.false.)
1153  call get_param(param_file, mdl, "MASK_SRESTORE_MARGINAL_SEAS", &
1154  cs%mask_srestore_marginal_seas, &
1155  "If true, disable SSS restoring in marginal seas. Only used when "//&
1156  "RESTORE_SALINITY is True.", default=.false.)
1157  call get_param(param_file, mdl, "BASIN_FILE", basin_file, &
1158  "A file in which to find the basin masks, in variable 'basin'.", &
1159  default="basin.nc")
1160  basin_file = trim(cs%inputdir) // trim(basin_file)
1161  call safe_alloc_ptr(cs%basin_mask,isd,ied,jsd,jed) ; cs%basin_mask(:,:) = 1.0
1162  if (cs%mask_srestore_marginal_seas) then
1163  call mom_read_data(basin_file,'basin',cs%basin_mask,g%domain, timelevel=1)
1164  do j=jsd,jed ; do i=isd,ied
1165  if (cs%basin_mask(i,j) >= 6.0) then ; cs%basin_mask(i,j) = 0.0
1166  else ; cs%basin_mask(i,j) = 1.0 ; endif
1167  enddo ; enddo
1168  endif
1169  call get_param(param_file, mdl, "MASK_SRESTORE", cs%mask_srestore, &
1170  "If true, read a file (salt_restore_mask) containing "//&
1171  "a mask for SSS restoring.", default=.false.)
1172  endif
1173 
1174  if (restore_temp) then
1175  call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
1176  "The constant that relates the restoring surface fluxes "//&
1177  "to the relative surface anomalies (akin to a piston "//&
1178  "velocity). Note the non-MKS units.", units="m day-1", &
1179  fail_if_missing=.true.)
1180  call get_param(param_file, mdl, "SST_RESTORE_FILE", cs%temp_restore_file, &
1181  "A file in which to find the surface temperature to use for restoring.", &
1182  default="temp_restore.nc")
1183  call get_param(param_file, mdl, "SST_RESTORE_VARIABLE", cs%temp_restore_var_name, &
1184  "The name of the surface temperature variable to read from "//&
1185  "SST_RESTORE_FILE for restoring sst.", &
1186  default="temp")
1187  ! Convert CS%Flux_const from m day-1 to m s-1.
1188  cs%Flux_const = cs%Flux_const / 86400.0
1189 
1190  call get_param(param_file, mdl, "MAX_DELTA_TRESTORE", cs%max_delta_trestore, &
1191  "The maximum sst difference used in restoring terms.", &
1192  units="degC ", default=999.0)
1193 
1194  call get_param(param_file, mdl, "MASK_TRESTORE", cs%mask_trestore, &
1195  "If true, read a file (temp_restore_mask) containing "//&
1196  "a mask for SST restoring.", default=.false.)
1197  endif
1198 
1199 ! Optionally read tidal amplitude from input file (m s-1) on model grid.
1200 ! Otherwise use default tidal amplitude for bottom frictionally-generated
1201 ! dissipation. Default cd_tides is chosen to yield approx 1 TWatt of
1202 ! work done against tides globally using OSU tidal amplitude.
1203  call get_param(param_file, mdl, "CD_TIDES", cs%cd_tides, &
1204  "The drag coefficient that applies to the tides.", &
1205  units="nondim", default=1.0e-4)
1206  call get_param(param_file, mdl, "READ_TIDEAMP", cs%read_TIDEAMP, &
1207  "If true, read a file (given by TIDEAMP_FILE) containing "//&
1208  "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1209  if (cs%read_TIDEAMP) then
1210  call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
1211  "The path to the file containing the spatially varying "//&
1212  "tidal amplitudes with INT_TIDE_DISSIPATION.", &
1213  default="tideamp.nc")
1214  cs%utide=0.0
1215  else
1216  call get_param(param_file, mdl, "UTIDE", cs%utide, &
1217  "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1218  units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1219  endif
1220 
1221  call safe_alloc_ptr(cs%TKE_tidal,isd,ied,jsd,jed)
1222  call safe_alloc_ptr(cs%ustar_tidal,isd,ied,jsd,jed)
1223 
1224  if (cs%read_TIDEAMP) then
1225  tideamp_file = trim(cs%inputdir) // trim(tideamp_file)
1226  call mom_read_data(tideamp_file,'tideamp',cs%TKE_tidal,g%domain,timelevel=1, scale=us%m_to_Z*us%T_to_s)
1227  do j=jsd, jed; do i=isd, ied
1228  utide = cs%TKE_tidal(i,j)
1229  cs%TKE_tidal(i,j) = g%mask2dT(i,j)*cs%Rho0*cs%cd_tides*(utide*utide*utide)
1230  cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1231  enddo ; enddo
1232  else
1233  do j=jsd,jed; do i=isd,ied
1234  utide = cs%utide
1235  cs%TKE_tidal(i,j) = cs%Rho0*cs%cd_tides*(utide*utide*utide)
1236  cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1237  enddo ; enddo
1238  endif
1239 
1240  call time_interp_external_init
1241 
1242 ! Optionally read a x-y gustiness field in place of a global
1243 ! constant.
1244 
1245  call get_param(param_file, mdl, "READ_GUST_2D", cs%read_gust_2d, &
1246  "If true, use a 2-dimensional gustiness supplied from "//&
1247  "an input file", default=.false.)
1248  call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
1249  "The background gustiness in the winds.", units="Pa", &
1250  default=0.02)
1251  if (cs%read_gust_2d) then
1252  call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
1253  "The file in which the wind gustiness is found in "//&
1254  "variable gustiness.")
1255 
1256  call safe_alloc_ptr(cs%gust,isd,ied,jsd,jed)
1257  gust_file = trim(cs%inputdir) // trim(gust_file)
1258  call mom_read_data(gust_file,'gustiness',cs%gust,g%domain, timelevel=1, &
1259  scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z) ! units in file should be Pa
1260  endif
1261 
1262 ! See whether sufficiently thick sea ice should be treated as rigid.
1263  call get_param(param_file, mdl, "USE_RIGID_SEA_ICE", cs%rigid_sea_ice, &
1264  "If true, sea-ice is rigid enough to exert a "//&
1265  "nonhydrostatic pressure that resist vertical motion.", &
1266  default=.false.)
1267  if (cs%rigid_sea_ice) then
1268  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
1269  "The gravitational acceleration of the Earth.", &
1270  units="m s-2", default = 9.80)
1271  call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", cs%density_sea_ice, &
1272  "A typical density of sea ice, used with the kinematic "//&
1273  "viscosity, when USE_RIGID_SEA_ICE is true.", units="kg m-3", &
1274  default=900.0)
1275  call get_param(param_file, mdl, "SEA_ICE_VISCOSITY", cs%Kv_sea_ice, &
1276  "The kinematic viscosity of sufficiently thick sea ice "//&
1277  "for use in calculating the rigidity of sea ice.", &
1278  units="m2 s-1", default=1.0e9)
1279  call get_param(param_file, mdl, "SEA_ICE_RIGID_MASS", cs%rigid_sea_ice_mass, &
1280  "The mass of sea-ice per unit area at which the sea-ice "//&
1281  "starts to exhibit rigidity", units="kg m-2", default=1000.0)
1282  endif
1283 
1284  call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
1285  "If true, makes available diagnostics of fluxes from icebergs "//&
1286  "as seen by MOM6.", default=.false.)
1287  call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles, &
1288  use_berg_fluxes=iceberg_flux_diags)
1289 
1290  call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", cs%allow_flux_adjustments, &
1291  "If true, allows flux adjustments to specified via the "//&
1292  "data_table using the component name 'OCN'.", default=.false.)
1293  if (cs%allow_flux_adjustments) then
1294  call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
1295  endif
1296 
1297  if (present(restore_salt)) then ; if (restore_salt) then
1298  salt_file = trim(cs%inputdir) // trim(cs%salt_restore_file)
1299  cs%id_srestore = init_external_field(salt_file, cs%salt_restore_var_name, domain=g%Domain%mpp_domain)
1300  call safe_alloc_ptr(cs%srestore_mask,isd,ied,jsd,jed); cs%srestore_mask(:,:) = 1.0
1301  if (cs%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes
1302  flnam = trim(cs%inputdir) // 'salt_restore_mask.nc'
1303  call mom_read_data(flnam,'mask', cs%srestore_mask, g%domain, timelevel=1)
1304  endif
1305  endif ; endif
1306 
1307  if (present(restore_temp)) then ; if (restore_temp) then
1308  temp_file = trim(cs%inputdir) // trim(cs%temp_restore_file)
1309  cs%id_trestore = init_external_field(temp_file, cs%temp_restore_var_name, domain=g%Domain%mpp_domain)
1310  call safe_alloc_ptr(cs%trestore_mask,isd,ied,jsd,jed); cs%trestore_mask(:,:) = 1.0
1311  if (cs%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes
1312  flnam = trim(cs%inputdir) // 'temp_restore_mask.nc'
1313  call mom_read_data(flnam, 'mask', cs%trestore_mask, g%domain, timelevel=1)
1314  endif
1315  endif ; endif
1316 
1317  ! Set up any restart fields associated with the forcing.
1318  call restart_init(param_file, cs%restart_CSp, "MOM_forcing.res")
1319  call restart_init_end(cs%restart_CSp)
1320 
1321  if (associated(cs%restart_CSp)) then
1322  call get_mom_input(dirs=dirs)
1323 
1324  new_sim = .false.
1325  if ((dirs%input_filename(1:1) == 'n') .and. &
1326  (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1327  if (.not.new_sim) then
1328  call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1329  g, cs%restart_CSp)
1330  endif
1331  endif
1332 
1333  call user_revise_forcing_init(param_file, cs%urf_CS)
1334 
1335  call cpu_clock_end(id_clock_forcing)
1336 end subroutine surface_forcing_init
1337 
1338 !> Clean up and deallocate any memory associated with this module and its children.
1339 subroutine surface_forcing_end(CS, fluxes)
1340  type(surface_forcing_cs), pointer :: cs !< A pointer to the control structure returned by
1341  !! a previous call to surface_forcing_init, it will
1342  !! be deallocated here.
1343  type(forcing), optional, intent(inout) :: fluxes !< A structure containing pointers to all
1344  !! possible mass, heat or salt flux forcing fields.
1345  !! If present, it will be deallocated here.
1346 
1347  if (present(fluxes)) call deallocate_forcing_type(fluxes)
1348 
1349  if (associated(cs)) deallocate(cs)
1350  cs => null()
1351 
1352 end subroutine surface_forcing_end
1353 
1354 !> Write out a set of messages with checksums of the fields in an ice_ocen_boundary type
1355 subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
1357  character(len=*), intent(in) :: id !< An identifying string for this call
1358  integer, intent(in) :: timestep !< The number of elapsed timesteps
1359  type(ice_ocean_boundary_type), &
1360  intent(in) :: iobt !< An ice-ocean boundary type with fluxes to drive the
1361  !! ocean in a coupled model whose checksums are reported
1362 
1363  ! local variables
1364  integer :: n,m, outunit
1365 
1366  outunit = stdout()
1367 
1368  write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
1369  write(outunit,100) 'iobt%u_flux ' , mpp_chksum( iobt%u_flux )
1370  write(outunit,100) 'iobt%v_flux ' , mpp_chksum( iobt%v_flux )
1371  write(outunit,100) 'iobt%t_flux ' , mpp_chksum( iobt%t_flux )
1372  write(outunit,100) 'iobt%q_flux ' , mpp_chksum( iobt%q_flux )
1373  write(outunit,100) 'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux )
1374  write(outunit,100) 'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat)
1375  write(outunit,100) 'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt )
1376  write(outunit,100) 'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux )
1377  write(outunit,100) 'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir)
1378  write(outunit,100) 'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif)
1379  write(outunit,100) 'iobt%sw_flux_nir_dir' , mpp_chksum( iobt%sw_flux_nir_dir)
1380  write(outunit,100) 'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif)
1381  write(outunit,100) 'iobt%lprec ' , mpp_chksum( iobt%lprec )
1382  write(outunit,100) 'iobt%fprec ' , mpp_chksum( iobt%fprec )
1383  write(outunit,100) 'iobt%runoff ' , mpp_chksum( iobt%runoff )
1384  write(outunit,100) 'iobt%calving ' , mpp_chksum( iobt%calving )
1385  write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p )
1386  if (associated(iobt%ustar_berg)) &
1387  write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg )
1388  if (associated(iobt%area_berg)) &
1389  write(outunit,100) 'iobt%area_berg ' , mpp_chksum( iobt%area_berg )
1390  if (associated(iobt%mass_berg)) &
1391  write(outunit,100) 'iobt%mass_berg ' , mpp_chksum( iobt%mass_berg )
1392 100 FORMAT(" CHECKSUM::",a20," = ",z20)
1393 
1394  call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')
1395 
1396 end subroutine ice_ocn_bnd_type_chksum
1397 
1398 end module mom_surface_forcing_mct
user_revise_forcing::user_revise_forcing_cs
Control structure for user_revise_forcing.
Definition: user_revise_forcing.F90:23
mom_restart::restore_state
subroutine, public restore_state(filename, directory, day, G, CS)
restore_state reads the model state from previously generated files. All restart variables are read f...
Definition: MOM_restart.F90:984
mom_forcing_type::allocate_forcing_type
subroutine, public allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, iceberg, salt)
Conditionally allocate fields within the forcing type.
Definition: MOM_forcing_type.F90:2811
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:187
mom_variables::surface
Pointers to various fields which may be used describe the surface state of MOM, and which will be ret...
Definition: MOM_variables.F90:38
mom_spatial_means
Functions and routines to take area, volume, mass-weighted, layerwise, zonal or meridional means.
Definition: MOM_spatial_means.F90:2
mom_spatial_means::adjust_area_mean_to_zero
subroutine, public adjust_area_mean_to_zero(array, G, scaling, unit_scale)
Adjust 2d array such that area mean is zero without moving the zero contour.
Definition: MOM_spatial_means.F90:357
mom_surface_forcing_mct::surface_forcing_end
subroutine, private surface_forcing_end(CS, fluxes)
Clean up and deallocate any memory associated with this module and its children.
Definition: mom_surface_forcing_mct.F90:1340
mom_surface_forcing_mct::ice_ocn_bnd_type_chksum
subroutine, public ice_ocn_bnd_type_chksum(id, timestep, iobt)
Write out a set of messages with checksums of the fields in an ice_ocen_boundary type.
Definition: mom_surface_forcing_mct.F90:1356
mom_surface_forcing_mct::apply_force_adjustments
subroutine, private apply_force_adjustments(G, US, CS, Time, forces)
Adds mechanical forcing adjustments obtained via data_override Component name is 'OCN' Available adju...
Definition: mom_surface_forcing_mct.F90:930
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_constants
Provides a few physical constants.
Definition: MOM_constants.F90:2
user_revise_forcing
Provides a template for users to code updating the forcing fluxes.
Definition: user_revise_forcing.F90:2
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_surface_forcing_mct::convert_iob_to_forces
subroutine, public convert_iob_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
This subroutine translates the Ice_ocean_boundary_type into a MOM mechanical forcing type,...
Definition: mom_surface_forcing_mct.F90:591
mom_forcing_type::deallocate_forcing_type
subroutine, public deallocate_forcing_type(fluxes)
Deallocate the forcing type.
Definition: MOM_forcing_type.F90:2931
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_forcing_type::allocate_mech_forcing
subroutine, public allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg)
Conditionally allocate fields within the mechanical forcing type.
Definition: MOM_forcing_type.F90:2879
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_surface_forcing_mct::surface_forcing_cs
Contains pointers to the forcing fields which may be used to drive MOM. All fluxes are positive downw...
Definition: mom_surface_forcing_mct.F90:58
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:72
mom_forcing_type::mech_forcing_diags
subroutine, public mech_forcing_diags(forces, dt, G, time_end, diag, handles)
Offer mechanical forcing fields for diagnostics for those fields registered as part of register_forci...
Definition: MOM_forcing_type.F90:2201
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_surface_forcing_mct::surface_forcing_init
subroutine, public surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp)
Initialize the surface forcing, including setting parameters and allocating permanent memory.
Definition: mom_surface_forcing_mct.F90:1007
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:54
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_domains::to_all
integer, parameter, public to_all
A flag for passing in all directions.
Definition: MOM_domains.F90:132
mom_surface_forcing_mct::forcing_save_restart
subroutine, public forcing_save_restart(CS, G, Time, directory, time_stamped, filename_suffix)
Save any restart files associated with the surface forcing.
Definition: mom_surface_forcing_mct.F90:988
mom_surface_forcing_mct::convert_iob_to_fluxes
subroutine, public convert_iob_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, sfc_state, restore_salt, restore_temp)
This subroutine translates the Ice_ocean_boundary_type into a MOM thermodynamic forcing type,...
Definition: mom_surface_forcing_mct.F90:197
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_surface_forcing_mct::id_clock_forcing
integer id_clock_forcing
Definition: mom_surface_forcing_mct.F90:188
user_revise_forcing::user_alter_forcing
subroutine, public user_alter_forcing(state, fluxes, day, G, CS)
This subroutine sets the surface wind stresses.
Definition: user_revise_forcing.F90:34
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_get_input::get_mom_input
subroutine, public get_mom_input(param_file, dirs, check_params, default_input_filename, ensemble_num)
Get the names of the I/O directories and initialization file. Also calls the subroutine that opens ru...
Definition: MOM_get_input.F90:34
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
user_revise_forcing::user_revise_forcing_init
subroutine, public user_revise_forcing_init(param_file, CS)
Initialize the user_revise_forcing control structure.
Definition: user_revise_forcing.F90:49
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:107
mom_surface_forcing_mct::ice_ocean_boundary_type
Structure corresponding to forcing, but with the elements, units, and conventions that exactly confor...
Definition: mom_surface_forcing_mct.F90:147
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_restart::restart_init
subroutine, public restart_init(param_file, CS, restart_root)
Initialize this module and set up a restart control structure.
Definition: MOM_restart.F90:1421
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:50
mom_surface_forcing_mct::apply_flux_adjustments
subroutine, private apply_flux_adjustments(G, US, CS, Time, fluxes)
Adds thermodynamic flux adjustments obtained via data_override Component name is 'OCN' Available adju...
Definition: mom_surface_forcing_mct.F90:883
mom_restart::restart_init_end
subroutine, public restart_init_end(CS)
Indicate that all variables have now been registered.
Definition: MOM_restart.F90:1479
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_restart::save_restart
subroutine, public save_restart(directory, time, G, CS, time_stamped, filename, GV)
save_restart saves all registered variables to restart files.
Definition: MOM_restart.F90:781
mom_forcing_type::forcing_diags
Structure that defines the id handles for the forcing type.
Definition: MOM_forcing_type.F90:237
mom_surface_forcing_mct
Definition: mom_surface_forcing_mct.F90:1
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54
mom_forcing_type::deallocate_mech_forcing
subroutine, public deallocate_mech_forcing(forces)
Deallocate the mechanical forcing type.
Definition: MOM_forcing_type.F90:2983
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_forcing_type::register_forcing_type_diags
subroutine, public register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes)
Register members of the forcing type for diagnostics.
Definition: MOM_forcing_type.F90:1221
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_domains::fill_symmetric_edges
Do a set of halo updates that fill in the values at the duplicated edges of a staggered symmetric mem...
Definition: MOM_domains.F90:88
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239