ocean_model_MOM.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Top-level module for the MOM6 ocean model in coupled mode.
7
8! This is the top level module for the MOM6 ocean model. It contains routines
9! for initialization, termination and update of ocean model state. This
10! particular version wraps all of the calls for MOM6 in the calls that had
11! been used for MOM4.
12!
13! This code is a stop-gap wrapper of the MOM6 code to enable it to be called
14! in the same way as MOM4.
15
17use mom, only : extract_surface_state, allocate_surface_state, finish_mom_initialization
20use mom, only : save_mom_restart
21use mom_coms, only : field_chksum
22use mom_constants, only : celsius_kelvin_offset, hlf
23use mom_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type
24use mom_coupler_types, only : coupler_type_spawn, coupler_type_write_chksums
25use mom_coupler_types, only : coupler_type_initialized, coupler_type_copy_data
29use mom_domains, only : mom_domain_type, domain2d, clone_mom_domain, get_domain_extent
30use mom_domains, only : pass_var, pass_vector, agrid, bgrid_ne, cgrid_ne, to_all, omit_corners
31use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
33use mom_eos, only : gsw_sp_from_sr, gsw_pt_from_ct
34use mom_file_parser, only : get_param, log_version, close_param_file, param_file_type
35use mom_forcing_type, only : forcing, mech_forcing, allocate_forcing_type
39use mom_grid, only : ocean_grid_type
47use mom_time_manager, only : time_type, operator(>), operator(+), operator(-)
48use mom_time_manager, only : operator(*), operator(/), operator(/=)
49use mom_time_manager, only : operator(<=), operator(>=), operator(<)
50use mom_time_manager, only : real_to_time, time_to_real
54use mom_variables, only : surface
62use iso_fortran_env, only : int64
63
64#include <MOM_memory.h>
65
66#ifdef _USE_GENERIC_TRACER
68#endif
69
70implicit none ; private
71
77public ice_ocn_bnd_type_chksum
79public ocean_model_data_get
80public get_ocean_grid
82
83!> This interface extracts a named scalar field or array from the ocean surface or public type
84interface ocean_model_data_get
85 module procedure ocean_model_data1d_get
86 module procedure ocean_model_data2d_get
87end interface
88
89
90!> This type is used for communication with other components via the FMS coupler.
91!! The element names and types can be changed only with great deliberation, hence
92!! the persistence of things like the cutesy element name "avg_kount".
93type, public :: ocean_public_type
94 type(domain2d) :: domain !< The domain for the surface fields.
95 logical :: is_ocean_pe !< .true. on processors that run the ocean model.
96 character(len=32) :: instance_name = '' !< A name that can be used to identify
97 !! this instance of an ocean model, for example
98 !! in ensembles when writing messages.
99 integer, pointer, dimension(:) :: pelist => null() !< The list of ocean PEs.
100 logical, pointer, dimension(:,:) :: maskmap =>null() !< A pointer to an array
101 !! indicating which logical processors are actually used for
102 !! the ocean code. The other logical processors would be all
103 !! land points and are not assigned to actual processors.
104 !! This need not be assigned if all logical processors are used.
105
106 integer :: stagger = -999 !< The staggering relative to the tracer points
107 !! points of the two velocity components. Valid entries
108 !! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW,
109 !! corresponding to the community-standard Arakawa notation.
110 !! (These are named integers taken from the MOM_domains module.)
111 !! Following MOM5, stagger is BGRID_NE by default when the
112 !! ocean is initialized, but here it is set to -999 so that
113 !! a global max across ocean and non-ocean processors can be
114 !! used to determine its value.
115 real, pointer, dimension(:,:) :: &
116 t_surf => null(), & !< SST on t-cell [degrees Kelvin]
117 s_surf => null(), & !< SSS on t-cell [ppt]
118 u_surf => null(), & !< i-velocity at the locations indicated by stagger [m s-1].
119 v_surf => null(), & !< j-velocity at the locations indicated by stagger [m s-1].
120 sea_lev => null(), & !< Sea level in m after correction for surface pressure,
121 !! i.e. dzt(1) + eta_t + patm/rho0/grav [m]
122 frazil =>null(), & !< Accumulated heating [J m-2] from frazil
123 !! formation in the ocean.
124 melt_potential => null(), & !< Instantaneous heat used to melt sea ice [J m-2].
125 obld => null(), & !< Ocean boundary layer depth [m].
126 area => null(), & !< cell area of the ocean surface [m2].
127 calving => null(), &!< The mass per unit area of the ice shelf to convert to
128 !! bergs [kg m-2].
129 calving_hflx => null() !< Calving heat flux [W m-2].
130 type(coupler_2d_bc_type) :: fields !< A structure that may contain named
131 !! arrays of tracer-related surface fields.
132 integer :: avg_kount !< A count of contributions to running
133 !! sums, used externally by the FMS coupler
134 !! for accumulating averages of this type.
135 integer, dimension(2) :: axes = 0 !< Axis numbers that are available
136 !! for I/O using this surface data.
137end type ocean_public_type
138
139
140!> The ocean_state_type contains all information about the state of the ocean,
141!! with a format that is private so it can be readily changed without disrupting
142!! other coupled components.
143type, public :: ocean_state_type ; private
144 ! This type is private, and can therefore vary between different ocean models.
145 logical :: is_ocean_pe = .false. !< True if this is an ocean PE.
146 type(time_type) :: time !< The ocean model's time and master clock.
147 type(time_type) :: time_dyn !< The ocean model's time for the dynamics. Time and Time_dyn
148 !! should be the same after a full time step.
149 integer :: restart_control !< An integer that is bit-tested to determine whether
150 !! incremental restart files are saved and whether they
151 !! have a time stamped name. +1 (bit 0) for generic
152 !! files and +2 (bit 1) for time-stamped files. A
153 !! restart file is saved at the end of a run segment
154 !! unless Restart_control is negative.
155
156 integer :: nstep = 0 !< The number of calls to update_ocean that update the dynamics.
157 integer :: nstep_thermo = 0 !< The number of calls to update_ocean that update the thermodynamics.
158 logical :: use_ice_shelf !< If true, the ice shelf model is enabled.
159 logical :: use_waves !< If true use wave coupling.
160
161 logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
162 !! ocean dynamics and forcing fluxes.
163 real :: press_to_z !< A conversion factor between pressure and ocean depth,
164 !! usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1].
165 logical :: calve_ice_shelf_bergs = .false. !< If true, bergs are initialized according to
166 !! ice shelf flux through the ice front
167 real :: c_p !< The heat capacity of seawater [J degC-1 kg-1].
168 logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode
169 !! with the barotropic and baroclinic dynamics, thermodynamics,
170 !! etc. stepped forward integrated in time.
171 !! If true, all of the above are bypassed with all
172 !! fields necessary to integrate only the tracer advection
173 !! and diffusion equation read in from files stored from
174 !! a previous integration of the prognostic model.
175
176 logical :: single_step_call !< If true, advance the state of MOM with a single
177 !! step including both dynamics and thermodynamics.
178 !! If false, the two phases are advanced with
179 !! separate calls. The default is true.
180 ! The following 3 variables are only used here if single_step_call is false.
181 real :: dt !< (baroclinic) dynamics time step [T ~> s]
182 real :: dt_therm !< thermodynamics time step [T ~> s]
183 logical :: thermo_spans_coupling !< If true, thermodynamic and tracer time
184 !! steps can span multiple coupled time steps.
185 logical :: diabatic_first !< If true, apply diabatic and thermodynamic
186 !! processes before time stepping the dynamics.
187
188 type(directories) :: dirs !< A structure containing several relevant directory paths.
189 type(mech_forcing) :: forces !< A structure with the driving mechanical surface forces
190 type(forcing) :: fluxes !< A structure containing pointers to
191 !! the thermodynamic ocean forcing fields.
192 type(forcing) :: flux_tmp !< A secondary structure containing pointers to the
193 !! ocean forcing fields for when multiple coupled
194 !! timesteps are taken per thermodynamic step.
195 type(surface) :: sfc_state !< A structure containing pointers to
196 !! the ocean surface state fields.
197 type(ocean_grid_type), pointer :: &
198 grid => null() !< A pointer to a grid structure containing metrics
199 !! and related information.
200 type(verticalgrid_type), pointer :: &
201 gv => null() !< A pointer to a structure containing information
202 !! about the vertical grid.
203 type(unit_scale_type), pointer :: &
204 us => null() !< A pointer to a structure containing dimensional
205 !! unit scaling factors.
206 type(mom_control_struct) :: mom_csp
207 !< MOM control structure
208 type(ice_shelf_cs), pointer :: &
209 ice_shelf_csp => null() !< A pointer to the control structure for the
210 !! ice shelf model that couples with MOM6. This
211 !! is null if there is no ice shelf.
212 type(marine_ice_cs), pointer :: &
213 marine_ice_csp => null() !< A pointer to the control structure for the
214 !! marine ice effects module.
215 type(wave_parameters_cs), pointer :: &
216 waves => null() !< A pointer to the surface wave control structure
217 type(surface_forcing_cs), pointer :: &
218 forcing_csp => null() !< A pointer to the MOM forcing control structure
219 type(diag_ctrl), pointer :: &
220 diag => null() !< A pointer to the diagnostic regulatory structure
221end type ocean_state_type
222
223contains
224
225!> ocean_model_init initializes the ocean model, including registering fields
226!! for restarts and reading restart files if appropriate.
227!!
228!! This subroutine initializes both the ocean state and the ocean surface type.
229!! Because of the way that indices and domains are handled, Ocean_sfc must have
230!! been used in a previous call to initialize_ocean_type.
231subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas_fields_ocn, calve_ice_shelf_bergs)
232 type(ocean_public_type), target, &
233 intent(inout) :: ocean_sfc !< A structure containing various publicly
234 !! visible ocean surface properties after initialization,
235 !! the data in this type is intent out.
236 type(ocean_state_type), pointer :: os !< A structure whose internal
237 !! contents are private to ocean_model_mod that may be used to
238 !! contain all information about the ocean's interior state.
239 type(time_type), intent(in) :: time_init !< The start time for the coupled model's calendar
240 type(time_type), intent(in) :: time_in !< The time at which to initialize the ocean model.
241 integer, optional, intent(in) :: wind_stagger !< If present, the staggering of the winds that are
242 !! being provided in calls to update_ocean_model
243 type(coupler_1d_bc_type), &
244 optional, intent(in) :: gas_fields_ocn !< If present, this type describes the
245 !! ocean and surface-ice fields that will participate
246 !! in the calculation of additional gas or other
247 !! tracer fluxes, and can be used to spawn related
248 !! internal variables in the ice model.
249 logical, optional, intent(in) :: calve_ice_shelf_bergs !< If true, track ice shelf flux through a
250 !! static ice shelf, so that it can be converted into icebergs
251 ! Local variables
252 real :: rho0 ! The Boussinesq ocean density [R ~> kg m-3]
253 real :: g_earth ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
254 real :: hfrz !< If HFrz > 0 [Z ~> m], melt potential will be computed.
255 !! The actual depth over which melt potential is computed will
256 !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
257 !! If HFrz <= 0 (default), melt potential will not be computed.
258 logical :: use_melt_pot !< If true, allocate melt_potential array
259 logical :: point_calving ! Equals calve_ice_shelf_bergs if calve_ice_shelf_bergs is present
260
261 ! This include declares and sets the variable "version".
262# include "version_variable.h"
263 character(len=40) :: mdl = "ocean_model_init" ! This module's name.
264 character(len=48) :: stagger ! A string indicating the staggering locations for the
265 ! surface velocities returned to the coupler.
266 type(param_file_type) :: param_file !< A structure to parse for run-time parameters
267 logical :: use_temperature ! If true, temperature and salinity are state variables.
268
269 call calltree_enter("ocean_model_init(), ocean_model_MOM.F90")
270 if (associated(os)) then
271 call mom_error(warning, "ocean_model_init called with an associated "// &
272 "ocean_state_type structure. Model is already initialized.")
273 return
274 endif
275 allocate(os)
276
277! allocate(OS%fluxes)
278! allocate(OS%forces)
279! allocate(OS%flux_tmp)
280
281 os%is_ocean_pe = ocean_sfc%is_ocean_pe
282 if (.not.os%is_ocean_pe) return
283
284 os%Time = time_in ; os%Time_dyn = time_in
285 ! Call initialize MOM with an optional Ice Shelf CS which, if present triggers
286 ! initialization of ice shelf parameters and arrays.
287 point_calving = .false. ; if (present(calve_ice_shelf_bergs)) point_calving = calve_ice_shelf_bergs
288 call initialize_mom(os%Time, time_init, param_file, os%dirs, os%MOM_CSp, &
289 time_in, offline_tracer_mode=os%offline_tracer_mode, &
290 diag_ptr=os%diag, count_calls=.true., ice_shelf_csp=os%ice_shelf_CSp, &
291 waves_csp=os%Waves, calve_ice_shelf_bergs=point_calving)
292 call get_mom_state_elements(os%MOM_CSp, g=os%grid, gv=os%GV, us=os%US, c_p=os%C_p, &
293 c_p_scaled=os%fluxes%C_p, use_temp=use_temperature)
294
295 ! Read all relevant parameters and write them to the model log.
296 call log_version(param_file, mdl, version, "")
297
298 call get_param(param_file, mdl, "SINGLE_STEPPING_CALL", os%single_step_call, &
299 "If true, advance the state of MOM with a single step "//&
300 "including both dynamics and thermodynamics. If false, "//&
301 "the two phases are advanced with separate calls.", default=.true.)
302 call get_param(param_file, mdl, "DT", os%dt, &
303 "The (baroclinic) dynamics time step. The time-step that is actually "//&
304 "used will be an integer fraction of the forcing time-step.", &
305 units="s", scale=os%US%s_to_T, fail_if_missing=.true.)
306 call get_param(param_file, mdl, "DT_THERM", os%dt_therm, &
307 "The thermodynamic and tracer advection time step. "//&
308 "Ideally DT_THERM should be an integer multiple of DT "//&
309 "and less than the forcing or coupling time-step, unless "//&
310 "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
311 "can be an integer multiple of the coupling timestep. By "//&
312 "default DT_THERM is set to DT.", &
313 units="s", scale=os%US%s_to_T, default=os%US%T_to_s*os%dt)
314 call get_param(param_file, "MOM", "THERMO_SPANS_COUPLING", os%thermo_spans_coupling, &
315 "If true, the MOM will take thermodynamic and tracer "//&
316 "timesteps that can be longer than the coupling timestep. "//&
317 "The actual thermodynamic timestep that is used in this "//&
318 "case is the largest integer multiple of the coupling "//&
319 "timestep that is less than or equal to DT_THERM.", default=.false.)
320 call get_param(param_file, mdl, "DIABATIC_FIRST", os%diabatic_first, &
321 "If true, apply diabatic and thermodynamic processes, "//&
322 "including buoyancy forcing and mass gain or loss, "//&
323 "before stepping the dynamics forward.", default=.false.)
324
325 call get_param(param_file, mdl, "RESTART_CONTROL", os%Restart_control, &
326 "An integer whose bits encode which restart files are "//&
327 "written. Add 2 (bit 1) for a time-stamped file, and odd "//&
328 "(bit 0) for a non-time-stamped file. A restart file "//&
329 "will be saved at the end of the run segment for any "//&
330 "non-negative value.", default=1)
331 call get_param(param_file, mdl, "OCEAN_SURFACE_STAGGER", stagger, &
332 "A case-insensitive character string to indicate the "//&
333 "staggering of the surface velocity field that is "//&
334 "returned to the coupler. Valid values include "//&
335 "'A', 'B', or 'C'.", default="C")
336 if (uppercase(stagger(1:1)) == 'A') then ; ocean_sfc%stagger = agrid
337 elseif (uppercase(stagger(1:1)) == 'B') then ; ocean_sfc%stagger = bgrid_ne
338 elseif (uppercase(stagger(1:1)) == 'C') then ; ocean_sfc%stagger = cgrid_ne
339 else ; call mom_error(fatal,"ocean_model_init: OCEAN_SURFACE_STAGGER = "// &
340 trim(stagger)//" is invalid.") ; endif
341
342 call get_param(param_file, mdl, "RHO_0", rho0, &
343 "The mean ocean density used with BOUSSINESQ true to "//&
344 "calculate accelerations and the mass for conservation "//&
345 "properties, or with BOUSSINESQ false to convert some "//&
346 "parameters from vertical units of m to kg m-2.", &
347 units="kg m-3", default=1035.0, scale=os%US%kg_m3_to_R)
348 call get_param(param_file, mdl, "G_EARTH", g_earth, &
349 "The gravitational acceleration of the Earth.", &
350 units="m s-2", default=9.80, scale=os%US%m_s_to_L_T**2*os%US%Z_to_m)
351
352 call get_param(param_file, mdl, "ICE_SHELF", os%use_ice_shelf, &
353 "If true, enables the ice shelf model.", default=.false.)
354
355 call get_param(param_file, mdl, "ICEBERGS_APPLY_RIGID_BOUNDARY", os%icebergs_alter_ocean, &
356 "If true, allows icebergs to change boundary condition felt by ocean", default=.false.)
357
358 os%press_to_z = 1.0 / (rho0*g_earth)
359
360 ! Consider using a run-time flag to determine whether to do the diagnostic
361 ! vertical integrals, since the related 3-d sums are not negligible in cost.
362 call get_param(param_file, mdl, "HFREEZE", hfrz, &
363 "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
364 "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
365 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
366 "melt potential will not be computed.", &
367 units="m", default=-1.0, scale=os%US%m_to_Z, do_not_log=.true.)
368
369 if (hfrz > 0.0) then
370 use_melt_pot=.true.
371 else
372 use_melt_pot=.false.
373 endif
374
375 !allocate(OS%sfc_state)
376 call allocate_surface_state(os%sfc_state, os%grid, use_temperature, do_integrals=.true., &
377 gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot, &
378 use_iceshelves=os%use_ice_shelf)
379
380 if (present(wind_stagger)) then
381 call surface_forcing_init(time_in, os%grid, os%US, param_file, os%diag, &
382 os%forcing_CSp, wind_stagger)
383 else
384 call surface_forcing_init(time_in, os%grid, os%US, param_file, os%diag, &
385 os%forcing_CSp)
386 endif
387
388 if (os%use_ice_shelf) then
389 call initialize_ice_shelf_fluxes(os%ice_shelf_CSp, os%grid, os%US, os%fluxes)
390 call initialize_ice_shelf_forces(os%ice_shelf_CSp, os%grid, os%US, os%forces)
391 endif
392
393 if (os%icebergs_alter_ocean) then
394 call marine_ice_init(os%Time, os%grid, param_file, os%diag, os%marine_ice_CSp)
395 if (.not. os%use_ice_shelf) &
396 call allocate_forcing_type(os%grid, os%fluxes, shelf=.true.)
397 endif
398
399 call get_param(param_file, mdl, "USE_WAVES", os%Use_Waves, &
400 "If true, enables surface wave modules.", default=.false.)
401 ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because
402 ! it also initializes statistical waves.
403 call mom_wave_interface_init(os%Time, os%grid, os%GV, os%US, param_file, os%Waves, os%diag)
404
405 call initialize_ocean_public_type(os%grid%Domain, ocean_sfc, os%diag, &
406 gas_fields_ocn=gas_fields_ocn)
407
408 ! This call can only occur here if the coupler_bc_type variables have been
409 ! initialized already using the information from gas_fields_ocn.
410 if (present(gas_fields_ocn)) then
411 call coupler_type_set_diags(ocean_sfc%fields, "ocean_sfc", &
412 ocean_sfc%axes(1:2), time_in)
413
414 call extract_surface_state(os%MOM_CSp, os%sfc_state)
415
416 if (os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
417 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
418
419 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
420
421 endif
422
423 if (present(calve_ice_shelf_bergs)) then
424 if (calve_ice_shelf_bergs) then
425 call convert_shelf_state_to_ocean_type(ocean_sfc, os%Ice_shelf_CSp, os%US)
426 os%calve_ice_shelf_bergs=.true.
427 endif
428 endif
429
430 call close_param_file(param_file)
432
433 call mom_mesg('==== Completed MOM6 Coupled Initialization ====', 2)
434
435 call calltree_leave("ocean_model_init(")
436end subroutine ocean_model_init
437
438!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
439!! ocean model's state from the input value of Ocean_state (which must be for
440!! time time_start_update) for a time interval of Ocean_coupling_time_step,
441!! returning the publicly visible ocean surface properties in Ocean_sfc and
442!! storing the new ocean properties in Ocean_state.
443subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_update, &
444 Ocean_coupling_time_step, update_dyn, update_thermo, &
445 Ocn_fluxes_used, start_cycle, end_cycle, cycle_length)
447 intent(in) :: ice_ocean_boundary !< A structure containing the various
448 !! forcing fields coming from the ice and atmosphere.
449 type(ocean_state_type), &
450 pointer :: os !< A pointer to a private structure containing the
451 !! internal ocean state.
452 type(ocean_public_type), &
453 intent(inout) :: ocean_sfc !< A structure containing all the publicly visible
454 !! ocean surface fields after a coupling time step.
455 !! The data in this type is intent out.
456 type(time_type), intent(in) :: time_start_update !< The time at the beginning of the update step.
457 type(time_type), intent(in) :: ocean_coupling_time_step !< The amount of time over which to
458 !! advance the ocean.
459 logical, optional, intent(in) :: update_dyn !< If present and false, do not do updates
460 !! due to the ocean dynamics.
461 logical, optional, intent(in) :: update_thermo !< If present and false, do not do updates
462 !! due to the ocean thermodynamics or remapping.
463 logical, optional, intent(in) :: ocn_fluxes_used !< If present, this indicates whether the
464 !! cumulative thermodynamic fluxes from the ocean,
465 !! like frazil, have been used and should be reset.
466 logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
467 !! treated as the first call to step_MOM in a
468 !! time-stepping cycle; missing is like true.
469 logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
470 !! treated as the last call to step_MOM in a
471 !! time-stepping cycle; missing is like true.
472 real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle [s].
473
474 ! Local variables
475 type(time_type) :: time_seg_start ! Stores the dynamic or thermodynamic ocean model time at the
476 ! start of this call to allow step_MOM to temporarily change the time
477 ! as seen by internal modules.
478 type(time_type) :: time_thermo_start ! Stores the ocean model thermodynamics time at the start of
479 ! this call to allow step_MOM to temporarily change the time as seen by
480 ! internal modules.
481 type(time_type) :: time1 ! The value of the ocean model's time at the start of a call to step_MOM.
482 integer :: index_bnds(4) ! The computational domain index bounds in the ice-ocean boundary type.
483 real :: weight ! Flux accumulation weight of the current fluxes [nondim].
484 real :: dt_coupling ! The coupling time step [T ~> s].
485 real :: dt_therm ! A limited and quantized version of OS%dt_therm [T ~> s].
486 real :: dt_dyn ! The dynamics time step [T ~> s].
487 real :: dtdia ! The diabatic time step [T ~> s].
488 real :: t_elapsed_seg ! The elapsed time in this update segment [T ~> s].
489 integer :: n ! The internal iteration counter.
490 integer :: nts ! The number of baroclinic dynamics time steps in a thermodynamic step.
491 integer :: n_max ! The number of calls to step_MOM dynamics in this call to update_ocean_model.
492 integer :: n_last_thermo ! The iteration number the last time thermodynamics were updated.
493 logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans multiple dynamic timesteps.
494 logical :: do_dyn ! If true, step the ocean dynamics and transport.
495 logical :: do_thermo ! If true, step the ocean thermodynamics.
496 logical :: step_thermo ! If true, take a thermodynamic step.
497 integer :: is, ie, js, je
498
499 call calltree_enter("update_ocean_model(), ocean_model_MOM.F90")
500 dt_coupling = time_to_real(ocean_coupling_time_step, scale=os%US%s_to_T)
501
502 if (.not.associated(os)) then
503 call mom_error(fatal, "update_ocean_model called with an unassociated "// &
504 "ocean_state_type structure. ocean_model_init must be "// &
505 "called first to allocate this structure.")
506 return
507 endif
508
509 do_dyn = .true. ; if (present(update_dyn)) do_dyn = update_dyn
510 do_thermo = .true. ; if (present(update_thermo)) do_thermo = update_thermo
511
512 if (do_thermo .and. (time_start_update /= os%Time)) &
513 call mom_error(warning, "update_ocean_model: internal clock does not "//&
514 "agree with time_start_update argument.")
515 if (do_dyn .and. (time_start_update /= os%Time_dyn)) &
516 call mom_error(warning, "update_ocean_model: internal dynamics clock does not "//&
517 "agree with time_start_update argument.")
518
519 if (.not.(do_dyn .or. do_thermo)) call mom_error(fatal, &
520 "update_ocean_model called without updating either dynamics or thermodynamics.")
521 if (do_dyn .and. do_thermo .and. (os%Time /= os%Time_dyn)) call mom_error(fatal, &
522 "update_ocean_model called to update both dynamics and thermodynamics with inconsistent clocks.")
523
524 ! This is benign but not necessary if ocean_model_init_sfc was called or if
525 ! OS%sfc_state%tr_fields was spawned in ocean_model_init. Consider removing it.
526 is = os%grid%isc ; ie = os%grid%iec ; js = os%grid%jsc ; je = os%grid%jec
527 call coupler_type_spawn(ocean_sfc%fields, os%sfc_state%tr_fields, &
528 (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
529
530 ! Translate Ice_ocean_boundary into fluxes and forces.
531 call get_domain_extent(ocean_sfc%Domain, index_bnds(1), index_bnds(2), index_bnds(3), index_bnds(4))
532
533 if (do_dyn) then
534 call convert_iob_to_forces(ice_ocean_boundary, os%forces, index_bnds, os%Time_dyn, os%grid, os%US, &
535 os%forcing_CSp, dt_forcing=dt_coupling, reset_avg=os%fluxes%fluxes_used)
536 if (os%use_ice_shelf) &
537 call add_shelf_forces(os%grid, os%US, os%Ice_shelf_CSp, os%forces)
538 if (os%icebergs_alter_ocean) &
539 call iceberg_forces(os%grid, os%forces, os%use_ice_shelf, &
540 os%sfc_state, dt_coupling, os%marine_ice_CSp)
541 endif
542
543 if (do_thermo) then
544 if (os%fluxes%fluxes_used) then
545 call convert_iob_to_fluxes(ice_ocean_boundary, os%fluxes, index_bnds, os%Time, dt_coupling, &
546 os%grid, os%US, os%forcing_CSp, os%sfc_state)
547
548 ! Add ice shelf fluxes
549 if (os%use_ice_shelf) &
550 call shelf_calc_flux(os%sfc_state, os%fluxes, os%Time, dt_coupling, os%Ice_shelf_CSp)
551 if (os%icebergs_alter_ocean) &
552 call iceberg_fluxes(os%grid, os%US, os%fluxes, os%use_ice_shelf, &
553 os%sfc_state, dt_coupling, os%marine_ice_CSp)
554
555#ifdef _USE_GENERIC_TRACER
556 call enable_averages(dt_coupling, os%Time + ocean_coupling_time_step, os%diag) !Is this needed?
557 call mom_generic_tracer_fluxes_accumulate(os%fluxes, 1.0) ! Here weight=1, so just store the current fluxes
558 call disable_averaging(os%diag)
559#endif
560 else
561 ! The previous fluxes have not been used yet, so translate the input fluxes
562 ! into a temporary type and then accumulate them in about 20 lines.
563 os%flux_tmp%C_p = os%fluxes%C_p
564 call convert_iob_to_fluxes(ice_ocean_boundary, os%flux_tmp, index_bnds, os%Time, dt_coupling, &
565 os%grid, os%US, os%forcing_CSp, os%sfc_state)
566
567 if (os%use_ice_shelf) &
568 call shelf_calc_flux(os%sfc_state, os%flux_tmp, os%Time,dt_coupling, os%Ice_shelf_CSp)
569 if (os%icebergs_alter_ocean) &
570 call iceberg_fluxes(os%grid, os%US, os%flux_tmp, os%use_ice_shelf, &
571 os%sfc_state, dt_coupling, os%marine_ice_CSp)
572
573 call fluxes_accumulate(os%flux_tmp, os%fluxes, os%grid, weight)
574#ifdef _USE_GENERIC_TRACER
575 ! Incorporate the current tracer fluxes into the running averages
576 call mom_generic_tracer_fluxes_accumulate(os%flux_tmp, weight)
577#endif
578 endif
579 endif
580
581 ! The net mass forcing is not currently used in the MOM6 dynamics solvers, so this is may be unnecessary.
582 if (do_dyn .and. associated(os%forces%net_mass_src) .and. .not.os%forces%net_mass_src_set) &
583 call get_net_mass_forcing(os%fluxes, os%grid, os%US, os%forces%net_mass_src)
584
585 if (os%use_waves .and. do_thermo) then
586 ! For now, the waves are only updated on the thermodynamics steps, because that is where
587 ! the wave intensities are actually used to drive mixing. At some point, the wave updates
588 ! might also need to become a part of the ocean dynamics, according to B. Reichl.
589 call update_surface_waves(os%grid, os%GV, os%US, os%time, ocean_coupling_time_step, os%waves)
590 endif
591
592 if ((os%nstep==0) .and. (os%nstep_thermo==0)) then ! This is the first call to update_ocean_model.
593 call finish_mom_initialization(os%Time, os%dirs, os%MOM_CSp)
594 endif
595
596 time_thermo_start = os%Time
597 time_seg_start = os%Time ; if (do_dyn) time_seg_start = os%Time_dyn
598 time1 = time_seg_start
599
600 if (os%offline_tracer_mode .and. do_thermo) then
601 call step_offline(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp)
602 elseif ((.not.do_thermo) .or. (.not.do_dyn)) then
603 ! The call sequence is being orchestrated from outside of update_ocean_model.
604 if (present(cycle_length)) then
605 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, &
606 waves=os%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, &
607 start_cycle=start_cycle, end_cycle=end_cycle, cycle_length=os%US%s_to_T*cycle_length, &
608 reset_therm=ocn_fluxes_used)
609 else
610 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, &
611 waves=os%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, &
612 start_cycle=start_cycle, end_cycle=end_cycle, reset_therm=ocn_fluxes_used)
613 endif
614 elseif (os%single_step_call) then
615 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, waves=os%Waves)
616 else ! Step both the dynamics and thermodynamics with separate calls.
617 n_max = 1 ; if (dt_coupling > os%dt) n_max = ceiling(dt_coupling/os%dt - 0.001)
618 dt_dyn = dt_coupling / real(n_max)
619 thermo_does_span_coupling = (os%thermo_spans_coupling .and. (os%dt_therm > 1.5*dt_coupling))
620
621 if (thermo_does_span_coupling) then
622 dt_therm = dt_coupling * floor(os%dt_therm / dt_coupling + 0.001)
623 nts = floor(dt_therm/dt_dyn + 0.001)
624 else
625 nts = max(1,min(n_max,floor(os%dt_therm/dt_dyn + 0.001)))
626 n_last_thermo = 0
627 endif
628
629 time1 = time_seg_start ; t_elapsed_seg = 0.0
630 do n=1,n_max
631 if (os%diabatic_first) then
632 if (thermo_does_span_coupling) call mom_error(fatal, &
633 "MOM is not yet set up to have restarts that work with "//&
634 "THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
635 if (modulo(n-1,nts)==0) then
636 dtdia = dt_dyn*min(nts,n_max-(n-1))
637 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dtdia, os%MOM_CSp, &
638 waves=os%Waves, do_dynamics=.false., do_thermodynamics=.true., &
639 start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
640 endif
641
642 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_dyn, os%MOM_CSp, &
643 waves=os%Waves, do_dynamics=.true., do_thermodynamics=.false., &
644 start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
645 else
646 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_dyn, os%MOM_CSp, &
647 waves=os%Waves, do_dynamics=.true., do_thermodynamics=.false., &
648 start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
649
650 step_thermo = .false.
651 if (thermo_does_span_coupling) then
652 dtdia = dt_therm
653 step_thermo = mom_state_is_synchronized(os%MOM_CSp, adv_dyn=.true.)
654 elseif ((modulo(n,nts)==0) .or. (n==n_max)) then
655 dtdia = dt_dyn*(n - n_last_thermo)
656 n_last_thermo = n
657 step_thermo = .true.
658 endif
659
660 if (step_thermo) then
661 ! Back up Time1 to the start of the thermodynamic segment.
662 time1 = time1 - real_to_time(dtdia - dt_dyn, unscale=os%US%T_to_s)
663 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dtdia, os%MOM_CSp, &
664 waves=os%Waves, do_dynamics=.false., do_thermodynamics=.true., &
665 start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
666 endif
667 endif
668
669 t_elapsed_seg = t_elapsed_seg + dt_dyn
670 time1 = time_seg_start + real_to_time(t_elapsed_seg, unscale=os%US%T_to_s)
671 enddo
672 endif
673
674 if (do_dyn) os%Time_dyn = time_seg_start + ocean_coupling_time_step
675 if (do_dyn) os%nstep = os%nstep + 1
676 os%Time = time_thermo_start ! Reset the clock to compensate for shared pointers.
677 if (do_thermo) os%Time = os%Time + ocean_coupling_time_step
678 if (do_thermo) os%nstep_thermo = os%nstep_thermo + 1
679
680 if (do_dyn) then
681 call mech_forcing_diags(os%forces, dt_coupling, os%grid, os%Time_dyn, os%diag, os%forcing_CSp%handles)
682 endif
683
684 if (os%fluxes%fluxes_used .and. do_thermo) then
685 call forcing_diagnostics(os%fluxes, os%sfc_state, os%grid, os%US, os%Time, os%diag, os%forcing_CSp%handles)
686 endif
687
688 !only ,ale ice-shelf frazil adjustments if sfc_state%frazil was updated (do_thermo=True)
689 if (do_thermo .and. os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
690 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
691
692! Translate state into Ocean.
693! call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
694! OS%fluxes%p_surf_full, OS%press_to_z)
695 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
696 if (os%calve_ice_shelf_bergs) call convert_shelf_state_to_ocean_type(ocean_sfc,os%Ice_shelf_CSp, os%US)
697 time1 = os%Time ; if (do_dyn) time1 = os%Time_dyn
698 call coupler_type_send_data(ocean_sfc%fields, time1)
699
700 call calltree_leave("update_ocean_model()")
701end subroutine update_ocean_model
702
703!> This subroutine writes out the ocean model restart file.
704subroutine ocean_model_restart(OS, timestamp)
705 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
706 !! internal ocean state being saved to a restart file
707 character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be
708 !! prepended to the file name. (Currently this is unused.)
709
710 if (.not.mom_state_is_synchronized(os%MOM_CSp)) &
711 call mom_error(warning, "End of MOM_main reached with inconsistent "//&
712 "dynamics and advective times. Additional restart fields "//&
713 "that have not been coded yet would be required for reproducibility.")
714 if (.not.os%fluxes%fluxes_used) call mom_error(fatal, "ocean_model_restart "//&
715 "was called with unused buoyancy fluxes. For conservation, the ocean "//&
716 "restart files can only be created after the buoyancy forcing is applied.")
717
718 if (btest(os%Restart_control,1)) then
719 call save_mom_restart(os%MOM_CSp, os%dirs%restart_output_dir, os%Time, &
720 os%grid, time_stamped=.true., gv=os%GV)
721 call forcing_save_restart(os%forcing_CSp, os%grid, os%Time, &
722 os%dirs%restart_output_dir, .true.)
723 if (os%use_ice_shelf) then
724 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir, .true.)
725 endif
726 endif
727 if (btest(os%Restart_control,0)) then
728 call save_mom_restart(os%MOM_CSp, os%dirs%restart_output_dir, os%Time, &
729 os%grid, gv=os%GV)
730 call forcing_save_restart(os%forcing_CSp, os%grid, os%Time, &
731 os%dirs%restart_output_dir)
732 if (os%use_ice_shelf) then
733 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir)
734 endif
735 endif
736
737end subroutine ocean_model_restart
738! </SUBROUTINE> NAME="ocean_model_restart"
739
740!> ocean_model_end terminates the model run, saving the ocean state in a restart
741!! and deallocating any data associated with the ocean.
742subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time)
743 type(ocean_public_type), intent(inout) :: ocean_sfc !< An ocean_public_type structure that is
744 !! to be deallocated upon termination.
745 type(ocean_state_type), pointer :: ocean_state !< A pointer to the structure containing
746 !! the internal ocean state to be deallocated
747 !! upon termination.
748 type(time_type), intent(in) :: time !< The model time, used for writing restarts.
749
750 call ocean_model_save_restart(ocean_state, time)
751 call diag_mediator_end(time, ocean_state%diag)
752 call mom_end(ocean_state%MOM_CSp)
753 if (ocean_state%use_ice_shelf) call ice_shelf_end(ocean_state%Ice_shelf_CSp)
754end subroutine ocean_model_end
755
756!> ocean_model_save_restart causes restart files associated with the ocean to be
757!! written out.
758subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix)
759 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
760 !! internal ocean state (in).
761 type(time_type), intent(in) :: time !< The model time at this call, needed for writing files.
762 character(len=*), optional, intent(in) :: directory !< An optional directory into which to
763 !! write these restart files.
764 character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-stamp)
765 !! to append to the restart file names.
766! Note: This is a new routine - it will need to exist for the new incremental
767! checkpointing. It will also be called by ocean_model_end, giving the same
768! restart behavior as now in FMS.
769 character(len=200) :: restart_dir
770
771 if (.not.mom_state_is_synchronized(os%MOM_CSp)) &
772 call mom_error(warning, "ocean_model_save_restart called with inconsistent "//&
773 "dynamics and advective times. Additional restart fields "//&
774 "that have not been coded yet would be required for reproducibility.")
775 if (.not.os%fluxes%fluxes_used) call mom_error(fatal, "ocean_model_save_restart "//&
776 "was called with unused buoyancy fluxes. For conservation, the ocean "//&
777 "restart files can only be created after the buoyancy forcing is applied.")
778
779 if (present(directory)) then ; restart_dir = directory
780 else ; restart_dir = os%dirs%restart_output_dir ; endif
781
782 call save_mom_restart(os%MOM_CSp, restart_dir, time, os%grid, gv=os%GV)
783
784 call forcing_save_restart(os%forcing_CSp, os%grid, time, restart_dir)
785
786 if (os%use_ice_shelf) then
787 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir)
788 endif
789end subroutine ocean_model_save_restart
790
791!> Initialize the public ocean type
792subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_fields_ocn)
793 type(mom_domain_type), intent(in) :: input_domain !< The ocean model domain description
794 type(ocean_public_type), intent(inout) :: Ocean_sfc !< A structure containing various publicly
795 !! visible ocean surface properties after
796 !! initialization, whose elements are allocated here.
797 type(diag_ctrl), intent(in) :: diag !< A structure that regulates diagnostic output
798 type(coupler_1d_bc_type), &
799 optional, intent(in) :: gas_fields_ocn !< If present, this type describes the
800 !! ocean and surface-ice fields that will participate
801 !! in the calculation of additional gas or other
802 !! tracer fluxes.
803
804 integer :: xsz, ysz, layout(2)
805 ! ice-ocean-boundary fields are always allocated using absolute indices
806 ! and have no halos.
807 integer :: isc, iec, jsc, jec
808
809 call clone_mom_domain(input_domain, ocean_sfc%Domain, halo_size=0, symmetric=.false.)
810
811 call get_domain_extent(ocean_sfc%Domain, isc, iec, jsc, jec)
812
813 allocate ( ocean_sfc%t_surf (isc:iec,jsc:jec), &
814 ocean_sfc%s_surf (isc:iec,jsc:jec), &
815 ocean_sfc%u_surf (isc:iec,jsc:jec), &
816 ocean_sfc%v_surf (isc:iec,jsc:jec), &
817 ocean_sfc%sea_lev(isc:iec,jsc:jec), &
818 ocean_sfc%calving(isc:iec,jsc:jec), &
819 ocean_sfc%calving_hflx(isc:iec,jsc:jec), &
820 ocean_sfc%area (isc:iec,jsc:jec), &
821 ocean_sfc%melt_potential(isc:iec,jsc:jec), &
822 ocean_sfc%OBLD (isc:iec,jsc:jec), &
823 ocean_sfc%frazil (isc:iec,jsc:jec))
824
825 ocean_sfc%t_surf(:,:) = 0.0 ! time averaged sst (Kelvin) passed to atmosphere/ice model
826 ocean_sfc%s_surf(:,:) = 0.0 ! time averaged sss (psu) passed to atmosphere/ice models
827 ocean_sfc%u_surf(:,:) = 0.0 ! time averaged u-current (m/sec) passed to atmosphere/ice models
828 ocean_sfc%v_surf(:,:) = 0.0 ! time averaged v-current (m/sec) passed to atmosphere/ice models
829 ocean_sfc%sea_lev(:,:) = 0.0 ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav
830 ocean_sfc%calving(:,:) = 0.0 ! time accumulated ice sheet calving (kg m-2) passed to ice model
831 ocean_sfc%calving_hflx(:,:) = 0.0 ! time accumulated ice sheet calving heat flux (W m-2) passed to ice model
832 ocean_sfc%frazil(:,:) = 0.0 ! time accumulated frazil (J/m^2) passed to ice model
833 ocean_sfc%melt_potential(:,:) = 0.0 ! time accumulated melt potential (J/m^2) passed to ice model
834 ocean_sfc%OBLD(:,:) = 0.0 ! ocean boundary layer depth (m)
835 ocean_sfc%area(:,:) = 0.0
836 ocean_sfc%axes = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics
837
838 if (present(gas_fields_ocn)) then
839 call coupler_type_spawn(gas_fields_ocn, ocean_sfc%fields, (/isc,isc,iec,iec/), &
840 (/jsc,jsc,jec,jec/), suffix = '_ocn', as_needed=.true.)
841 endif
842
843end subroutine initialize_ocean_public_type
844
845!> This subroutine translates the coupler's ocean_data_type into MOM's
846!! surface state variable. This may eventually be folded into the MOM
847!! code that calculates the surface state in the first place.
848!! Note the offset in the arrays because the ocean_data_type has no
849!! halo points in its arrays and always uses absolute indices.
850subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_to_z)
851 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
852 !! describe the surface state of the ocean.
853 type(ocean_public_type), &
854 target, intent(inout) :: Ocean_sfc !< A structure containing various publicly
855 !! visible ocean surface fields, whose elements
856 !! have their data set here.
857 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
858 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
859 real, optional, intent(in) :: patm(:,:) !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
860 real, optional, intent(in) :: press_to_z !< A conversion factor between pressure and ocean
861 !! depth, usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1]
862 ! Local variables
863 character(len=48) :: val_str
864 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
865 integer :: i, j, i0, j0, is, ie, js, je
866
867 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
868 call pass_vector(sfc_state%u, sfc_state%v, g%Domain)
869
870 call get_domain_extent(ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
871 if (present(patm)) then
872 ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd).
873 if (.not.present(press_to_z)) call mom_error(fatal, &
874 'convert_state_to_ocean_type: press_to_z must be present if patm is.')
875 endif
876
877 i0 = is - isc_bnd ; j0 = js - jsc_bnd
878 if (sfc_state%T_is_conT) then
879 ! Convert the surface T from conservative T to potential T.
880 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
881 ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(us%S_to_ppt*sfc_state%SSS(i+i0,j+j0), &
882 us%C_to_degC*sfc_state%SST(i+i0,j+j0)) + celsius_kelvin_offset
883 enddo ; enddo
884 else
885 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
886 ocean_sfc%t_surf(i,j) = us%C_to_degC*sfc_state%SST(i+i0,j+j0) + celsius_kelvin_offset
887 enddo ; enddo
888 endif
889 if (sfc_state%S_is_absS) then
890 ! Convert the surface S from absolute salinity to practical salinity.
891 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
892 ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(us%S_to_ppt*sfc_state%SSS(i+i0,j+j0))
893 enddo ; enddo
894 else
895 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
896 ocean_sfc%s_surf(i,j) = us%S_to_ppt*sfc_state%SSS(i+i0,j+j0)
897 enddo ; enddo
898 endif
899
900 if (present(patm)) then
901 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
902 ocean_sfc%sea_lev(i,j) = us%Z_to_m * (sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z)
903 ocean_sfc%area(i,j) = us%L_to_m**2 * g%areaT(i+i0,j+j0)
904 enddo ; enddo
905 else
906 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
907 ocean_sfc%sea_lev(i,j) = us%Z_to_m * sfc_state%sea_lev(i+i0,j+j0)
908 ocean_sfc%area(i,j) = us%L_to_m**2 * g%areaT(i+i0,j+j0)
909 enddo ; enddo
910 endif
911
912 if (allocated(sfc_state%frazil)) then
913 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
914 ocean_sfc%frazil(i,j) = us%Q_to_J_kg*us%RZ_to_kg_m2 * sfc_state%frazil(i+i0,j+j0)
915 enddo ; enddo
916 endif
917
918 if (allocated(sfc_state%melt_potential)) then
919 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
920 ocean_sfc%melt_potential(i,j) = us%Q_to_J_kg*us%RZ_to_kg_m2 * sfc_state%melt_potential(i+i0,j+j0)
921 enddo ; enddo
922 endif
923
924 if (allocated(sfc_state%Hml)) then
925 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
926 ocean_sfc%OBLD(i,j) = us%Z_to_m * sfc_state%Hml(i+i0,j+j0)
927 enddo ; enddo
928 endif
929
930 if (ocean_sfc%stagger == agrid) then
931 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
932 ocean_sfc%u_surf(i,j) = g%mask2dT(i+i0,j+j0) * us%L_T_to_m_s * &
933 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i-1+i0,j+j0))
934 ocean_sfc%v_surf(i,j) = g%mask2dT(i+i0,j+j0) * us%L_T_to_m_s * &
935 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0,j-1+j0))
936 enddo ; enddo
937 elseif (ocean_sfc%stagger == bgrid_ne) then
938 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
939 ocean_sfc%u_surf(i,j) = g%mask2dBu(i+i0,j+j0) * us%L_T_to_m_s * &
940 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i+i0,j+j0+1))
941 ocean_sfc%v_surf(i,j) = g%mask2dBu(i+i0,j+j0) * us%L_T_to_m_s * &
942 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0+1,j+j0))
943 enddo ; enddo
944 elseif (ocean_sfc%stagger == cgrid_ne) then
945 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
946 ocean_sfc%u_surf(i,j) = g%mask2dCu(i+i0,j+j0)*us%L_T_to_m_s * sfc_state%u(i+i0,j+j0)
947 ocean_sfc%v_surf(i,j) = g%mask2dCv(i+i0,j+j0)*us%L_T_to_m_s * sfc_state%v(i+i0,j+j0)
948 enddo ; enddo
949 else
950 write(val_str, '(I8)') ocean_sfc%stagger
951 call mom_error(fatal, "convert_state_to_ocean_type: "//&
952 "Ocean_sfc%stagger has the unrecognized value of "//trim(val_str))
953 endif
954
955 if (coupler_type_initialized(sfc_state%tr_fields)) then
956 if (.not.coupler_type_initialized(ocean_sfc%fields)) then
957 call mom_error(fatal, "convert_state_to_ocean_type: "//&
958 "Ocean_sfc%fields has not been initialized.")
959 endif
960 call coupler_type_copy_data(sfc_state%tr_fields, ocean_sfc%fields)
961 endif
962
963end subroutine convert_state_to_ocean_type
964
965!> Converts the ice-shelf-to-ocean calving and calving_hflx variables from the ice-shelf state (ISS) type
966!! to the ocean public type
967subroutine convert_shelf_state_to_ocean_type(Ocean_sfc, CS, US)
968 type(ocean_public_type), &
969 target, intent(inout) :: Ocean_sfc !< A structure containing various publicly
970 !! visible ocean surface fields, whose elements
971 !! have their data set here.
972 type(ice_shelf_cs), pointer :: CS !< A pointer to the ice shelf control structure
973 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
974 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd, i, j
975
976 call get_domain_extent(ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
977
978 call ice_sheet_calving_to_ocean_sfc(cs,us,ocean_sfc%calving(isc_bnd:iec_bnd,jsc_bnd:jec_bnd),&
979 ocean_sfc%calving_hflx(isc_bnd:iec_bnd,jsc_bnd:jec_bnd))
980
982
983!> This subroutine extracts the surface properties from the ocean's internal
984!! state and stores them in the ocean type returned to the calling ice model.
985!! It has to be separate from the ocean_initialization call because the coupler
986!! module allocates the space for some of these variables.
987subroutine ocean_model_init_sfc(OS, Ocean_sfc)
988 type(ocean_state_type), pointer :: os !< The structure with the complete ocean state
989 type(ocean_public_type), intent(inout) :: ocean_sfc !< A structure containing various publicly
990 !! visible ocean surface properties after initialization, whose
991 !! elements have their data set here.
992 integer :: is, ie, js, je
993
994 is = os%grid%isc ; ie = os%grid%iec ; js = os%grid%jsc ; je = os%grid%jec
995 call coupler_type_spawn(ocean_sfc%fields, os%sfc_state%tr_fields, &
996 (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
997
998 call extract_surface_state(os%MOM_CSp, os%sfc_state)
999
1000 if (os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
1001 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
1002
1003 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
1004
1005end subroutine ocean_model_init_sfc
1006
1007!> ocean_model_flux_init is used to initialize properties of the air-sea fluxes
1008!! as determined by various run-time parameters. It can be called from
1009!! non-ocean PEs, or PEs that have not yet been initialized, and it can safely
1010!! be called multiple times.
1011subroutine ocean_model_flux_init(OS, verbosity)
1012 type(ocean_state_type), optional, pointer :: os !< An optional pointer to the ocean state,
1013 !! used to figure out if this is an ocean PE that
1014 !! has already been initialized.
1015 integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity.
1016
1017 logical :: os_is_set
1018 integer :: verbose
1019
1020 os_is_set = .false. ; if (present(os)) os_is_set = associated(os)
1021
1022 ! Use this to control the verbosity of output; consider rethinking this logic later.
1023 verbose = 5 ; if (os_is_set) verbose = 3
1024 if (present(verbosity)) verbose = verbosity
1025
1026 call call_tracer_flux_init(verbosity=verbose)
1027
1028end subroutine ocean_model_flux_init
1029
1030!> Ocean_stock_pe - returns the integrated stocks of heat, water, etc. for conservation checks.
1031!! Because of the way FMS is coded, only the root PE has the integrated amount,
1032!! while all other PEs get 0.
1033subroutine ocean_stock_pe(OS, index, value, time_index)
1034 use stock_constants_mod, only : istock_water, istock_heat,istock_salt
1035 type(ocean_state_type), pointer :: os !< A structure containing the internal ocean state.
1036 !! The data in OS is intent in.
1037 integer, intent(in) :: index !< The stock index for the quantity of interest.
1038 real, intent(out) :: value !< Sum returned for the conservation quantity of interest [various]
1039 integer, optional, intent(in) :: time_index !< An unused optional argument, present only for
1040 !! interfacial compatibility with other models.
1041! Arguments: OS - A structure containing the internal ocean state.
1042! (in) index - Index of conservation quantity of interest.
1043! (in) value - Sum returned for the conservation quantity of interest.
1044! (in,opt) time_index - Index for time level to use if this is necessary.
1045
1046 real :: salt ! The total salt in the ocean [kg]
1047
1048 value = 0.0
1049 if (.not.associated(os)) return
1050 if (.not.os%is_ocean_pe) return
1051
1052 select case (index)
1053 case (istock_water) ! Return the mass of fresh water in the ocean in [kg].
1054 if (os%GV%Boussinesq) then
1055 call get_ocean_stocks(os%MOM_CSp, mass=value, on_pe_only=.true.)
1056 else ! In non-Boussinesq mode, the mass of salt needs to be subtracted.
1057 call get_ocean_stocks(os%MOM_CSp, mass=value, salt=salt, on_pe_only=.true.)
1058 value = value - salt
1059 endif
1060 case (istock_heat) ! Return the heat content of the ocean in [J].
1061 call get_ocean_stocks(os%MOM_CSp, heat=value, on_pe_only=.true.)
1062 case (istock_salt) ! Return the mass of the salt in the ocean in [kg].
1063 call get_ocean_stocks(os%MOM_CSp, salt=value, on_pe_only=.true.)
1064 case default ; value = 0.0
1065 end select
1066 ! If the FMS coupler is changed so that Ocean_stock_PE is only called on
1067 ! ocean PEs, uncomment the following and eliminate the on_PE_only flags above.
1068 ! if (.not.is_root_pe()) value = 0.0
1069
1070end subroutine ocean_stock_pe
1071
1072!> This subroutine extracts a named 2-D field from the ocean surface or public type
1073subroutine ocean_model_data2d_get(OS, Ocean, name, array2D, isc, jsc)
1074 use mom_constants, only : celsius_kelvin_offset
1075 type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the
1076 !! internal ocean state (intent in).
1077 type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly
1078 !! visible ocean surface fields.
1079 character(len=*) , intent(in) :: name !< The name of the field to extract
1080 integer , intent(in) :: isc !< The starting i-index of array2D
1081 integer , intent(in) :: jsc !< The starting j-index of array2D
1082 real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must
1083 !! cover only the computational domain [various]
1084
1085 integer :: g_isc, g_iec, g_jsc, g_jec, g_isd, g_ied, g_jsd, g_jed, i, j
1086
1087 if (.not.associated(os)) return
1088 if (.not.os%is_ocean_pe) return
1089
1090 ! The problem is that %areaT is on MOM domain but Ice_Ocean_Boundary%... is on a haloless domain.
1091 ! We want to return the MOM data on the haloless (compute) domain
1092 call get_domain_extent(os%grid%Domain, g_isc, g_iec, g_jsc, g_jec, g_isd, g_ied, g_jsd, g_jed)
1093
1094 g_isc = g_isc-g_isd+1 ; g_iec = g_iec-g_isd+1 ; g_jsc = g_jsc-g_jsd+1 ; g_jec = g_jec-g_jsd+1
1095
1096 select case(name)
1097 case('area')
1098 array2d(isc:,jsc:) = os%US%L_to_m**2*os%grid%areaT(g_isc:g_iec,g_jsc:g_jec)
1099 case('mask')
1100 array2d(isc:,jsc:) = os%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec)
1101!OR same result
1102! do j=g_jsc,g_jec ; do i=g_isc,g_iec
1103! array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j)
1104! enddo ; enddo
1105 case('t_surf')
1106 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1107 case('t_pme')
1108 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1109 case('t_runoff')
1110 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1111 case('t_calving')
1112 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1113 case('btfHeat')
1114 array2d(isc:,jsc:) = 0
1115 case('cos_rot')
1116 array2d(isc:,jsc:) = os%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1
1117 case('sin_rot')
1118 array2d(isc:,jsc:) = os%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0
1119 case('s_surf')
1120 array2d(isc:,jsc:) = ocean%s_surf(isc:,jsc:)
1121 case('sea_lev')
1122 array2d(isc:,jsc:) = ocean%sea_lev(isc:,jsc:)
1123 case('frazil')
1124 array2d(isc:,jsc:) = ocean%frazil(isc:,jsc:)
1125 case('melt_pot')
1126 array2d(isc:,jsc:) = ocean%melt_potential(isc:,jsc:)
1127 case('obld')
1128 array2d(isc:,jsc:) = ocean%OBLD(isc:,jsc:)
1129 case default
1130 call mom_error(fatal,'get_ocean_grid_data2D: unknown argument name='//name)
1131 end select
1132
1133end subroutine ocean_model_data2d_get
1134
1135!> This subroutine extracts a named scalar field from the ocean surface or public type
1136subroutine ocean_model_data1d_get(OS, Ocean, name, value)
1137 type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the
1138 !! internal ocean state (intent in).
1139 type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly
1140 !! visible ocean surface fields.
1141 character(len=*), intent(in) :: name !< The name of the field to extract
1142 real, intent(out):: value !< The value of the named field [various]
1143
1144 if (.not.associated(os)) return
1145 if (.not.os%is_ocean_pe) return
1146
1147 select case(name)
1148 case('c_p')
1149 value = os%C_p
1150 case default
1151 call mom_error(fatal,'get_ocean_grid_data1D: unknown argument name='//name)
1152 end select
1153
1154end subroutine ocean_model_data1d_get
1155
1156!> Write out checksums for fields from the ocean surface state
1157subroutine ocean_public_type_chksum(id, timestep, ocn)
1158
1159 character(len=*), intent(in) :: id !< An identifying string for this call
1160 integer, intent(in) :: timestep !< The number of elapsed timesteps
1161 type(ocean_public_type), intent(in) :: ocn !< A structure containing various publicly
1162 !! visible ocean surface fields.
1163 ! Local variables
1164 integer(kind=int64) :: chks ! A checksum for the field
1165 logical :: root ! True only on the root PE
1166 integer :: outunit ! The output unit to write to
1167
1168 root = is_root_pe()
1169 outunit = stdout_if_root()
1170
1171 if (root) write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep
1172 chks = field_chksum(ocn%t_surf ) ; if (root) write(outunit,100) 'ocean%t_surf ', chks
1173 chks = field_chksum(ocn%s_surf ) ; if (root) write(outunit,100) 'ocean%s_surf ', chks
1174 chks = field_chksum(ocn%u_surf ) ; if (root) write(outunit,100) 'ocean%u_surf ', chks
1175 chks = field_chksum(ocn%v_surf ) ; if (root) write(outunit,100) 'ocean%v_surf ', chks
1176 chks = field_chksum(ocn%sea_lev) ; if (root) write(outunit,100) 'ocean%sea_lev ', chks
1177 chks = field_chksum(ocn%frazil ) ; if (root) write(outunit,100) 'ocean%frazil ', chks
1178 chks = field_chksum(ocn%melt_potential) ; if (root) write(outunit,100) 'ocean%melt_potential ', chks
1179 call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%')
1180100 FORMAT(" CHECKSUM::",a20," = ",z20)
1181
1182end subroutine ocean_public_type_chksum
1183
1184!> This subroutine gives a handle to the grid from ocean state
1185subroutine get_ocean_grid(OS, Gridp)
1186 ! Obtain the ocean grid.
1187 type(ocean_state_type) :: os !< A structure containing the
1188 !! internal ocean state
1189 type(ocean_grid_type) , pointer :: gridp !< The ocean's grid structure
1190
1191 gridp => os%grid
1192 return
1193end subroutine get_ocean_grid
1194
1195!> This subroutine extracts a named (u- or v-) 2-D surface current from ocean internal state
1196subroutine ocean_model_get_uv_surf(OS, Ocean, name, array2D, isc, jsc)
1197
1198 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
1199 !! internal ocean state (intent in).
1200 type(ocean_public_type), intent(in) :: ocean !< A structure containing various publicly
1201 !! visible ocean surface fields.
1202 character(len=*) , intent(in) :: name !< The name of the current (ua or va) to extract
1203 integer , intent(in) :: isc !< The starting i-index of array2D
1204 integer , intent(in) :: jsc !< The starting j-index of array2D
1205 real, dimension(isc:,jsc:), intent(out):: array2d !< The values of the named field, it must
1206 !! cover only the computational domain [L T-1 ~> m s-1]
1207
1208 type(ocean_grid_type) , pointer :: g !< The ocean's grid structure
1209 type(surface), pointer :: sfc_state !< A structure containing fields that
1210 !! describe the surface state of the ocean.
1211
1212 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
1213 integer :: i, j, i0, j0
1214 integer :: is, ie, js, je
1215
1216 if (.not.associated(os)) return
1217 if (.not.os%is_ocean_pe) return
1218
1219 g => os%grid
1220 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1221
1222 call get_domain_extent(ocean%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
1223
1224 i0 = is - isc_bnd ; j0 = js - jsc_bnd
1225
1226 sfc_state => os%sfc_state
1227
1228 call pass_vector(sfc_state%u, sfc_state%v, g%Domain)
1229
1230 select case(name)
1231 case('ua')
1232 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1233 array2d(i,j) = g%mask2dT(i+i0,j+j0) * &
1234 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i-1+i0,j+j0))
1235 enddo ; enddo
1236 case('va')
1237 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1238 array2d(i,j) = g%mask2dT(i+i0,j+j0) * &
1239 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0,j-1+j0))
1240 enddo ; enddo
1241 case('ub')
1242 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1243 array2d(i,j) = g%mask2dBu(i+i0,j+j0) * &
1244 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i+i0,j+j0+1))
1245 enddo ; enddo
1246 case('vb')
1247 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1248 array2d(i,j) = g%mask2dBu(i+i0,j+j0) * &
1249 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0+1,j+j0))
1250 enddo ; enddo
1251 case('uc')
1252 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1253 array2d(i,j) = g%mask2dCu(i+i0,j+j0) * sfc_state%u(i+i0,j+j0)
1254 enddo ; enddo
1255 case('vc')
1256 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1257 array2d(i,j) = g%mask2dCv(i+i0,j+j0) * sfc_state%v(i+i0,j+j0)
1258 enddo ; enddo
1259 case default
1260 call mom_error(fatal,'ocean_model_get_UV_surf: unknown argument name='//name)
1261 end select
1262
1263end subroutine ocean_model_get_uv_surf
1264
1265end module ocean_model_mod