diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 5cd929c901..99e2feeabe 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -29,244 +29,338 @@ module ocn_comp_mct shr_file_setLogUnit, shr_file_setLogLevel ! MOM6 modules -use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end -use MOM, only: calculate_surface_state, allocate_surface_state -use MOM, only: finish_MOM_initialization, step_offline -use MOM_forcing_type, only: forcing, forcing_diagnostics, mech_forcing_diags, forcing_accumulate -use MOM_forcing_type, only: allocate_forcing_type -use MOM_surface_forcing, only: surface_forcing_init, convert_IOB_to_fluxes -use MOM_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS -use MOM_surface_forcing, only: forcing_save_restart -use MOM_restart, only: save_restart -use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here -use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE -use MOM_domains, only: pass_var, AGRID -use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_verticalGrid, only: verticalGrid_type -use MOM_variables, only: surface -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING -use MOM_error_handler, only: callTree_enter, callTree_leave -use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date -use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only: operator(/=), operator(>), get_time -use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file -use MOM_get_input, only: Get_MOM_Input, directories -use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging -use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end -use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS -use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart -use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS -use MOM_sum_output, only: write_energy, accumulate_net_input +use MOM_coms, only : reproducing_sum +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT +use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end +use MOM, only: calculate_surface_state, allocate_surface_state +use MOM, only: finish_MOM_initialization, step_offline +use MOM_forcing_type, only: forcing, forcing_diags, register_forcing_type_diags +use MOM_forcing_type, only: allocate_forcing_type, deallocate_forcing_type +use MOM_forcing_type, only: mech_forcing_diags, forcing_accumulate, forcing_diagnostics +use MOM_restart, only: save_restart +use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here +use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE, To_All +use MOM_domains, only: pass_var, AGRID, fill_symmetric_edges +use MOM_grid, only: ocean_grid_type, get_global_grid_size +use MOM_verticalGrid, only: verticalGrid_type +use MOM_variables, only: surface +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING +use MOM_error_handler, only: callTree_enter, callTree_leave +use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date +use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only: operator(/=), operator(>), get_time +use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file +use MOM_get_input, only: Get_MOM_Input, directories +use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging +use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end +use MOM_diag_mediator, only: safe_alloc_ptr +use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart +use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS +use MOM_sum_output, only: write_energy, accumulate_net_input use MOM_string_functions, only: uppercase -use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf +use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf, hlv use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct - -! FMS -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain -use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain - -! GFDL coupler +use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init +use user_revise_forcing, only : user_revise_forcing_CS +use MOM_restart, only : restart_init, MOM_restart_CS +use MOM_restart, only : restart_init_end, save_restart, restore_state +use data_override_mod, only : data_override_init, data_override +use MOM_io, only : slasher, write_version_number +use MOM_spatial_means, only : adjust_area_mean_to_zero + +! FMS modules +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain +use time_interp_external_mod, only : init_external_field, time_interp_external +use time_interp_external_mod, only : time_interp_external_init +use fms_mod, only : read_data + +! GFDL coupler modules use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type -use coupler_types_mod, only : coupler_type_spawn -use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data +use coupler_types_mod, only : coupler_type_spawn +use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data ! By default make data private implicit none; private - ! Public member functions - public :: ocn_init_mct - public :: ocn_run_mct - public :: ocn_final_mct - ! Flag for debugging - logical, parameter :: debug=.true. - - !> Structure with MCT attribute vectors and indices - type cpl_indices - - ! ocean to coupler - integer :: o2x_So_t !< Surface potential temperature (deg C) - integer :: o2x_So_u !< Surface zonal velocity (m/s) - integer :: o2x_So_v !< Surface meridional velocity (m/s) - integer :: o2x_So_s !< Surface salinity (PSU) - integer :: o2x_So_dhdx !< Zonal slope in the sea surface height - integer :: o2x_So_dhdy !< Meridional lope in the sea surface height - integer :: o2x_So_bldepth !< Boundary layer depth (m) - integer :: o2x_Fioo_q !< Heat flux? - integer :: o2x_Faoo_fco2_ocn!< CO2 flux - integer :: o2x_Faoo_fdms_ocn!< DMS flux - - ! coupler to ocean - integer :: x2o_Si_ifrac !< Fractional ice wrt ocean - integer :: x2o_So_duu10n !< 10m wind speed squared (m^2/s^2) - integer :: x2o_Sa_pslv !< Sea-level pressure (Pa) - integer :: x2o_Sa_co2prog !< Bottom atm level prognostic CO2 - integer :: x2o_Sa_co2diag !< Bottom atm level diagnostic CO2 - integer :: x2o_Sw_lamult !< Wave model langmuir multiplier - integer :: x2o_Sw_ustokes !< Surface Stokes drift, x-component - integer :: x2o_Sw_vstokes !< Surface Stokes drift, y-component - integer :: x2o_Foxx_taux !< Zonal wind stress (W/m2) - integer :: x2o_Foxx_tauy !< Meridonal wind stress (W/m2) - integer :: x2o_Foxx_swnet !< Net short-wave heat flux (W/m2) - integer :: x2o_Foxx_sen !< Sensible heat flux (W/m2) - integer :: x2o_Foxx_lat !< Latent heat flux (W/m2) - integer :: x2o_Foxx_lwup !< Longwave radiation, up (W/m2) - integer :: x2o_Faxa_lwdn !< Longwave radiation, down (W/m2) - integer :: x2o_Fioi_melth !< Heat flux from snow & ice melt (W/m2) - integer :: x2o_Fioi_meltw !< Snow melt flux (kg/m2/s) - integer :: x2o_Fioi_bcpho !< Black Carbon hydrophobic release - !! from sea ice component - integer :: x2o_Fioi_bcphi !< Black Carbon hydrophilic release from - !! sea ice component - integer :: x2o_Fioi_flxdst !< Dust release from sea ice component - integer :: x2o_Fioi_salt !< Salt flux (kg(salt)/m2/s) - integer :: x2o_Foxx_evap !< Evaporation flux (kg/m2/s) - integer :: x2o_Faxa_prec !< Total precipitation flux (kg/m2/s) - integer :: x2o_Faxa_snow !< Water flux due to snow (kg/m2/s) - integer :: x2o_Faxa_rain !< Water flux due to rain (kg/m2/s) - integer :: x2o_Faxa_bcphidry !< Black Carbon hydrophilic dry deposition - integer :: x2o_Faxa_bcphodry !< Black Carbon hydrophobic dry deposition - integer :: x2o_Faxa_bcphiwet !< Black Carbon hydrophilic wet deposition - integer :: x2o_Faxa_ocphidry !< Organic Carbon hydrophilic dry deposition - integer :: x2o_Faxa_ocphodry !< Organic Carbon hydrophobic dry deposition - integer :: x2o_Faxa_ocphiwet !< Organic Carbon hydrophilic dry deposition - integer :: x2o_Faxa_dstwet1 !< Size 1 dust -- wet deposition - integer :: x2o_Faxa_dstwet2 !< Size 2 dust -- wet deposition - integer :: x2o_Faxa_dstwet3 !< Size 3 dust -- wet deposition - integer :: x2o_Faxa_dstwet4 !< Size 4 dust -- wet deposition - integer :: x2o_Faxa_dstdry1 !< Size 1 dust -- dry deposition - integer :: x2o_Faxa_dstdry2 !< Size 2 dust -- dry deposition - integer :: x2o_Faxa_dstdry3 !< Size 3 dust -- dry deposition - integer :: x2o_Faxa_dstdry4 !< Size 4 dust -- dry deposition - integer :: x2o_Foxx_rofl !< River runoff flux (kg/m2/s) - integer :: x2o_Foxx_rofi !< Ice runoff flux (kg/m2/s) - - ! optional per thickness category fields - - integer, dimension(:), allocatable :: x2o_frac_col !< Fraction of ocean cell, - !! per column - integer, dimension(:), allocatable :: x2o_fracr_col!< Fraction of ocean cell used - !! in radiation computations, - !! per column - integer, dimension(:), allocatable :: x2o_qsw_fracr_col !< qsw * fracr, per column +#include + +! Public member functions +public :: ocn_init_mct +public :: ocn_run_mct +public :: ocn_final_mct +! Flag for debugging +logical, parameter :: debug=.true. + +!> Structure with MCT attribute vectors and indices +type cpl_indices + + ! ocean to coupler + integer :: o2x_So_t !< Surface potential temperature (deg C) + integer :: o2x_So_u !< Surface zonal velocity (m/s) + integer :: o2x_So_v !< Surface meridional velocity (m/s) + integer :: o2x_So_s !< Surface salinity (PSU) + integer :: o2x_So_dhdx !< Zonal slope in the sea surface height + integer :: o2x_So_dhdy !< Meridional lope in the sea surface height + integer :: o2x_So_bldepth !< Boundary layer depth (m) + integer :: o2x_Fioo_q !< Heat flux? + integer :: o2x_Faoo_fco2_ocn!< CO2 flux + integer :: o2x_Faoo_fdms_ocn!< DMS flux + + ! coupler to ocean + integer :: x2o_Si_ifrac !< Fractional ice wrt ocean + integer :: x2o_So_duu10n !< 10m wind speed squared (m^2/s^2) + integer :: x2o_Sa_pslv !< Sea-level pressure (Pa) + integer :: x2o_Sa_co2prog !< Bottom atm level prognostic CO2 + integer :: x2o_Sa_co2diag !< Bottom atm level diagnostic CO2 + integer :: x2o_Sw_lamult !< Wave model langmuir multiplier + integer :: x2o_Sw_ustokes !< Surface Stokes drift, x-component + integer :: x2o_Sw_vstokes !< Surface Stokes drift, y-component + integer :: x2o_Foxx_taux !< Zonal wind stress (W/m2) + integer :: x2o_Foxx_tauy !< Meridonal wind stress (W/m2) + integer :: x2o_Foxx_swnet !< Net short-wave heat flux (W/m2) + integer :: x2o_Foxx_sen !< Sensible heat flux (W/m2) + integer :: x2o_Foxx_lat !< Latent heat flux (W/m2) + integer :: x2o_Foxx_lwup !< Longwave radiation, up (W/m2) + integer :: x2o_Faxa_lwdn !< Longwave radiation, down (W/m2) + integer :: x2o_Fioi_melth !< Heat flux from snow & ice melt (W/m2) + integer :: x2o_Fioi_meltw !< Snow melt flux (kg/m2/s) + integer :: x2o_Fioi_bcpho !< Black Carbon hydrophobic release + !! from sea ice component + integer :: x2o_Fioi_bcphi !< Black Carbon hydrophilic release from + !! sea ice component + integer :: x2o_Fioi_flxdst !< Dust release from sea ice component + integer :: x2o_Fioi_salt !< Salt flux (kg(salt)/m2/s) + integer :: x2o_Foxx_evap !< Evaporation flux (kg/m2/s) + integer :: x2o_Faxa_prec !< Total precipitation flux (kg/m2/s) + integer :: x2o_Faxa_snow !< Water flux due to snow (kg/m2/s) + integer :: x2o_Faxa_rain !< Water flux due to rain (kg/m2/s) + integer :: x2o_Faxa_bcphidry !< Black Carbon hydrophilic dry deposition + integer :: x2o_Faxa_bcphodry !< Black Carbon hydrophobic dry deposition + integer :: x2o_Faxa_bcphiwet !< Black Carbon hydrophilic wet deposition + integer :: x2o_Faxa_ocphidry !< Organic Carbon hydrophilic dry deposition + integer :: x2o_Faxa_ocphodry !< Organic Carbon hydrophobic dry deposition + integer :: x2o_Faxa_ocphiwet !< Organic Carbon hydrophilic dry deposition + integer :: x2o_Faxa_dstwet1 !< Size 1 dust -- wet deposition + integer :: x2o_Faxa_dstwet2 !< Size 2 dust -- wet deposition + integer :: x2o_Faxa_dstwet3 !< Size 3 dust -- wet deposition + integer :: x2o_Faxa_dstwet4 !< Size 4 dust -- wet deposition + integer :: x2o_Faxa_dstdry1 !< Size 1 dust -- dry deposition + integer :: x2o_Faxa_dstdry2 !< Size 2 dust -- dry deposition + integer :: x2o_Faxa_dstdry3 !< Size 3 dust -- dry deposition + integer :: x2o_Faxa_dstdry4 !< Size 4 dust -- dry deposition + integer :: x2o_Foxx_rofl !< River runoff flux (kg/m2/s) + integer :: x2o_Foxx_rofi !< Ice runoff flux (kg/m2/s) + + ! optional per thickness category fields + integer, dimension(:), allocatable :: x2o_frac_col !< Fraction of ocean cell, + !! per column + integer, dimension(:), allocatable :: x2o_fracr_col!< Fraction of ocean cell used + !! in radiation computations, + !! per column + integer, dimension(:), allocatable :: x2o_qsw_fracr_col !< qsw * fracr, per column end type cpl_indices - !> This type is used for communication with other components via the FMS coupler. - ! The element names and types can be changed only with great deliberation, hence - ! the persistnce of things like the cutsy element name "avg_kount". - type, public :: ocean_public_type - type(domain2d) :: Domain ! The domain for the surface fields. - logical :: is_ocean_pe ! .true. on processors that run the ocean model. - character(len=32) :: instance_name = '' ! A name that can be used to identify - ! this instance of an ocean model, for example - ! in ensembles when writing messages. - integer, pointer, dimension(:) :: pelist => NULL() ! The list of ocean PEs. - logical, pointer, dimension(:,:) :: maskmap =>NULL() ! A pointer to an array - ! indicating which logical processors are actually used for - ! the ocean code. The other logical processors would be all - ! land points and are not assigned to actual processors. - ! This need not be assigned if all logical processors are used. - - integer :: stagger = -999 ! The staggering relative to the tracer points - ! points of the two velocity components. Valid entries - ! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, - ! corresponding to the community-standard Arakawa notation. - ! (These are named integers taken from mpp_parameter_mod.) - ! Following MOM, this is BGRID_NE by default when the ocean - ! is initialized, but here it is set to -999 so that a - ! global max across ocean and non-ocean processors can be - ! used to determine its value. - real, pointer, dimension(:,:) :: & - t_surf => NULL(), & ! SST on t-cell (degrees Kelvin) - s_surf => NULL(), & ! SSS on t-cell (psu) - u_surf => NULL(), & ! i-velocity at the locations indicated by stagger, m/s. - v_surf => NULL(), & ! j-velocity at the locations indicated by stagger, m/s. - sea_lev => NULL(), & ! Sea level in m after correction for surface pressure, - ! i.e. dzt(1) + eta_t + patm/rho0/grav (m) - frazil =>NULL(), & ! Accumulated heating (in Joules/m^2) from frazil - ! formation in the ocean. - area => NULL() ! cell area of the ocean surface, in m2. - type(coupler_2d_bc_type) :: fields ! A structure that may contain an - ! array of named tracer-related fields. - integer :: avg_kount ! Used for accumulating averages of this type. - integer, dimension(2) :: axes = 0 ! Axis numbers that are available +!> This type is used for communication with other components via the FMS coupler. +! The element names and types can be changed only with great deliberation, hence +! the persistnce of things like the cutsy element name "avg_kount". +type, public :: ocean_public_type + type(domain2d) :: Domain !< The domain for the surface fields. + logical :: is_ocean_pe !! .true. on processors that run the ocean model. + character(len=32) :: instance_name = '' !< A name that can be used to identify + !! this instance of an ocean model, for example + !! in ensembles when writing messages. + integer, pointer, dimension(:) :: pelist => NULL() !< The list of ocean PEs. + logical, pointer, dimension(:,:) :: maskmap =>NULL() !< A pointer to an array + !! indicating which logical processors are actually + !! used for the ocean code. The other logical + !! processors would be all land points and are not + !! assigned to actual processors. This need not be + !! assigned if all logical processors are used. + + integer :: stagger = -999 !< The staggering relative to the tracer points + !! of the two velocity components. Valid entries + !! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, + !! corresponding to the community-standard Arakawa notation. + !! (These are named integers taken from mpp_parameter_mod.) + !! Following MOM, this is BGRID_NE by default when the ocean + !! is initialized, but here it is set to -999 so that a + !! global max across ocean and non-ocean processors can be + !! used to determine its value. + real, pointer, dimension(:,:) :: & + t_surf => NULL(), & !< SST on t-cell (degrees Kelvin) + s_surf => NULL(), & !< SSS on t-cell (psu) + u_surf => NULL(), & !< i-velocity at the locations indicated by stagger, m/s. + v_surf => NULL(), & !< j-velocity at the locations indicated by stagger, m/s. + sea_lev => NULL(), & !< Sea level in m after correction for surface pressure, + !! i.e. dzt(1) + eta_t + patm/rho0/grav (m) + frazil =>NULL(), & !< Accumulated heating (in Joules/m^2) from frazil + !! formation in the ocean. + area => NULL() !< cell area of the ocean surface, in m2. + type(coupler_2d_bc_type) :: fields !< A structure that may contain an + !! array of named tracer-related fields. + integer :: avg_kount !< Used for accumulating averages of this type. + integer, dimension(2) :: axes = 0 !< Axis numbers that are available ! for I/O using this surface data. - end type ocean_public_type - - !> Contains information about the ocean state, although it is not necessary that - !! this is implemented with all models. This type is private, and can therefore vary - !! between different ocean models. - type, public :: ocean_state_type ; private - logical :: is_ocean_PE = .false. !< True if this is an ocean PE. - type(time_type) :: Time !< The ocean model's time and master clock. - integer :: Restart_control !< An integer that is bit-tested to determine whether - !! incremental restart files are saved and whether they - !! have a time stamped name. +1 (bit 0) for generic - !! files and +2 (bit 1) for time-stamped files. A - !! restart file is saved at the end of a run segment - !! unless Restart_control is negative. - type(time_type) :: energysavedays ! The interval between writing the energies - ! and other integral quantities of the run. - type(time_type) :: write_energy_time ! The next time to write to the energy file. - - integer :: nstep = 0 ! The number of calls to update_ocean. - logical :: use_ice_shelf ! If true, the ice shelf model is enabled. - logical :: icebergs_apply_rigid_boundary ! If true, the icebergs can change ocean bd condition. - real :: kv_iceberg ! The viscosity of the icebergs in m2/s (for ice rigidity) - real :: berg_area_threshold ! Fraction of grid cell which iceberg must occupy - !so that fluxes below are set to zero. (0.5 is a - !good value to use. Not applied for negative values. - real :: latent_heat_fusion ! Latent heat of fusion - real :: density_iceberg ! A typical density of icebergs in kg/m3 (for ice rigidity) - type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() - logical :: restore_salinity ! If true, the coupled MOM driver adds a term to - ! restore salinity to a specified value. - logical :: restore_temp ! If true, the coupled MOM driver adds a term to - ! restore sst to a specified value. - real :: press_to_z ! A conversion factor between pressure and ocean - ! depth in m, usually 1/(rho_0*g), in m Pa-1. - real :: C_p ! The heat capacity of seawater, in J K-1 kg-1. - - type(directories) :: dirs ! A structure containing several relevant directory paths. - type(forcing) :: fluxes ! A structure containing pointers to - ! the ocean forcing fields. - type(forcing) :: flux_tmp ! A secondary structure containing pointers to the - ! ocean forcing fields for when multiple coupled - ! timesteps are taken per thermodynamic step. - type(surface) :: state ! A structure containing pointers to - ! the ocean surface state fields. - type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure - ! containing metrics and related information. - type(verticalGrid_type), pointer :: GV => NULL() ! A pointer to a vertical grid - ! structure containing metrics and related information. - type(MOM_control_struct), pointer :: MOM_CSp => NULL() - type(surface_forcing_CS), pointer :: forcing_CSp => NULL() - type(sum_output_CS), pointer :: sum_output_CSp => NULL() - end type ocean_state_type - - !> Control structure for this module - type MCT_MOM_Data - - type(ocean_state_type), pointer :: ocn_state => NULL() !< The private state of ocean - type(ocean_public_type), pointer :: ocn_public => NULL() !< The public state of ocean - type(ocean_grid_type), pointer :: grid => NULL() !< The grid structure - type(surface), pointer :: ocn_surface => NULL() !< The ocean surface state - type(ice_ocean_boundary_type) :: ice_ocean_boundary !< The ice ocean boundary type - type(seq_infodata_type), pointer :: infodata !< The input info type - type(cpl_indices), public :: ind !< Variable IDs - ! runtime params - logical :: sw_decomp !< Controls whether shortwave is decomposed into four components - real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition - ! i/o - character(len=384) :: pointer_filename !< Name of the ascii file that contains the path - !! and filename of the latest restart file. - integer :: stdout !< standard output unit. (by default, it should point to ocn.log.* file) - end type MCT_MOM_Data - - type(MCT_MOM_Data) :: glb !< global structure +end type ocean_public_type + +!> Contains pointers to the forcing fields which may be used to drive MOM. +!! All fluxes are positive downward. +type, public :: surface_forcing_CS ; private + integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integer values + !! from MOM_domains) to indicate the staggering of + !! the winds that are being provided in calls to + !! update_ocean_model. CIME uses AGRID, so this option + !! is being hard coded for now. + logical :: use_temperature !< If true, temp and saln used as state variables + real :: wind_stress_multiplier!< A multiplier applied to incoming wind stress (nondim). + ! smg: remove when have A=B code reconciled + logical :: bulkmixedlayer !< If true, model based on bulk mixed layer code + real :: Rho0 !< Boussinesq reference density (kg/m^3) + real :: area_surf = -1.0 !< total ocean surface area (m^2) + real :: latent_heat_fusion !< latent heat of fusion (J/kg) + real :: latent_heat_vapor !< latent heat of vaporization (J/kg) + real :: max_p_surf !< maximum surface pressure that can be + !! exerted by the atmosphere and floating sea-ice, + !! in Pa. This is needed because the FMS coupling + !! structure does not limit the water that can be + !! frozen out of the ocean and the ice-ocean heat + !! fluxes are treated explicitly. + logical :: use_limited_P_SSH !< If true, return the sea surface height with + !! the correction for the atmospheric (and sea-ice) + !! pressure limited by max_p_surf instead of the + !! full atmospheric pressure. The default is true. + real :: gust_const !< constant unresolved background gustiness for ustar (Pa) + logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied + !! from an input file. + real, pointer, dimension(:,:) :: & + TKE_tidal => NULL(), & !< turbulent kinetic energy introduced to the + !! bottom boundary layer by drag on the tidal flows, + !! in W m-2. + gust => NULL(), & !< spatially varying unresolved background + !! gustiness that contributes to ustar (Pa). + !! gust is used when read_gust_2d is true. + ustar_tidal => NULL() !< tidal contribution to the bottom friction velocity (m/s) + real :: cd_tides !< drag coefficient that applies to the tides (nondimensional) + real :: utide !< constant tidal velocity to use if read_tideamp + !! is false, in m s-1. + logical :: read_tideamp !< If true, spatially varying tidal amplitude read from a file. + logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts + !! to damp surface deflections (especially surface + !! gravity waves). The default is false. + real :: Kv_sea_ice !< viscosity in sea-ice that resists sheared vertical motions (m^2/s) + real :: density_sea_ice !< typical density of sea-ice (kg/m^3). The value is + !! only used to convert the ice pressure into + !! appropriate units for use with Kv_sea_ice. + real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which + !! sea-ice viscosity becomes effective, in kg m-2, + !! typically of order 1000 kg m-2. + logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments + real :: Flux_const !< piston velocity for surface restoring (m/s) + logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux + logical :: adjust_net_srestore_to_zero !< adjust srestore to zero (for both salt_flux or vprec) + logical :: adjust_net_srestore_by_scaling !< adjust srestore w/o moving zero contour + logical :: adjust_net_fresh_water_to_zero !< adjust net surface fresh-water (w/ restoring) to zero + logical :: adjust_net_fresh_water_by_scaling !< adjust net surface fresh-water w/o moving zero contour + logical :: mask_srestore_under_ice !< If true, use an ice mask defined by frazil + !! criteria for salinity restoring. + real :: ice_salt_concentration !< salt concentration for sea ice (kg/kg) + logical :: mask_srestore_marginal_seas !< if true, then mask SSS restoring in marginal seas + real :: max_delta_srestore !< maximum delta salinity used for restoring + real :: max_delta_trestore !< maximum delta sst used for restoring + real, pointer, dimension(:,:) :: basin_mask => NULL() !< mask for SSS restoring + type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing + character(len=200) :: inputdir !< directory where NetCDF input files are + character(len=200) :: salt_restore_file !< filename for salt restoring data + character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file + character(len=200) :: temp_restore_file !< filename for sst restoring data + character(len=30) :: temp_restore_var_name !< name of surface temperature in temp_restore_file + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. + type(forcing_diags), public :: handles !< diagnostics handles + !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL() + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< restart pointer + type(user_revise_forcing_CS), pointer :: urf_CS => NULL()!< user revise pointer +end type surface_forcing_CS + +!> Contains information about the ocean state, although it is not necessary that +!! this is implemented with all models. This type is private, and can therefore vary +!! between different ocean models. +type, public :: ocean_state_type ; private + logical :: is_ocean_PE = .false. !< True if this is an ocean PE. + type(time_type) :: Time !< The ocean model's time and master clock. + integer :: Restart_control !< An integer that is bit-tested to determine whether + !! incremental restart files are saved and whether they + !! have a time stamped name. +1 (bit 0) for generic + !! files and +2 (bit 1) for time-stamped files. A + !! restart file is saved at the end of a run segment + !! unless Restart_control is negative. + type(time_type) :: energysavedays !< The interval between writing the energies + !! and other integral quantities of the run. + type(time_type) :: write_energy_time !< The next time to write to the energy file. + integer :: nstep = 0 !< The number of calls to update_ocean. + logical :: use_ice_shelf !< If true, the ice shelf model is enabled. + logical :: icebergs_apply_rigid_boundary !< If true, the icebergs can change ocean bd condition. + real :: kv_iceberg !< The viscosity of the icebergs in m2/s (for ice rigidity) + real :: berg_area_threshold !< Fraction of grid cell which iceberg must occupy + !! so that fluxes below are set to zero. (0.5 is a + !! good value to use. Not applied for negative values. + real :: latent_heat_fusion !< Latent heat of fusion + real :: density_iceberg !< A typical density of icebergs in kg/m3 (for ice rigidity) + type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() !< ice shelf structure. + logical :: restore_salinity !< If true, the coupled MOM driver adds a term to + !! restore salinity to a specified value. + logical :: restore_temp !< If true, the coupled MOM driver adds a term to + !! restore sst to a specified value. + real :: press_to_z !< A conversion factor between pressure and ocean + !! depth in m, usually 1/(rho_0*g), in m Pa-1. + real :: C_p !< The heat capacity of seawater, in J K-1 kg-1. + type(directories) :: dirs !< A structure containing several relevant directory paths. + type(forcing) :: fluxes !< A structure containing pointers to + !! the ocean forcing fields. + type(forcing) :: flux_tmp !< A secondary structure containing pointers to the + !! ocean forcing fields for when multiple coupled + !! timesteps are taken per thermodynamic step. + type(surface) :: state !< A structure containing pointers to + !! the ocean surface state fields. + type(ocean_grid_type), pointer :: grid => NULL() !< A pointer to a grid structure + !! containing metrics and related information. + type(verticalGrid_type), pointer :: GV => NULL() !< A pointer to a vertical grid + !! structure containing metrics and related information. + type(MOM_control_struct), pointer :: MOM_CSp => NULL() + type(surface_forcing_CS), pointer :: forcing_CSp => NULL() + type(sum_output_CS), pointer :: sum_output_CSp => NULL() +end type ocean_state_type + +!> Control structure for this module +type MCT_MOM_Data + + type(ocean_state_type), pointer :: ocn_state => NULL() !< The private state of ocean + type(ocean_public_type), pointer :: ocn_public => NULL() !< The public state of ocean + type(ocean_grid_type), pointer :: grid => NULL() !< The grid structure + type(surface), pointer :: ocn_surface => NULL() !< The ocean surface state + type(forcing) :: fluxes !< Structure that contains pointers to the + !! boundary forcing used to drive the liquid + !! ocean simulated by MOM. + type(seq_infodata_type), pointer :: infodata !< The input info type + type(cpl_indices), public :: ind !< Variable IDs + ! runtime params + logical :: sw_decomp !< Controls whether shortwave is decomposed into four components + real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + ! i/o + character(len=384) :: pointer_filename !< Name of the ascii file that contains the path + !! and filename of the latest restart file. + integer :: stdout !< standard output unit. (by default, it should point to ocn.log.* file) +end type MCT_MOM_Data + +type(MCT_MOM_Data) :: glb !< global structure +integer :: id_clock_forcing contains -!> Initializes MOM6 +!> This subroutine initializes MOM6. subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) type(ESMF_Clock), intent(inout) :: EClock !< Time and time step ? \todo Why must this !! be intent(inout)? @@ -317,13 +411,11 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) seconds_in_day = 86400.0d0, & minutes_in_hour = 60.0d0 - ! ocn model input namelist filename - character(len=99) :: ocn_modelio_name - integer :: shrlogunit ! original log file unit - integer :: shrloglev ! original log level + character(len=99) :: ocn_modelio_name !< ocn model input namelist filename + integer :: shrlogunit !< original log file unit + integer :: shrloglev !< original log level - ! instance control vars (these are local for now) - integer(kind=4) :: inst_index + integer(kind=4) :: inst_index !< instance control vars (these are local for now) character(len=16) :: inst_name character(len=16) :: inst_suffix @@ -560,32 +652,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! Size of global domain call get_global_grid_size(glb%grid, ni, nj) - ! Allocate ice_ocean_boundary using global indexing without halos - isc = glb%grid%isc + glb%grid%idg_offset - iec = glb%grid%iec + glb%grid%idg_offset - jsc = glb%grid%jsc + glb%grid%jdg_offset - jec = glb%grid%jec + glb%grid%jdg_offset - allocate(glb%ice_ocean_boundary%u_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%u_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%v_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%v_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%t_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%t_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%q_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%q_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%latent(isc:iec,jsc:jec)); glb%ice_ocean_boundary%latent(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%salt_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%salt_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%lw_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%lw_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_vis_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dir(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_vis_dif(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dif(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_nir_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_nir_dir(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_nir_dif(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_nir_dif(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%lprec(isc:iec,jsc:jec)); glb%ice_ocean_boundary%lprec(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%fprec(isc:iec,jsc:jec)); glb%ice_ocean_boundary%fprec(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%runoff(isc:iec,jsc:jec)); glb%ice_ocean_boundary%runoff(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%calving(isc:iec,jsc:jec)); glb%ice_ocean_boundary%calving(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%runoff_hflx(isc:iec,jsc:jec)); glb%ice_ocean_boundary%runoff_hflx(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%calving_hflx(isc:iec,jsc:jec)); glb%ice_ocean_boundary%calving_hflx(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%p(isc:iec,jsc:jec)); glb%ice_ocean_boundary%p(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%mi(isc:iec,jsc:jec)); glb%ice_ocean_boundary%mi(:,:) = 0.0 - - if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" call seq_infodata_PutData( glb%infodata, & @@ -593,7 +659,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call seq_infodata_PutData( glb%infodata, & ocn_prognostic=.true., ocnrof_prognostic=.true.) - if (debug .and. root_pe().eq.pe_here()) print *, "leaving ocean_init_mct" ! Reset shr logging to original values @@ -733,23 +798,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! internal variables in the ice model. character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read -! This subroutine initializes both the ocean state and the ocean surface type. +! This subroutine initializes both the ocean state and the ocean surface type. ! Because of the way that indicies and domains are handled, Ocean_sfc must have ! been used in a previous call to initialize_ocean_type. -! Arguments: Ocean_sfc - A structure containing various publicly visible ocean -! surface properties after initialization, this is intent(out). -! (out,private) OS - A structure whose internal contents are private -! to ocean_model_mod that may be used to contain all -! information about the ocean's interior state. -! (in) Time_init - The start time for the coupled model's calendar. -! (in) Time_in - The time at which to initialize the ocean model. - real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS. - real :: Rho0 ! The Boussinesq ocean density, in kg m-3. - real :: G_Earth ! The gravitational acceleration in m s-2. -! This include declares and sets the variable "version". + real :: Time_unit !< The time unit in seconds for ENERGYSAVEDAYS. + real :: Rho0 !< The Boussinesq ocean density, in kg m-3. + real :: G_Earth !< The gravitational acceleration in m s-2. + !! This include declares and sets the variable "version". #include "version_variable.h" - character(len=40) :: mdl = "ocean_model_init" ! This module's name. + character(len=40) :: mdl = "ocean_model_init" !< This module's name. character(len=48) :: stagger integer :: secs, days type(param_file_type) :: param_file !< A structure to parse for run-time parameters @@ -897,7 +955,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call callTree_leave("ocean_model_init(") end subroutine ocean_model_init -!> This subroutine extracts the surface properties from the ocean's internal +!> Extracts the surface properties from the ocean's internal !! state and stores them in the ocean type returned to the calling ice model. !! It has to be separate from the ocean_initialization call because the coupler !! module allocates the space for some of these variables. @@ -920,21 +978,334 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) end subroutine ocean_model_init_sfc +!> Initializes surface forcing: get relevant parameters and allocate arrays. +subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, restore_temp) + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output + type(surface_forcing_CS), pointer :: CS !< A pointer that is set to point to the + !! control structure for this module + logical, optional, intent(in) :: restore_salt, restore_temp !< If present and true, + !! temp/salt restoring will be applied + + ! local variables + real :: utide !< The RMS tidal velocity, in m s-1. + type(directories) :: dirs + logical :: new_sim, iceberg_flux_diags + type(time_type) :: Time_frc + character(len=200) :: TideAmp_file, gust_file, salt_file, temp_file ! Input file names. +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "ocn_comp_mct" ! This module's name. + character(len=48) :: stagger + character(len=240) :: basin_file + integer :: i, j, isd, ied, jsd, jed + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + if (associated(CS)) then + call MOM_error(WARNING, "surface_forcing_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + id_clock_forcing=cpu_clock_id('Ocean surface forcing', grain=CLOCK_SUBCOMPONENT) + call cpu_clock_begin(id_clock_forcing) + + CS%diag => diag + + call write_version_number (version) + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + + call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, & + "The directory in which all input files are found.", & + default=".") + CS%inputdir = slasher(CS%inputdir) + call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & + "If true, Temperature and salinity are used as state \n"//& + "variables.", default=.true.) + call get_param(param_file, mdl, "RHO_0", CS%Rho0, & + "The mean ocean density used with BOUSSINESQ true to \n"//& + "calculate accelerations and the mass for conservation \n"//& + "properties, or with BOUSSINSEQ false to convert some \n"//& + "parameters from vertical units of m to kg m-2.", & + units="kg m-3", default=1035.0) + call get_param(param_file, mdl, "LATENT_HEAT_FUSION", CS%latent_heat_fusion, & + "The latent heat of fusion.", units="J/kg", default=hlf) + call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", CS%latent_heat_vapor, & + "The latent heat of fusion.", units="J/kg", default=hlv) + call get_param(param_file, mdl, "MAX_P_SURF", CS%max_p_surf, & + "The maximum surface pressure that can be exerted by the \n"//& + "atmosphere and floating sea-ice or ice shelves. This is \n"//& + "needed because the FMS coupling structure does not \n"//& + "limit the water that can be frozen out of the ocean and \n"//& + "the ice-ocean heat fluxes are treated explicitly. No \n"//& + "limit is applied if a negative value is used.", units="Pa", & + default=-1.0) + call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", & + CS%adjust_net_srestore_to_zero, & + "If true, adjusts the salinity restoring seen to zero\n"//& + "whether restoring is via a salt flux or virtual precip.",& + default=restore_salt) + call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_BY_SCALING", & + CS%adjust_net_srestore_by_scaling, & + "If true, adjustments to salt restoring to achieve zero net are\n"//& + "made by scaling values without moving the zero contour.",& + default=.false.) + call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_TO_ZERO", & + CS%adjust_net_fresh_water_to_zero, & + "If true, adjusts the net fresh-water forcing seen \n"//& + "by the ocean (including restoring) to zero.", default=.false.) + call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", & + CS%adjust_net_fresh_water_by_scaling, & + "If true, adjustments to net fresh water to achieve zero net are\n"//& + "made by scaling values without moving the zero contour.",& + default=.false.) + call get_param(param_file, mdl, "ICE_SALT_CONCENTRATION", & + CS%ice_salt_concentration, & + "The assumed sea-ice salinity needed to reverse engineer the \n"//& + "melt flux (or ice-ocean fresh-water flux).", & + units="kg/kg", default=0.005) + call get_param(param_file, mdl, "USE_LIMITED_PATM_SSH", CS%use_limited_P_SSH, & + "If true, return the sea surface height with the \n"//& + "correction for the atmospheric (and sea-ice) pressure \n"//& + "limited by max_p_surf instead of the full atmospheric \n"//& + "pressure.", default=.true.) + +! smg: should get_param call should be removed when have A=B code reconciled. +! this param is used to distinguish how to diagnose surface heat content from water. + call get_param(param_file, mdl, "BULKMIXEDLAYER", CS%bulkmixedlayer, & + default=CS%use_temperature,do_not_log=.true.) + + call get_param(param_file, mdl, "WIND_STAGGER", stagger, & + "A case-insensitive character string to indicate the \n"//& + "staggering of the input wind stress field. Valid \n"//& + "values are 'A', 'B', or 'C'.", default="C") + if (uppercase(stagger(1:1)) == 'A') then ; CS%wind_stagger = AGRID + elseif (uppercase(stagger(1:1)) == 'B') then ; CS%wind_stagger = BGRID_NE + elseif (uppercase(stagger(1:1)) == 'C') then ; CS%wind_stagger = CGRID_NE + else ; call MOM_error(FATAL,"surface_forcing_init: WIND_STAGGER = "// & + trim(stagger)//" is invalid.") ; endif + call get_param(param_file, mdl, "WIND_STRESS_MULTIPLIER", CS%wind_stress_multiplier, & + "A factor multiplying the wind-stress given to the ocean by the\n"//& + "coupler. This is used for testing and should be =1.0 for any\n"//& + "production runs.", default=1.0) + + if (restore_salt) then + call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & + "The constant that relates the restoring surface fluxes \n"//& + "to the relative surface anomalies (akin to a piston \n"//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "SALT_RESTORE_FILE", CS%salt_restore_file, & + "A file in which to find the surface salinity to use for restoring.", & + default="salt_restore.nc") + call get_param(param_file, mdl, "SALT_RESTORE_VARIABLE", CS%salt_restore_var_name, & + "The name of the surface salinity variable to read from "//& + "SALT_RESTORE_FILE for restoring salinity.", & + default="salt") +! Convert CS%Flux_const from m day-1 to m s-1. + CS%Flux_const = CS%Flux_const / 86400.0 + + call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & + "If true, the restoring of salinity is applied as a salt \n"//& + "flux instead of as a freshwater flux.", default=.false.) + call get_param(param_file, mdl, "MAX_DELTA_SRESTORE", CS%max_delta_srestore, & + "The maximum salinity difference used in restoring terms.", & + units="PSU or g kg-1", default=999.0) + call get_param(param_file, mdl, "MASK_SRESTORE_UNDER_ICE", & + CS%mask_srestore_under_ice, & + "If true, disables SSS restoring under sea-ice based on a frazil\n"//& + "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", & + default=.false.) + call get_param(param_file, mdl, "MASK_SRESTORE_MARGINAL_SEAS", & + CS%mask_srestore_marginal_seas, & + "If true, disable SSS restoring in marginal seas. Only used when\n"//& + "RESTORE_SALINITY is True.", default=.false.) + call get_param(param_file, mdl, "BASIN_FILE", basin_file, & + "A file in which to find the basin masks, in variable 'basin'.", & + default="basin.nc") + basin_file = trim(CS%inputdir) // trim(basin_file) + call safe_alloc_ptr(CS%basin_mask,isd,ied,jsd,jed) ; CS%basin_mask(:,:) = 1.0 + if (CS%mask_srestore_marginal_seas) then + call read_data(basin_file,'basin',CS%basin_mask,domain=G%domain%mpp_domain,timelevel=1) + do j=jsd,jed ; do i=isd,ied + if (CS%basin_mask(i,j) >= 6.0) then ; CS%basin_mask(i,j) = 0.0 + else ; CS%basin_mask(i,j) = 1.0 ; endif + enddo ; enddo + endif + endif + + if (restore_temp) then + call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & + "The constant that relates the restoring surface fluxes \n"//& + "to the relative surface anomalies (akin to a piston \n"//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "SST_RESTORE_FILE", CS%temp_restore_file, & + "A file in which to find the surface temperature to use for restoring.", & + default="temp_restore.nc") + call get_param(param_file, mdl, "SST_RESTORE_VARIABLE", CS%temp_restore_var_name, & + "The name of the surface temperature variable to read from "//& + "SST_RESTORE_FILE for restoring sst.", & + default="temp") +! Convert CS%Flux_const from m day-1 to m s-1. + CS%Flux_const = CS%Flux_const / 86400.0 + + call get_param(param_file, mdl, "MAX_DELTA_TRESTORE", CS%max_delta_trestore, & + "The maximum sst difference used in restoring terms.", & + units="degC ", default=999.0) + + endif + +! Optionally read tidal amplitude from input file (m s-1) on model grid. +! Otherwise use default tidal amplitude for bottom frictionally-generated +! dissipation. Default cd_tides is chosen to yield approx 1 TWatt of +! work done against tides globally using OSU tidal amplitude. + call get_param(param_file, mdl, "CD_TIDES", CS%cd_tides, & + "The drag coefficient that applies to the tides.", & + units="nondim", default=1.0e-4) + call get_param(param_file, mdl, "READ_TIDEAMP", CS%read_TIDEAMP, & + "If true, read a file (given by TIDEAMP_FILE) containing \n"//& + "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + if (CS%read_TIDEAMP) then + call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & + "The path to the file containing the spatially varying \n"//& + "tidal amplitudes with INT_TIDE_DISSIPATION.", & + default="tideamp.nc") + CS%utide=0.0 + else + call get_param(param_file, mdl, "UTIDE", CS%utide, & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + units="m s-1", default=0.0) + endif + + call safe_alloc_ptr(CS%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%ustar_tidal,isd,ied,jsd,jed) + + if (CS%read_TIDEAMP) then + TideAmp_file = trim(CS%inputdir) // trim(TideAmp_file) + call read_data(TideAmp_file,'tideamp',CS%TKE_tidal,domain=G%domain%mpp_domain,timelevel=1) + do j=jsd, jed; do i=isd, ied + utide = CS%TKE_tidal(i,j) + CS%TKE_tidal(i,j) = G%mask2dT(i,j)*CS%Rho0*CS%cd_tides*(utide*utide*utide) + CS%ustar_tidal(i,j)=sqrt(CS%cd_tides)*utide + enddo ; enddo + else + do j=jsd,jed; do i=isd,ied + utide=CS%utide + CS%TKE_tidal(i,j) = CS%Rho0*CS%cd_tides*(utide*utide*utide) + CS%ustar_tidal(i,j)=sqrt(CS%cd_tides)*utide + enddo ; enddo + endif + + call time_interp_external_init + +! Optionally read a x-y gustiness field in place of a global +! constant. + + call get_param(param_file, mdl, "READ_GUST_2D", CS%read_gust_2d, & + "If true, use a 2-dimensional gustiness supplied from \n"//& + "an input file", default=.false.) + call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & + "The background gustiness in the winds.", units="Pa", & + default=0.02) + if (CS%read_gust_2d) then + call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & + "The file in which the wind gustiness is found in \n"//& + "variable gustiness.") + + call safe_alloc_ptr(CS%gust,isd,ied,jsd,jed) + gust_file = trim(CS%inputdir) // trim(gust_file) + call read_data(gust_file,'gustiness',CS%gust,domain=G%domain%mpp_domain, & + timelevel=1) ! units should be Pa + endif + +! See whether sufficiently thick sea ice should be treated as rigid. + call get_param(param_file, mdl, "USE_RIGID_SEA_ICE", CS%rigid_sea_ice, & + "If true, sea-ice is rigid enough to exert a \n"//& + "nonhydrostatic pressure that resist vertical motion.", & + default=.false.) + if (CS%rigid_sea_ice) then + call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", CS%density_sea_ice, & + "A typical density of sea ice, used with the kinematic \n"//& + "viscosity, when USE_RIGID_SEA_ICE is true.", units="kg m-3", & + default=900.0) + call get_param(param_file, mdl, "SEA_ICE_VISCOSITY", CS%Kv_sea_ice, & + "The kinematic viscosity of sufficiently thick sea ice \n"//& + "for use in calculating the rigidity of sea ice.", & + units="m2 s-1", default=1.0e9) + call get_param(param_file, mdl, "SEA_ICE_RIGID_MASS", CS%rigid_sea_ice_mass, & + "The mass of sea-ice per unit area at which the sea-ice \n"//& + "starts to exhibit rigidity", units="kg m-2", default=1000.0) + endif + + call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, & + "If true, makes available diagnostics of fluxes from icebergs\n"//& + "as seen by MOM6.", default=.false.) + call register_forcing_type_diags(Time, diag, CS%use_temperature, CS%handles, & + use_berg_fluxes=iceberg_flux_diags) + + call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & + "If true, allows flux adjustments to specified via the \n"//& + "data_table using the component name 'OCN'.", default=.false.) + if (CS%allow_flux_adjustments) then + call data_override_init(Ocean_domain_in=G%Domain%mpp_domain) + endif + + if (present(restore_salt)) then ; if (restore_salt) then + salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) + CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + endif ; endif + + if (present(restore_temp)) then ; if (restore_temp) then + temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) + CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + endif ; endif + + ! Set up any restart fields associated with the forcing. + call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") +!### call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, & +!### CS%restart_CSp) + call restart_init_end(CS%restart_CSp) + + if (associated(CS%restart_CSp)) then + call Get_MOM_Input(dirs=dirs) + + new_sim = .false. + if ((dirs%input_filename(1:1) == 'n') .and. & + (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. + if (.not.new_sim) then + call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & + G, CS%restart_CSp) + endif + endif + +!### call controlled_forcing_init(Time, G, param_file, diag, CS%ctrl_forcing_CSp) + + call user_revise_forcing_init(param_file, CS%urf_CS) + + call cpu_clock_end(id_clock_forcing) +end subroutine surface_forcing_init + +!> Initializes domain and state variables contained in the ocean public type. subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, & gas_fields_ocn) - type(domain2D), intent(in) :: input_domain - type(ocean_public_type), intent(inout) :: Ocean_sfc - type(diag_ctrl), intent(in) :: diag - logical, intent(in), optional :: maskmap(:,:) + type(domain2D), intent(in) :: input_domain !< The FMS domain for the input structure + type(ocean_public_type), intent(inout) :: Ocean_sfc !< Ocean surface state + type(diag_ctrl), intent(in) :: diag !< A structure used to control diagnostics. + logical, intent(in), optional :: maskmap(:,:) !< A pointer to an array indicating which + !! logical processors are actually used for the ocean code. type(coupler_1d_bc_type), & optional, intent(in) :: gas_fields_ocn !< If present, this type describes the - !! ocean and surface-ice fields that will participate - !! in the calculation of additional gas or other - !! tracer fluxes. - + !! ocean and surface-ice fields that will participate + !! in the calculation of additional gas or other + !! tracer fluxes. + ! local variables integer :: xsz, ysz, layout(2) - ! ice-ocean-boundary fields are always allocated using absolute indicies - ! and have no halos. integer :: isc, iec, jsc, jec call mpp_get_layout(input_domain,layout) @@ -972,17 +1343,17 @@ end subroutine initialize_ocean_public_type !> Translates the coupler's ocean_data_type into MOM6's surface state variable. !! This may eventually be folded into the MOM6's code that calculates the -!! surface state in the first place. Note the offset in the arrays because -!! the ocean_data_type has no halo points in its arrays and always uses -!! absolute indicies. +!! surface state in the first place. subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, & patm, press_to_z) type(surface), intent(inout) :: state - type(ocean_public_type), target, intent(inout) :: Ocean_sfc - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - logical, intent(in) :: use_conT_absS - real, optional, intent(in) :: patm(:,:) - real, optional, intent(in) :: press_to_z + type(ocean_public_type), target, intent(inout) :: Ocean_sfc !< Ocean surface state + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + logical, intent(in) :: use_conT_absS !< If true, , the prognostics + !! T&S are the conservative temperature + real, optional, intent(in) :: patm(:,:) !< Atmospheric pressure. + real, optional, intent(in) :: press_to_z !< Factor to tranform atmospheric + !! pressure to z? ! local variables real :: IgR0 @@ -1069,7 +1440,9 @@ subroutine get_state_pointers(OS, grid, surf) end subroutine get_state_pointers -!> Maps outgoing ocean data to MCT buffer +!> Maps outgoing ocean data to MCT buffer. +!! See \ref section_ocn_export for a summary of the data +!! that is transferred from MOM6 to MCT. subroutine ocn_export(ind, ocn_public, grid, o2x) type(cpl_indices), intent(inout) :: ind !< Structure with coupler !! indices and vectors @@ -1161,123 +1534,6 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) end subroutine ocn_export -!> Fills the ice ocean boundary type -subroutine fill_ice_ocean_bnd(ice_ocean_boundary, grid, x2o_o, ind, & - sw_decomp, c1, c2, c3, c4) - - type(ice_ocean_boundary_type), intent(inout) :: ice_ocean_boundary !< A type for the ice ocean boundary - type(ocean_grid_type), intent(in) :: grid !< - !type(mct_aVect), intent(in) :: x2o_o - real(kind=8), intent(in) :: x2o_o(:,:) - type(cpl_indices), intent(inout) :: ind - logical, intent(in) :: sw_decomp !< controls if shortwave is - !!decomposed into four components - real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition - - ! local variables - integer :: i, j, k, ig, jg !< grid indices - - ! /TODO The following comments summarizes the mismatches between MCT and MOM6 in terms - ! of ice ocean fluxes. - - ! Redundancies: - ! x2o_Foxx_lat, x2o_Faxa_prec, x2o_Faxa_evap - **latent from MOM6 and coupler are different** - ! x2o_Faxa_prec = x2o_Faxa_rain + x2o_Faxa_snow - - ! Variables that **could not** be verified so far: - ! x2o_Foxx_rofl - ! x2o_Foxx_rof - - ! Variables in IOB that are **NOT** filled by the coupler: - ! ustar_berg, frictional velocity beneath icebergs (m/s) - ! area_berg, area covered by icebergs(m2/m2) - ! mass_berg, mass of icebergs(kg/m2) - ! runoff_hflx, heat content of liquid runoff (W/m2) - ! calving_hflx, heat content of frozen runoff (W/m2) - ! mi, mass of ice (kg/m2) - - ! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in IOB): - ! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2) - ! x2o_Fioi_meltw, snow melt flux (kg/m2/s) - ! x2o_Si_ifrac, fractional ice wrt ocean - ! x2o_So_duu10n, 10m wind speed squared (m^2/s^2) - ! x2o_Sa_co2prog, bottom atm level prognostic CO2 - ! x2o_Sa_co2diag, bottom atm level diagnostic CO2 - ! x2o_Foxx_lat, latent heat flux (W/m2) - - ! Langmuir related fields: - ! surface Stokes drift, x-comp. (x2o_Sw_ustokes) - ! surface Stokes drift, y-comp. (x2o_Sw_vstokes) - ! wave model langmuir multiplier (x2o_Sw_lamult) - - ! Biogeochemistry: - ! x2o_Fioi_bcpho, Black Carbon hydrophobic release from sea ice component - ! x2o_Fioi_bcphi, Black Carbon hydrophilic release from sea ice component - ! x2o_Fioi_flxdst, Dust release from sea ice component - ! x2o_Faxa_bcphidry, Black Carbon hydrophilic dry deposition - ! x2o_Faxa_bcphodry, Black Carbon hydrophobic dry deposition - ! x2o_Faxa_bcphiwet, Black Carbon hydrophobic wet deposition - ! x2o_Faxa_ocphidry, Organic Carbon hydrophilic dry deposition - ! x2o_Faxa_ocphodry, Organic Carbon hydrophobic dry deposition - ! x2o_Faxa_ocphiwet, Organic Carbon hydrophilic dry deposition - ! x2o_Faxa_dstwet, Sizes 1 to 4 dust - wet deposition - ! x2o_Faxa_dstdry, Sizes 1 to 4 dust - dry depositi - - k = 0 - do j = grid%jsc, grid%jec - jg = j + grid%jdg_offset - do i = grid%isc, grid%iec - k = k + 1 ! Increment position within gindex - ig = i + grid%idg_offset - ! sea-level pressure (Pa) - ice_ocean_boundary%p(ig,jg) = x2o_o(ind%x2o_Sa_pslv,k) - ! zonal wind stress (taux) - ice_ocean_boundary%u_flux(ig,jg) = x2o_o(ind%x2o_Foxx_taux,k) - ! meridional wind stress (tauy) - ice_ocean_boundary%v_flux(ig,jg) = x2o_o(ind%x2o_Foxx_tauy,k) - ! sensible heat flux (W/m2) - ice_ocean_boundary%t_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_sen,k) - ! latent heat flux (W/m^2) - ice_ocean_boundary%latent(ig,jg) = x2o_o(ind%x2o_Foxx_lat,k) - ! salt flux - ice_ocean_boundary%salt_flux(ig,jg) = -x2o_o(ind%x2o_Fioi_salt,k) - ! heat content from frozen runoff - ice_ocean_boundary%calving_hflx(ig,jg) = 0.0 - ! snow melt flux - !ice_ocean_boundary%fprec(ig,jg) = x2o_o(ind%x2o_Fioi_meltw,k) - ! river runoff flux - ice_ocean_boundary%runoff(ig,jg) = x2o_o(ind%x2o_Foxx_rofl,k) - ! ice runoff flux - ice_ocean_boundary%calving(ig,jg) = x2o_o(ind%x2o_Foxx_rofi,k) - ! liquid precipitation (rain) - ice_ocean_boundary%lprec(ig,jg) = x2o_o(ind%x2o_Faxa_rain,k) - ! frozen precipitation (snow) - ice_ocean_boundary%fprec(ig,jg) = x2o_o(ind%x2o_Faxa_snow,k) - ! evaporation flux, MOM6 calls q_flux specific humidity (kg/m2/s) - ice_ocean_boundary%q_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_evap,k) - ! longwave radiation, sum up and down (W/m2) - ice_ocean_boundary%lw_flux(ig,jg) = x2o_o(ind%x2o_Faxa_lwdn,k) + x2o_o(ind%x2o_Foxx_lwup,k) - if (sw_decomp) then - ! Use runtime coefficients to decompose net short-wave heat flux into 4 components - ! 1) visible, direct shortwave (W/m2) - ice_ocean_boundary%sw_flux_vis_dir(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c1 - ! 2) visible, diffuse shortwave (W/m2) - ice_ocean_boundary%sw_flux_vis_dif(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c2 - ! 3) near-IR, direct shortwave (W/m2) - ice_ocean_boundary%sw_flux_nir_dir(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c3 - ! 4) near-IR, diffuse shortwave (W/m2) - ice_ocean_boundary%sw_flux_nir_dif(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c4 - else - call MOM_error(FATAL,"fill_ice_ocean_bnd: this option has not been implemented yet."// & - "Shortwave must be decomposed using coeffs. c1, c2, c3, c4."); - endif - enddo - enddo - - ice_ocean_boundary%wind_stagger = AGRID - -end subroutine fill_ice_ocean_bnd - !> Step forward ocean model for coupling interval subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(ESMF_Clock), intent(inout) :: EClock !< Time and time step ? \todo Why must this be intent(inout)? @@ -1340,14 +1596,11 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) ! set (actually, get from mct) the cdata pointers: ! \todo this was done in _init_, is it needed again. Does this infodata need to be in glb%? + ! GMM, check if this is needed! call seq_cdata_setptrs(cdata_o, infodata=glb%infodata) - ! fill ice ocean boundary - call fill_ice_ocean_bnd(glb%ice_ocean_boundary, glb%grid, x2o_o%rattr, glb%ind, glb%sw_decomp, & - glb%c1, glb%c2, glb%c3, glb%c4) - - call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, & - time_start, coupling_timestep) + call update_ocean_model(glb%ocn_state, glb%ocn_public, time_start, coupling_timestep, & + x2o_o%rattr, glb%ind, glb%sw_decomp, glb%c1, glb%c2, glb%c3, glb%c4) ! return export state to driver call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr) @@ -1398,40 +1651,55 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) end subroutine ocn_run_mct +!> Saves restart fields associated with the forcing +subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & + filename_suffix) + type(surface_forcing_CS), pointer :: CS !< pointer to the control structure + !! returned by a previous call to + !! surface_forcing_init + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(time_type), intent(in) :: Time !< model time at this call + character(len=*), intent(in) :: directory !< optional directory into which + !! to write these restart files + logical, optional, intent(in) :: time_stamped !< If true, the restart file + !! names include a unique time + !! stamp + character(len=*), optional, intent(in) :: filename_suffix !< optional suffix + !! (e.g., a time-stamp) to append to the + !! restart file names + if (.not.associated(CS)) return + if (.not.associated(CS%restart_CSp)) return + call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + +end subroutine forcing_save_restart + !> Updates the ocean model fields. This code wraps the call to step_MOM with MOM6's call. -!! It uses the forcing in Ice_ocean_boundary to advance the ocean model's state from the +!! It uses the forcing to advance the ocean model's state from the !! input value of Ocean_state (which must be for time time_start_update) for a time interval -!! of Ocean_coupling_time_step, eturning the publicly visible ocean surface properties in +!! of Ocean_coupling_time_step, returning the publicly visible ocean surface properties in !! Ocean_sfc and storing the new ocean properties in Ocean_state. -subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & - time_start_update, Ocean_coupling_time_step) - type(ice_ocean_boundary_type), intent(inout) :: Ice_ocean_boundary - type(ocean_state_type), pointer :: OS - type(ocean_public_type), intent(inout) :: Ocean_sfc - type(time_type), intent(in) :: time_start_update - type(time_type), intent(in) :: Ocean_coupling_time_step - -! Arguments: Ice_ocean_boundary - A structure containing the various forcing -! fields coming from the ice. It is intent in. -! (inout) Ocean_state - A structure containing the internal ocean state. -! (out) Ocean_sfc - A structure containing all the publicly visible ocean -! surface fields after a coupling time step. -! (in) time_start_update - The time at the beginning of the update step. -! (in) Ocean_coupling_time_step - The amount of time over which to advance -! the ocean. - -! Note: although several types are declared intent(inout), this is to allow for -! the possibility of halo updates and to keep previously allocated memory. -! In practice, Ice_ocean_boundary is intent in, Ocean_state is private to -! this module and intent inout, and Ocean_sfc is intent out. - type(time_type) :: Master_time ! This allows step_MOM to temporarily change - ! the time that is seen by internal modules. - type(time_type) :: Time1 ! The value of the ocean model's time at the - ! start of a call to step_MOM. - integer :: index_bnds(4) ! The computational domain index bounds in the - ! ice-ocean boundary type. - real :: weight ! Flux accumulation weight - real :: time_step ! The time step of a call to step_MOM in seconds. +subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & + Ocean_coupling_time_step, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4) + type(ocean_state_type), pointer :: OS !< Structure containing the internal ocean state + type(ocean_public_type), intent(inout) :: Ocean_sfc !< Structure containing all the publicly + !! visible ocean surface fields after a coupling time step + type(time_type), intent(in) :: time_start_update !< Time at the beginning of the update step + type(time_type), intent(in) :: Ocean_coupling_time_step !< Amount of time over which to + !! advance the ocean + real(kind=8), intent(in) :: x2o_o(:,:) !< Fluxes from coupler to ocean, computed by ocean + type(cpl_indices), intent(inout) :: ind !< Structure with MCT attribute vectors and indices + logical, intent(in) :: sw_decomp !< controls if shortwave is + !!decomposed into four components + real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + + ! local variables + type(time_type) :: Master_time !< This allows step_MOM to temporarily change + !! the time that is seen by internal modules. + type(time_type) :: Time1 !< The value of the ocean model's time at the + !! start of a call to step_MOM. + real :: weight !< Flux accumulation weight + real :: time_step !< The time step of a call to step_MOM in seconds. integer :: secs, days integer :: is, ie, js, je @@ -1457,16 +1725,13 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) -! Translate Ice_ocean_boundary into fluxes. - call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), & - index_bnds(3), index_bnds(4)) - weight = 1.0 if (OS%fluxes%fluxes_used) then - call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) ! Needed to allow diagnostics in convert_IOB - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) + ! GMM, is enable_averaging needed now? + call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) + call ocn_import(OS%fluxes, OS%Time, OS%grid, OS%forcing_CSp, OS%state, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4, OS%restore_salinity,OS%restore_temp) #ifdef _USE_GENERIC_TRACER call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes #endif @@ -1487,8 +1752,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%fluxes%dt_buoy_accum = time_step else OS%flux_tmp%C_p = OS%fluxes%C_p - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) + ! Import fluxes from coupler to ocean. Also, perform do SST and SSS restoring, if needed. + call ocn_import(OS%flux_tmp, OS%Time, OS%grid, OS%forcing_CSp, & + OS%state, x2o_o, ind, sw_decomp, c1, c2, c3, c4, & + OS%restore_salinity,OS%restore_temp) + if (OS%use_ice_shelf) then call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) endif @@ -1499,6 +1767,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%State, time_step, OS%berg_area_threshold) !endif + ! Accumulate the forcing over time steps call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight) #ifdef _USE_GENERIC_TRACER call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average @@ -1557,6 +1826,609 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call callTree_leave("update_ocean_model()") end subroutine update_ocean_model +!> This function has a few purposes: 1) it allocates and initializes the data +!! in the fluxes structure; 2) it imports surface fluxes using data from +!! the coupler; and 3) it can apply restoring in SST and SSS. +!! See \ref section_ocn_import for a summary of the surface fluxes that are +!! passed from MCT to MOM6, including fluxes that need to be included in +!! the future. +subroutine ocn_import(fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4, restore_salt, restore_temp) + type(forcing), intent(inout) :: fluxes !< Surface fluxes structure + type(time_type), intent(in) :: Time !< Model time structure + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(surface_forcing_CS), pointer :: CS !< control structure returned by + !! a previous call to surface_forcing_init + type(surface), intent(in) :: state !< control structure to ocean + !! surface state fields. + real(kind=8), intent(in) :: x2o_o(:,:)!< Fluxes from coupler to ocean, computed by ocean + type(cpl_indices), intent(inout) :: ind !< Structure with MCT attribute vectors and indices + logical, intent(in) :: sw_decomp !< controls if shortwave is + !!decomposed into four components + real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + logical, optional, intent(in) :: restore_salt, restore_temp !< Controls if salt and temp are + !! restored + + ! local variables + real, dimension(SZIB_(G),SZJB_(G)) :: & + taux_at_q, & ! Zonal wind stresses at q points (Pa) + tauy_at_q ! Meridional wind stresses at q points (Pa) + + real, dimension(SZI_(G),SZJ_(G)) :: & + taux_at_h, & ! Zonal wind stresses at h points (Pa) + tauy_at_h, & ! Meridional wind stresses at h points (Pa) + data_restore, & ! The surface value toward which to restore (g/kg or degC) + SST_anom, & ! Instantaneous sea surface temperature anomalies from a target value (deg C) + SSS_anom, & ! Instantaneous sea surface salinity anomalies from a target value (g/kg) + SSS_mean, & ! A (mean?) salinity about which to normalize local salinity + ! anomalies when calculating restorative precipitation anomalies (g/kg) + PmE_adj, & ! The adjustment to PminusE that will cause the salinity + ! to be restored toward its target value (kg/(m^2 * s)) + net_FW, & ! The area integrated net freshwater flux into the ocean (kg/s) + net_FW2, & ! The area integrated net freshwater flux into the ocean (kg/s) + work_sum, & ! A 2-d array that is used as the work space for a global + ! sum, used with units of m2 or (kg/s) + open_ocn_mask ! a binary field indicating where ice is present based on frazil criteria + + real :: gustiness ! unresolved gustiness that contributes to ustar (Pa) + real :: Irho0 ! inverse of the mean density in (m^3/kg) + real :: taux2, tauy2 ! squared wind stresses (Pa^2) + real :: tau_mag ! magnitude of the wind stress (Pa) + real :: I_GEarth ! 1.0 / G%G_Earth (s^2/m) + real :: Kv_rho_ice ! (CS%kv_sea_ice / CS%density_sea_ice) ( m^5/(s*kg) ) + real :: mass_ice ! mass of sea ice at a face (kg/m^2) + real :: mass_eff ! effective mass of sea ice for rigidity (kg/m^2) + + integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer + integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd + + logical :: restore_salinity ! local copy of the argument restore_salt, if it + ! is present, or false (no restoring) otherwise. + logical :: restore_sst ! local copy of the argument restore_temp, if it + ! is present, or false (no restoring) otherwise. + real :: delta_sss ! temporary storage for sss diff from restoring value + real :: delta_sst ! temporary storage for sst diff from restoring value + + real :: C_p ! heat capacity of seawater ( J/(K kg) ) + + call cpu_clock_begin(id_clock_forcing) + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1 + + C_p = fluxes%C_p + Irho0 = 1.0/CS%Rho0 + open_ocn_mask(:,:) = 1.0 + pme_adj(:,:) = 0.0 + fluxes%vPrecGlobalAdj = 0.0 + fluxes%vPrecGlobalScl = 0.0 + fluxes%saltFluxGlobalAdj = 0.0 + fluxes%saltFluxGlobalScl = 0.0 + fluxes%netFWGlobalAdj = 0.0 + fluxes%netFWGlobalScl = 0.0 + + restore_salinity = .false. + if (present(restore_salt)) restore_salinity = restore_salt + restore_sst = .false. + if (present(restore_temp)) restore_sst = restore_temp + + ! if true, allocation and initialization + if (fluxes%dt_buoy_accum < 0) then + call allocate_forcing_type(G, fluxes, stress=.true., ustar=.true., & + water=.true., heat=.true.) + + call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%p_surf_SSH,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed) + + if (CS%allow_flux_adjustments) then + call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + endif + + do j=js-2,je+2 ; do i=is-2,ie+2 + fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) + fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) + enddo; enddo + + if (CS%rigid_sea_ice) then + call safe_alloc_ptr(fluxes%rigidity_ice_u,IsdB,IedB,jsd,jed) + call safe_alloc_ptr(fluxes%rigidity_ice_v,isd,ied,JsdB,JedB) + endif + + if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) + + fluxes%dt_buoy_accum = 0.0 + endif ! endif for allocation and initialization + + if (CS%allow_flux_adjustments) then + fluxes%heat_added(:,:)=0.0 + fluxes%salt_flux_added(:,:)=0.0 + endif + + if (CS%area_surf < 0.0) then + do j=js,je ; do i=is,ie + work_sum(i,j) = G%areaT(i,j) * G%mask2dT(i,j) + enddo ; enddo + CS%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer) + endif ! endif for allocation and initialization + + do j=js,je ; do i=is,ie + fluxes%salt_flux(i,j) = 0.0 + fluxes%vprec(i,j) = 0.0 + enddo ; enddo + + ! Salinity restoring logic + if (restore_salinity) then + call time_interp_external(CS%id_srestore,Time,data_restore) + ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) + open_ocn_mask(:,:) = 1.0 + if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice + do j=js,je ; do i=is,ie + if (state%SST(i,j) .le. -0.0539*state%SSS(i,j)) open_ocn_mask(i,j)=0.0 + enddo; enddo + endif + if (CS%salt_restore_as_sflux) then + do j=js,je ; do i=is,ie + delta_sss = data_restore(i,j)- state%SSS(i,j) + delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) + fluxes%salt_flux(i,j) = 1.e-3*G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)* & + (CS%basin_mask(i,j)*open_ocn_mask(i,j)) *delta_sss ! kg Salt m-2 s-1 + enddo; enddo + if (CS%adjust_net_srestore_to_zero) then + if (CS%adjust_net_srestore_by_scaling) then + call adjust_area_mean_to_zero(fluxes%salt_flux, G, fluxes%saltFluxGlobalScl) + fluxes%saltFluxGlobalAdj = 0. + else + work_sum(is:ie,js:je) = G%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je) + fluxes%saltFluxGlobalAdj = reproducing_sum(work_sum(:,:), isr,ier, jsr,jer)/CS%area_surf + fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - fluxes%saltFluxGlobalAdj + endif + endif + fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic + else + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.5) then + delta_sss = state%SSS(i,j) - data_restore(i,j) + delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) + fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j))* & + (CS%Rho0*CS%Flux_const) * & + delta_sss / (0.5*(state%SSS(i,j) + data_restore(i,j))) + endif + enddo; enddo + if (CS%adjust_net_srestore_to_zero) then + if (CS%adjust_net_srestore_by_scaling) then + call adjust_area_mean_to_zero(fluxes%vprec, G, fluxes%vPrecGlobalScl) + fluxes%vPrecGlobalAdj = 0. + else + work_sum(is:ie,js:je) = G%areaT(is:ie,js:je)*fluxes%vprec(is:ie,js:je) + fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / CS%area_surf + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%vPrecGlobalAdj ) * G%mask2dT(i,j) + enddo; enddo + endif + endif + endif + endif + + ! SST restoring logic + if (restore_sst) then + call time_interp_external(CS%id_trestore,Time,data_restore) + do j=js,je ; do i=is,ie + delta_sst = data_restore(i,j)- state%SST(i,j) + delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) + fluxes%heat_added(i,j) = G%mask2dT(i,j) * (CS%Rho0*fluxes%C_p) * delta_sst * CS%Flux_const ! W m-2 + enddo; enddo + endif + + ! GMM, CIME uses AGRID. All the BGRID_NE code can be cleaned later + wind_stagger = AGRID + + if (wind_stagger == BGRID_NE) then + ! This is necessary to fill in the halo points. + taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0 + endif + if (wind_stagger == AGRID) then + ! This is necessary to fill in the halo points. + taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0 + endif + + k = 0 + do j=js,je ; do i=is,ie + k = k + 1 ! Increment position within gindex + + if (wind_stagger == BGRID_NE) then + taux_at_q(I,J) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + tauy_at_q(I,J) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier + ! GMM, cime uses AGRID + elseif (wind_stagger == AGRID) then + taux_at_h(i,j) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + tauy_at_h(i,j) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier + else ! C-grid wind stresses. + fluxes%taux(I,j) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + fluxes%tauy(i,J) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier + endif + + ! liquid precipitation (rain) + if (ASSOCIATED(fluxes%lprec)) & + fluxes%lprec(i,j) = x2o_o(ind%x2o_Faxa_rain,k) * G%mask2dT(i,j) + + ! frozen precipitation (snow) + if (ASSOCIATED(fluxes%fprec)) & + fluxes%fprec(i,j) = x2o_o(ind%x2o_Faxa_snow,k) * G%mask2dT(i,j) + + ! evaporation + if (ASSOCIATED(fluxes%evap)) & + fluxes%evap(i,j) = x2o_o(ind%x2o_Foxx_evap,k) * G%mask2dT(i,j) + + ! river runoff flux + if (ASSOCIATED(fluxes%lrunoff)) & + fluxes%lrunoff(i,j) = x2o_o(ind%x2o_Foxx_rofl,k) * G%mask2dT(i,j) + + ! ice runoff flux + if (ASSOCIATED(fluxes%frunoff)) & + fluxes%frunoff(i,j) = x2o_o(ind%x2o_Foxx_rofi,k) * G%mask2dT(i,j) + + ! GMM, we don't have an icebergs yet so the following is not needed + !if (((ASSOCIATED(IOB%ustar_berg) .and. (.not. ASSOCIATED(fluxes%ustar_berg))) & + ! .or. (ASSOCIATED(IOB%area_berg) .and. (.not. ASSOCIATED(fluxes%area_berg)))) & + ! .or. (ASSOCIATED(IOB%mass_berg) .and. (.not. ASSOCIATED(fluxes%mass_berg)))) & + ! call allocate_forcing_type(G, fluxes, iceberg=.true.) + !if (ASSOCIATED(IOB%ustar_berg)) & + ! fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) + !if (ASSOCIATED(IOB%area_berg)) & + ! fluxes%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) + !if (ASSOCIATED(IOB%mass_berg)) & + ! fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) + + ! GMM, cime does not not have an equivalent for heat_content_lrunoff and + ! heat_content_frunoff. I am seeting these to zero for now. + if (ASSOCIATED(fluxes%heat_content_lrunoff)) & + fluxes%heat_content_lrunoff(i,j) = 0.0 * G%mask2dT(i,j) + + if (ASSOCIATED(fluxes%heat_content_frunoff)) & + fluxes%heat_content_frunoff(i,j) = 0.0 * G%mask2dT(i,j) + + ! longwave radiation, sum up and down (W/m2) + if (ASSOCIATED(fluxes%LW)) & + fluxes%LW(i,j) = (x2o_o(ind%x2o_Faxa_lwdn,k) + x2o_o(ind%x2o_Foxx_lwup,k)) * G%mask2dT(i,j) + + ! sensible heat flux (W/m2) + if (ASSOCIATED(fluxes%sens)) & + fluxes%sens(i,j) = x2o_o(ind%x2o_Foxx_sen,k) * G%mask2dT(i,j) + + ! latent heat flux (W/m^2) + if (ASSOCIATED(fluxes%latent)) & + fluxes%latent(i,j) = x2o_o(ind%x2o_Foxx_lat,k) * G%mask2dT(i,j) + + if (sw_decomp) then + ! Use runtime coefficients to decompose net short-wave heat flux into 4 components + ! 1) visible, direct shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_vis_dir)) & + fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c1 + ! 2) visible, diffuse shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_vis_dif)) & + fluxes%sw_vis_dif(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c2 + ! 3) near-IR, direct shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_nir_dir)) & + fluxes%sw_nir_dir(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c3 + ! 4) near-IR, diffuse shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_nir_dif)) & + fluxes%sw_nir_dif(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c4 + + fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + & + fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j) + else + call MOM_error(FATAL,"fill_ice_ocean_bnd: this option has not been implemented yet."// & + "Shortwave must be decomposed using coeffs. c1, c2, c3, c4."); + endif + + ! applied surface pressure from atmosphere and cryosphere + ! sea-level pressure (Pa) + if (ASSOCIATED(fluxes%p_surf_full) .and. ASSOCIATED(fluxes%p_surf)) then + fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Sa_pslv,k) + if (CS%max_p_surf >= 0.0) then + fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf) + else + fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j) + endif + + if (CS%use_limited_P_SSH) then + fluxes%p_surf_SSH(i,j) = fluxes%p_surf(i,j) + else + fluxes%p_surf_SSH(i,j) = fluxes%p_surf_full(i,j) + endif + + endif + + ! salt flux + ! more salt restoring logic + if (ASSOCIATED(fluxes%salt_flux)) & + fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(x2o_o(ind%x2o_Fioi_salt,k) + fluxes%salt_flux(i,j)) + + if (ASSOCIATED(fluxes%salt_flux_in)) & + fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*x2o_o(ind%x2o_Fioi_salt,k) + + enddo ; enddo + ! ############################ END OF MCT to MOM ############################## + + ! adjust the NET fresh-water flux to zero, if flagged + if (CS%adjust_net_fresh_water_to_zero) then + do j=js,je ; do i=is,ie + net_FW(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + & + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + & + (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * G%areaT(i,j) + ! The following contribution appears to be calculating the volume flux of sea-ice + ! melt. This calculation is clearly WRONG if either sea-ice has variable + ! salinity or the sea-ice is completely fresh. + ! Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system + ! is constant. + ! To do this correctly we will need a sea-ice melt field added to IOB. -AJA + if (ASSOCIATED(fluxes%salt_flux) .and. (CS%ice_salt_concentration>0.0)) & + net_FW(i,j) = net_FW(i,j) - G%areaT(i,j) * & + (fluxes%salt_flux(i,j) / CS%ice_salt_concentration) + net_FW2(i,j) = net_FW(i,j) + enddo ; enddo + + if (CS%adjust_net_fresh_water_by_scaling) then + call adjust_area_mean_to_zero(net_FW2, G, fluxes%netFWGlobalScl) + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = fluxes%vprec(i,j) + (net_FW2(i,j) - net_FW(i,j)) * G%mask2dT(i,j) + enddo; enddo + else + fluxes%netFWGlobalAdj = reproducing_sum(net_FW(:,:), isr, ier, jsr, jer) / CS%area_surf + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%netFWGlobalAdj ) * G%mask2dT(i,j) + enddo; enddo + endif + + endif + + ! surface momentum stress related fields as function of staggering + if (wind_stagger == BGRID_NE) then + if (G%symmetric) & + call fill_symmetric_edges(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE) + call pass_vector(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE) + + do j=js,je ; do I=Isq,Ieq + fluxes%taux(I,j) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) & + fluxes%taux(I,j) = (G%mask2dBu(I,J)*taux_at_q(I,J) + & + G%mask2dBu(I,J-1)*taux_at_q(I,J-1)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) + enddo ; enddo + + do J=Jsq,Jeq ; do i=is,ie + fluxes%tauy(i,J) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) & + fluxes%tauy(i,J) = (G%mask2dBu(I,J)*tauy_at_q(I,J) + & + G%mask2dBu(I-1,J)*tauy_at_q(I-1,J)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) + enddo ; enddo + + ! ustar is required for the bulk mixed layer formulation. The background value + ! of 0.02 Pa is a relatively small value intended to give reasonable behavior + ! in regions of very weak winds. + + do j=js,je ; do i=is,ie + tau_mag = 0.0 ; gustiness = CS%gust_const + if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then + tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_at_q(I,J)**2 + tauy_at_q(I,J)**2) + & + G%mask2dBu(I-1,J-1)*(taux_at_q(I-1,J-1)**2 + tauy_at_q(I-1,J-1)**2)) + & + (G%mask2dBu(I,J-1)*(taux_at_q(I,J-1)**2 + tauy_at_q(I,J-1)**2) + & + G%mask2dBu(I-1,J)*(taux_at_q(I-1,J)**2 + tauy_at_q(I-1,J)**2)) ) / & + ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) + if (CS%read_gust_2d) gustiness = CS%gust(i,j) + endif + fluxes%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) + enddo ; enddo + + elseif (wind_stagger == AGRID) then + call pass_vector(taux_at_h, tauy_at_h, G%Domain,stagger=AGRID) + + do j=js,je ; do I=Isq,Ieq + fluxes%taux(I,j) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) & + fluxes%taux(I,j) = (G%mask2dT(i,j)*taux_at_h(i,j) + & + G%mask2dT(i+1,j)*taux_at_h(i+1,j)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + enddo ; enddo + + do J=Jsq,Jeq ; do i=is,ie + fluxes%tauy(i,J) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) & + fluxes%tauy(i,J) = (G%mask2dT(i,j)*tauy_at_h(i,j) + & + G%mask2dT(i,J+1)*tauy_at_h(i,j+1)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + enddo ; enddo + + do j=js,je ; do i=is,ie + gustiness = CS%gust_const + if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) + fluxes%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & + sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) + enddo ; enddo + + else ! C-grid wind stresses. + if (G%symmetric) & + call fill_symmetric_edges(fluxes%taux, fluxes%tauy, G%Domain) + call pass_vector(fluxes%taux, fluxes%tauy, G%Domain) + + do j=js,je ; do i=is,ie + taux2 = 0.0 + if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) & + taux2 = (G%mask2dCu(I-1,j)*fluxes%taux(I-1,j)**2 + & + G%mask2dCu(I,j)*fluxes%taux(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) + + tauy2 = 0.0 + if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) & + tauy2 = (G%mask2dCv(i,J-1)*fluxes%tauy(i,J-1)**2 + & + G%mask2dCv(i,J)*fluxes%tauy(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) + + if (CS%read_gust_2d) then + fluxes%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) + else + fluxes%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2)) + endif + enddo ; enddo + + endif ! endif for wind related fields + + + ! sea ice related fields + if (CS%rigid_sea_ice) then + ! The commented out code here and in the following lines is the correct + ! version, but the incorrect version is being retained temporarily to avoid + ! changing answers. + call pass_var(fluxes%p_surf_full, G%Domain) + I_GEarth = 1.0 / G%G_Earth + Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice) + do I=isd,ied-1 ; do j=jsd,jed + mass_ice = min(fluxes%p_surf_full(i,j), fluxes%p_surf_full(i+1,j)) * I_GEarth + mass_eff = 0.0 + if (mass_ice > CS%rigid_sea_ice_mass) then + mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / & + (mass_ice + CS%rigid_sea_ice_mass) + endif + ! CAUTION: with both rigid_sea_ice and ice shelves, we will need to make this + ! a maximum for the second call. + fluxes%rigidity_ice_u(I,j) = Kv_rho_ice * mass_eff + enddo ; enddo + do i=isd,ied ; do J=jsd,jed-1 + mass_ice = min(fluxes%p_surf_full(i,j), fluxes%p_surf_full(i,j+1)) * I_GEarth + mass_eff = 0.0 + if (mass_ice > CS%rigid_sea_ice_mass) then + mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / & + (mass_ice + CS%rigid_sea_ice_mass) + endif + fluxes%rigidity_ice_v(i,J) = Kv_rho_ice * mass_eff + enddo ; enddo + endif + + if (CS%allow_flux_adjustments) then + ! Apply adjustments to fluxes + call apply_flux_adjustments(G, CS, Time, fluxes) + endif + + ! Allow for user-written code to alter fluxes after all the above + call user_alter_forcing(state, fluxes, Time, G, CS%urf_CS) + + call cpu_clock_end(id_clock_forcing) + +end subroutine ocn_import + +!> Adds flux adjustments obtained via data_override +!! Component name is 'OCN' +!! Available adjustments are: +!! - taux_adj (Zonal wind stress delta, positive to the east, in Pa) +!! - tauy_adj (Meridional wind stress delta, positive to the north, in Pa) +subroutine apply_flux_adjustments(G, CS, Time, fluxes) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(surface_forcing_CS), pointer :: CS !< Surface forcing control structure + type(time_type), intent(in) :: Time !< Model time structure + type(forcing), optional, intent(inout) :: fluxes !< Surface fluxes structure + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points (Pa) + real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points (Pa) + real, dimension(SZI_(G),SZJ_(G)) :: temp_at_h ! Fluxes at h points (W m-2 or kg m-2 s-1) + + integer :: isc, iec, jsc, jec, i, j + real :: dLonDx, dLonDy, rDlon, cosA, sinA, zonal_tau, merid_tau + logical :: overrode_x, overrode_y, overrode_h + + isc = G%isc; iec = G%iec + jsc = G%jsc; jec = G%jec + + overrode_h = .false. + call data_override('OCN', 'hflx_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%heat_added, G%Domain) + + overrode_h = .false. + call data_override('OCN', 'sflx_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%salt_flux_added, G%Domain) + overrode_h = .false. + + call data_override('OCN', 'prcme_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%vprec(i,j) = fluxes%vprec(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%vprec, G%Domain) + + + tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0 + ! Either reads data or leaves contents unchanged + overrode_x = .false. ; overrode_y = .false. + call data_override('OCN', 'taux_adj', tempx_at_h(isc:iec,jsc:jec), Time, override=overrode_x) + call data_override('OCN', 'tauy_adj', tempy_at_h(isc:iec,jsc:jec), Time, override=overrode_y) + + if (overrode_x .or. overrode_y) then + if (.not. (overrode_x .and. overrode_y)) call MOM_error(FATAL,"apply_flux_adjustments: "//& + "Both taux_adj and tauy_adj must be specified, or neither, in data_table") + + ! Rotate winds + call pass_vector(tempx_at_h, tempy_at_h, G%Domain, To_All, AGRID) + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + dLonDx = G%geoLonCu(I,j) - G%geoLonCu(I-1,j) + dLonDy = G%geoLonCv(i,J) - G%geoLonCv(i,J-1) + rDlon = sqrt( dLonDx * dLonDx + dLonDy * dLonDy ) + if (rDlon > 0.) rDlon = 1. / rDlon + cosA = dLonDx * rDlon + sinA = dLonDy * rDlon + zonal_tau = tempx_at_h(i,j) + merid_tau = tempy_at_h(i,j) + tempx_at_h(i,j) = cosA * zonal_tau - sinA * merid_tau + tempy_at_h(i,j) = sinA * zonal_tau + cosA * merid_tau + enddo ; enddo + + ! Average to C-grid locations + do j=G%jsc,G%jec ; do I=G%isc-1,G%iec + fluxes%taux(I,j) = fluxes%taux(I,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) ) + enddo ; enddo + + do J=G%jsc-1,G%jec ; do i=G%isc,G%iec + fluxes%tauy(i,J) = fluxes%tauy(i,J) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) ) + enddo ; enddo + endif ! overrode_x .or. overrode_y + +end subroutine apply_flux_adjustments + !> Finalizes MOM6 !! !! \todo This needs to be done here. @@ -1579,12 +2451,12 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! !! ocean state to be deallocated upon termination. type(time_type), intent(in) :: Time !< The model time, used for writing restarts. - if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 1' + !if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 1' !GMM call save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) - if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 2' + !if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 2' end subroutine ocean_model_end @@ -1716,4 +2588,69 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) end subroutine ocn_domain_mct +!> \namespace ocn_comp_mct +!! +!! \section section_ocn_import Fluxes imported from the coupler (MCT) to MOM6 +!! The following summarizes the mismatches between MCT and MOM6 in terms +!! of ice ocean fluxes. +!! +!! Redundancies: +!! x2o_Faxa_prec = x2o_Faxa_rain + x2o_Faxa_snow +!! +!! Variables whose units and sign **could not** be verified so far: +!! x2o_Foxx_rofl +!! x2o_Foxx_rof +!! +!! Variables in MOM6 fluxes that are **NOT** filled by the coupler: +!! ustar_berg, frictional velocity beneath icebergs (m/s) +!! area_berg, area covered by icebergs(m2/m2) +!! mass_berg, mass of icebergs(kg/m2) +!! runoff_hflx, heat content of liquid runoff (W/m2) +!! calving_hflx, heat content of frozen runoff (W/m2) +!! mi, mass of ice (kg/m2) +!! +!! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in fluxes): +!! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2) +!! x2o_Fioi_meltw, snow melt flux (kg/m2/s) +!! x2o_Si_ifrac, fractional ice wrt ocean +!! x2o_So_duu10n, 10m wind speed squared (m^2/s^2) +!! x2o_Sa_co2prog, bottom atm level prognostic CO2 +!! x2o_Sa_co2diag, bottom atm level diagnostic CO2 +!! +!! \TODO Langmuir related fields: +!! surface Stokes drift, x-comp. (x2o_Sw_ustokes) +!! surface Stokes drift, y-comp. (x2o_Sw_vstokes) +!! wave model langmuir multiplier (x2o_Sw_lamult) +!! +!! \TODO Biogeochemistry: +!! x2o_Fioi_bcpho, Black Carbon hydrophobic release from sea ice component +!! x2o_Fioi_bcphi, Black Carbon hydrophilic release from sea ice component +!! x2o_Fioi_flxdst, Dust release from sea ice component +!! x2o_Faxa_bcphidry, Black Carbon hydrophilic dry deposition +!! x2o_Faxa_bcphodry, Black Carbon hydrophobic dry deposition +!! x2o_Faxa_bcphiwet, Black Carbon hydrophobic wet deposition +!! x2o_Faxa_ocphidry, Organic Carbon hydrophilic dry deposition +!! x2o_Faxa_ocphodry, Organic Carbon hydrophobic dry deposition +!! x2o_Faxa_ocphiwet, Organic Carbon hydrophilic dry deposition +!! x2o_Faxa_dstwet, Sizes 1 to 4 dust - wet deposition +!! x2o_Faxa_dstdry, Sizes 1 to 4 dust - dry deposition +!! +!! \section section_ocn_export Fluxes exported from MOM6 to the coupler (MCT) +!! +!! Variables that are currently being exported: +!! +!! Surface temperature (Kelvin) +!! Surface salinity (psu) +!! Surface eastward velocity (m/s) +!! Surface northward velocity (m/s) +!! Zonal slope in the sea surface height +!! Meridional slope in the sea surface height +!! +!! \TODO Variables that **are not** currently being exported: +!! +!! Boundary layer depth +!! CO2 +!! DMS +!! o2x_Fioo_q !< Heat flux? + end module ocn_comp_mct