From 69f5ac70d77d7a03277bb19cc0b09aa87a170d1d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Aug 2017 14:23:06 -0600 Subject: [PATCH] Separates ocn_comp_mct from ocean_model_mod This PR is the first attempt to separates ocn_comp_mct from ocean_model_mod. --- config_src/mct_driver/ocn_comp_mct.F90 | 712 ++++++++++++++++++++++++- 1 file changed, 690 insertions(+), 22 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 847c3daba9..b750264589 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -14,7 +14,7 @@ module ocn_comp_mct mct_aVect_nRattr use mct_mod, only: mct_gGrid, mct_gGrid_init, mct_gGrid_importRAttr, & mct_gGrid_importIAttr -use mct_mod, only : mct_avect_indexra, mct_aVect_clean +use mct_mod, only: mct_avect_indexra, mct_aVect_clean use seq_flds_mod, only: seq_flds_x2o_fields, seq_flds_o2x_fields, seq_flds_dom_coord, & seq_flds_dom_other use seq_infodata_mod, only: seq_infodata_type, seq_infodata_GetData, & @@ -26,17 +26,44 @@ module ocn_comp_mct use shr_kind_mod, only: shr_kind_r8 ! MOM6 modules -use ocean_model_mod, only: ocean_state_type, ocean_public_type, ocean_model_init_sfc -use ocean_model_mod, only: ocean_model_init, get_state_pointers, ocean_model_restart -use ocean_model_mod, only: ice_ocean_boundary_type, update_ocean_model -use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here -use MOM_domains, only : pass_var, AGRID -use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_variables, only: surface -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe -use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP -use MOM_file_parser, only: get_param, log_version, param_file_type -use MOM_get_input, only: Get_MOM_Input, directories +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_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS +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_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_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_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 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 ! By default make data private implicit none; private @@ -118,6 +145,99 @@ module ocn_comp_mct 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 + ! 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 @@ -156,6 +276,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc character(len=384) :: runid !< Run ID character(len=384) :: runtype !< Run type + character(len=384) :: restartfile !< Path/Name of restart file character(len=32) :: starttype !< infodata start type integer :: mpicom_ocn !< MPI ocn communicator integer :: npes, pe0 !< # of processors and current processor @@ -303,8 +424,22 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) glb%c1 = 0.0; glb%c2 = 0.0; glb%c3 = 0.0; glb%c4 = 0.0 endif + restartfile = '/glade/scratch/gmarques/mom_test/run/RESTART/MOM.res_Y0001_D002_S00000.nc' ! Initialize the MOM6 model - call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) + if (runtype == "initial") then ! startup (new run) - 'n' is needed below since we don't + ! specify input_filename in input.nml + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = 'n') + else ! hybrid or branch or continuos runs + !call seq_infodata_GetData(glb%infodata, restart_file=restartfile) + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = restartfile) + endif + + if (debug .and. root_pe().eq.pe_here()) then + write(6,*)'runtype, restartfile',runtype, restartfile + endif + + ! Initialize the MOM6 model + !call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) ! Initialize ocn_state%state out of sight call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) @@ -528,6 +663,362 @@ subroutine coupler_indices_init(ind) end subroutine coupler_indices_init +!> Initializes the ocean model, including registering fields +!! for restarts and reading restart files if appropriate. +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file) + type(ocean_public_type), target, & + intent(inout) :: Ocean_sfc !< A structure containing various + !! publicly visible ocean surface properties after initialization, + !! the data in this type is intent(out). + type(ocean_state_type), pointer :: 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. + type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar + type(time_type), intent(in) :: Time_in !< The time at which to initialize the ocean model. + 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, and can be used to spawn related + !! 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. +! 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". +#include "version_variable.h" + 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 + logical :: offline_tracer_mode + + call callTree_enter("ocean_model_init(), ocn_comp_mct.F90") + if (associated(OS)) then + call MOM_error(WARNING, "ocean_model_init called with an associated "// & + "ocean_state_type structure. Model is already initialized.") + return + endif + allocate(OS) + + OS%is_ocean_pe = Ocean_sfc%is_ocean_pe + if (.not.OS%is_ocean_pe) return + + OS%Time = Time_in + call initialize_MOM(OS%Time, param_file, OS%dirs, OS%MOM_CSp, Time_in, & + offline_tracer_mode=offline_tracer_mode, input_restart_file=input_restart_file) + OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV + OS%C_p = OS%MOM_CSp%tv%C_p + OS%fluxes%C_p = OS%MOM_CSp%tv%C_p + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "RESTART_CONTROL", OS%Restart_control, & + "An integer whose bits encode which restart files are \n"//& + "written. Add 2 (bit 1) for a time-stamped file, and odd \n"//& + "(bit 0) for a non-time-stamped file. A restart file \n"//& + "will be saved at the end of the run segment for any \n"//& + "non-negative value.", default=1) + call get_param(param_file, mdl, "TIMEUNIT", Time_unit, & + "The time unit for ENERGYSAVEDAYS.", & + units="s", default=86400.0) + call get_param(param_file, mdl, "ENERGYSAVEDAYS",OS%energysavedays, & + "The interval in units of TIMEUNIT between saves of the \n"//& + "energies of the run and other globally summed diagnostics.", & + default=set_time(0,days=1), timeunit=Time_unit) + + call get_param(param_file, mdl, "OCEAN_SURFACE_STAGGER", stagger, & + "A case-insensitive character string to indicate the \n"//& + "staggering of the surface velocity field that is \n"//& + "returned to the coupler. Valid values include \n"//& + "'A', 'B', or 'C'.", default="C") + if (uppercase(stagger(1:1)) == 'A') then ; Ocean_sfc%stagger = AGRID + elseif (uppercase(stagger(1:1)) == 'B') then ; Ocean_sfc%stagger = BGRID_NE + elseif (uppercase(stagger(1:1)) == 'C') then ; Ocean_sfc%stagger = CGRID_NE + else ; call MOM_error(FATAL,"ocean_model_init: OCEAN_SURFACE_STAGGER = "// & + trim(stagger)//" is invalid.") ; endif + + call get_param(param_file, mdl, "RESTORE_SALINITY",OS%restore_salinity, & + "If true, the coupled driver will add a globally-balanced \n"//& + "fresh-water flux that drives sea-surface salinity \n"//& + "toward specified values.", default=.false.) + call get_param(param_file, mdl, "RESTORE_TEMPERATURE",OS%restore_temp, & + "If true, the coupled driver will add a \n"//& + "heat flux that drives sea-surface temperauture \n"//& + "toward specified values.", default=.false.) + call get_param(param_file, mdl, "RHO_0", 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, "G_EARTH", G_Earth, & + "The gravitational acceleration of the Earth.", & + units="m s-2", default = 9.80) + + call get_param(param_file, mdl, "ICE_SHELF", OS%use_ice_shelf, & + "If true, enables the ice shelf model.", default=.false.) + + call get_param(param_file, mdl, "ICEBERGS_APPLY_RIGID_BOUNDARY", OS%icebergs_apply_rigid_boundary, & + "If true, allows icebergs to change boundary condition felt by ocean", default=.false.) + + if (OS%icebergs_apply_rigid_boundary) then + call get_param(param_file, mdl, "KV_ICEBERG", OS%kv_iceberg, & + "The viscosity of the icebergs", units="m2 s-1",default=1.0e10) + call get_param(param_file, mdl, "DENSITY_ICEBERGS", OS%density_iceberg, & + "A typical density of icebergs.", units="kg m-3", default=917.0) + call get_param(param_file, mdl, "LATENT_HEAT_FUSION", OS%latent_heat_fusion, & + "The latent heat of fusion.", units="J/kg", default=hlf) + call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", OS%berg_area_threshold, & + "Fraction of grid cell which iceberg must occupy, so that fluxes \n"//& + "below berg are set to zero. Not applied for negative \n"//& + " values.", units="non-dim", default=-1.0) + endif + + OS%press_to_z = 1.0/(Rho0*G_Earth) + + ! Consider using a run-time flag to determine whether to do the diagnostic + ! vertical integrals, since the related 3-d sums are not negligible in cost. + call allocate_surface_state(OS%state, OS%grid, OS%MOM_CSp%use_temperature, & + do_integrals=.true., gas_fields_ocn=gas_fields_ocn) + + call surface_forcing_init(Time_in, OS%grid, param_file, OS%MOM_CSp%diag, & + OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) + + if (OS%use_ice_shelf) then + call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & + OS%MOM_CSp%diag, OS%fluxes) + endif + if (OS%icebergs_apply_rigid_boundary) then + !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.) + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + if (.not. OS%use_ice_shelf) call allocate_forcing_type(OS%grid, OS%fluxes, ustar=.true., shelf=.true.) + endif + + call MOM_sum_output_init(OS%grid, param_file, OS%dirs%output_directory, & + OS%MOM_CSp%ntrunc, Time_init, OS%sum_output_CSp) + + ! This call has been moved into the first call to update_ocean_model. +! call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & +! OS%Time, 0, OS%grid, OS%GV, OS%sum_output_CSp, OS%MOM_CSp%tracer_flow_CSp) + + ! write_energy_time is the next integral multiple of energysavedays. + OS%write_energy_time = Time_init + OS%energysavedays * & + (1 + (OS%Time - Time_init) / OS%energysavedays) + + if (ASSOCIATED(OS%grid%Domain%maskmap)) then + call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & + OS%MOM_CSp%diag, maskmap=OS%grid%Domain%maskmap, & + gas_fields_ocn=gas_fields_ocn) + else + call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & + OS%MOM_CSp%diag, gas_fields_ocn=gas_fields_ocn) + endif + + ! This call can only occur here if the coupler_bc_type variables have been + ! initialized already using the information from gas_fields_ocn. + if (present(gas_fields_ocn)) then + call calculate_surface_state(OS%state, OS%MOM_CSp%u, & + OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& + OS%grid, OS%GV, OS%MOM_CSp) + + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + endif + + call close_param_file(param_file) + call diag_mediator_close_registration(OS%MOM_CSp%diag) + + if (is_root_pe()) & + write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' + + call callTree_leave("ocean_model_init(") +end subroutine ocean_model_init + +!> This subroutine 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. +subroutine ocean_model_init_sfc(OS, Ocean_sfc) + type(ocean_state_type), pointer :: OS + type(ocean_public_type), intent(inout) :: Ocean_sfc + + integer :: is, ie, js, je + + is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec + call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, & + (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) + + call calculate_surface_state(OS%state, OS%MOM_CSp%u, & + OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& + OS%grid, OS%GV, OS%MOM_CSp) + + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + +end subroutine ocean_model_init_sfc + +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(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. + + 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) + call mpp_get_global_domain(input_domain, xsize=xsz, ysize=ysz) + if(PRESENT(maskmap)) then + call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain, maskmap=maskmap) + else + call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain) + endif + call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec) + + allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), & + Ocean_sfc%s_surf (isc:iec,jsc:jec), & + Ocean_sfc%u_surf (isc:iec,jsc:jec), & + Ocean_sfc%v_surf (isc:iec,jsc:jec), & + Ocean_sfc%sea_lev(isc:iec,jsc:jec), & + Ocean_sfc%area (isc:iec,jsc:jec), & + Ocean_sfc%frazil (isc:iec,jsc:jec)) + + Ocean_sfc%t_surf = 0.0 ! time averaged sst (Kelvin) passed to atmosphere/ice model + Ocean_sfc%s_surf = 0.0 ! time averaged sss (psu) passed to atmosphere/ice models + Ocean_sfc%u_surf = 0.0 ! time averaged u-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%v_surf = 0.0 ! time averaged v-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%sea_lev = 0.0 ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav + Ocean_sfc%frazil = 0.0 ! time accumulated frazil (J/m^2) passed to ice model + Ocean_sfc%area = 0.0 + Ocean_sfc%axes = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics + + if (present(gas_fields_ocn)) then + call coupler_type_spawn(gas_fields_ocn, Ocean_sfc%fields, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix = '_ocn', as_needed=.true.) + endif + +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. +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 + + ! local variables + real :: IgR0 + character(len=48) :: val_str + integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd + integer :: i, j, i0, j0, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + call pass_vector(state%u,state%v,G%Domain) + + call mpp_get_compute_domain(Ocean_sfc%Domain, isc_bnd, iec_bnd, & + jsc_bnd, jec_bnd) + if (present(patm)) then + ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd). + if (.not.present(press_to_z)) call MOM_error(FATAL, & + 'convert_state_to_ocean_type: press_to_z must be present if patm is.') + endif + + i0 = is - isc_bnd ; j0 = js - jsc_bnd + !If directed convert the surface T&S + !from conservative T to potential T and + !from absolute (reference) salinity to practical salinity + ! + if(use_conT_absS) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(state%SSS(i+i0,j+j0)) + Ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(state%SSS(i+i0,j+j0),state%SST(i+i0,j+j0)) + CELSIUS_KELVIN_OFFSET + enddo ; enddo + else + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%t_surf(i,j) = state%SST(i+i0,j+j0) + CELSIUS_KELVIN_OFFSET + Ocean_sfc%s_surf(i,j) = state%SSS(i+i0,j+j0) + enddo ; enddo + endif + + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%sea_lev(i,j) = state%sea_lev(i+i0,j+j0) + if (present(patm)) & + Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + patm(i,j) * press_to_z + if (associated(state%frazil)) & + Ocean_sfc%frazil(i,j) = state%frazil(i+i0,j+j0) + Ocean_sfc%area(i,j) = G%areaT(i+i0,j+j0) + enddo ; enddo + + if (Ocean_sfc%stagger == AGRID) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0)*0.5*(state%u(I+i0,j+j0)+state%u(I-1+i0,j+j0)) + Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0)*0.5*(state%v(i+i0,J+j0)+state%v(i+i0,J-1+j0)) + enddo ; enddo + elseif (Ocean_sfc%stagger == BGRID_NE) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0)*0.5*(state%u(I+i0,j+j0)+state%u(I+i0,j+j0+1)) + Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0)*0.5*(state%v(i+i0,J+j0)+state%v(i+i0+1,J+j0)) + enddo ; enddo + elseif (Ocean_sfc%stagger == CGRID_NE) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*state%u(I+i0,j+j0) + Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*state%v(i+i0,J+j0) + enddo ; enddo + else + write(val_str, '(I8)') Ocean_sfc%stagger + call MOM_error(FATAL, "convert_state_to_ocean_type: "//& + "Ocean_sfc%stagger has the unrecognized value of "//trim(val_str)) + endif + + if (coupler_type_initialized(state%tr_fields)) then + if (.not.coupler_type_initialized(Ocean_sfc%fields)) then + call MOM_error(FATAL, "convert_state_to_ocean_type: "//& + "Ocean_sfc%fields has not been initialized.") + endif + call coupler_type_copy_data(state%tr_fields, Ocean_sfc%fields) + endif + +end subroutine convert_state_to_ocean_type + +!> Returns pointers to objects within ocean_state_type +subroutine get_state_pointers(OS, grid, surf) + type(ocean_state_type), pointer :: OS !< Ocean state type + type(ocean_grid_type), optional, pointer :: grid !< Ocean grid + type(surface), optional, pointer :: surf !< Ocean surface state + + if (present(grid)) grid => OS%grid + if (present(surf)) surf=> OS%state + +end subroutine get_state_pointers + !> Maps outgoing ocean data to MCT buffer subroutine ocn_export(ind, ocn_public, grid, o2x) type(cpl_indices), intent(inout) :: ind !< Structure with coupler @@ -741,12 +1232,13 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(ESMF_time) :: time_start_ESMF !< Time at the start of the coupling interval type(ESMF_timeInterval) :: ocn_cpl_interval !< The length of one ocean coupling interval integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc - logical :: write_restart_at_eod + logical :: write_restart_at_eod !< Controls if restart files must be written type(time_type) :: time_start !< Start of coupled time interval to pass to MOM6 type(time_type) :: coupling_timestep !< Coupled time interval to pass to MOM6 - character(len=128) :: err_msg + character(len=128) :: err_msg !< Error message character(len=32) :: timestamp !< Name of intermediate restart file - + character(len=80) :: restartname !< The restart file name (no dir) + character(len=384) :: runid !< Run ID ! Compute the time at the start of this coupling interval call ESMF_ClockGet(EClock, PrevTime=time_start_ESMF, rc=rc) call ESMF_TimeGet(time_start_ESMF, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) @@ -783,7 +1275,6 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) ! 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) - if (debug .and. is_root_pe()) write(6,*) 'fill_ice_ocean_bnd' call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, & time_start, coupling_timestep) @@ -793,15 +1284,192 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) if (debug .and. is_root_pe()) write(6,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod - if (write_restart_at_eod) then - !timestamp = date_to_string(EClock) - ! \todo add time stamp to ocean_model_restart + !!if (write_restart_at_eod) then + ! case name + !!call seq_infodata_GetData( glb%infodata, case_name=runid ) + ! add time stamp to restart filename + ! \todo use MCT call to get year, month etc + !!call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) + !call get_date(glb%ocn_state%Time,year,month,days,hour,minute,seconds) + !!seconds = seconds + hour*3600 + minute*60 + !!write(restartname,'(A,".mom6.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') runid, year, month, day, seconds !call ocean_model_restart(glb%ocn_state, timestamp) - call ocean_model_restart(glb%ocn_state) - endif + !!call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & + !! glb%ocn_state%MOM_CSp%restart_CSp, filename=restartname,GV=glb%ocn_state%GV) + + !! GM call via ocean_model_MOM.F90 call ocean_model_restart(glb%ocn_state) + + ! Is this needed? + !call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & + ! OS%dirs%restart_output_dir, .true.) + ! Once we start using the ice shelf module, the following will be needed + !if (glb%ocn_state%use_ice_shelf) then + ! call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir, .true.) + !endif + + !!endif end subroutine ocn_run_mct +!> 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 +!! 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 +!! 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. + integer :: secs, days + integer :: is, ie, js, je + + call callTree_enter("update_ocean_model(), ocn_comp_mct.F90") + call get_time(Ocean_coupling_time_step, secs, days) + time_step = 86400.0*real(days) + real(secs) + + if (time_start_update /= OS%Time) then + call MOM_error(WARNING, "update_ocean_model: internal clock does not "//& + "agree with time_start_update argument.") + endif + + if (.not.associated(OS)) then + call MOM_error(FATAL, "update_ocean_model called with an unassociated "// & + "ocean_state_type structure. ocean_model_init must be "// & + "called first to allocate this structure.") + return + endif + + ! This is benign but not necessary if ocean_model_init_sfc was called or if + ! OS%state%tr_fields was spawnded in ocean_model_init. Consider removing it. + is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec + 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) +#ifdef _USE_GENERIC_TRACER + call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes +#endif + + ! Add ice shelf fluxes + if (OS%use_ice_shelf) then + call shelf_calc_flux(OS%State, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp) + endif + + ! GMM, check ocean_model_MOM.F90 to enable the following option + !if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + ! call add_berg_flux_to_shelf(OS%grid, OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%State, time_step, OS%berg_area_threshold) + !endif + + ! Indicate that there are new unused fluxes. + OS%fluxes%fluxes_used = .false. + 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) + if (OS%use_ice_shelf) then + call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) + endif + + ! GMM, check ocean_model_MOM.F90 to enable the following option + !if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + ! 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 + + 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 +#endif + endif + + if (OS%nstep==0) then + call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%fluxes) + + call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & + OS%Time, 0, OS%grid, OS%GV, OS%sum_output_CSp, & + OS%MOM_CSp%tracer_flow_CSp) + endif + + call disable_averaging(OS%MOM_CSp%diag) + Master_time = OS%Time ; Time1 = OS%Time + + if(OS%MOM_Csp%offline_tracer_mode) then + call step_offline(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + else + call step_MOM(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + endif + + OS%Time = Master_time + Ocean_coupling_time_step + OS%nstep = OS%nstep + 1 + + call enable_averaging(time_step, OS%Time, OS%MOM_CSp%diag) + call mech_forcing_diags(OS%fluxes, time_step, OS%grid, & + OS%MOM_CSp%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%MOM_CSp%diag) + + if (OS%fluxes%fluxes_used) then + call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%MOM_CSp%diag) + call forcing_diagnostics(OS%fluxes, OS%state, OS%fluxes%dt_buoy_accum, & + OS%grid, OS%MOM_CSp%diag, OS%forcing_CSp%handles) + call accumulate_net_input(OS%fluxes, OS%state, OS%fluxes%dt_buoy_accum, & + OS%grid, OS%sum_output_CSp) + call disable_averaging(OS%MOM_CSp%diag) + endif + +! See if it is time to write out the energy. + if ((OS%Time + ((Ocean_coupling_time_step)/2) > OS%write_energy_time) .and. & + (OS%MOM_CSp%t_dyn_rel_adv==0.0)) then + call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & + OS%Time, OS%nstep, OS%grid, OS%GV, OS%sum_output_CSp, & + OS%MOM_CSp%tracer_flow_CSp) + OS%write_energy_time = OS%write_energy_time + OS%energysavedays + endif + +! Translate state into Ocean. +! call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & +! Ice_ocean_boundary%p, OS%press_to_z) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + + call callTree_leave("update_ocean_model()") +end subroutine update_ocean_model + !> Finalizes MOM6 !! !! \todo This needs to be done here.