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_fluxes(os%ice_shelf_CSp, os%grid, os%US, os%flux_tmp)
391 call initialize_ice_shelf_forces(os%ice_shelf_CSp, os%grid, os%US, os%forces)
392 endif
393
394 if (os%icebergs_alter_ocean) then
395 call marine_ice_init(os%Time, os%grid, param_file, os%diag, os%marine_ice_CSp)
396 if (.not. os%use_ice_shelf) &
397 call allocate_forcing_type(os%grid, os%fluxes, shelf=.true.)
398 endif
399
400 call get_param(param_file, mdl, "USE_WAVES", os%Use_Waves, &
401 "If true, enables surface wave modules.", default=.false.)
402 ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because
403 ! it also initializes statistical waves.
404 call mom_wave_interface_init(os%Time, os%grid, os%GV, os%US, param_file, os%Waves, os%diag)
405
406 call initialize_ocean_public_type(os%grid%Domain, ocean_sfc, os%diag, &
407 gas_fields_ocn=gas_fields_ocn)
408
409 ! This call can only occur here if the coupler_bc_type variables have been
410 ! initialized already using the information from gas_fields_ocn.
411 if (present(gas_fields_ocn)) then
412 call coupler_type_set_diags(ocean_sfc%fields, "ocean_sfc", &
413 ocean_sfc%axes(1:2), time_in)
414
415 call extract_surface_state(os%MOM_CSp, os%sfc_state)
416
417 if (os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
418 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
419
420 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
421
422 endif
423
424 if (present(calve_ice_shelf_bergs)) then
425 if (calve_ice_shelf_bergs) then
426 call convert_shelf_state_to_ocean_type(ocean_sfc, os%Ice_shelf_CSp, os%US)
427 os%calve_ice_shelf_bergs=.true.
428 endif
429 endif
430
431 call close_param_file(param_file)
433
434 call mom_mesg('==== Completed MOM6 Coupled Initialization ====', 2)
435
436 call calltree_leave("ocean_model_init(")
437end subroutine ocean_model_init
438
439!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
440!! ocean model's state from the input value of Ocean_state (which must be for
441!! time time_start_update) for a time interval of Ocean_coupling_time_step,
442!! returning the publicly visible ocean surface properties in Ocean_sfc and
443!! storing the new ocean properties in Ocean_state.
444subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_update, &
445 Ocean_coupling_time_step, update_dyn, update_thermo, &
446 Ocn_fluxes_used, start_cycle, end_cycle, cycle_length)
448 intent(in) :: ice_ocean_boundary !< A structure containing the various
449 !! forcing fields coming from the ice and atmosphere.
450 type(ocean_state_type), &
451 pointer :: os !< A pointer to a private structure containing the
452 !! internal ocean state.
453 type(ocean_public_type), &
454 intent(inout) :: ocean_sfc !< A structure containing all the publicly visible
455 !! ocean surface fields after a coupling time step.
456 !! The data in this type is intent out.
457 type(time_type), intent(in) :: time_start_update !< The time at the beginning of the update step.
458 type(time_type), intent(in) :: ocean_coupling_time_step !< The amount of time over which to
459 !! advance the ocean.
460 logical, optional, intent(in) :: update_dyn !< If present and false, do not do updates
461 !! due to the ocean dynamics.
462 logical, optional, intent(in) :: update_thermo !< If present and false, do not do updates
463 !! due to the ocean thermodynamics or remapping.
464 logical, optional, intent(in) :: ocn_fluxes_used !< If present, this indicates whether the
465 !! cumulative thermodynamic fluxes from the ocean,
466 !! like frazil, have been used and should be reset.
467 logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
468 !! treated as the first call to step_MOM in a
469 !! time-stepping cycle; missing is like true.
470 logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
471 !! treated as the last call to step_MOM in a
472 !! time-stepping cycle; missing is like true.
473 real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle [s].
474
475 ! Local variables
476 type(time_type) :: time_seg_start ! Stores the dynamic or thermodynamic ocean model time at the
477 ! start of this call to allow step_MOM to temporarily change the time
478 ! as seen by internal modules.
479 type(time_type) :: time_thermo_start ! Stores the ocean model thermodynamics time at the start of
480 ! this call to allow step_MOM to temporarily change the time as seen by
481 ! internal modules.
482 type(time_type) :: time1 ! The value of the ocean model's time at the start of a call to step_MOM.
483 integer :: index_bnds(4) ! The computational domain index bounds in the ice-ocean boundary type.
484 real :: weight ! Flux accumulation weight of the current fluxes [nondim].
485 real :: dt_coupling ! The coupling time step [T ~> s].
486 real :: dt_therm ! A limited and quantized version of OS%dt_therm [T ~> s].
487 real :: dt_dyn ! The dynamics time step [T ~> s].
488 real :: dtdia ! The diabatic time step [T ~> s].
489 real :: t_elapsed_seg ! The elapsed time in this update segment [T ~> s].
490 integer :: n ! The internal iteration counter.
491 integer :: nts ! The number of baroclinic dynamics time steps in a thermodynamic step.
492 integer :: n_max ! The number of calls to step_MOM dynamics in this call to update_ocean_model.
493 integer :: n_last_thermo ! The iteration number the last time thermodynamics were updated.
494 logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans multiple dynamic timesteps.
495 logical :: do_dyn ! If true, step the ocean dynamics and transport.
496 logical :: do_thermo ! If true, step the ocean thermodynamics.
497 logical :: step_thermo ! If true, take a thermodynamic step.
498 integer :: is, ie, js, je
499
500 call calltree_enter("update_ocean_model(), ocean_model_MOM.F90")
501 dt_coupling = time_to_real(ocean_coupling_time_step, scale=os%US%s_to_T)
502
503 if (.not.associated(os)) then
504 call mom_error(fatal, "update_ocean_model called with an unassociated "// &
505 "ocean_state_type structure. ocean_model_init must be "// &
506 "called first to allocate this structure.")
507 return
508 endif
509
510 do_dyn = .true. ; if (present(update_dyn)) do_dyn = update_dyn
511 do_thermo = .true. ; if (present(update_thermo)) do_thermo = update_thermo
512
513 if (do_thermo .and. (time_start_update /= os%Time)) &
514 call mom_error(warning, "update_ocean_model: internal clock does not "//&
515 "agree with time_start_update argument.")
516 if (do_dyn .and. (time_start_update /= os%Time_dyn)) &
517 call mom_error(warning, "update_ocean_model: internal dynamics clock does not "//&
518 "agree with time_start_update argument.")
519
520 if (.not.(do_dyn .or. do_thermo)) call mom_error(fatal, &
521 "update_ocean_model called without updating either dynamics or thermodynamics.")
522 if (do_dyn .and. do_thermo .and. (os%Time /= os%Time_dyn)) call mom_error(fatal, &
523 "update_ocean_model called to update both dynamics and thermodynamics with inconsistent clocks.")
524
525 ! This is benign but not necessary if ocean_model_init_sfc was called or if
526 ! OS%sfc_state%tr_fields was spawned in ocean_model_init. Consider removing it.
527 is = os%grid%isc ; ie = os%grid%iec ; js = os%grid%jsc ; je = os%grid%jec
528 call coupler_type_spawn(ocean_sfc%fields, os%sfc_state%tr_fields, &
529 (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
530
531 ! Translate Ice_ocean_boundary into fluxes and forces.
532 call get_domain_extent(ocean_sfc%Domain, index_bnds(1), index_bnds(2), index_bnds(3), index_bnds(4))
533
534 if (do_dyn) then
535 call convert_iob_to_forces(ice_ocean_boundary, os%forces, index_bnds, os%Time_dyn, os%grid, os%US, &
536 os%forcing_CSp, dt_forcing=dt_coupling, reset_avg=os%fluxes%fluxes_used)
537 if (os%use_ice_shelf) &
538 call add_shelf_forces(os%grid, os%US, os%Ice_shelf_CSp, os%forces)
539 if (os%icebergs_alter_ocean) &
540 call iceberg_forces(os%grid, os%forces, os%use_ice_shelf, &
541 os%sfc_state, dt_coupling, os%marine_ice_CSp)
542 endif
543
544 if (do_thermo) then
545 if (os%fluxes%fluxes_used) then
546 call convert_iob_to_fluxes(ice_ocean_boundary, os%fluxes, index_bnds, os%Time, dt_coupling, &
547 os%grid, os%US, os%forcing_CSp, os%sfc_state)
548
549 ! Add ice shelf fluxes
550 if (os%use_ice_shelf) &
551 call shelf_calc_flux(os%sfc_state, os%fluxes, os%Time, dt_coupling, os%Ice_shelf_CSp)
552 if (os%icebergs_alter_ocean) &
553 call iceberg_fluxes(os%grid, os%US, os%fluxes, os%use_ice_shelf, &
554 os%sfc_state, dt_coupling, os%marine_ice_CSp)
555
556#ifdef _USE_GENERIC_TRACER
557 call enable_averages(dt_coupling, os%Time + ocean_coupling_time_step, os%diag) !Is this needed?
558 call mom_generic_tracer_fluxes_accumulate(os%fluxes, 1.0) ! Here weight=1, so just store the current fluxes
559 call disable_averaging(os%diag)
560#endif
561 else
562 ! The previous fluxes have not been used yet, so translate the input fluxes
563 ! into a temporary type and then accumulate them in about 20 lines.
564 os%flux_tmp%C_p = os%fluxes%C_p
565 call convert_iob_to_fluxes(ice_ocean_boundary, os%flux_tmp, index_bnds, os%Time, dt_coupling, &
566 os%grid, os%US, os%forcing_CSp, os%sfc_state)
567
568 if (os%use_ice_shelf) &
569 call shelf_calc_flux(os%sfc_state, os%flux_tmp, os%Time,dt_coupling, os%Ice_shelf_CSp)
570 if (os%icebergs_alter_ocean) &
571 call iceberg_fluxes(os%grid, os%US, os%flux_tmp, os%use_ice_shelf, &
572 os%sfc_state, dt_coupling, os%marine_ice_CSp)
573
574 call fluxes_accumulate(os%flux_tmp, os%fluxes, os%grid, weight)
575#ifdef _USE_GENERIC_TRACER
576 ! Incorporate the current tracer fluxes into the running averages
577 call mom_generic_tracer_fluxes_accumulate(os%flux_tmp, weight)
578#endif
579 endif
580 endif
581
582 ! The net mass forcing is not currently used in the MOM6 dynamics solvers, so this is may be unnecessary.
583 if (do_dyn .and. associated(os%forces%net_mass_src) .and. .not.os%forces%net_mass_src_set) &
584 call get_net_mass_forcing(os%fluxes, os%grid, os%US, os%forces%net_mass_src)
585
586 if (os%use_waves .and. do_thermo) then
587 ! For now, the waves are only updated on the thermodynamics steps, because that is where
588 ! the wave intensities are actually used to drive mixing. At some point, the wave updates
589 ! might also need to become a part of the ocean dynamics, according to B. Reichl.
590 call update_surface_waves(os%grid, os%GV, os%US, os%time, ocean_coupling_time_step, os%waves)
591 endif
592
593 if ((os%nstep==0) .and. (os%nstep_thermo==0)) then ! This is the first call to update_ocean_model.
594 call finish_mom_initialization(os%Time, os%dirs, os%MOM_CSp)
595 endif
596
597 time_thermo_start = os%Time
598 time_seg_start = os%Time ; if (do_dyn) time_seg_start = os%Time_dyn
599 time1 = time_seg_start
600
601 if (os%offline_tracer_mode .and. do_thermo) then
602 call step_offline(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp)
603 elseif ((.not.do_thermo) .or. (.not.do_dyn)) then
604 ! The call sequence is being orchestrated from outside of update_ocean_model.
605 if (present(cycle_length)) then
606 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, &
607 waves=os%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, &
608 start_cycle=start_cycle, end_cycle=end_cycle, cycle_length=os%US%s_to_T*cycle_length, &
609 reset_therm=ocn_fluxes_used)
610 else
611 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, &
612 waves=os%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, &
613 start_cycle=start_cycle, end_cycle=end_cycle, reset_therm=ocn_fluxes_used)
614 endif
615 elseif (os%single_step_call) then
616 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_coupling, os%MOM_CSp, waves=os%Waves)
617 else ! Step both the dynamics and thermodynamics with separate calls.
618 n_max = 1 ; if (dt_coupling > os%dt) n_max = ceiling(dt_coupling/os%dt - 0.001)
619 dt_dyn = dt_coupling / real(n_max)
620 thermo_does_span_coupling = (os%thermo_spans_coupling .and. (os%dt_therm > 1.5*dt_coupling))
621
622 if (thermo_does_span_coupling) then
623 dt_therm = dt_coupling * floor(os%dt_therm / dt_coupling + 0.001)
624 nts = floor(dt_therm/dt_dyn + 0.001)
625 else
626 nts = max(1,min(n_max,floor(os%dt_therm/dt_dyn + 0.001)))
627 n_last_thermo = 0
628 endif
629
630 time1 = time_seg_start ; t_elapsed_seg = 0.0
631 do n=1,n_max
632 if (os%diabatic_first) then
633 if (thermo_does_span_coupling) call mom_error(fatal, &
634 "MOM is not yet set up to have restarts that work with "//&
635 "THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
636 if (modulo(n-1,nts)==0) then
637 dtdia = dt_dyn*min(nts,n_max-(n-1))
638 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dtdia, os%MOM_CSp, &
639 waves=os%Waves, do_dynamics=.false., do_thermodynamics=.true., &
640 start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
641 endif
642
643 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_dyn, os%MOM_CSp, &
644 waves=os%Waves, do_dynamics=.true., do_thermodynamics=.false., &
645 start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
646 else
647 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dt_dyn, os%MOM_CSp, &
648 waves=os%Waves, do_dynamics=.true., do_thermodynamics=.false., &
649 start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
650
651 step_thermo = .false.
652 if (thermo_does_span_coupling) then
653 dtdia = dt_therm
654 step_thermo = mom_state_is_synchronized(os%MOM_CSp, adv_dyn=.true.)
655 elseif ((modulo(n,nts)==0) .or. (n==n_max)) then
656 dtdia = dt_dyn*(n - n_last_thermo)
657 n_last_thermo = n
658 step_thermo = .true.
659 endif
660
661 if (step_thermo) then
662 ! Back up Time1 to the start of the thermodynamic segment.
663 time1 = time1 - real_to_time(dtdia - dt_dyn, unscale=os%US%T_to_s)
664 call step_mom(os%forces, os%fluxes, os%sfc_state, time1, dtdia, os%MOM_CSp, &
665 waves=os%Waves, do_dynamics=.false., do_thermodynamics=.true., &
666 start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
667 endif
668 endif
669
670 t_elapsed_seg = t_elapsed_seg + dt_dyn
671 time1 = time_seg_start + real_to_time(t_elapsed_seg, unscale=os%US%T_to_s)
672 enddo
673 endif
674
675 if (do_dyn) os%Time_dyn = time_seg_start + ocean_coupling_time_step
676 if (do_dyn) os%nstep = os%nstep + 1
677 os%Time = time_thermo_start ! Reset the clock to compensate for shared pointers.
678 if (do_thermo) os%Time = os%Time + ocean_coupling_time_step
679 if (do_thermo) os%nstep_thermo = os%nstep_thermo + 1
680
681 if (do_dyn) then
682 call mech_forcing_diags(os%forces, dt_coupling, os%grid, os%Time_dyn, os%diag, os%forcing_CSp%handles)
683 endif
684
685 if (os%fluxes%fluxes_used .and. do_thermo) then
686 call forcing_diagnostics(os%fluxes, os%sfc_state, os%grid, os%US, os%Time, os%diag, os%forcing_CSp%handles)
687 endif
688
689 !only ,ale ice-shelf frazil adjustments if sfc_state%frazil was updated (do_thermo=True)
690 if (do_thermo .and. os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
691 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
692
693! Translate state into Ocean.
694! call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
695! OS%fluxes%p_surf_full, OS%press_to_z)
696 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
697 if (os%calve_ice_shelf_bergs) call convert_shelf_state_to_ocean_type(ocean_sfc,os%Ice_shelf_CSp, os%US)
698 time1 = os%Time ; if (do_dyn) time1 = os%Time_dyn
699 call coupler_type_send_data(ocean_sfc%fields, time1)
700
701 call calltree_leave("update_ocean_model()")
702end subroutine update_ocean_model
703
704!> This subroutine writes out the ocean model restart file.
705subroutine ocean_model_restart(OS, timestamp)
706 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
707 !! internal ocean state being saved to a restart file
708 character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be
709 !! prepended to the file name. (Currently this is unused.)
710
711 if (.not.mom_state_is_synchronized(os%MOM_CSp)) &
712 call mom_error(warning, "End of MOM_main reached with inconsistent "//&
713 "dynamics and advective times. Additional restart fields "//&
714 "that have not been coded yet would be required for reproducibility.")
715 if (.not.os%fluxes%fluxes_used) call mom_error(fatal, "ocean_model_restart "//&
716 "was called with unused buoyancy fluxes. For conservation, the ocean "//&
717 "restart files can only be created after the buoyancy forcing is applied.")
718
719 if (btest(os%Restart_control,1)) then
720 call save_mom_restart(os%MOM_CSp, os%dirs%restart_output_dir, os%Time, &
721 os%grid, time_stamped=.true., gv=os%GV)
722 call forcing_save_restart(os%forcing_CSp, os%grid, os%Time, &
723 os%dirs%restart_output_dir, .true.)
724 if (os%use_ice_shelf) then
725 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir, .true.)
726 endif
727 endif
728 if (btest(os%Restart_control,0)) then
729 call save_mom_restart(os%MOM_CSp, os%dirs%restart_output_dir, os%Time, &
730 os%grid, gv=os%GV)
731 call forcing_save_restart(os%forcing_CSp, os%grid, os%Time, &
732 os%dirs%restart_output_dir)
733 if (os%use_ice_shelf) then
734 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir)
735 endif
736 endif
737
738end subroutine ocean_model_restart
739! </SUBROUTINE> NAME="ocean_model_restart"
740
741!> ocean_model_end terminates the model run, saving the ocean state in a restart
742!! and deallocating any data associated with the ocean.
743subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time)
744 type(ocean_public_type), intent(inout) :: ocean_sfc !< An ocean_public_type structure that is
745 !! to be deallocated upon termination.
746 type(ocean_state_type), pointer :: ocean_state !< A pointer to the structure containing
747 !! the internal ocean state to be deallocated
748 !! upon termination.
749 type(time_type), intent(in) :: time !< The model time, used for writing restarts.
750
751 call ocean_model_save_restart(ocean_state, time)
752 call diag_mediator_end(time, ocean_state%diag)
753 call mom_end(ocean_state%MOM_CSp)
754 if (ocean_state%use_ice_shelf) call ice_shelf_end(ocean_state%Ice_shelf_CSp)
755end subroutine ocean_model_end
756
757!> ocean_model_save_restart causes restart files associated with the ocean to be
758!! written out.
759subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix)
760 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
761 !! internal ocean state (in).
762 type(time_type), intent(in) :: time !< The model time at this call, needed for writing files.
763 character(len=*), optional, intent(in) :: directory !< An optional directory into which to
764 !! write these restart files.
765 character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-stamp)
766 !! to append to the restart file names.
767! Note: This is a new routine - it will need to exist for the new incremental
768! checkpointing. It will also be called by ocean_model_end, giving the same
769! restart behavior as now in FMS.
770 character(len=200) :: restart_dir
771
772 if (.not.mom_state_is_synchronized(os%MOM_CSp)) &
773 call mom_error(warning, "ocean_model_save_restart called with inconsistent "//&
774 "dynamics and advective times. Additional restart fields "//&
775 "that have not been coded yet would be required for reproducibility.")
776 if (.not.os%fluxes%fluxes_used) call mom_error(fatal, "ocean_model_save_restart "//&
777 "was called with unused buoyancy fluxes. For conservation, the ocean "//&
778 "restart files can only be created after the buoyancy forcing is applied.")
779
780 if (present(directory)) then ; restart_dir = directory
781 else ; restart_dir = os%dirs%restart_output_dir ; endif
782
783 call save_mom_restart(os%MOM_CSp, restart_dir, time, os%grid, gv=os%GV)
784
785 call forcing_save_restart(os%forcing_CSp, os%grid, time, restart_dir)
786
787 if (os%use_ice_shelf) then
788 call ice_shelf_save_restart(os%Ice_shelf_CSp, os%Time, os%dirs%restart_output_dir)
789 endif
790end subroutine ocean_model_save_restart
791
792!> Initialize the public ocean type
793subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_fields_ocn)
794 type(mom_domain_type), intent(in) :: input_domain !< The ocean model domain description
795 type(ocean_public_type), intent(inout) :: Ocean_sfc !< A structure containing various publicly
796 !! visible ocean surface properties after
797 !! initialization, whose elements are allocated here.
798 type(diag_ctrl), intent(in) :: diag !< A structure that regulates diagnostic output
799 type(coupler_1d_bc_type), &
800 optional, intent(in) :: gas_fields_ocn !< If present, this type describes the
801 !! ocean and surface-ice fields that will participate
802 !! in the calculation of additional gas or other
803 !! tracer fluxes.
804
805 integer :: xsz, ysz, layout(2)
806 ! ice-ocean-boundary fields are always allocated using absolute indices
807 ! and have no halos.
808 integer :: isc, iec, jsc, jec
809
810 call clone_mom_domain(input_domain, ocean_sfc%Domain, halo_size=0, symmetric=.false.)
811
812 call get_domain_extent(ocean_sfc%Domain, isc, iec, jsc, jec)
813
814 allocate ( ocean_sfc%t_surf (isc:iec,jsc:jec), &
815 ocean_sfc%s_surf (isc:iec,jsc:jec), &
816 ocean_sfc%u_surf (isc:iec,jsc:jec), &
817 ocean_sfc%v_surf (isc:iec,jsc:jec), &
818 ocean_sfc%sea_lev(isc:iec,jsc:jec), &
819 ocean_sfc%calving(isc:iec,jsc:jec), &
820 ocean_sfc%calving_hflx(isc:iec,jsc:jec), &
821 ocean_sfc%area (isc:iec,jsc:jec), &
822 ocean_sfc%melt_potential(isc:iec,jsc:jec), &
823 ocean_sfc%OBLD (isc:iec,jsc:jec), &
824 ocean_sfc%frazil (isc:iec,jsc:jec))
825
826 ocean_sfc%t_surf(:,:) = 0.0 ! time averaged sst (Kelvin) passed to atmosphere/ice model
827 ocean_sfc%s_surf(:,:) = 0.0 ! time averaged sss (psu) passed to atmosphere/ice models
828 ocean_sfc%u_surf(:,:) = 0.0 ! time averaged u-current (m/sec) passed to atmosphere/ice models
829 ocean_sfc%v_surf(:,:) = 0.0 ! time averaged v-current (m/sec) passed to atmosphere/ice models
830 ocean_sfc%sea_lev(:,:) = 0.0 ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav
831 ocean_sfc%calving(:,:) = 0.0 ! time accumulated ice sheet calving (kg m-2) passed to ice model
832 ocean_sfc%calving_hflx(:,:) = 0.0 ! time accumulated ice sheet calving heat flux (W m-2) passed to ice model
833 ocean_sfc%frazil(:,:) = 0.0 ! time accumulated frazil (J/m^2) passed to ice model
834 ocean_sfc%melt_potential(:,:) = 0.0 ! time accumulated melt potential (J/m^2) passed to ice model
835 ocean_sfc%OBLD(:,:) = 0.0 ! ocean boundary layer depth (m)
836 ocean_sfc%area(:,:) = 0.0
837 ocean_sfc%axes = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics
838
839 if (present(gas_fields_ocn)) then
840 call coupler_type_spawn(gas_fields_ocn, ocean_sfc%fields, (/isc,isc,iec,iec/), &
841 (/jsc,jsc,jec,jec/), suffix = '_ocn', as_needed=.true.)
842 endif
843
844end subroutine initialize_ocean_public_type
845
846!> This subroutine translates the coupler's ocean_data_type into MOM's
847!! surface state variable. This may eventually be folded into the MOM
848!! code that calculates the surface state in the first place.
849!! Note the offset in the arrays because the ocean_data_type has no
850!! halo points in its arrays and always uses absolute indices.
851subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_to_z)
852 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
853 !! describe the surface state of the ocean.
854 type(ocean_public_type), &
855 target, intent(inout) :: Ocean_sfc !< A structure containing various publicly
856 !! visible ocean surface fields, whose elements
857 !! have their data set here.
858 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
859 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
860 real, optional, intent(in) :: patm(:,:) !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
861 real, optional, intent(in) :: press_to_z !< A conversion factor between pressure and ocean
862 !! depth, usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1]
863 ! Local variables
864 character(len=48) :: val_str
865 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
866 integer :: i, j, i0, j0, is, ie, js, je
867
868 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
869 call pass_vector(sfc_state%u, sfc_state%v, g%Domain)
870
871 call get_domain_extent(ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
872 if (present(patm)) then
873 ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd).
874 if (.not.present(press_to_z)) call mom_error(fatal, &
875 'convert_state_to_ocean_type: press_to_z must be present if patm is.')
876 endif
877
878 i0 = is - isc_bnd ; j0 = js - jsc_bnd
879 if (sfc_state%T_is_conT) then
880 ! Convert the surface T from conservative T to potential T.
881 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
882 ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(us%S_to_ppt*sfc_state%SSS(i+i0,j+j0), &
883 us%C_to_degC*sfc_state%SST(i+i0,j+j0)) + celsius_kelvin_offset
884 enddo ; enddo
885 else
886 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
887 ocean_sfc%t_surf(i,j) = us%C_to_degC*sfc_state%SST(i+i0,j+j0) + celsius_kelvin_offset
888 enddo ; enddo
889 endif
890 if (sfc_state%S_is_absS) then
891 ! Convert the surface S from absolute salinity to practical salinity.
892 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
893 ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(us%S_to_ppt*sfc_state%SSS(i+i0,j+j0))
894 enddo ; enddo
895 else
896 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
897 ocean_sfc%s_surf(i,j) = us%S_to_ppt*sfc_state%SSS(i+i0,j+j0)
898 enddo ; enddo
899 endif
900
901 if (present(patm)) then
902 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
903 ocean_sfc%sea_lev(i,j) = us%Z_to_m * (sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z)
904 ocean_sfc%area(i,j) = us%L_to_m**2 * g%areaT(i+i0,j+j0)
905 enddo ; enddo
906 else
907 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
908 ocean_sfc%sea_lev(i,j) = us%Z_to_m * sfc_state%sea_lev(i+i0,j+j0)
909 ocean_sfc%area(i,j) = us%L_to_m**2 * g%areaT(i+i0,j+j0)
910 enddo ; enddo
911 endif
912
913 if (allocated(sfc_state%frazil)) then
914 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
915 ocean_sfc%frazil(i,j) = us%Q_to_J_kg*us%RZ_to_kg_m2 * sfc_state%frazil(i+i0,j+j0)
916 enddo ; enddo
917 endif
918
919 if (allocated(sfc_state%melt_potential)) then
920 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
921 ocean_sfc%melt_potential(i,j) = us%Q_to_J_kg*us%RZ_to_kg_m2 * sfc_state%melt_potential(i+i0,j+j0)
922 enddo ; enddo
923 endif
924
925 if (allocated(sfc_state%Hml)) then
926 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
927 ocean_sfc%OBLD(i,j) = us%Z_to_m * sfc_state%Hml(i+i0,j+j0)
928 enddo ; enddo
929 endif
930
931 if (ocean_sfc%stagger == agrid) then
932 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
933 ocean_sfc%u_surf(i,j) = g%mask2dT(i+i0,j+j0) * us%L_T_to_m_s * &
934 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i-1+i0,j+j0))
935 ocean_sfc%v_surf(i,j) = g%mask2dT(i+i0,j+j0) * us%L_T_to_m_s * &
936 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0,j-1+j0))
937 enddo ; enddo
938 elseif (ocean_sfc%stagger == bgrid_ne) then
939 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
940 ocean_sfc%u_surf(i,j) = g%mask2dBu(i+i0,j+j0) * us%L_T_to_m_s * &
941 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i+i0,j+j0+1))
942 ocean_sfc%v_surf(i,j) = g%mask2dBu(i+i0,j+j0) * us%L_T_to_m_s * &
943 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0+1,j+j0))
944 enddo ; enddo
945 elseif (ocean_sfc%stagger == cgrid_ne) then
946 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
947 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)
948 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)
949 enddo ; enddo
950 else
951 write(val_str, '(I8)') ocean_sfc%stagger
952 call mom_error(fatal, "convert_state_to_ocean_type: "//&
953 "Ocean_sfc%stagger has the unrecognized value of "//trim(val_str))
954 endif
955
956 if (coupler_type_initialized(sfc_state%tr_fields)) then
957 if (.not.coupler_type_initialized(ocean_sfc%fields)) then
958 call mom_error(fatal, "convert_state_to_ocean_type: "//&
959 "Ocean_sfc%fields has not been initialized.")
960 endif
961 call coupler_type_copy_data(sfc_state%tr_fields, ocean_sfc%fields)
962 endif
963
964end subroutine convert_state_to_ocean_type
965
966!> Converts the ice-shelf-to-ocean calving and calving_hflx variables from the ice-shelf state (ISS) type
967!! to the ocean public type
968subroutine convert_shelf_state_to_ocean_type(Ocean_sfc, CS, US)
969 type(ocean_public_type), &
970 target, intent(inout) :: Ocean_sfc !< A structure containing various publicly
971 !! visible ocean surface fields, whose elements
972 !! have their data set here.
973 type(ice_shelf_cs), pointer :: CS !< A pointer to the ice shelf control structure
974 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
975 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd, i, j
976
977 call get_domain_extent(ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
978
979 call ice_sheet_calving_to_ocean_sfc(cs,us,ocean_sfc%calving(isc_bnd:iec_bnd,jsc_bnd:jec_bnd),&
980 ocean_sfc%calving_hflx(isc_bnd:iec_bnd,jsc_bnd:jec_bnd))
981
983
984!> This subroutine extracts the surface properties from the ocean's internal
985!! state and stores them in the ocean type returned to the calling ice model.
986!! It has to be separate from the ocean_initialization call because the coupler
987!! module allocates the space for some of these variables.
988subroutine ocean_model_init_sfc(OS, Ocean_sfc)
989 type(ocean_state_type), pointer :: os !< The structure with the complete ocean state
990 type(ocean_public_type), intent(inout) :: ocean_sfc !< A structure containing various publicly
991 !! visible ocean surface properties after initialization, whose
992 !! elements have their data set here.
993 integer :: is, ie, js, je
994
995 is = os%grid%isc ; ie = os%grid%iec ; js = os%grid%jsc ; je = os%grid%jec
996 call coupler_type_spawn(ocean_sfc%fields, os%sfc_state%tr_fields, &
997 (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
998
999 call extract_surface_state(os%MOM_CSp, os%sfc_state)
1000
1001 if (os%use_ice_shelf .and. allocated(os%sfc_state%frazil)) &
1002 call adjust_ice_sheet_frazil(os%sfc_state, os%fluxes, os%Ice_shelf_CSp)
1003
1004 call convert_state_to_ocean_type(os%sfc_state, ocean_sfc, os%grid, os%US)
1005
1006end subroutine ocean_model_init_sfc
1007
1008!> ocean_model_flux_init is used to initialize properties of the air-sea fluxes
1009!! as determined by various run-time parameters. It can be called from
1010!! non-ocean PEs, or PEs that have not yet been initialized, and it can safely
1011!! be called multiple times.
1012subroutine ocean_model_flux_init(OS, verbosity)
1013 type(ocean_state_type), optional, pointer :: os !< An optional pointer to the ocean state,
1014 !! used to figure out if this is an ocean PE that
1015 !! has already been initialized.
1016 integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity.
1017
1018 logical :: os_is_set
1019 integer :: verbose
1020
1021 os_is_set = .false. ; if (present(os)) os_is_set = associated(os)
1022
1023 ! Use this to control the verbosity of output; consider rethinking this logic later.
1024 verbose = 5 ; if (os_is_set) verbose = 3
1025 if (present(verbosity)) verbose = verbosity
1026
1027 call call_tracer_flux_init(verbosity=verbose)
1028
1029end subroutine ocean_model_flux_init
1030
1031!> Ocean_stock_pe - returns the integrated stocks of heat, water, etc. for conservation checks.
1032!! Because of the way FMS is coded, only the root PE has the integrated amount,
1033!! while all other PEs get 0.
1034subroutine ocean_stock_pe(OS, index, value, time_index)
1035 use stock_constants_mod, only : istock_water, istock_heat,istock_salt
1036 type(ocean_state_type), pointer :: os !< A structure containing the internal ocean state.
1037 !! The data in OS is intent in.
1038 integer, intent(in) :: index !< The stock index for the quantity of interest.
1039 real, intent(out) :: value !< Sum returned for the conservation quantity of interest [various]
1040 integer, optional, intent(in) :: time_index !< An unused optional argument, present only for
1041 !! interfacial compatibility with other models.
1042! Arguments: OS - A structure containing the internal ocean state.
1043! (in) index - Index of conservation quantity of interest.
1044! (in) value - Sum returned for the conservation quantity of interest.
1045! (in,opt) time_index - Index for time level to use if this is necessary.
1046
1047 real :: salt ! The total salt in the ocean [kg]
1048
1049 value = 0.0
1050 if (.not.associated(os)) return
1051 if (.not.os%is_ocean_pe) return
1052
1053 select case (index)
1054 case (istock_water) ! Return the mass of fresh water in the ocean in [kg].
1055 if (os%GV%Boussinesq) then
1056 call get_ocean_stocks(os%MOM_CSp, mass=value, on_pe_only=.true.)
1057 else ! In non-Boussinesq mode, the mass of salt needs to be subtracted.
1058 call get_ocean_stocks(os%MOM_CSp, mass=value, salt=salt, on_pe_only=.true.)
1059 value = value - salt
1060 endif
1061 case (istock_heat) ! Return the heat content of the ocean in [J].
1062 call get_ocean_stocks(os%MOM_CSp, heat=value, on_pe_only=.true.)
1063 case (istock_salt) ! Return the mass of the salt in the ocean in [kg].
1064 call get_ocean_stocks(os%MOM_CSp, salt=value, on_pe_only=.true.)
1065 case default ; value = 0.0
1066 end select
1067 ! If the FMS coupler is changed so that Ocean_stock_PE is only called on
1068 ! ocean PEs, uncomment the following and eliminate the on_PE_only flags above.
1069 ! if (.not.is_root_pe()) value = 0.0
1070
1071end subroutine ocean_stock_pe
1072
1073!> This subroutine extracts a named 2-D field from the ocean surface or public type
1074subroutine ocean_model_data2d_get(OS, Ocean, name, array2D, isc, jsc)
1075 use mom_constants, only : celsius_kelvin_offset
1076 type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the
1077 !! internal ocean state (intent in).
1078 type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly
1079 !! visible ocean surface fields.
1080 character(len=*) , intent(in) :: name !< The name of the field to extract
1081 integer , intent(in) :: isc !< The starting i-index of array2D
1082 integer , intent(in) :: jsc !< The starting j-index of array2D
1083 real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must
1084 !! cover only the computational domain [various]
1085
1086 integer :: g_isc, g_iec, g_jsc, g_jec, g_isd, g_ied, g_jsd, g_jed, i, j
1087
1088 if (.not.associated(os)) return
1089 if (.not.os%is_ocean_pe) return
1090
1091 ! The problem is that %areaT is on MOM domain but Ice_Ocean_Boundary%... is on a haloless domain.
1092 ! We want to return the MOM data on the haloless (compute) domain
1093 call get_domain_extent(os%grid%Domain, g_isc, g_iec, g_jsc, g_jec, g_isd, g_ied, g_jsd, g_jed)
1094
1095 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
1096
1097 select case(name)
1098 case('area')
1099 array2d(isc:,jsc:) = os%US%L_to_m**2*os%grid%areaT(g_isc:g_iec,g_jsc:g_jec)
1100 case('mask')
1101 array2d(isc:,jsc:) = os%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec)
1102!OR same result
1103! do j=g_jsc,g_jec ; do i=g_isc,g_iec
1104! array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j)
1105! enddo ; enddo
1106 case('t_surf')
1107 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1108 case('t_pme')
1109 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1110 case('t_runoff')
1111 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1112 case('t_calving')
1113 array2d(isc:,jsc:) = ocean%t_surf(isc:,jsc:)-celsius_kelvin_offset
1114 case('btfHeat')
1115 array2d(isc:,jsc:) = 0
1116 case('cos_rot')
1117 array2d(isc:,jsc:) = os%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1
1118 case('sin_rot')
1119 array2d(isc:,jsc:) = os%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0
1120 case('s_surf')
1121 array2d(isc:,jsc:) = ocean%s_surf(isc:,jsc:)
1122 case('sea_lev')
1123 array2d(isc:,jsc:) = ocean%sea_lev(isc:,jsc:)
1124 case('frazil')
1125 array2d(isc:,jsc:) = ocean%frazil(isc:,jsc:)
1126 case('melt_pot')
1127 array2d(isc:,jsc:) = ocean%melt_potential(isc:,jsc:)
1128 case('obld')
1129 array2d(isc:,jsc:) = ocean%OBLD(isc:,jsc:)
1130 case default
1131 call mom_error(fatal,'get_ocean_grid_data2D: unknown argument name='//name)
1132 end select
1133
1134end subroutine ocean_model_data2d_get
1135
1136!> This subroutine extracts a named scalar field from the ocean surface or public type
1137subroutine ocean_model_data1d_get(OS, Ocean, name, value)
1138 type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the
1139 !! internal ocean state (intent in).
1140 type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly
1141 !! visible ocean surface fields.
1142 character(len=*), intent(in) :: name !< The name of the field to extract
1143 real, intent(out):: value !< The value of the named field [various]
1144
1145 if (.not.associated(os)) return
1146 if (.not.os%is_ocean_pe) return
1147
1148 select case(name)
1149 case('c_p')
1150 value = os%C_p
1151 case default
1152 call mom_error(fatal,'get_ocean_grid_data1D: unknown argument name='//name)
1153 end select
1154
1155end subroutine ocean_model_data1d_get
1156
1157!> Write out checksums for fields from the ocean surface state
1158subroutine ocean_public_type_chksum(id, timestep, ocn)
1159
1160 character(len=*), intent(in) :: id !< An identifying string for this call
1161 integer, intent(in) :: timestep !< The number of elapsed timesteps
1162 type(ocean_public_type), intent(in) :: ocn !< A structure containing various publicly
1163 !! visible ocean surface fields.
1164 ! Local variables
1165 integer(kind=int64) :: chks ! A checksum for the field
1166 logical :: root ! True only on the root PE
1167 integer :: outunit ! The output unit to write to
1168
1169 root = is_root_pe()
1170 outunit = stdout_if_root()
1171
1172 if (root) write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep
1173 chks = field_chksum(ocn%t_surf ) ; if (root) write(outunit,100) 'ocean%t_surf ', chks
1174 chks = field_chksum(ocn%s_surf ) ; if (root) write(outunit,100) 'ocean%s_surf ', chks
1175 chks = field_chksum(ocn%u_surf ) ; if (root) write(outunit,100) 'ocean%u_surf ', chks
1176 chks = field_chksum(ocn%v_surf ) ; if (root) write(outunit,100) 'ocean%v_surf ', chks
1177 chks = field_chksum(ocn%sea_lev) ; if (root) write(outunit,100) 'ocean%sea_lev ', chks
1178 chks = field_chksum(ocn%frazil ) ; if (root) write(outunit,100) 'ocean%frazil ', chks
1179 chks = field_chksum(ocn%melt_potential) ; if (root) write(outunit,100) 'ocean%melt_potential ', chks
1180 call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%')
1181100 FORMAT(" CHECKSUM::",a20," = ",z20)
1182
1183end subroutine ocean_public_type_chksum
1184
1185!> This subroutine gives a handle to the grid from ocean state
1186subroutine get_ocean_grid(OS, Gridp)
1187 ! Obtain the ocean grid.
1188 type(ocean_state_type) :: os !< A structure containing the
1189 !! internal ocean state
1190 type(ocean_grid_type) , pointer :: gridp !< The ocean's grid structure
1191
1192 gridp => os%grid
1193 return
1194end subroutine get_ocean_grid
1195
1196!> This subroutine extracts a named (u- or v-) 2-D surface current from ocean internal state
1197subroutine ocean_model_get_uv_surf(OS, Ocean, name, array2D, isc, jsc)
1198
1199 type(ocean_state_type), pointer :: os !< A pointer to the structure containing the
1200 !! internal ocean state (intent in).
1201 type(ocean_public_type), intent(in) :: ocean !< A structure containing various publicly
1202 !! visible ocean surface fields.
1203 character(len=*) , intent(in) :: name !< The name of the current (ua or va) to extract
1204 integer , intent(in) :: isc !< The starting i-index of array2D
1205 integer , intent(in) :: jsc !< The starting j-index of array2D
1206 real, dimension(isc:,jsc:), intent(out):: array2d !< The values of the named field, it must
1207 !! cover only the computational domain [L T-1 ~> m s-1]
1208
1209 type(ocean_grid_type) , pointer :: g !< The ocean's grid structure
1210 type(surface), pointer :: sfc_state !< A structure containing fields that
1211 !! describe the surface state of the ocean.
1212
1213 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
1214 integer :: i, j, i0, j0
1215 integer :: is, ie, js, je
1216
1217 if (.not.associated(os)) return
1218 if (.not.os%is_ocean_pe) return
1219
1220 g => os%grid
1221 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1222
1223 call get_domain_extent(ocean%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
1224
1225 i0 = is - isc_bnd ; j0 = js - jsc_bnd
1226
1227 sfc_state => os%sfc_state
1228
1229 call pass_vector(sfc_state%u, sfc_state%v, g%Domain)
1230
1231 select case(name)
1232 case('ua')
1233 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1234 array2d(i,j) = g%mask2dT(i+i0,j+j0) * &
1235 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i-1+i0,j+j0))
1236 enddo ; enddo
1237 case('va')
1238 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1239 array2d(i,j) = g%mask2dT(i+i0,j+j0) * &
1240 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0,j-1+j0))
1241 enddo ; enddo
1242 case('ub')
1243 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1244 array2d(i,j) = g%mask2dBu(i+i0,j+j0) * &
1245 0.5*(sfc_state%u(i+i0,j+j0)+sfc_state%u(i+i0,j+j0+1))
1246 enddo ; enddo
1247 case('vb')
1248 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1249 array2d(i,j) = g%mask2dBu(i+i0,j+j0) * &
1250 0.5*(sfc_state%v(i+i0,j+j0)+sfc_state%v(i+i0+1,j+j0))
1251 enddo ; enddo
1252 case('uc')
1253 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1254 array2d(i,j) = g%mask2dCu(i+i0,j+j0) * sfc_state%u(i+i0,j+j0)
1255 enddo ; enddo
1256 case('vc')
1257 do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
1258 array2d(i,j) = g%mask2dCv(i+i0,j+j0) * sfc_state%v(i+i0,j+j0)
1259 enddo ; enddo
1260 case default
1261 call mom_error(fatal,'ocean_model_get_UV_surf: unknown argument name='//name)
1262 end select
1263
1264end subroutine ocean_model_get_uv_surf
1265
1266end module ocean_model_mod