diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 756c648120..250d554361 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -17,7 +17,7 @@ module ocean_model_mod ! in the same way as MOM4. ! -use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end +use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_state_type, MOM_end use MOM, only : calculate_surface_state, allocate_surface_state, finish_MOM_initialization use MOM, only : step_offline use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf @@ -35,9 +35,7 @@ module ocean_model_mod use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type use MOM_io, only : close_file, file_exists, read_data, write_version_number -use MOM_restart, only : save_restart -use MOM_sum_output, only : write_energy, accumulate_net_input -use MOM_sum_output, only : MOM_sum_output_init, sum_output_CS +use MOM_restart, only : MOM_restart_CS, save_restart use MOM_string_functions, only : uppercase use MOM_surface_forcing, only : surface_forcing_init, convert_IOB_to_fluxes use MOM_surface_forcing, only : ice_ocn_bnd_type_chksum @@ -140,29 +138,16 @@ module ocean_model_mod type, public :: ocean_state_type ; private ! This type is private, and can therefore vary between different ocean models. 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) :: energysavedays_geometric ! The starting interval for computing a geometric - ! progression of time deltas between calls to - ! write_energy. This interval will increase by a factor of 2. - ! after each call to write_energy. - logical :: energysave_geometric ! Logical to control whether calls to write_energy should - ! follow a geometric progression - type(time_type) :: write_energy_time ! The next time to write to the energy file. - type(time_type) :: geometric_end_time ! Time at which to stop the geometric progression - ! of calls to write_energy and revert to the standard - ! energysavedays interval - - integer :: nstep = 0 ! The number of calls to update_ocean. - logical :: use_ice_shelf ! If true, the ice shelf model is enabled. + 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. + + integer :: nstep = 0 !< The number of calls to update_ocean. + logical :: use_ice_shelf !< If true, the ice shelf model is enabled. ! Many of the following variables do not appear to belong here. -RWH logical :: icebergs_apply_rigid_boundary ! If true, the icebergs can change ocean bd condition. @@ -173,31 +158,51 @@ module ocean_model_mod 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. + 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. + logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode + !! with the barotropic and baroclinic dynamics, thermodynamics, + !! etc. stepped forward integrated in time. + !! If true, all of the above are bypassed with all + !! fields necessary to integrate only the tracer advection + !! and diffusion equation read in from files stored from + !! a previous integration of the prognostic model. + + type(directories) :: dirs !< A structure containing several relevant directory paths. type(mech_forcing) :: forces !< A structure with the driving mechanical surface forces - type(forcing) :: fluxes ! A structure containing pointers to - ! the thermodynamic 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) :: sfc_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() + type(forcing) :: fluxes !< A structure containing pointers to + !! the thermodynamic 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) :: sfc_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 structure containing information + !! about the vertical grid. + type(MOM_control_struct), pointer :: & + MOM_CSp => NULL() !< A pointer to the MOM control structure + type(MOM_state_type), pointer :: & + MSp => NULL() !< A pointer to the MOM_state_type + type(ice_shelf_CS), pointer :: & + Ice_shelf_CSp => NULL() !< A pointer to the control structure for the + !! ice shelf model that couples with MOM6. This + !! is null if there is no ice shelf. + type(surface_forcing_CS), pointer :: & + forcing_CSp => NULL() !< A pointer to the MOM forcing control structure + type(MOM_restart_CS), pointer :: & + restart_CSp => NULL() !< A pointer set to the restart control structure + !! that will be used for MOM restart files. + type(diag_ctrl), pointer :: & + diag => NULL() !< A pointer to the diagnostic regulatory structure end type ocean_state_type contains @@ -239,7 +244,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) ! 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". @@ -248,7 +252,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) 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 + logical :: use_temperature type(time_type) :: dt_geometric, dt_savedays, dt_from_base call callTree_enter("ocean_model_init(), ocean_model_MOM.F90") @@ -263,11 +267,13 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) 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) - 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 + call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MSp, OS%MOM_CSp, & + OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & + diag_ptr=OS%diag, count_calls=.true.) + OS%grid => OS%MSp%G ; OS%GV => OS%MSp%GV + OS%C_p = OS%MSp%tv%C_p + OS%fluxes%C_p = OS%MSp%tv%C_p + use_temperature = ASSOCIATED(OS%MSp%tv%T) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") @@ -277,27 +283,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) "(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, "ENERGYSAVEDAYS_GEOMETRIC",OS%energysavedays_geometric, & - "The starting interval in units of TIMEUNIT for the first call \n"//& - "to save the energies of the run and other globally summed diagnostics. \n"//& - "The interval increases by a factor of 2. after each call to write_energy.",& - default=set_time(seconds=0), timeunit=Time_unit) - - if ((time_type_to_real(OS%energysavedays_geometric) > 0.) .and. & - (OS%energysavedays_geometric < OS%energysavedays)) then - OS%energysave_geometric = .true. - else - OS%energysave_geometric = .false. - endif - 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"//& @@ -350,15 +335,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) ! 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%sfc_state, OS%grid, OS%MOM_CSp%use_temperature, & + call allocate_surface_state(OS%sfc_state, OS%grid, 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, & + call surface_forcing_init(Time_in, OS%grid, param_file, OS%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%forces, OS%fluxes) + OS%diag, OS%forces, OS%fluxes) endif if (OS%icebergs_apply_rigid_boundary) then !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.) @@ -367,35 +352,13 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) call allocate_forcing_type(OS%grid, OS%fluxes, 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. - if (OS%energysave_geometric) then - if (OS%energysavedays_geometric < OS%energysavedays) then - OS%write_energy_time = OS%Time + OS%energysavedays_geometric - OS%geometric_end_time = Time_init + OS%energysavedays * & - (1 + (OS%Time - Time_init) / OS%energysavedays) - else - OS%write_energy_time = Time_init + OS%energysavedays * & - (1 + (OS%Time - Time_init) / OS%energysavedays) - endif - else - OS%write_energy_time = Time_init + OS%energysavedays * & - (1 + (OS%Time - Time_init) / OS%energysavedays) - endif - 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, & + OS%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) + OS%diag, gas_fields_ocn=gas_fields_ocn) endif ! This call can only occur here if the coupler_bc_type variables have been @@ -404,16 +367,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) call coupler_type_set_diags(Ocean_sfc%fields, "ocean_sfc", & Ocean_sfc%axes(1:2), Time_in) - call calculate_surface_state(OS%sfc_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 calculate_surface_state(OS%sfc_state, OS%MSp%u, & + OS%MSp%v, OS%MSp%h, OS%MSp%ave_ssh,& + OS%grid, OS%GV, OS%MSp, OS%MOM_CSp) call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) endif call close_param_file(param_file) - call diag_mediator_close_registration(OS%MOM_CSp%diag) + call diag_mediator_close_registration(OS%diag) if (is_root_pe()) & write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' @@ -480,7 +443,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & real :: time_step ! The time step of a call to step_MOM in seconds. integer :: secs, days integer :: is, ie, js, je - type(time_type) :: write_energy_time_geometric call callTree_enter("update_ocean_model(), ocean_model_MOM.F90") call get_time(Ocean_coupling_time_step, secs, days) @@ -510,7 +472,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & 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 enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%diag) ! Needed to allow diagnostics in convert_IOB call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%forces, OS%fluxes, index_bnds, OS%Time, & OS%grid, OS%forcing_CSp, OS%sfc_state, OS%restore_salinity,OS%restore_temp) @@ -560,67 +522,34 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) 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) + call finish_MOM_initialization(OS%Time, OS%dirs, OS%MSp, OS%MOM_CSp, OS%fluxes, & + OS%restart_CSp) endif - call disable_averaging(OS%MOM_CSp%diag) + call disable_averaging(OS%diag) Master_time = OS%Time ; Time1 = OS%Time - if(OS%MOM_Csp%offline_tracer_mode) then - call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MOM_CSp) + if(OS%offline_tracer_mode) then + call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MSp, OS%MOM_CSp) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MOM_CSp) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MSp, 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 enable_averaging(time_step, OS%Time, OS%diag) call mech_forcing_diags(OS%forces, OS%fluxes, time_step, OS%grid, & - OS%MOM_CSp%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%MOM_CSp%diag) + OS%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%diag) if (OS%fluxes%fluxes_used) then - call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%MOM_CSp%diag) + call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & - OS%grid, OS%MOM_CSp%diag, OS%forcing_CSp%handles) - call accumulate_net_input(OS%fluxes, OS%sfc_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%energysave_geometric) then - if ((OS%Time + ((Ocean_coupling_time_step)/2) > OS%geometric_end_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%geometric_end_time + OS%energysavedays - OS%energysave_geometric = .false. ! stop geometric progression - endif + OS%grid, OS%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%diag) endif - 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) - if (OS%energysave_geometric) then - OS%energysavedays_geometric = OS%energysavedays_geometric*2 - OS%write_energy_time = OS%write_energy_time + OS%energysavedays_geometric - else - OS%write_energy_time = OS%write_energy_time + OS%energysavedays - endif - endif - - - ! Translate state into Ocean. ! call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, & ! Ice_ocean_boundary%p, OS%press_to_z) @@ -746,7 +675,7 @@ subroutine ocean_model_restart(OS, timestamp) type(ocean_state_type), pointer :: OS character(len=*), intent(in), optional :: timestamp - if (OS%MOM_CSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& + if (OS%MSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& "with inconsistent dynamics and advective times. Additional restart fields "//& "that have not been coded yet would be required for reproducibility.") if (.not.OS%fluxes%fluxes_used) call MOM_error(FATAL, "ocean_model_restart "//& @@ -755,7 +684,7 @@ subroutine ocean_model_restart(OS, timestamp) if (BTEST(OS%Restart_control,1)) then call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%MOM_CSp%restart_CSp, .true., GV=OS%GV) + OS%restart_CSp, .true., GV=OS%GV) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, .true.) if (OS%use_ice_shelf) then @@ -764,7 +693,7 @@ subroutine ocean_model_restart(OS, timestamp) endif if (BTEST(OS%Restart_control,0)) then call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%MOM_CSp%restart_CSp, GV=OS%GV) + OS%restart_CSp, GV=OS%GV) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then @@ -802,8 +731,8 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! (in) Time - The model time, used for writing restarts. call ocean_model_save_restart(Ocean_state, Time) - call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) - call MOM_end(Ocean_state%MOM_CSp) + call diag_mediator_end(Time, Ocean_state%diag) + call MOM_end(Ocean_state%MSp, Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end ! NAME="ocean_model_end" @@ -830,7 +759,7 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) ! restart behavior as now in FMS. character(len=200) :: restart_dir - if (OS%MOM_CSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& + if (OS%MSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& "with inconsistent dynamics and advective times. Additional restart fields "//& "that have not been coded yet would be required for reproducibility.") if (.not.OS%fluxes%fluxes_used) call MOM_error(FATAL, "ocean_model_save_restart "//& @@ -840,7 +769,7 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) if (present(directory)) then ; restart_dir = directory else ; restart_dir = OS%dirs%restart_output_dir ; endif - call save_restart(restart_dir, Time, OS%grid, OS%MOM_CSp%restart_CSp, GV=OS%GV) + call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) @@ -1029,9 +958,9 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) call coupler_type_spawn(Ocean_sfc%fields, OS%sfc_state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) - call calculate_surface_state(OS%sfc_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 calculate_surface_state(OS%sfc_state, OS%MSp%u, & + OS%MSp%v, OS%MSp%h, OS%MSp%ave_ssh,& + OS%grid, OS%GV, OS%MSp, OS%MOM_CSp) call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) @@ -1094,30 +1023,30 @@ subroutine Ocean_stock_pe(OS, index, value, time_index) to_mass = OS%GV%H_to_kg_m2 if (OS%GV%Boussinesq) then do k=1,nz ; do j=js,je ; do i=is,ie ; if (OS%grid%mask2dT(i,j) > 0.5) then - value = value + to_mass*(OS%MOM_CSp%h(i,j,k) * OS%grid%areaT(i,j)) + value = value + to_mass*(OS%MSp%h(i,j,k) * OS%grid%areaT(i,j)) endif ; enddo ; enddo ; enddo else ! In non-Boussinesq mode, the mass of salt needs to be subtracted. PSU_to_kg = 1.0e-3 do k=1,nz ; do j=js,je ; do i=is,ie ; if (OS%grid%mask2dT(i,j) > 0.5) then - value = value + to_mass * ((1.0 - PSU_to_kg*OS%MOM_CSp%tv%S(i,j,k))*& - (OS%MOM_CSp%h(i,j,k) * OS%grid%areaT(i,j))) + value = value + to_mass * ((1.0 - PSU_to_kg*OS%MSp%tv%S(i,j,k))*& + (OS%MSp%h(i,j,k) * OS%grid%areaT(i,j))) endif ; enddo ; enddo ; enddo endif case (ISTOCK_HEAT) ! Return the heat content of the ocean on this PE in J. to_heat = OS%GV%H_to_kg_m2 * OS%C_p do k=1,nz ; do j=js,je ; do i=is,ie ; if (OS%grid%mask2dT(i,j) > 0.5) then - value = value + (to_heat * OS%MOM_CSp%tv%T(i,j,k)) * & - (OS%MOM_CSp%h(i,j,k)*OS%grid%areaT(i,j)) + value = value + (to_heat * OS%MSp%tv%T(i,j,k)) * & + (OS%MSp%h(i,j,k)*OS%grid%areaT(i,j)) endif ; enddo ; enddo ; enddo case (ISTOCK_SALT) ! Return the mass of the salt in the ocean on this PE in kg. ! The 1000 converts salinity in PSU to salt in kg kg-1. to_salt = OS%GV%H_to_kg_m2 / 1000.0 do k=1,nz ; do j=js,je ; do i=is,ie ; if (OS%grid%mask2dT(i,j) > 0.5) then - value = value + (to_salt * OS%MOM_CSp%tv%S(i,j,k)) * & - (OS%MOM_CSp%h(i,j,k)*OS%grid%areaT(i,j)) + value = value + (to_salt * OS%MSp%tv%S(i,j,k)) * & + (OS%MSp%h(i,j,k)*OS%grid%areaT(i,j)) endif ; enddo ; enddo ; enddo case default ; value = 0.0 end select diff --git a/config_src/ice_solo_driver/ice_shelf_driver.F90 b/config_src/ice_solo_driver/ice_shelf_driver.F90 index f879807916..628b138639 100644 --- a/config_src/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/ice_solo_driver/ice_shelf_driver.F90 @@ -283,8 +283,6 @@ program SHELF_main call log_param(param_file, mdl, "ELAPSED TIME AS MASTER", elapsed_time_master) ! i don't think we'll use this... -! call MOM_sum_output_init(grid, param_file, dirs%output_directory, & -! MOM_CSp%ntrunc, Start_time, sum_output_CSp) call MOM_write_cputime_init(param_file, dirs%output_directory, Start_time, & write_CPU_CSp) call MOM_mesg("Done MOM_write_cputime_init.", 5) @@ -292,7 +290,7 @@ program SHELF_main ! Close the param_file. No further parsing of input is possible after this. call close_param_file(param_file) - call diag_mediator_close_registration(MOM_CSp%diag) +! call diag_mediator_close_registration(diag) ! Write out a time stamp file. call open_file(unit, 'time_stamp.out', form=ASCII_FILE, action=APPEND_FILE, & diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 52aeebb9de..39e477c2eb 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -33,7 +33,7 @@ module ocn_comp_mct 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: initialize_MOM, step_MOM, MOM_control_struct, MOM_state_type, 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 @@ -61,8 +61,6 @@ module ocn_comp_mct 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, hlv use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct @@ -165,7 +163,7 @@ module ocn_comp_mct !! in radiation computations, !! per column integer, dimension(:), allocatable :: x2o_qsw_fracr_col !< qsw * fracr, per column - end type cpl_indices +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 @@ -301,9 +299,6 @@ module ocn_comp_mct !! 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. @@ -321,22 +316,34 @@ module ocn_comp_mct 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. + logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode + !! with the barotropic and baroclinic dynamics, thermodynamics, + !! etc. stepped forward integrated in time. + !! If true, all of the above are bypassed with all + !! fields necessary to integrate only the tracer advection + !! and diffusion equation read in from files stored from + !! a previous integration of the prognostic model. type(directories) :: dirs !< A structure containing several relevant directory paths. - type(mech_forcing) :: forces!< A structure with the driving mechanical forces - type(forcing) :: fluxes !< A structure containing pointers to + type(mech_forcing) :: forces!< A structure with the driving mechanical surface forces + type(forcing) :: fluxes !< A structure containing pointers to !! the ocean forcing fields. - type(forcing) :: flux_tmp !< A secondary structure containing pointers to the + 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 + type(surface) :: sfc_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(MOM_state_type), pointer :: MSp => NULL() type(surface_forcing_CS), pointer :: forcing_CSp => NULL() - type(sum_output_CS), pointer :: sum_output_CSp => NULL() + type(MOM_restart_CS), pointer :: & + restart_CSp => NULL() !< A pointer set to the restart control structure + !! that will be used for MOM restart files. + type(diag_ctrl), pointer :: & + diag => NULL() !< A pointer to the diagnostic regulatory structure end type ocean_state_type !> Control structure for this module @@ -551,7 +558,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call ocean_model_init(glb%ocn_public, glb%ocn_state, time0, time0, input_restart_file=trim(restartfile)) endif - ! Initialize ocn_state%state out of sight + ! Initialize ocn_state%sfc_state out of sight call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) ! store pointers to components inside MOM @@ -777,16 +784,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! Because of the way that indicies and domains are handled, Ocean_sfc must have ! been used in a previous call to initialize_ocean_type. - 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 + logical :: use_temperature 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 @@ -800,11 +806,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i 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 + call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MSp, OS%MOM_CSp, & + OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & + input_restart_file=input_restart_file, diag_ptr=OS%diag, & + count_calls=.true.) + OS%grid => OS%MSp%G ; OS%GV => OS%MSp%GV + OS%C_p = OS%MSp%tv%C_p + OS%fluxes%C_p = OS%MSp%tv%C_p + use_temperature = ASSOCIATED(OS%MSp%tv%T) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") @@ -814,14 +823,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i "(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"//& @@ -874,15 +875,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! 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, & + call allocate_surface_state(OS%sfc_state, OS%grid, 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, & + call surface_forcing_init(Time_in, OS%grid, param_file, OS%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%forces, OS%fluxes) + OS%diag, OS%forces, OS%fluxes) endif if (OS%icebergs_apply_rigid_boundary) then !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.) @@ -890,38 +891,27 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i 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, & + OS%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) + OS%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 calculate_surface_state(OS%sfc_state, OS%MSp%u, & + OS%MSp%v, OS%MSp%h, OS%MSp%ave_ssh,& + OS%grid, OS%GV, OS%MSp, OS%MOM_CSp) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid) + call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) endif call close_param_file(param_file) - call diag_mediator_close_registration(OS%MOM_CSp%diag) + call diag_mediator_close_registration(OS%diag) if (is_root_pe()) & write(glb%stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' @@ -940,14 +930,14 @@ subroutine ocean_model_init_sfc(OS, 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, & + call coupler_type_spawn(Ocean_sfc%fields, OS%sfc_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 calculate_surface_state(OS%sfc_state, OS%MSp%u, & + OS%MSp%v, OS%MSp%h, OS%MSp%ave_ssh,& + OS%grid, OS%GV, OS%Msp, OS%MOM_CSp) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid) + call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) end subroutine ocean_model_init_sfc @@ -1412,7 +1402,7 @@ subroutine get_state_pointers(OS, grid, surf) type(surface), optional, pointer :: surf !< Ocean surface state if (present(grid)) grid => OS%grid - if (present(surf)) surf=> OS%state + if (present(surf)) surf=> OS%sfc_state end subroutine get_state_pointers @@ -1619,7 +1609,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') trim(runid), year, month, day, seconds call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & - glb%ocn_state%MOM_CSp%restart_CSp, .false., filename=restartname,GV=glb%ocn_state%GV) + glb%ocn_state%restart_CSp, .false., filename=restartname,GV=glb%ocn_state%GV) ! write name of restart file in the rpointer file nu = shr_file_getUnit() @@ -1719,17 +1709,17 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & 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. + ! OS%sfc_state%tr_fields was spawned 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, & + call coupler_type_spawn(Ocean_sfc%fields, OS%sfc_state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) weight = 1.0 if (OS%fluxes%fluxes_used) then ! 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%forces, OS%fluxes, OS%Time, OS%grid, OS%forcing_CSp, OS%state, x2o_o, ind, sw_decomp, & + call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%diag) + call ocn_import(OS%forces, OS%fluxes, OS%Time, OS%grid, OS%forcing_CSp, OS%sfc_state, x2o_o, ind, sw_decomp, & c1, c2, c3, c4, OS%restore_salinity,OS%restore_temp) ! Fields that exist in both the forcing and mech_forcing types must be copied. @@ -1741,13 +1731,13 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & ! Add ice shelf fluxes if (OS%use_ice_shelf) then - call shelf_calc_flux(OS%State, OS%forces, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp) + call shelf_calc_flux(OS%sfc_state, OS%forces, 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) + ! 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%forces,OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%sfc_state, time_step, OS%berg_area_threshold) !endif ! Indicate that there are new unused fluxes. @@ -1757,17 +1747,17 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & OS%flux_tmp%C_p = OS%fluxes%C_p ! Import fluxes from coupler to ocean. Also, perform do SST and SSS restoring, if needed. call ocn_import(OS%forces, OS%flux_tmp, OS%Time, OS%grid, OS%forcing_CSp, & - OS%state, x2o_o, ind, sw_decomp, c1, c2, c3, c4, & + OS%sfc_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%forces, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) + call shelf_calc_flux(OS%sfc_state, OS%forces, 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%forces, 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) + ! call add_berg_flux_to_shelf(OS%grid, OS%forces, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%sfc_state, time_step, OS%berg_area_threshold) !endif ! Accumulate the forcing over time steps @@ -1784,52 +1774,38 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) 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) + call finish_MOM_initialization(OS%Time, OS%dirs, OS%MSp, OS%MOM_CSp, OS%fluxes, & + OS%restart_CSp) endif - call disable_averaging(OS%MOM_CSp%diag) + call disable_averaging(OS%diag) Master_time = OS%Time ; Time1 = OS%Time - if(OS%MOM_Csp%offline_tracer_mode) then - call step_offline(OS%forces, OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + if(OS%offline_tracer_mode) then + call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MSp, OS%MOM_CSp) else - call step_MOM(OS%forces, OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, time_step, OS%MSp, 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 enable_averaging(time_step, OS%Time, OS%diag) call mech_forcing_diags(OS%forces, OS%fluxes, time_step, OS%grid, & - OS%MOM_CSp%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%MOM_CSp%diag) + OS%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%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 + call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) + call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & + OS%grid, OS%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%diag) endif ! Translate state into Ocean. -! call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & +! call convert_state_to_ocean_type(OS%sfc_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) + call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) call callTree_leave("update_ocean_model()") end subroutine update_ocean_model @@ -1937,9 +1913,8 @@ subroutine ocn_import(forces, fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp, 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(forces%p_surf,isd,ied,jsd,jed) - call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed) - call safe_alloc_ptr(forces%p_surf_SSH,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%salt_flux,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) @@ -1958,6 +1933,9 @@ subroutine ocn_import(forces, fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp, fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) enddo; enddo + call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed) + call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed) + if (CS%rigid_sea_ice) then call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,JsdB,JedB) @@ -2161,9 +2139,9 @@ subroutine ocn_import(forces, fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp, endif if (CS%use_limited_P_SSH) then - forces%p_surf_SSH(i,j) = forces%p_surf(i,j) + forces%p_surf_SSH => forces%p_surf else - forces%p_surf_SSH(i,j) = forces%p_surf_full(i,j) + forces%p_surf_SSH => forces%p_surf_full endif endif @@ -2462,10 +2440,10 @@ 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. - call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag, end_diag_manager=.true.) + call diag_mediator_end(Time, Ocean_state%diag, end_diag_manager=.true.) ! print time stats call MOM_infra_end - call MOM_end(Ocean_state%MOM_CSp) + call MOM_end(Ocean_state%MSp, Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 5e727ed250..62a0dc12bf 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -27,10 +27,10 @@ program MOM_main use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_diag_mediator, only : enable_averaging, disable_averaging, diag_mediator_end - use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end + use MOM_diag_mediator, only : diag_ctrl, diag_mediator_close_registration use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : calculate_surface_state, finish_MOM_initialization - use MOM, only : step_offline + use MOM, only : MOM_state_type, step_offline use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_error_handler, only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -43,10 +43,8 @@ program MOM_main use MOM_io, only : file_exists, open_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE - use MOM_restart, only : save_restart + use MOM_restart, only : MOM_restart_CS, save_restart use MOM_string_functions,only : uppercase - use MOM_sum_output, only : write_energy, accumulate_net_input - use MOM_sum_output, only : MOM_sum_output_init, sum_output_CS use MOM_surface_forcing, only : set_forcing, forcing_save_restart use MOM_surface_forcing, only : surface_forcing_init, surface_forcing_CS use MOM_time_manager, only : time_type, set_date, set_time, get_date, time_type_to_real @@ -55,6 +53,7 @@ program MOM_main use MOM_time_manager, only : increment_date, set_calendar_type, month_name use MOM_time_manager, only : JULIAN, GREGORIAN, NOLEAP, THIRTY_DAY_MONTHS use MOM_time_manager, only : NO_CALENDAR + use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use MOM_write_cputime, only : write_cputime, MOM_write_cputime_init @@ -91,11 +90,11 @@ program MOM_main ! This is .true. if incremental restart files may be saved. logical :: permit_incr_restart = .true. - integer :: n + integer :: ns ! nmax is the number of iterations after which to stop so that the ! simulation does not exceed its CPU time limit. nmax is determined by - ! evaluating the CPU time used between successive calls to write_energy. + ! evaluating the CPU time used between successive calls to write_cputime. ! Initially it is set to be very large. integer :: nmax=2000000000; @@ -113,7 +112,6 @@ program MOM_main type(time_type) :: Start_time ! The start time of the simulation. type(time_type) :: segment_start_time ! The start time of this run segment. type(time_type) :: Time_end ! End time for the segment or experiment. - type(time_type) :: write_energy_time ! The next time to write to the energy file. type(time_type) :: restart_time ! The next time to write restart files. type(time_type) :: Time_step_ocean ! A time_type version of time_step. @@ -138,9 +136,8 @@ program MOM_main real :: Time_unit ! The time unit in seconds for the following input fields. type(time_type) :: restint ! The time between saves of the restart file. type(time_type) :: daymax ! The final day of the simulation. - type(time_type) :: energysavedays ! The interval between writing the energies - ! and other integral quantities of the run. + integer :: CPU_steps ! The number of steps between writing CPU time. integer :: date_init(6)=0 ! The start date of the whole simulation. integer :: date(6)=-1 ! Possibly the start date of this run segment. integer :: years=0, months=0, days=0 ! These may determine the segment run @@ -160,6 +157,7 @@ program MOM_main logical :: unit_in_use integer :: initClock, mainClock, termClock + logical :: debug ! If true, write verbose checksums for debugging purposes. logical :: offline_tracer_mode ! If false, use the model in prognostic mode where ! the barotropic and baroclinic dynamics, thermodynamics, ! etc. are stepped forward integrated in time. @@ -169,10 +167,18 @@ program MOM_main ! a previous integration of the prognostic model type(MOM_control_struct), pointer :: MOM_CSp => NULL() + type(MOM_state_type), pointer :: MSp => NULL() + !> A pointer to the tracer flow control structure. + type(tracer_flow_control_CS), pointer :: & + tracer_flow_CSp => NULL() !< A pointer to the tracer flow control structure type(surface_forcing_CS), pointer :: surface_forcing_CSp => NULL() - type(sum_output_CS), pointer :: sum_output_CSp => NULL() type(write_cputime_CS), pointer :: write_CPU_CSp => NULL() type(ice_shelf_CS), pointer :: ice_shelf_CSp => NULL() + type(MOM_restart_CS), pointer :: & + restart_CSp => NULL() !< A pointer to the restart control structure + !! that will be used for MOM restart files. + type(diag_ctrl), pointer :: & + diag => NULL() !< A pointer to the diagnostic regulatory structure !----------------------------------------------------------------------- character(len=4), parameter :: vers_num = 'v2.0' @@ -278,25 +284,29 @@ program MOM_main ! In this case, the segment starts at a time fixed by ocean_solo.res segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) Time = segment_start_time - ! Note the not before CS%d - call initialize_MOM(Time, param_file, dirs, MOM_CSp, segment_start_time, offline_tracer_mode = offline_tracer_mode) + call initialize_MOM(Time, Start_time, param_file, dirs, MSp, MOM_CSp, restart_CSp, & + segment_start_time, offline_tracer_mode=offline_tracer_mode, & + diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp) else ! In this case, the segment starts at a time read from the MOM restart file ! or left as Start_time by MOM_initialize. Time = Start_time - call initialize_MOM(Time, param_file, dirs, MOM_CSp, offline_tracer_mode=offline_tracer_mode) + call initialize_MOM(Time, Start_time, param_file, dirs, MSp, MOM_CSp, restart_CSp, & + offline_tracer_mode=offline_tracer_mode, diag_ptr=diag, & + tracer_flow_CSp=tracer_flow_CSp) endif - fluxes%C_p = MOM_CSp%tv%C_p ! Copy the heat capacity for consistency. + fluxes%C_p = MSp%tv%C_p ! Copy the heat capacity for consistency. Master_Time = Time - grid => MOM_CSp%G - GV => MOM_CSp%GV - call calculate_surface_state(sfc_state, MOM_CSp%u, MOM_CSp%v, MOM_CSp%h, & - MOM_CSp%ave_ssh, grid, GV, MOM_CSp) + grid => MSp%G + GV => MSp%GV + call callTree_waypoint("done initialize_MOM") + call calculate_surface_state(sfc_state, MSp%u, MSp%v, MSp%h, & + MSp%ave_ssh, grid, GV, MSp, MOM_CSp) - call surface_forcing_init(Time, grid, param_file, MOM_CSp%diag, & - surface_forcing_CSp, MOM_CSp%tracer_flow_CSp) + call surface_forcing_init(Time, grid, param_file, diag, & + surface_forcing_CSp, tracer_flow_CSp) call callTree_waypoint("done surface_forcing_init") call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & @@ -305,15 +315,9 @@ program MOM_main ! These arrays are not initialized in most solo cases, but are needed ! when using an ice shelf call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & - MOM_CSp%diag, forces, fluxes) + diag, forces, fluxes) endif - call MOM_sum_output_init(grid, param_file, dirs%output_directory, & - MOM_CSp%ntrunc, Start_time, sum_output_CSp) - call MOM_write_cputime_init(param_file, dirs%output_directory, Start_time, & - write_CPU_CSp) - call callTree_waypoint("done MOM_sum_output_init") - segment_start_time = Time elapsed_time = 0.0 @@ -375,16 +379,23 @@ program MOM_main "of TIMEUNIT. Use 0 (the default) to not save \n"//& "incremental restart files at all.", default=set_time(0), & timeunit=Time_unit) - call get_param(param_file, mod_name, "ENERGYSAVEDAYS", energysavedays, & - "The interval in units of TIMEUNIT between saves of the \n"//& - "energies of the run and other globally summed diagnostics.", & - default=set_time(int(time_step+0.5)), timeunit=Time_unit) + call get_param(param_file, mod_name, "WRITE_CPU_STEPS", cpu_steps, & + "The number of coupled timesteps between writing the cpu \n"//& + "time. If this is not positive, do not check cpu time, and \n"//& + "the segment run-length can not be set via an elapsed CPU time.", & + default=1000) + call get_param(param_file, "MOM", "DEBUG", debug, & + "If true, write out verbose debugging data.", default=.false.) call log_param(param_file, mod_name, "ELAPSED TIME AS MASTER", elapsed_time_master) + if (cpu_steps > 0) & + call MOM_write_cputime_init(param_file, dirs%output_directory, Start_time, & + write_CPU_CSp) + ! Close the param_file. No further parsing of input is possible after this. call close_param_file(param_file) - call diag_mediator_close_registration(MOM_CSp%diag) + call diag_mediator_close_registration(diag) ! Write out a time stamp file. if (calendar_type /= NO_CALENDAR) then @@ -399,13 +410,7 @@ program MOM_main call close_file(unit) endif -! This has been moved inside the loop to be applied when n=1. -! call write_energy(MOM_CSp%u, MOM_CSp%v, MOM_CSp%h, & -! MOM_CSp%tv, Time, 0, grid, GV, sum_output_CSp, MOM_CSp%tracer_flow_CSp) - call write_cputime(Time, 0, nmax, write_CPU_CSp) - - write_energy_time = Start_time + energysavedays * & - (1 + (Time - Start_time) / energysavedays) + if (cpu_steps > 0) call write_cputime(Time, 0, nmax, write_CPU_CSp) if (((.not.BTEST(Restart_control,1)) .and. (.not.BTEST(Restart_control,0))) & .or. (Restart_control < 0)) permit_incr_restart = .false. @@ -424,16 +429,16 @@ program MOM_main call cpu_clock_begin(mainClock) !begin main loop - n = 1 - do while ((n < nmax) .and. (Time < Time_end)) - call callTree_enter("Main loop, MOM_driver.F90",n) + ns = 1 + do while ((ns < nmax) .and. (Time < Time_end)) + call callTree_enter("Main loop, MOM_driver.F90",ns) ! Set the forcing for the next steps. if (.not. offline_tracer_mode) then call set_forcing(sfc_state, forces, fluxes, Time, Time_step_ocean, grid, & surface_forcing_CSp) endif - if (MOM_CSp%debug) then + if (debug) then call MOM_mech_forcing_chksum("After set forcing", forces, grid, haloshift=0) call MOM_forcing_chksum("After set forcing", fluxes, grid, haloshift=0) endif @@ -447,20 +452,16 @@ program MOM_main fluxes%fluxes_used = .false. fluxes%dt_buoy_accum = time_step - if (n==1) then - call finish_MOM_initialization(Time, dirs, MOM_CSp, fluxes) - - call write_energy(MOM_CSp%u, MOM_CSp%v, MOM_CSp%h, MOM_CSp%tv, & - Time, 0, grid, GV, sum_output_CSp, MOM_CSp%tracer_flow_CSp, & - MOM_CSp%OBC) + if (ns==1) then + call finish_MOM_initialization(Time, dirs, MSp, MOM_CSp, fluxes, restart_CSp) endif ! This call steps the model over a time time_step. Time1 = Master_Time ; Time = Master_Time if (offline_tracer_mode) then - call step_offline(forces, fluxes, sfc_state, Time1, time_step, MOM_CSp) + call step_offline(forces, fluxes, sfc_state, Time1, time_step, MSp, MOM_CSp) else - call step_MOM(forces, fluxes, sfc_state, Time1, time_step, MOM_CSp) + call step_MOM(forces, fluxes, sfc_state, Time1, time_step, MSp, MOM_CSp) endif ! Time = Time + Time_step_ocean @@ -469,7 +470,7 @@ program MOM_main if (elapsed_time > 2e9) then ! This is here to ensure that the conversion from a real to an integer ! can be accurately represented in long runs (longer than ~63 years). - ! It will also ensure that elapsed time does not loose resolution of order + ! It will also ensure that elapsed time does not lose resolution of order ! the timetype's resolution, provided that the timestep and tick are ! larger than 10-5 seconds. If a clock with a finer resolution is used, ! a smaller value would be required. @@ -483,40 +484,33 @@ program MOM_main endif Time = Master_Time - call enable_averaging(time_step, Time, MOM_CSp%diag) - call mech_forcing_diags(forces, fluxes, time_step, grid, MOM_CSp%diag, & + if (cpu_steps > 0) then ; if (MOD(ns, cpu_steps) == 0) then + call write_cputime(Time, ns+ntstep-1, nmax, write_CPU_CSp) + endif ; endif + + call enable_averaging(time_step, Time, diag) + call mech_forcing_diags(forces, fluxes, time_step, grid, diag, & surface_forcing_CSp%handles) - call disable_averaging(MOM_CSp%diag) + call disable_averaging(diag) if (.not. offline_tracer_mode) then if (fluxes%fluxes_used) then - call enable_averaging(fluxes%dt_buoy_accum, Time, MOM_CSp%diag) + call enable_averaging(fluxes%dt_buoy_accum, Time, diag) call forcing_diagnostics(fluxes, sfc_state, fluxes%dt_buoy_accum, grid, & - MOM_CSp%diag, surface_forcing_CSp%handles) - call accumulate_net_input(fluxes, sfc_state, fluxes%dt_buoy_accum, grid, sum_output_CSp) - call disable_averaging(MOM_CSp%diag) + diag, surface_forcing_CSp%handles) + call disable_averaging(diag) else call MOM_error(FATAL, "The solo MOM_driver is not yet set up to handle "//& "thermodynamic time steps that are longer than the coupling timestep.") endif endif -! See if it is time to write out the energy. - if ((Time + (Time_step_ocean/2) > write_energy_time) .and. & - (MOM_CSp%t_dyn_rel_adv == 0.0)) then - call write_energy(MOM_CSp%u, MOM_CSp%v, MOM_CSp%h, & - MOM_CSp%tv, Time, n+ntstep-1, grid, GV, sum_output_CSp, & - MOM_CSp%tracer_flow_CSp) - call write_cputime(Time, n+ntstep-1, nmax, write_CPU_CSp) - write_energy_time = write_energy_time + energysavedays - endif - ! See if it is time to write out a restart file - timestamped or not. if ((permit_incr_restart) .and. (fluxes%fluxes_used) .and. & (Time + (Time_step_ocean/2) > restart_time)) then if (BTEST(Restart_control,1)) then call save_restart(dirs%restart_output_dir, Time, grid, & - MOM_CSp%restart_CSp, .true., GV=GV) + restart_CSp, .true., GV=GV) call forcing_save_restart(surface_forcing_CSp, grid, Time, & dirs%restart_output_dir, .true.) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & @@ -524,7 +518,7 @@ program MOM_main endif if (BTEST(Restart_control,0)) then call save_restart(dirs%restart_output_dir, Time, grid, & - MOM_CSp%restart_CSp, GV=GV) + restart_CSp, GV=GV) call forcing_save_restart(surface_forcing_CSp, grid, Time, & dirs%restart_output_dir) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & @@ -533,14 +527,14 @@ program MOM_main restart_time = restart_time + restint endif - n = n + ntstep + ns = ns + ntstep call callTree_leave("Main loop") enddo call cpu_clock_end(mainClock) call cpu_clock_begin(termClock) if (Restart_control>=0) then - if (MOM_CSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& + if (MSp%t_dyn_rel_adv > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& "with inconsistent dynamics and advective times. Additional restart fields "//& "that have not been coded yet would be required for reproducibility.") if (.not.fluxes%fluxes_used .and. .not.offline_tracer_mode) call MOM_error(FATAL, & @@ -548,7 +542,7 @@ program MOM_main "For conservation, the ocean restart files can only be "//& "created after the buoyancy forcing is applied.") - call save_restart(dirs%restart_output_dir, Time, grid, MOM_CSp%restart_CSp, GV=GV) + call save_restart(dirs%restart_output_dir, Time, grid, restart_CSp, GV=GV) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & dirs%restart_output_dir) ! Write ocean solo restart file. @@ -582,12 +576,12 @@ program MOM_main endif call callTree_waypoint("End MOM_main") - call diag_mediator_end(Time, MOM_CSp%diag, end_diag_manager=.true.) + call diag_mediator_end(Time, diag, end_diag_manager=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end - call MOM_end(MOM_CSp) + call MOM_end(MSp, MOM_CSp) if (use_ice_shelf) call ice_shelf_end(ice_shelf_CSp) end program MOM_main diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3244416468..0b499e75cf 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -113,6 +113,8 @@ module MOM use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_init use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS use MOM_sponge, only : init_sponge_diags, sponge_CS +use MOM_sum_output, only : write_energy, accumulate_net_input +use MOM_sum_output, only : MOM_sum_output_init, sum_output_CS use MOM_ALE_sponge, only : init_ALE_sponge_diags, ALE_sponge_CS use MOM_thickness_diffuse, only : thickness_diffuse, thickness_diffuse_init, thickness_diffuse_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -162,8 +164,8 @@ module MOM integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = -1 end type transport_diag_IDs -!> Control structure for this module -type, public :: MOM_control_struct +!> Structure describing the state of the ocean. +type, public :: MOM_state_type real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: & h, & !< layer thickness (m or kg/m2 (H)) T, & !< potential temperature (degrees C) @@ -178,15 +180,21 @@ module MOM vhtr !< accumulated meridional thickness fluxes to advect tracers (m3 or kg) real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & ave_ssh !< time-averaged (ave over baroclinic time steps) sea surface height (meter) - real, pointer, dimension(:,:) :: Hml => NULL() !< active mixed layer depth, in m - real, pointer, dimension(:,:,:) :: & - u_prev => NULL(), & !< previous value of u stored for diagnostics - v_prev => NULL() !< previous value of v stored for diagnostics type(ocean_grid_type) :: G !< structure containing metrics and grid info type(verticalGrid_type), pointer :: GV => NULL() !< structure containing vertical grid info type(thermo_var_ptrs) :: tv !< structure containing pointers to available !! thermodynamic fields + real :: t_dyn_rel_adv !< The time of the dynamics relative to tracer + !! advection and lateral mixing (in seconds), or + !! equivalently the elapsed time since advectively + !! updating the tracers. t_dyn_rel_adv is invariably + !! positive and may span multiple coupling timesteps. +end type MOM_state_type + + +!> Control structure for this module +type, public :: MOM_control_struct ; private type(diag_ctrl) :: diag !< structure to regulate diagnostic output timing type(vertvisc_type) :: visc !< structure containing vertical viscosities, !! bottom drag viscosities, and related fields @@ -196,6 +204,10 @@ module MOM !! for derived diagnostics (e.g., energy budgets) type(cont_diag_ptrs) :: CDp !< structure containing pointers continuity equation !! terms, for derived diagnostics (e.g., energy budgets) + real, pointer, dimension(:,:) :: Hml => NULL() !< active mixed layer depth, in m + real, pointer, dimension(:,:,:) :: & + u_prev => NULL(), & !< previous value of u stored for diagnostics + v_prev => NULL() !< previous value of v stored for diagnostics logical :: split !< If true, use the split time stepping scheme. logical :: use_RK2 !< If true, use RK2 instead of RK3 in unsplit mode @@ -233,11 +245,6 @@ module MOM real :: dt_therm !< thermodynamics time step (seconds) logical :: thermo_spans_coupling !< If true, thermodynamic and tracer time !! steps can span multiple coupled time steps. - real :: t_dyn_rel_adv !< The time of the dynamics relative to tracer - !! advection and lateral mixing (in seconds), or - !! equivalently the elapsed time since advectively - !! updating the tracers. t_dyn_rel_adv is invariably - !! positive and may span multiple coupling timesteps. real :: t_dyn_rel_thermo !< The time of the dynamics relative to diabatic !! processes and remapping (in seconds). t_dyn_rel_thermo !! can be negative or positive depending on whether @@ -274,6 +281,10 @@ module MOM character(len=120) :: IC_file !< A file into which the initial conditions are !! written in a new run if SAVE_INITIAL_CONDS is true. + integer :: nstep_tot = 0 !< The total number of dynamic timesteps taken + !! so far in this run segment + logical :: count_calls = .false. !< If true, count the calls to step_MOM, rather than the + !! number of dynamics steps in nstep_tot integer :: ntrunc !< number u,v truncations since last call to write_energy logical :: check_bad_surface_vals !< If true, scan surface state for ridiculous values. real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message @@ -315,13 +326,13 @@ module MOM type(tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL() type(diagnostics_CS), pointer :: diagnostics_CSp => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() - type(MOM_restart_CS), pointer :: restart_CSp => NULL() type(update_OBC_CS), pointer :: update_OBC_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() type(sponge_CS), pointer :: sponge_CSp => NULL() type(ALE_sponge_CS), pointer :: ALE_sponge_CSp => NULL() type(ALE_CS), pointer :: ALE_CSp => NULL() type(offline_transport_CS), pointer :: offline_CSp => NULL() + type(sum_output_CS), pointer :: sum_output_CSp => NULL() ! These are used for group halo updates. type(group_pass_type) :: pass_tau_ustar_psurf @@ -368,12 +379,13 @@ module MOM !! The action of lateral processes on tracers occur in calls to !! advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping !! occur inside of diabatic. -subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) +subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, MS, CS) type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(forcing), intent(inout) :: fluxes !< pointers to forcing fields type(surface), intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_interval !< time interval covered by this run segment, in s. + type(MOM_state_type), pointer :: MS !< structure describing the MOM state type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM ! local @@ -411,7 +423,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) logical :: do_calc_bbl ! If true, calculate the boundary layer properties. logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans ! multiple dynamic timesteps. - real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: & + real, dimension(SZI_(MS%G),SZJ_(MS%G)) :: & eta_av, & ! average sea surface height or column mass over a timestep (meter or kg/m2) ssh ! sea surface height based on eta_av (meter or kg/m2) @@ -422,9 +434,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! Store the layer thicknesses, temperature, and salinity before any changes by the dynamics. ! This is necessary for remapped mass transport diagnostics - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: h_pre_dyn - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: T_pre_dyn - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: S_pre_dyn + real, dimension(SZI_(MS%G),SZJ_(MS%G),SZK_(MS%G)) :: h_pre_dyn + real, dimension(SZI_(MS%G),SZJ_(MS%G),SZK_(MS%G)) :: T_pre_dyn + real, dimension(SZI_(MS%G),SZJ_(MS%G),SZK_(MS%G)) :: S_pre_dyn real :: tot_wt_ssh, Itot_wt_ssh type(time_type) :: Time_local, end_time_thermo @@ -432,19 +444,19 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! These are used for group halo passes. logical :: do_pass_Ray, do_pass_kv_bbl_thick - G => CS%G ; GV => CS%GV ; IDs => CS%IDs + G => MS%G ; GV => MS%GV ; IDs => CS%IDs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke 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 - u => CS%u ; v => CS%v ; h => CS%h + u => MS%u ; v => MS%v ; h => MS%h call cpu_clock_begin(id_clock_ocean) call cpu_clock_begin(id_clock_other) if (CS%debug) then - call MOM_state_chksum("Beginning of step_MOM ", u, v, h, CS%uh, CS%vh, G, GV) - call hchksum(CS%h,"CS%h beginning of step_MOM",G%HI, scale=GV%H_to_m) + call MOM_state_chksum("Beginning of step_MOM ", u, v, h, MS%uh, MS%vh, G, GV) + call hchksum(MS%h,"MS%h beginning of step_MOM",G%HI, scale=GV%H_to_m) endif showCallTree = callTree_showQuery() @@ -507,15 +519,15 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) if (.not.CS%adiabatic .AND. CS%use_ALE_algorithm ) then if (CS%use_temperature) then - call create_group_pass(CS%pass_T_S_h, CS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) - call create_group_pass(CS%pass_T_S_h, CS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(CS%pass_T_S_h, MS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(CS%pass_T_S_h, MS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) endif call create_group_pass(CS%pass_T_S_h, h, G%Domain, To_All+Omit_Corners, halo=1) endif if ((CS%adiabatic .OR. CS%diabatic_first) .AND. CS%use_temperature) then - call create_group_pass(CS%pass_T_S, CS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) - call create_group_pass(CS%pass_T_S, CS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(CS%pass_T_S, MS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(CS%pass_T_S, MS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) endif !---------- End setup for group halo pass @@ -527,20 +539,20 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) endif call cpu_clock_end(id_clock_pass) - if (ASSOCIATED(CS%tv%frazil)) CS%tv%frazil(:,:) = 0.0 - if (ASSOCIATED(CS%tv%salt_deficit)) CS%tv%salt_deficit(:,:) = 0.0 - if (ASSOCIATED(CS%tv%TempxPmE)) CS%tv%TempxPmE(:,:) = 0.0 - if (ASSOCIATED(CS%tv%internal_heat)) CS%tv%internal_heat(:,:) = 0.0 + if (ASSOCIATED(MS%tv%frazil)) MS%tv%frazil(:,:) = 0.0 + if (ASSOCIATED(MS%tv%salt_deficit)) MS%tv%salt_deficit(:,:) = 0.0 + if (ASSOCIATED(MS%tv%TempxPmE)) MS%tv%TempxPmE(:,:) = 0.0 + if (ASSOCIATED(MS%tv%internal_heat)) MS%tv%internal_heat(:,:) = 0.0 CS%rel_time = 0.0 tot_wt_ssh = 0.0 - do j=js,je ; do i=is,ie ; CS%ave_ssh(i,j) = 0.0 ; ssh(i,j) = CS%missing; enddo ; enddo + do j=js,je ; do i=is,ie ; MS%ave_ssh(i,j) = 0.0 ; ssh(i,j) = CS%missing; enddo ; enddo if (associated(CS%VarMix)) then call enable_averaging(time_interval, Time_start+set_time(int(time_interval)), & CS%diag) - call calc_resoln_function(h, CS%tv, G, GV, CS%VarMix) + call calc_resoln_function(h, MS%tv, G, GV, CS%VarMix) call disable_averaging(CS%diag) endif @@ -561,7 +573,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) endif if (CS%debug) then - call MOM_state_chksum("Before steps ", u, v, h, CS%uh, CS%vh, G, GV) + call MOM_state_chksum("Before steps ", u, v, h, MS%uh, MS%vh, G, GV) call MOM_forcing_chksum("Before steps", fluxes, G, haloshift=0) call MOM_mech_forcing_chksum("Before steps", forces, G, haloshift=0) call check_redundant("Before steps ", u, v, G) @@ -583,7 +595,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) !=========================================================================== ! This is the first place where the diabatic processes and remapping could occur. - if (CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0)) then ! do thermodynamics. + if (CS%diabatic_first .and. (MS%t_dyn_rel_adv==0.0)) then ! do thermodynamics. if (thermo_does_span_coupling) then dtdia = dt_therm @@ -603,13 +615,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) end_time_thermo = Time_local + set_time(int(floor(dtdia-dt+0.5))) ! Apply diabatic forcing, do mixing, and regrid. - call step_MOM_thermo(CS, G, GV, u, v, h, CS%tv, fluxes, dtdia, end_time_thermo, .true.) + call step_MOM_thermo(MS, CS, G, GV, u, v, h, MS%tv, fluxes, dtdia, end_time_thermo, .true.) ! The diabatic processes are now ahead of the dynamics by dtdia. CS%t_dyn_rel_thermo = -dtdia if (showCallTree) call callTree_waypoint("finished diabatic_first (step_MOM)") - endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))" + endif ! end of block "(CS%diabatic_first .and. (MS%t_dyn_rel_adv==0.0))" !=========================================================================== ! This is the start of the dynamics stepping part of the algorithm. @@ -617,7 +629,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) call cpu_clock_begin(id_clock_dynamics) call disable_averaging(CS%diag) - if ((CS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then + if ((MS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then if (thermo_does_span_coupling) then dtth = dt_therm else @@ -627,8 +639,8 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) call enable_averaging(dtth,Time_local+set_time(int(floor(dtth-dt+0.5))), CS%diag) call cpu_clock_begin(id_clock_thick_diff) if (associated(CS%VarMix)) & - call calc_slope_functions(h, CS%tv, dt, G, GV, CS%VarMix) - call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dtth, G, GV, & + call calc_slope_functions(h, MS%tv, dt, G, GV, CS%VarMix) + call thickness_diffuse(h, MS%uhtr, MS%vhtr, MS%tv, dtth, G, GV, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) call cpu_clock_end(id_clock_thick_diff) call pass_var(h, G%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil)) @@ -643,14 +655,14 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! The bottom boundary layer properties are out-of-date and need to be ! recalculated. This always occurs at the start of a coupling time ! step because the externally prescribed stresses may have changed. - do_calc_bbl = ((CS%t_dyn_rel_adv == 0.0) .or. (n==1)) + do_calc_bbl = ((MS%t_dyn_rel_adv == 0.0) .or. (n==1)) if (do_calc_bbl) then ! Calculate the BBL properties and store them inside visc (u,h). call cpu_clock_begin(id_clock_BBL_visc) - bbl_time_int = max(dt, min(dt_therm - CS%t_dyn_rel_adv, dt*(1+n_max-n)) ) + bbl_time_int = max(dt, min(dt_therm - MS%t_dyn_rel_adv, dt*(1+n_max-n)) ) call enable_averaging(bbl_time_int, & Time_local+set_time(int(bbl_time_int-dt+0.5)), CS%diag) - call set_viscous_BBL(u, v, h, CS%tv, CS%visc, G, GV, CS%set_visc_CSp) + call set_viscous_BBL(u, v, h, MS%tv, CS%visc, G, GV, CS%set_visc_CSp) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)") @@ -696,8 +708,8 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) if (transport_remap_grid_needed(CS%transport_IDs)) then do k=1,nz ; do j=jsd,jed ; do i=isd,ied h_pre_dyn(i,j,k) = h(i,j,k) - if (associated(CS%tv%T)) T_pre_dyn(i,j,k) = CS%tv%T(i,j,k) - if (associated(CS%tv%S)) S_pre_dyn(i,j,k) = CS%tv%S(i,j,k) + if (associated(MS%tv%T)) T_pre_dyn(i,j,k) = MS%tv%T(i,j,k) + if (associated(MS%tv%S)) S_pre_dyn(i,j,k) = MS%tv%S(i,j,k) enddo ; enddo ; enddo endif @@ -721,9 +733,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) endif mass_src_time = CS%t_dyn_rel_thermo - call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, & + call step_MOM_dyn_split_RK2(u, v, h, MS%tv, CS%visc, & Time_local, dt, forces, CS%p_surf_begin, CS%p_surf_end, & - mass_src_time, dt_therm, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & + mass_src_time, dt_therm, MS%uh, MS%vh, MS%uhtr, MS%vhtr, & eta_av, G, GV, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, CS%MEKE) if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)") @@ -736,12 +748,12 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! useful for debugging purposes. if (CS%use_RK2) then - call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & - CS%p_surf_begin, CS%p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & + call step_MOM_dyn_unsplit_RK2(u, v, h, MS%tv, CS%visc, Time_local, dt, forces, & + CS%p_surf_begin, CS%p_surf_end, MS%uh, MS%vh, MS%uhtr, MS%vhtr, & eta_av, G, GV, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE) else - call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & - CS%p_surf_begin, CS%p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & + call step_MOM_dyn_unsplit(u, v, h, MS%tv, CS%visc, Time_local, dt, forces, & + CS%p_surf_begin, CS%p_surf_end, MS%uh, MS%vh, MS%uhtr, MS%vhtr, & eta_av, G, GV, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE) endif if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)") @@ -755,8 +767,8 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) if (CS%debug) call hchksum(h,"Pre-thickness_diffuse h", G%HI, haloshift=0, scale=GV%H_to_m) if (associated(CS%VarMix)) & - call calc_slope_functions(h, CS%tv, dt, G, GV, CS%VarMix) - call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, & + call calc_slope_functions(h, MS%tv, dt, G, GV, CS%VarMix) + call thickness_diffuse(h, MS%uhtr, MS%vhtr, MS%tv, dt, G, GV, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_m) @@ -770,17 +782,17 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) if (CS%debug) then call hchksum(h,"Pre-mixedlayer_restrat h", G%HI, haloshift=1, scale=GV%H_to_m) call uvchksum("Pre-mixedlayer_restrat uhtr", & - CS%uhtr, CS%vhtr, G%HI, haloshift=0, scale=GV%H_to_m) + MS%uhtr, MS%vhtr, G%HI, haloshift=0, scale=GV%H_to_m) endif call cpu_clock_begin(id_clock_ml_restrat) - call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, & + call mixedlayer_restrat(h, MS%uhtr, MS%vhtr, MS%tv, forces, dt, CS%visc%MLD, & CS%VarMix, G, GV, CS%mixedlayer_restrat_CSp) call cpu_clock_end(id_clock_ml_restrat) call pass_var(h, G%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil)) if (CS%debug) then call hchksum(h,"Post-mixedlayer_restrat h", G%HI, haloshift=1, scale=GV%H_to_m) call uvchksum("Post-mixedlayer_restrat [uv]htr", & - CS%uhtr, CS%vhtr, G%HI, haloshift=0, scale=GV%H_to_m) + MS%uhtr, MS%vhtr, G%HI, haloshift=0, scale=GV%H_to_m) endif endif @@ -789,11 +801,11 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) call diag_update_remap_grids(CS%diag) if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, & - CS%visc, dt, G, GV, CS%MEKE_CSp, CS%uhtr, CS%vhtr) + CS%visc, dt, G, GV, CS%MEKE_CSp, MS%uhtr, MS%vhtr) call disable_averaging(CS%diag) ! Advance the dynamics time by dt. - CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt + MS%t_dyn_rel_adv = MS%t_dyn_rel_adv + dt CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt CS%t_dyn_rel_diag = CS%t_dyn_rel_diag + dt @@ -803,19 +815,19 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! This is the start of the tracer advection part of the algorithm. if (thermo_does_span_coupling) then - do_advection = (CS%t_dyn_rel_adv + 0.5*dt > dt_therm) + do_advection = (MS%t_dyn_rel_adv + 0.5*dt > dt_therm) else do_advection = ((MOD(n,ntstep) == 0) .or. (n==n_max)) endif if (do_advection) then ! Do advective transport and lateral tracer mixing. - call step_MOM_tracer_dyn(CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & - CS%tv, Time_local) + call step_MOM_tracer_dyn(MS, CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & + MS%tv, Time_local) endif !=========================================================================== ! This is the second place where the diabatic processes and remapping could occur. - if (CS%t_dyn_rel_adv == 0.0) then + if (MS%t_dyn_rel_adv == 0.0) then if (.not.CS%diabatic_first) then dtdia = CS%t_dyn_rel_thermo if (thermo_does_span_coupling .and. (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then @@ -823,7 +835,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) "before call to diabatic.") endif ! Apply diabatic forcing, do mixing, and regrid. - call step_MOM_thermo(CS, G, GV, u, v, h, CS%tv, fluxes, dtdia, Time_local, .false.) + call step_MOM_thermo(MS, CS, G, GV, u, v, h, MS%tv, fluxes, dtdia, Time_local, .false.) endif if (CS%diabatic_first .and. abs(CS%t_dyn_rel_thermo) > 1e-6*dt) call MOM_error(FATAL, & @@ -838,9 +850,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) ! Determining the time-average sea surface height is part of the algorithm. ! This may be eta_av if Boussinesq, or need to be diagnosed if not. tot_wt_ssh = tot_wt_ssh + dt - call find_eta(h, CS%tv, GV%g_Earth, G, GV, ssh, eta_av) + call find_eta(h, MS%tv, GV%g_Earth, G, GV, ssh, eta_av) do j=js,je ; do i=is,ie - CS%ave_ssh(i,j) = CS%ave_ssh(i,j) + dt*ssh(i,j) + MS%ave_ssh(i,j) = MS%ave_ssh(i,j) + dt*ssh(i,j) enddo ; enddo call cpu_clock_end(id_clock_dynamics) @@ -856,11 +868,11 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) if (IDs%id_ssh_inst > 0) call post_data(IDs%id_ssh_inst, ssh, CS%diag) call disable_averaging(CS%diag) - if (CS%t_dyn_rel_adv == 0.0) then + if (MS%t_dyn_rel_adv == 0.0) then ! Diagnostics that require the complete state to be up-to-date can be calculated. call enable_averaging(CS%t_dyn_rel_diag, Time_local, CS%diag) - call calculate_diagnostic_fields(u, v, h, CS%uh, CS%vh, CS%tv, CS%ADp, & + call calculate_diagnostic_fields(u, v, h, MS%uh, MS%vh, MS%tv, CS%ADp, & CS%CDp, fluxes, CS%t_dyn_rel_diag, G, GV, CS%diagnostics_CSp) call post_tracer_diagnostics(CS%Tracer_reg, h, CS%diag, G, GV, CS%t_dyn_rel_diag) if (showCallTree) call callTree_waypoint("finished calculate_diagnostic_fields (step_MOM)") @@ -881,33 +893,48 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS) endif call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other) + if (.not.CS%count_calls) CS%nstep_tot = CS%nstep_tot + 1 + if (showCallTree) call callTree_leave("DT cycles (step_MOM)") enddo ! complete the n loop + if (CS%count_calls) CS%nstep_tot = CS%nstep_tot + 1 + call cpu_clock_begin(id_clock_other) Itot_wt_ssh = 1.0/tot_wt_ssh do j=js,je ; do i=is,ie - CS%ave_ssh(i,j) = CS%ave_ssh(i,j)*Itot_wt_ssh - ssh(i,j) = CS%ave_ssh(i,j) + MS%ave_ssh(i,j) = MS%ave_ssh(i,j)*Itot_wt_ssh + ssh(i,j) = MS%ave_ssh(i,j) enddo ; enddo - call adjust_ssh_for_p_atm(CS%tv, G, GV, CS%ave_ssh, forces%p_surf_SSH, CS%calc_rho_for_sea_lev) + call adjust_ssh_for_p_atm(MS%tv, G, GV, MS%ave_ssh, forces%p_surf_SSH, CS%calc_rho_for_sea_lev) if (CS%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied CS%p_surf_prev(i,j) = forces%p_surf(i,j) enddo ; enddo ; endif if (showCallTree) call callTree_waypoint("calling calculate_surface_state (step_MOM)") - call calculate_surface_state(sfc_state, u, v, h, CS%ave_ssh, G, GV, CS) + call calculate_surface_state(sfc_state, u, v, h, MS%ave_ssh, G, GV, MS, CS) ! Do diagnostics that only occur at the end of a complete forcing step. call cpu_clock_begin(id_clock_diagnostics) call enable_averaging(dt*n_max, Time_local, CS%diag) - call post_surface_diagnostics(CS%sfc_IDs, G, GV, CS%diag, dt*n_max, sfc_state, CS%tv, ssh, fluxes) + call post_surface_diagnostics(CS%sfc_IDs, G, GV, CS%diag, dt*n_max, sfc_state, MS%tv, ssh, fluxes) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) + ! Accumulate the surface fluxes for assessing conservation + if (fluxes%fluxes_used) & + call accumulate_net_input(fluxes, sfc_state, fluxes%dt_buoy_accum, & + G, CS%sum_output_CSp) + + if (MS%t_dyn_rel_adv==0.0) & + call write_energy(MS%u, MS%v, MS%h, MS%tv, Time_local, & + CS%nstep_tot, G, GV, CS%sum_output_CSp, & + CS%tracer_flow_CSp, & + dt_forcing=set_time(int(floor(time_interval+0.5))) ) + call cpu_clock_end(id_clock_other) if (showCallTree) call callTree_leave("step_MOM()") @@ -918,8 +945,9 @@ end subroutine step_MOM !> step_MOM_tracer_dyn does tracer advection and lateral diffusion, bringing the !! tracers up to date with the changes in state due to the dynamics. Surface !! sources and sinks and remapping are handled via step_MOM_thermo. -subroutine step_MOM_tracer_dyn(CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & +subroutine step_MOM_tracer_dyn(MS, CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & tv, Time_local) + type(MOM_state_type), intent(inout) :: MS !< structure describing the MOM state type(MOM_control_struct), intent(inout) :: CS !< control structure type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -940,31 +968,31 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & if (CS%debug) then call cpu_clock_begin(id_clock_other) call hchksum(h,"Pre-advection h", G%HI, haloshift=1, scale=GV%H_to_m) - call uvchksum("Pre-advection uhtr", CS%uhtr, CS%vhtr, G%HI, & + call uvchksum("Pre-advection uhtr", MS%uhtr, MS%vhtr, G%HI, & haloshift=0, scale=GV%H_to_m) - if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre-advection T", G%HI, haloshift=1) - if (associated(CS%tv%S)) call hchksum(CS%tv%S, "Pre-advection S", G%HI, haloshift=1) - if (associated(CS%tv%frazil)) call hchksum(CS%tv%frazil, & + if (associated(MS%tv%T)) call hchksum(MS%tv%T, "Pre-advection T", G%HI, haloshift=1) + if (associated(MS%tv%S)) call hchksum(MS%tv%S, "Pre-advection S", G%HI, haloshift=1) + if (associated(MS%tv%frazil)) call hchksum(MS%tv%frazil, & "Pre-advection frazil", G%HI, haloshift=0) - if (associated(CS%tv%salt_deficit)) call hchksum(CS%tv%salt_deficit, & + if (associated(MS%tv%salt_deficit)) call hchksum(MS%tv%salt_deficit, & "Pre-advection salt deficit", G%HI, haloshift=0) - ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G) + ! call MOM_thermo_chksum("Pre-advection ", MS%tv, G) call cpu_clock_end(id_clock_other) endif call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer) - call enable_averaging(CS%t_dyn_rel_adv, Time_local, CS%diag) + call enable_averaging(MS%t_dyn_rel_adv, Time_local, CS%diag) - call advect_tracer(h, CS%uhtr, CS%vhtr, CS%OBC, CS%t_dyn_rel_adv, G, GV, & + call advect_tracer(h, MS%uhtr, MS%vhtr, CS%OBC, MS%t_dyn_rel_adv, G, GV, & CS%tracer_adv_CSp, CS%tracer_Reg) - call tracer_hordiff(h, CS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, G, GV, & - CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) + call tracer_hordiff(h, MS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, MS%tv) if (showCallTree) call callTree_waypoint("finished tracer advection/diffusion (step_MOM)") call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo) call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics) - call post_transport_diagnostics(G, GV, CS%uhtr, CS%vhtr, h, CS%transport_IDs, & - CS%diag, CS%t_dyn_rel_adv, CS%diag_to_Z_CSp, h_pre_dyn, T_pre_dyn, S_pre_dyn) + call post_transport_diagnostics(G, GV, MS%uhtr, MS%vhtr, h, CS%transport_IDs, & + CS%diag, MS%t_dyn_rel_adv, CS%diag_to_Z_CSp, h_pre_dyn, T_pre_dyn, S_pre_dyn) ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses ! from before the dynamics calls call diag_update_remap_grids(CS%diag) @@ -975,9 +1003,9 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, h, h_pre_dyn, T_pre_dyn, S_pre_dyn, & ! Reset the accumulated transports to 0 and record that the dynamics ! and advective times now agree. call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer) - CS%uhtr(:,:,:) = 0.0 - CS%vhtr(:,:,:) = 0.0 - CS%t_dyn_rel_adv = 0.0 + MS%uhtr(:,:,:) = 0.0 + MS%vhtr(:,:,:) = 0.0 + MS%t_dyn_rel_adv = 0.0 call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo) if (CS%diabatic_first .and. CS%use_temperature) then @@ -990,7 +1018,8 @@ end subroutine step_MOM_tracer_dyn !> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical !! remapping, via calls to diabatic (or adiabatic) and ALE_main. -subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_thermo, update_BBL) +subroutine step_MOM_thermo(MS, CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_thermo, update_BBL) + type(MOM_state_type), intent(inout) :: MS !< structure describing the MOM state type(MOM_control_struct), intent(inout) :: CS !< control structure type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure @@ -1024,8 +1053,8 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm if (CS%debug) then call uvchksum("Pre set_viscous_BBL [uv]", u, v, G%HI, haloshift=1) call hchksum(h,"Pre set_viscous_BBL h", G%HI, haloshift=1, scale=GV%H_to_m) - if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre set_viscous_BBL T", G%HI, haloshift=1) - if (associated(CS%tv%S)) call hchksum(CS%tv%S, "Pre set_viscous_BBL S", G%HI, haloshift=1) + if (associated(MS%tv%T)) call hchksum(MS%tv%T, "Pre set_viscous_BBL T", G%HI, haloshift=1) + if (associated(MS%tv%S)) call hchksum(MS%tv%S, "Pre set_viscous_BBL S", G%HI, haloshift=1) endif ! Calculate the BBL properties and store them inside visc (u,h). @@ -1033,7 +1062,7 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) - call set_viscous_BBL(u, v, h, CS%tv, CS%visc, G, GV, CS%set_visc_CSp) + call set_viscous_BBL(u, v, h, MS%tv, CS%visc, G, GV, CS%set_visc_CSp) call cpu_clock_end(id_clock_BBL_visc) if (.not.G%Domain%symmetric) then @@ -1051,10 +1080,10 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm if (CS%debug) then call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2) call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m) - call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, & + call uvchksum("Pre-diabatic [uv]h", MS%uhtr, MS%vhtr, G%HI, & haloshift=0, scale=GV%H_to_m) - ! call MOM_state_chksum("Pre-diabatic ",u, v, h, CS%uhtr, CS%vhtr, G, GV) - call MOM_thermo_chksum("Pre-diabatic ", CS%tv, G,haloshift=0) + ! call MOM_state_chksum("Pre-diabatic ",u, v, h, MS%uhtr, MS%vhtr, G, GV) + call MOM_thermo_chksum("Pre-diabatic ", MS%tv, G,haloshift=0) call check_redundant("Pre-diabatic ", u, v, G) call MOM_forcing_chksum("Pre-diabatic", fluxes, G, haloshift=0) endif @@ -1078,7 +1107,7 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm call preAle_tracer_diagnostics(CS%tracer_Reg, G, GV) if (CS%debug) then - call MOM_state_chksum("Pre-ALE ", u, v, h, CS%uh, CS%vh, G, GV) + call MOM_state_chksum("Pre-ALE ", u, v, h, MS%uh, MS%vh, G, GV) call hchksum(tv%T,"Pre-ALE T", G%HI, haloshift=1) call hchksum(tv%S,"Pre-ALE S", G%HI, haloshift=1) call check_redundant("Pre-ALE ", u, v, G) @@ -1098,7 +1127,7 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm call do_group_pass(CS%pass_uv_T_S_h, G%Domain, clock=id_clock_pass) if (CS%debug .and. CS%use_ALE_algorithm) then - call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G, GV) + call MOM_state_chksum("Post-ALE ", u, v, h, MS%uh, MS%vh, G, GV) call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1) call hchksum(tv%S, "Post-ALE S", G%HI, haloshift=1) call check_redundant("Post-ALE ", u, v, G) @@ -1115,10 +1144,10 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, Time_end_therm if (CS%debug) then call uvchksum("Post-diabatic u", u, v, G%HI, haloshift=2) call hchksum(h, "Post-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m) - call uvchksum("Post-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, & + call uvchksum("Post-diabatic [uv]h", MS%uhtr, MS%vhtr, G%HI, & haloshift=0, scale=GV%H_to_m) ! call MOM_state_chksum("Post-diabatic ", u, v, & - ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1) + ! h, MS%uhtr, MS%vhtr, G, GV, haloshift=1) if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1) if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1) if (associated(tv%frazil)) call hchksum(tv%frazil, & @@ -1158,12 +1187,13 @@ end subroutine step_MOM_thermo !! developed with ALE configurations in mind. Some work has been done in isopycnal configuration, but !! the work is very preliminary. Some more detail about this capability along with some of the subroutines !! called here can be found in tracers/MOM_offline_control.F90 -subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS) +subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, MS, CS) type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces type(forcing), intent(inout) :: fluxes !< pointers to forcing fields type(surface), intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_interval !< time interval + type(MOM_state_type), pointer :: MS !< structure describing the MOM state type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM ! Local pointers @@ -1192,12 +1222,12 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS h_end ! 2D Array for diagnostics - real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end + real, dimension(SZI_(MS%G),SZJ_(MS%G)) :: eta_pre, eta_end type(time_type) :: Time_end ! End time of a segment, as a time type ! Grid-related pointer assignments - G => CS%G - GV => CS%GV + G => MS%G + GV => MS%GV is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -1237,71 +1267,71 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if(is_root_pe()) print *, "Reading in new offline fields" ! Read in new transport and other fields ! call update_transport_from_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, & - ! CS%tv%T, CS%tv%S, fluxes, CS%use_ALE_algorithm) + ! MS%tv%T, MS%tv%S, fluxes, CS%use_ALE_algorithm) ! call update_transport_from_arrays(CS%offline_CSp) - call update_offline_fields(CS%offline_CSp, CS%h, fluxes, CS%use_ALE_algorithm) + call update_offline_fields(CS%offline_CSp, MS%h, fluxes, CS%use_ALE_algorithm) ! Apply any fluxes into the ocean - call offline_fw_fluxes_into_ocean(G, GV, CS%offline_CSp, fluxes, CS%h) + call offline_fw_fluxes_into_ocean(G, GV, CS%offline_CSp, fluxes, MS%h) if (.not.CS%diabatic_first) then call offline_advection_ale(fluxes, Time_start, time_interval, CS%offline_CSp, id_clock_ALE, & - CS%h, uhtr, vhtr, converged=adv_converged) + MS%h, uhtr, vhtr, converged=adv_converged) ! Redistribute any remaining transport - call offline_redistribute_residual(CS%offline_CSp, CS%h, uhtr, vhtr, adv_converged) + call offline_redistribute_residual(CS%offline_CSp, MS%h, uhtr, vhtr, adv_converged) ! Perform offline diffusion if requested if (.not. skip_diffusion) then if (associated(CS%VarMix)) then - call pass_var(CS%h,G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, CS%VarMix) - call calc_slope_functions(CS%h, CS%tv, REAL(dt_offline), G, GV, CS%VarMix) + call pass_var(MS%h,G%Domain) + call calc_resoln_function(MS%h, MS%tv, G, GV, CS%VarMix) + call calc_slope_functions(MS%h, MS%tv, REAL(dt_offline), G, GV, CS%VarMix) endif - call tracer_hordiff(CS%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, & - CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) + call tracer_hordiff(MS%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, MS%tv) endif endif endif ! The functions related to column physics of tracers is performed separately in ALE mode if (do_vertical) then - call offline_diabatic_ale(fluxes, Time_start, Time_end, CS%offline_CSp, CS%h, eatr, ebtr) + call offline_diabatic_ale(fluxes, Time_start, Time_end, CS%offline_CSp, MS%h, eatr, ebtr) endif ! Last thing that needs to be done is the final ALE remapping if(last_iter) then if (CS%diabatic_first) then call offline_advection_ale(fluxes, Time_start, time_interval, CS%offline_CSp, id_clock_ALE, & - CS%h, uhtr, vhtr, converged=adv_converged) + MS%h, uhtr, vhtr, converged=adv_converged) ! Redistribute any remaining transport and perform the remaining advection - call offline_redistribute_residual(CS%offline_CSp, CS%h, uhtr, vhtr, adv_converged) + call offline_redistribute_residual(CS%offline_CSp, MS%h, uhtr, vhtr, adv_converged) ! Perform offline diffusion if requested if (.not. skip_diffusion) then if (associated(CS%VarMix)) then - call pass_var(CS%h,G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, CS%VarMix) - call calc_slope_functions(CS%h, CS%tv, REAL(dt_offline), G, GV, CS%VarMix) + call pass_var(MS%h,G%Domain) + call calc_resoln_function(MS%h, MS%tv, G, GV, CS%VarMix) + call calc_slope_functions(MS%h, MS%tv, REAL(dt_offline), G, GV, CS%VarMix) endif - call tracer_hordiff(CS%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, & - CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) + call tracer_hordiff(MS%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, MS%tv) endif endif if(is_root_pe()) print *, "Last iteration of offline interval" ! Apply freshwater fluxes out of the ocean - call offline_fw_fluxes_out_ocean(G, GV, CS%offline_CSp, fluxes, CS%h) + call offline_fw_fluxes_out_ocean(G, GV, CS%offline_CSp, fluxes, MS%h) ! These diagnostic can be used to identify which grid points did not converge within ! the specified number of advection sub iterations - call post_offline_convergence_diags(CS%offline_CSp, CS%h, h_end, uhtr, vhtr) + call post_offline_convergence_diags(CS%offline_CSp, MS%h, h_end, uhtr, vhtr) ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses ! stored from the forward run call cpu_clock_begin(id_clock_ALE) - call ALE_offline_tracer_final( G, GV, CS%h, CS%tv, h_end, CS%tracer_Reg, CS%ALE_CSp) + call ALE_offline_tracer_final( G, GV, MS%h, MS%tv, h_end, CS%tracer_Reg, CS%ALE_CSp) call cpu_clock_end(id_clock_ALE) - call pass_var(CS%h, G%Domain) + call pass_var(MS%h, G%Domain) endif else ! NON-ALE MODE...NOT WELL TESTED call MOM_error(WARNING, & @@ -1313,30 +1343,30 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call MOM_error(FATAL, & "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval") endif - call update_offline_fields(CS%offline_CSp, CS%h, fluxes, CS%use_ALE_algorithm) + call update_offline_fields(CS%offline_CSp, MS%h, fluxes, CS%use_ALE_algorithm) call offline_advection_layer(fluxes, Time_start, time_interval, CS%offline_CSp, & - CS%h, eatr, ebtr, uhtr, vhtr) + MS%h, eatr, ebtr, uhtr, vhtr) ! Perform offline diffusion if requested if (.not. skip_diffusion) then call tracer_hordiff(h_end, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, & - CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) + CS%tracer_diff_CSp, CS%tracer_Reg, MS%tv) endif - CS%h = h_end + MS%h = h_end - call pass_var(CS%tv%T, G%Domain) - call pass_var(CS%tv%S, G%Domain) - call pass_var(CS%h, G%Domain) + call pass_var(MS%tv%T, G%Domain) + call pass_var(MS%tv%S, G%Domain) + call pass_var(MS%h, G%Domain) endif - call adjust_ssh_for_p_atm(CS%tv, G, GV, CS%ave_ssh, forces%p_surf_SSH, CS%calc_rho_for_sea_lev) - call calculate_surface_state(sfc_state, CS%u, CS%v, CS%h, CS%ave_ssh, G, GV, CS) + call adjust_ssh_for_p_atm(MS%tv, G, GV, MS%ave_ssh, forces%p_surf_SSH, CS%calc_rho_for_sea_lev) + call calculate_surface_state(sfc_state, MS%u, MS%v, MS%h, MS%ave_ssh, G, GV, MS, CS) call disable_averaging(CS%diag) - call pass_var(CS%tv%T,G%Domain) - call pass_var(CS%tv%S,G%Domain) - call pass_var(CS%h,G%Domain) + call pass_var(MS%tv%T,G%Domain) + call pass_var(MS%tv%S,G%Domain) + call pass_var(MS%h,G%Domain) fluxes%fluxes_used = .true. @@ -1345,16 +1375,30 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS end subroutine step_offline !> This subroutine initializes MOM. -subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mode, input_restart_file) +subroutine initialize_MOM(Time, Time_init, param_file, dirs, MS, CS, restart_CSp, & + Time_in, offline_tracer_mode, input_restart_file, diag_ptr, & + count_calls, tracer_flow_CSp) type(time_type), target, intent(inout) :: Time !< model time, set in this routine + type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(param_file_type), intent(out) :: param_file !< structure indicating paramater file to parse type(directories), intent(out) :: dirs !< structure with directory paths + type(MOM_state_type), pointer :: MS !< pointer set in this routine to structure describing the MOM state type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure + type(MOM_restart_CS), pointer :: restart_CSp !< pointer set in this routine to the + !! restart control structure that will + !! be used for MOM. type(time_type), optional, intent(in) :: Time_in !< time passed to MOM_initialize_state when !! model is not being started from a restart file - logical, optional, intent(out) :: offline_tracer_mode !< True if tracers are being run offline + logical, optional, intent(out) :: offline_tracer_mode !< True is returned if tracers are being run offline character(len=*),optional, intent(in) :: input_restart_file !< If present, name of restart file to read - + type(diag_ctrl), optional, pointer :: diag_ptr !< A pointer set in this routine to the diagnostic + !! regulatory structure + type(tracer_flow_control_CS), & + optional, pointer :: tracer_flow_CSp !< A pointer set in this routine to + !! the tracer flow control structure. + logical, optional, intent(in) :: count_calls !< If true, nstep_tot counts the number of + !! calls to step_MOM instead of the number of + !! dynamics timesteps. ! local type(ocean_grid_type), pointer :: G => NULL() ! A pointer to a structure with metrics and related type(hor_index_type) :: HI ! A hor_index_type for array extents @@ -1420,8 +1464,14 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo endif allocate(CS) + if (associated(MS)) then + call MOM_error(WARNING, "initialize_MOM called with an MOM state structure.") + return + endif + allocate(MS) + if (test_grid_copy) then ; allocate(G) - else ; G => CS%G ; endif + else ; G => MS%G ; endif CS%Time => Time @@ -1486,7 +1536,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo "to potential temperature and practical salinity before \n"//& "exchanging them with the coupler and/or reporting T&S diagnostics.\n", & default=.false.) - CS%tv%T_is_conT = use_conT_absS ; CS%tv%S_is_absS = use_conT_absS + MS%tv%T_is_conT = use_conT_absS ; MS%tv%S_is_absS = use_conT_absS call get_param(param_file, "MOM", "ADIABATIC", CS%adiabatic, & "There are no diapycnal mass fluxes if ADIABATIC is \n"//& "true. This assumes that KD = KDML = 0.0 and that \n"//& @@ -1497,19 +1547,19 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo "the gravity wave adjustment to h. This is a fragile feature and\n"//& "thus undocumented.", default=.true., do_not_log=.true. ) call get_param(param_file, "MOM", "ADVECT_TS", CS%advect_TS , & - "If True, advect temperature and salinity horizontally\n"//& + "If True, advect temperature and salinity horizontally \n"//& "If False, T/S are registered for advection.\n"//& - "This is intended only to be used in offline tracer mode.", & - "and is by default false in that case", & + "This is intended only to be used in offline tracer mode \n"//& + "and is by default false in that case.", & do_not_log = .true., default=.true. ) - if (present(offline_tracer_mode)) then ! Only read this parameter in solo mode + if (present(offline_tracer_mode)) then ! Only read this parameter in enabled modes call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", CS%offline_tracer_mode, & "If true, barotropic and baroclinic dynamics, thermodynamics\n"//& "are all bypassed with all the fields necessary to integrate\n"//& "the tracer advection and diffusion equation are read in from\n"//& "files stored from a previous integration of the prognostic model.\n"//& "NOTE: This option only used in the ocean_solo_driver.", default=.false.) - if(CS%offline_tracer_mode) then + if (CS%offline_tracer_mode) then call get_param(param_file, "MOM", "ADVECT_TS", CS%advect_TS , & "If True, advect temperature and salinity horizontally\n"//& "If False, T/S are registered for advection.\n"//& @@ -1607,7 +1657,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo endif ! This is here in case these values are used inappropriately. - CS%use_frazil = .false. ; CS%bound_salinity = .false. ; CS%tv%P_Ref = 2.0e7 + CS%use_frazil = .false. ; CS%bound_salinity = .false. ; MS%tv%P_Ref = 2.0e7 if (CS%use_temperature) then call get_param(param_file, "MOM", "FRAZIL", CS%use_frazil, & "If true, water freezes if it gets too cold, and the \n"//& @@ -1620,14 +1670,14 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo "If true, limit salinity to being positive. (The sea-ice \n"//& "model may ask for more salt than is available and \n"//& "drive the salinity negative otherwise.)", default=.false.) - call get_param(param_file, "MOM", "C_P", CS%tv%C_p, & + call get_param(param_file, "MOM", "C_P", MS%tv%C_p, & "The heat capacity of sea water, approximated as a \n"//& "constant. This is only used if ENABLE_THERMODYNAMICS is \n"//& "true. The default value is from the TEOS-10 definition \n"//& "of conservative temperature.", units="J kg-1 K-1", & default=3991.86795711963) endif - if (use_EOS) call get_param(param_file, "MOM", "P_REF", CS%tv%P_Ref, & + if (use_EOS) call get_param(param_file, "MOM", "P_REF", MS%tv%P_Ref, & "The pressure that is used for calculating the coordinate \n"//& "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) \n"//& "This is only used if USE_EOS and ENABLE_THERMODYNAMICS \n"//& @@ -1702,7 +1752,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo ((dirs%input_filename(1:1)=='n') .and. (LEN_TRIM(dirs%input_filename)==1)))) ! If the restart file type had been initialized, this could become: ! write_geom_files = ((write_geom==2) .or. & -! ((write_geom==1) .and. is_new_run(CS%restart_CSp))) +! ((write_geom==1) .and. is_new_run(restart_CSp))) ! Check for inconsistent parameter settings. if (CS%use_ALE_algorithm .and. CS%bulkmixedlayer) call MOM_error(FATAL, & @@ -1760,8 +1810,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call create_dyn_horgrid(dG, HI, bathymetry_at_vel=bathy_at_vel) call clone_MOM_domain(G%Domain, dG%Domain) - call verticalGridInit( param_file, CS%GV ) - GV => CS%GV + call verticalGridInit( param_file, MS%GV ) + GV => MS%GV ! dG%g_Earth = GV%g_Earth ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. @@ -1784,25 +1834,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo is = dG%isc ; ie = dG%iec ; js = dG%jsc ; je = dG%jec ; nz = GV%ke isd = dG%isd ; ied = dG%ied ; jsd = dG%jsd ; jed = dG%jed IsdB = dG%IsdB ; IedB = dG%IedB ; JsdB = dG%JsdB ; JedB = dG%JedB - ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz)) ; CS%u(:,:,:) = 0.0 - ALLOC_(CS%v(isd:ied,JsdB:JedB,nz)) ; CS%v(:,:,:) = 0.0 - ALLOC_(CS%h(isd:ied,jsd:jed,nz)) ; CS%h(:,:,:) = GV%Angstrom - ALLOC_(CS%uh(IsdB:IedB,jsd:jed,nz)) ; CS%uh(:,:,:) = 0.0 - ALLOC_(CS%vh(isd:ied,JsdB:JedB,nz)) ; CS%vh(:,:,:) = 0.0 + ALLOC_(MS%u(IsdB:IedB,jsd:jed,nz)) ; MS%u(:,:,:) = 0.0 + ALLOC_(MS%v(isd:ied,JsdB:JedB,nz)) ; MS%v(:,:,:) = 0.0 + ALLOC_(MS%h(isd:ied,jsd:jed,nz)) ; MS%h(:,:,:) = GV%Angstrom + ALLOC_(MS%uh(IsdB:IedB,jsd:jed,nz)) ; MS%uh(:,:,:) = 0.0 + ALLOC_(MS%vh(isd:ied,JsdB:JedB,nz)) ; MS%vh(:,:,:) = 0.0 if (CS%use_temperature) then - ALLOC_(CS%T(isd:ied,jsd:jed,nz)) ; CS%T(:,:,:) = 0.0 - ALLOC_(CS%S(isd:ied,jsd:jed,nz)) ; CS%S(:,:,:) = 0.0 - CS%tv%T => CS%T ; CS%tv%S => CS%S - if (CS%tv%T_is_conT) then + ALLOC_(MS%T(isd:ied,jsd:jed,nz)) ; MS%T(:,:,:) = 0.0 + ALLOC_(MS%S(isd:ied,jsd:jed,nz)) ; MS%S(:,:,:) = 0.0 + MS%tv%T => MS%T ; MS%tv%S => MS%S + if (MS%tv%T_is_conT) then CS%vd_T = var_desc(name="contemp", units="Celsius", longname="Conservative Temperature", & cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", & - conversion=CS%tv%C_p) + conversion=MS%tv%C_p) else CS%vd_T = var_desc(name="temp", units="degC", longname="Potential Temperature", & cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", & - conversion=CS%tv%C_p) + conversion=MS%tv%C_p) endif - if (CS%tv%S_is_absS) then + if (MS%tv%S_is_absS) then CS%vd_S = var_desc(name="abssalt",units="g kg-1",longname="Absolute Salinity", & cmor_field_name="so", cmor_longname="Sea Water Salinity", & conversion=0.001) @@ -1814,7 +1864,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo if(CS%advect_TS) then S_flux_units = get_tr_flux_units(GV, "psu") ! Could change to "kg m-2 s-1"? - conv2watt = GV%H_to_kg_m2 * CS%tv%C_p + conv2watt = GV%H_to_kg_m2 * MS%tv%C_p if (GV%Boussinesq) then conv2salt = GV%H_to_m ! Could change to GV%H_to_kg_m2 * 0.001? H_convert = GV%H_to_m @@ -1822,25 +1872,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo conv2salt = GV%H_to_kg_m2 H_convert = GV%H_to_kg_m2 endif - call register_tracer(CS%tv%T, CS%tracer_Reg, param_file, dG%HI, GV, & + call register_tracer(MS%tv%T, CS%tracer_Reg, param_file, dG%HI, GV, & tr_desc=CS%vd_T, registry_diags=.true., flux_nameroot='T', & flux_units='W m-2', flux_longname='Heat', & flux_scale=conv2watt, convergence_units='W m-2', & convergence_scale=conv2watt, CMOR_tendname="opottemp", diag_form=2) - call register_tracer(CS%tv%S, CS%tracer_Reg, param_file, dG%HI, GV, & + call register_tracer(MS%tv%S, CS%tracer_Reg, param_file, dG%HI, GV, & tr_desc=CS%vd_S, registry_diags=.true., flux_nameroot='S', & flux_units=S_flux_units, flux_longname='Salt', & flux_scale=conv2salt, convergence_units='kg m-2 s-1', & convergence_scale=0.001*GV%H_to_kg_m2, CMOR_tendname="osalt", diag_form=2) endif if (associated(CS%OBC)) & - call register_temp_salt_segments(GV, CS%OBC, CS%tv, CS%vd_T, CS%vd_S, param_file) + call register_temp_salt_segments(GV, CS%OBC, MS%tv, CS%vd_T, CS%vd_S, param_file) endif if (CS%use_frazil) then - allocate(CS%tv%frazil(isd:ied,jsd:jed)) ; CS%tv%frazil(:,:) = 0.0 + allocate(MS%tv%frazil(isd:ied,jsd:jed)) ; MS%tv%frazil(:,:) = 0.0 endif if (CS%bound_salinity) then - allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:)=0.0 + allocate(MS%tv%salt_deficit(isd:ied,jsd:jed)) ; MS%tv%salt_deficit(:,:)=0.0 endif if (CS%bulkmixedlayer .or. CS%use_temperature) then @@ -1856,9 +1906,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call get_param(param_file, "MOM", "NK_RHO_VARIES", GV%nk_rho_varies, default=0) ! Will default to nz later... -AJA endif - ALLOC_(CS%uhtr(IsdB:IedB,jsd:jed,nz)) ; CS%uhtr(:,:,:) = 0.0 - ALLOC_(CS%vhtr(isd:ied,JsdB:JedB,nz)) ; CS%vhtr(:,:,:) = 0.0 - CS%t_dyn_rel_adv = 0.0 ; CS%t_dyn_rel_thermo = 0.0 ; CS%t_dyn_rel_diag = 0.0 + ALLOC_(MS%uhtr(IsdB:IedB,jsd:jed,nz)) ; MS%uhtr(:,:,:) = 0.0 + ALLOC_(MS%vhtr(isd:ied,JsdB:JedB,nz)) ; MS%vhtr(:,:,:) = 0.0 + MS%t_dyn_rel_adv = 0.0 ; CS%t_dyn_rel_thermo = 0.0 ; CS%t_dyn_rel_diag = 0.0 if (CS%debug_truncations) then allocate(CS%u_prev(IsdB:IedB,jsd:jed,nz)) ; CS%u_prev(:,:,:) = 0.0 @@ -1871,68 +1921,68 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo endif endif - MOM_internal_state%u => CS%u ; MOM_internal_state%v => CS%v - MOM_internal_state%h => CS%h - MOM_internal_state%uh => CS%uh ; MOM_internal_state%vh => CS%vh + MOM_internal_state%u => MS%u ; MOM_internal_state%v => MS%v + MOM_internal_state%h => MS%h + MOM_internal_state%uh => MS%uh ; MOM_internal_state%vh => MS%vh if (CS%use_temperature) then - MOM_internal_state%T => CS%T ; MOM_internal_state%S => CS%S + MOM_internal_state%T => MS%T ; MOM_internal_state%S => MS%S endif - CS%CDp%uh => CS%uh ; CS%CDp%vh => CS%vh + CS%CDp%uh => MS%uh ; CS%CDp%vh => MS%vh if (CS%interp_p_surf) then allocate(CS%p_surf_prev(isd:ied,jsd:jed)) ; CS%p_surf_prev(:,:) = 0.0 endif - ALLOC_(CS%ave_ssh(isd:ied,jsd:jed)) ; CS%ave_ssh(:,:) = 0.0 + ALLOC_(MS%ave_ssh(isd:ied,jsd:jed)) ; MS%ave_ssh(:,:) = 0.0 ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate ! initialization routine for tv. - if (use_EOS) call EOS_init(param_file, CS%tv%eqn_of_state) + if (use_EOS) call EOS_init(param_file, MS%tv%eqn_of_state) if (CS%use_temperature) then - allocate(CS%tv%TempxPmE(isd:ied,jsd:jed)) - CS%tv%TempxPmE(:,:) = 0.0 + allocate(MS%tv%TempxPmE(isd:ied,jsd:jed)) + MS%tv%TempxPmE(:,:) = 0.0 if (use_geothermal) then - allocate(CS%tv%internal_heat(isd:ied,jsd:jed)) - CS%tv%internal_heat(:,:) = 0.0 + allocate(MS%tv%internal_heat(isd:ied,jsd:jed)) + MS%tv%internal_heat(:,:) = 0.0 endif endif call callTree_waypoint("state variables allocated (initialize_MOM)") ! Set the fields that are needed for bitwise identical restarting ! the time stepping scheme. - call restart_init(param_file, CS%restart_CSp) - call set_restart_fields(GV, param_file, CS) + call restart_init(param_file, restart_CSp) + call set_restart_fields(GV, param_file, MS, CS, restart_CSp) if (CS%split) then call register_restarts_dyn_split_RK2(dG%HI, GV, param_file, & - CS%dyn_split_RK2_CSp, CS%restart_CSp, CS%uh, CS%vh) + CS%dyn_split_RK2_CSp, restart_CSp, MS%uh, MS%vh) elseif (CS%use_RK2) then call register_restarts_dyn_unsplit_RK2(dG%HI, GV, param_file, & - CS%dyn_unsplit_RK2_CSp, CS%restart_CSp) + CS%dyn_unsplit_RK2_CSp, restart_CSp) else call register_restarts_dyn_unsplit(dG%HI, GV, param_file, & - CS%dyn_unsplit_CSp, CS%restart_CSp) + CS%dyn_unsplit_CSp, restart_CSp) endif ! This subroutine calls user-specified tracer registration routines. ! Additional calls can be added to MOM_tracer_flow_control.F90. call call_tracer_register(dG%HI, GV, param_file, CS%tracer_flow_CSp, & - CS%tracer_Reg, CS%restart_CSp) + CS%tracer_Reg, restart_CSp) - call MEKE_alloc_register_restart(dG%HI, param_file, CS%MEKE, CS%restart_CSp) - call set_visc_register_restarts(dG%HI, GV, param_file, CS%visc, CS%restart_CSp) - call mixedlayer_restrat_register_restarts(dG%HI, param_file, CS%mixedlayer_restrat_CSp, CS%restart_CSp) + call MEKE_alloc_register_restart(dG%HI, param_file, CS%MEKE, restart_CSp) + call set_visc_register_restarts(dG%HI, GV, param_file, CS%visc, restart_CSp) + call mixedlayer_restrat_register_restarts(dG%HI, param_file, CS%mixedlayer_restrat_CSp, restart_CSp) if (associated(CS%OBC)) & - call open_boundary_register_restarts(dg%HI, GV, CS%OBC, CS%restart_CSp) + call open_boundary_register_restarts(dg%HI, GV, CS%OBC, restart_CSp) call callTree_waypoint("restart registration complete (initialize_MOM)") ! Initialize dynamically evolving fields, perhaps from restart files. call cpu_clock_begin(id_clock_MOM_init) call MOM_initialize_coord(GV, param_file, write_geom_files, & - dirs%output_directory, CS%tv, dG%max_depth) + dirs%output_directory, MS%tv, dG%max_depth) call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)") if (CS%use_ALE_algorithm) then @@ -1957,8 +2007,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo ! Consider removing this later? G%ke = GV%ke ; G%g_Earth = GV%g_Earth - call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, & - dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + call MOM_initialize_state(MS%u, MS%v, MS%h, MS%tv, Time, G, GV, param_file, & + dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") @@ -1967,22 +2017,22 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo ! that will persist throughout the run has to be used. if (test_grid_copy) then - ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G. + ! Copy the data from the temporary grid to the dyn_hor_grid to MS%G. call create_dyn_horgrid(dG, G%HI) call clone_MOM_domain(G%Domain, dG%Domain) - call clone_MOM_domain(G%Domain, CS%G%Domain) - call MOM_grid_init(CS%G, param_file) + call clone_MOM_domain(G%Domain, MS%G%Domain) + call MOM_grid_init(MS%G, param_file) call copy_MOM_grid_to_dyngrid(G, dg) - call copy_dyngrid_to_MOM_grid(dg, CS%G) + call copy_dyngrid_to_MOM_grid(dg, MS%G) call destroy_dyn_horgrid(dG) call MOM_grid_end(G) ; deallocate(G) - G => CS%G - if (CS%debug .or. CS%G%symmetric) & - call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) + G => MS%G + if (CS%debug .or. MS%G%symmetric) & + call clone_MOM_domain(MS%G%Domain, MS%G%Domain_aux, symmetric=.false.) G%ke = GV%ke ; G%g_Earth = GV%g_Earth endif @@ -1991,17 +2041,17 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo ! remainder of this subroutine is controlled by the parameters that have ! have already been set. - if (ALE_remap_init_conds(CS%ALE_CSp) .and. .not. query_initialized(CS%h,"h",CS%restart_CSp)) then + if (ALE_remap_init_conds(CS%ALE_CSp) .and. .not. query_initialized(MS%h,"h",restart_CSp)) then ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION. ! \todo This block exists for legacy reasons and we should phase it out of ! all examples. !### if (CS%debug) then call uvchksum("Pre ALE adjust init cond [uv]", & - CS%u, CS%v, G%HI, haloshift=1) - call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m) + MS%u, MS%v, G%HI, haloshift=1) + call hchksum(MS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m) endif call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)") - call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h ) + call adjustGridForIntegrity(CS%ALE_CSp, G, GV, MS%h ) call callTree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)") if (use_ice_shelf) then filename = trim(inputdir)//trim(ice_shelf_file) @@ -2020,25 +2070,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo enddo ; enddo ! pass to the pointer shelf_area => frac_shelf_h - call ALE_main(G, GV, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, & + call ALE_main(G, GV, MS%h, MS%u, MS%v, MS%tv, CS%tracer_Reg, CS%ALE_CSp, & frac_shelf_h = shelf_area) else - call ALE_main( G, GV, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp ) + call ALE_main( G, GV, MS%h, MS%u, MS%v, MS%tv, CS%tracer_Reg, CS%ALE_CSp ) endif call cpu_clock_begin(id_clock_pass_init) - call create_group_pass(tmp_pass_uv_T_S_h, CS%u, CS%v, G%Domain) + call create_group_pass(tmp_pass_uv_T_S_h, MS%u, MS%v, G%Domain) if (CS%use_temperature) then - call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain, halo=1) - call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain, halo=1) + call create_group_pass(tmp_pass_uv_T_S_h, MS%tv%T, G%Domain, halo=1) + call create_group_pass(tmp_pass_uv_T_S_h, MS%tv%S, G%Domain, halo=1) endif - call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain, halo=1) + call create_group_pass(tmp_pass_uv_T_S_h, MS%h, G%Domain, halo=1) call do_group_pass(tmp_pass_uv_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass_init) if (CS%debug) then - call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1) - call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m) + call uvchksum("Post ALE adjust init cond [uv]", MS%u, MS%v, G%HI, haloshift=1) + call hchksum(MS%h, "Post ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m) endif endif if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV ) @@ -2046,6 +2096,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo diag => CS%diag ! Initialize the diag mediator. call diag_mediator_init(G, GV, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) + if (present(diag_ptr)) diag_ptr => CS%diag ! Initialize the diagnostics masks for native arrays. ! This step has to be done after call to MOM_initialize_state @@ -2053,8 +2104,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call diag_masks_set(G, GV%ke, diag) ! Set up pointers within diag mediator control structure, - ! this needs to occur _after_ CS%h etc. have been allocated. - call diag_set_state_ptrs(CS%h, CS%T, CS%S, CS%tv%eqn_of_state, diag) + ! this needs to occur _after_ MS%h etc. have been allocated. + call diag_set_state_ptrs(MS%h, MS%T, MS%S, MS%tv%eqn_of_state, diag) ! This call sets up the diagnostic axes. These are needed, ! e.g. to generate the target grids below. @@ -2072,7 +2123,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call set_masks_for_axes(G, diag) ! Diagnose static fields AND associate areas/volumes with axes - call write_static_fields(G, GV, CS%tv, CS%diag) + call write_static_fields(G, GV, MS%tv, CS%diag) call callTree_waypoint("static fields written (initialize_MOM)") ! Register the volume cell measure (must be one of first diagnostics) @@ -2085,25 +2136,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("ALE initialized (initialize_MOM)") - CS%useMEKE = MEKE_init(Time, G, param_file, diag, CS%MEKE_CSp, CS%MEKE, CS%restart_CSp) + CS%useMEKE = MEKE_init(Time, G, param_file, diag, CS%MEKE_CSp, CS%MEKE, restart_CSp) call VarMix_init(Time, G, param_file, diag, CS%VarMix) call set_visc_init(Time, G, GV, param_file, diag, CS%visc, CS%set_visc_CSp,CS%OBC) if (CS%split) then allocate(eta(SZI_(G),SZJ_(G))) ; eta(:,:) = 0.0 - call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & - G, GV, param_file, diag, CS%dyn_split_RK2_CSp, CS%restart_CSp, & + call initialize_dyn_split_RK2(MS%u, MS%v, MS%h, MS%uh, MS%vh, eta, Time, & + G, GV, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc) elseif (CS%use_RK2) then - call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, Time, G, GV, & - param_file, diag, CS%dyn_unsplit_RK2_CSp, CS%restart_CSp, & + call initialize_dyn_unsplit_RK2(MS%u, MS%v, MS%h, Time, G, GV, & + param_file, diag, CS%dyn_unsplit_RK2_CSp, restart_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%OBC, CS%update_OBC_CSp, & CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, CS%ntrunc) else - call initialize_dyn_unsplit(CS%u, CS%v, CS%h, Time, G, GV, & - param_file, diag, CS%dyn_unsplit_CSp, CS%restart_CSp, & + call initialize_dyn_unsplit(MS%u, MS%v, MS%h, Time, G, GV, & + param_file, diag, CS%dyn_unsplit_CSp, restart_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%OBC, CS%update_OBC_CSp, & CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, CS%ntrunc) endif @@ -2121,7 +2172,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo endif call MOM_diagnostics_init(MOM_internal_state, CS%ADp, CS%CDp, Time, G, GV, & - param_file, diag, CS%diagnostics_CSp, CS%tv) + param_file, diag, CS%diagnostics_CSp, MS%tv) CS%Z_diag_interval = set_time(int((CS%dt_therm) * & max(1,floor(0.01 + Z_diag_int/(CS%dt_therm))))) @@ -2139,57 +2190,58 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call adiabatic_driver_init(Time, G, param_file, diag, CS%diabatic_CSp, & CS%tracer_flow_CSp, CS%diag_to_Z_CSp) else - call diabatic_driver_init(Time, G, GV, param_file, CS%use_ALE_algorithm, diag, & + call diabatic_driver_init(Time, G, GV, param_file, CS%use_ALE_algorithm, diag, & CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%diag_to_Z_CSp) endif call tracer_advect_init(Time, G, param_file, diag, CS%tracer_adv_CSp) - call tracer_hor_diff_init(Time, G, param_file, diag, CS%tv%eqn_of_state, CS%tracer_diff_CSp, CS%neutral_diffusion_CSp) + call tracer_hor_diff_init(Time, G, param_file, diag, MS%tv%eqn_of_state, & + CS%tracer_diff_CSp, CS%neutral_diffusion_CSp) call lock_tracer_registry(CS%tracer_Reg) call callTree_waypoint("tracer registry now locked (initialize_MOM)") ! now register some diagnostics since the tracer registry is now locked - call register_surface_diags(Time, G, CS%sfc_IDs, CS%diag, CS%missing, CS%tv) + call register_surface_diags(Time, G, CS%sfc_IDs, CS%diag, CS%missing, MS%tv) call register_diags(Time, G, GV, CS%IDs, CS%diag, CS%missing) call register_transport_diags(Time, G, GV, CS%transport_IDs, CS%diag, CS%missing) - call register_tracer_diagnostics(CS%tracer_Reg, CS%h, Time, diag, G, GV, & + call register_tracer_diagnostics(CS%tracer_Reg, MS%h, Time, diag, G, GV, & CS%use_ALE_algorithm, CS%diag_to_Z_CSp) if (CS%use_ALE_algorithm) then call ALE_register_diags(Time, G, GV, diag, CS%ALE_CSp) endif ! This subroutine initializes any tracer packages. - new_sim = is_new_run(CS%restart_CSp) - call tracer_flow_control_init(.not.new_sim, Time, G, GV, CS%h, param_file, & + new_sim = is_new_run(restart_CSp) + call tracer_flow_control_init(.not.new_sim, Time, G, GV, MS%h, param_file, & CS%diag, CS%OBC, CS%tracer_flow_CSp, CS%sponge_CSp, & - CS%ALE_sponge_CSp, CS%diag_to_Z_CSp, CS%tv) - + CS%ALE_sponge_CSp, CS%diag_to_Z_CSp, MS%tv) + if (present(tracer_flow_CSp)) tracer_flow_CSp => CS%tracer_flow_CSp ! If running in offline tracer mode, initialize the necessary control structure and ! parameters - if(present(offline_tracer_mode)) offline_tracer_mode=CS%offline_tracer_mode + if (present(offline_tracer_mode)) offline_tracer_mode=CS%offline_tracer_mode - if(CS%offline_tracer_mode) then + if (CS%offline_tracer_mode) then ! Setup some initial parameterizations and also assign some of the subtypes call offline_transport_init(param_file, CS%offline_CSp, CS%diabatic_CSp, G, GV) call insert_offline_main( CS=CS%offline_CSp, ALE_CSp=CS%ALE_CSp, diabatic_CSp=CS%diabatic_CSp, & diag=CS%diag, OBC=CS%OBC, tracer_adv_CSp=CS%tracer_adv_CSp, & tracer_flow_CSp=CS%tracer_flow_CSp, tracer_Reg=CS%tracer_Reg, & - tv=CS%tv, x_before_y = (MOD(first_direction,2)==0), debug=CS%debug ) + tv=MS%tv, x_before_y = (MOD(first_direction,2)==0), debug=CS%debug ) call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp) endif !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM call cpu_clock_begin(id_clock_pass_init) dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) - call create_group_pass(CS%pass_uv_T_S_h, CS%u, CS%v, G%Domain, halo=dynamics_stencil) + call create_group_pass(CS%pass_uv_T_S_h, MS%u, MS%v, G%Domain, halo=dynamics_stencil) if (CS%use_temperature) then - call create_group_pass(CS%pass_uv_T_S_h, CS%tv%T, G%Domain, halo=dynamics_stencil) - call create_group_pass(CS%pass_uv_T_S_h, CS%tv%S, G%Domain, halo=dynamics_stencil) + call create_group_pass(CS%pass_uv_T_S_h, MS%tv%T, G%Domain, halo=dynamics_stencil) + call create_group_pass(CS%pass_uv_T_S_h, MS%tv%S, G%Domain, halo=dynamics_stencil) endif - call create_group_pass(CS%pass_uv_T_S_h, CS%h, G%Domain, halo=dynamics_stencil) + call create_group_pass(CS%pass_uv_T_S_h, MS%h, G%Domain, halo=dynamics_stencil) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) @@ -2201,26 +2253,31 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo call register_obsolete_diagnostics(param_file, CS%diag) if (CS%use_frazil) then - if (.not.query_initialized(CS%tv%frazil,"frazil",CS%restart_CSp)) & - CS%tv%frazil(:,:) = 0.0 + if (.not.query_initialized(MS%tv%frazil,"frazil",restart_CSp)) & + MS%tv%frazil(:,:) = 0.0 endif if (CS%interp_p_surf) then CS%p_surf_prev_set = & - query_initialized(CS%p_surf_prev,"p_surf_prev",CS%restart_CSp) + query_initialized(CS%p_surf_prev,"p_surf_prev",restart_CSp) if (CS%p_surf_prev_set) call pass_var(CS%p_surf_prev, G%domain) endif - if (.not.query_initialized(CS%ave_ssh,"ave_ssh",CS%restart_CSp)) then + if (.not.query_initialized(MS%ave_ssh,"ave_ssh",restart_CSp)) then if (CS%split) then - call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, CS%ave_ssh, eta) + call find_eta(MS%h, MS%tv, GV%g_Earth, G, GV, MS%ave_ssh, eta) else - call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, CS%ave_ssh) + call find_eta(MS%h, MS%tv, GV%g_Earth, G, GV, MS%ave_ssh) endif endif if (CS%split) deallocate(eta) + CS%nstep_tot = 0 + if (present(count_calls)) CS%count_calls = count_calls + call MOM_sum_output_init(G, param_file, dirs%output_directory, & + CS%ntrunc, Time_init, CS%sum_output_CSp) + ! Flag whether to save initial conditions in finish_MOM_initialization() or not. CS%write_IC = save_IC .and. & .not.((dirs%input_filename(1:1) == 'r') .and. & @@ -2232,11 +2289,14 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo end subroutine initialize_MOM !> This subroutine finishes initializing MOM and writes out the initial conditions. -subroutine finish_MOM_initialization(Time, dirs, CS, fluxes) - type(time_type), intent(in) :: Time !< model time, used in this routine - type(directories), intent(in) :: dirs !< structure with directory paths - type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure - type(forcing), intent(inout) :: fluxes !< pointers to forcing fields +subroutine finish_MOM_initialization(Time, dirs, MS, CS, fluxes, restart_CSp) + type(time_type), intent(in) :: Time !< model time, used in this routine + type(directories), intent(in) :: dirs !< structure with directory paths + type(MOM_state_type), pointer :: MS !< pointer to structure describing the MOM state + type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure + type(forcing), intent(inout) :: fluxes !< pointers to forcing fields + type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control + !! structure that will be used for MOM. ! Local variables type(ocean_grid_type), pointer :: G => NULL() type(verticalGrid_type), pointer :: GV => NULL() @@ -2249,14 +2309,14 @@ subroutine finish_MOM_initialization(Time, dirs, CS, fluxes) call callTree_enter("finish_MOM_initialization()") ! Pointers for convenience - G => CS%G ; GV => CS%GV + G => MS%G ; GV => MS%GV ! Write initial conditions if (CS%write_IC) then allocate(restart_CSp_tmp) - restart_CSp_tmp = CS%restart_CSp + restart_CSp_tmp = restart_CSp allocate(z_interface(SZI_(G),SZJ_(G),SZK_(G)+1)) - call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, z_interface) + call find_eta(MS%h, MS%tv, GV%g_Earth, G, GV, z_interface) call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, & "Interface heights", "meter", z_grid='i') @@ -2266,6 +2326,9 @@ subroutine finish_MOM_initialization(Time, dirs, CS, fluxes) deallocate(restart_CSp_tmp) endif + call write_energy(MS%u, MS%v, MS%h, MS%tv, Time, 0, G, GV, & + CS%sum_output_CSp, CS%tracer_flow_CSp) + call callTree_leave("finish_MOM_initialization()") call cpu_clock_end(id_clock_init) @@ -2368,7 +2431,7 @@ subroutine MOM_timing_init(CS) id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=CLOCK_MODULE) id_clock_Z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=CLOCK_MODULE) id_clock_ALE = cpu_clock_id('(Ocean ALE)', grain=CLOCK_MODULE) - if(CS%offline_tracer_mode) then + if (CS%offline_tracer_mode) then id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=CLOCK_SUBCOMPONENT) endif @@ -2643,10 +2706,13 @@ end subroutine write_static_fields !! This routine should be altered if there are any changes to the !! time stepping scheme. The CHECK_RESTART facility may be used to !! confirm that all needed restart fields have been included. -subroutine set_restart_fields(GV, param_file, CS) +subroutine set_restart_fields(GV, param_file, MS, CS, restart_CSp) type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters + type(MOM_state_type), intent(in) :: MS !< structure describing the MOM state type(MOM_control_struct), intent(in) :: CS !< control structure set up by inialize_MOM + type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control + !! structure that will be used for MOM. ! Local variables logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts character(len=48) :: thickness_units, flux_units @@ -2657,37 +2723,37 @@ subroutine set_restart_fields(GV, param_file, CS) flux_units = get_flux_units(GV) if (CS%use_temperature) then - call register_restart_field(CS%tv%T, "Temp", .true., CS%restart_CSp, & + call register_restart_field(MS%tv%T, "Temp", .true., restart_CSp, & "Potential Temperature", "degC") - call register_restart_field(CS%tv%S, "Salt", .true., CS%restart_CSp, & + call register_restart_field(MS%tv%S, "Salt", .true., restart_CSp, & "Salinity", "PPT") endif - call register_restart_field(CS%h, "h", .true., CS%restart_CSp, & + call register_restart_field(MS%h, "h", .true., restart_CSp, & "Layer Thickness", thickness_units) - call register_restart_field(CS%u, "u", .true., CS%restart_CSp, & + call register_restart_field(MS%u, "u", .true., restart_CSp, & "Zonal velocity", "m s-1", hor_grid='Cu') - call register_restart_field(CS%v, "v", .true., CS%restart_CSp, & + call register_restart_field(MS%v, "v", .true., restart_CSp, & "Meridional velocity", "m s-1", hor_grid='Cv') if (CS%use_frazil) then - call register_restart_field(CS%tv%frazil, "frazil", .false., CS%restart_CSp, & + call register_restart_field(MS%tv%frazil, "frazil", .false., restart_CSp, & "Frazil heat flux into ocean", "J m-2") endif if (CS%interp_p_surf) then - call register_restart_field(CS%p_surf_prev, "p_surf_prev", .false., CS%restart_CSp, & + call register_restart_field(CS%p_surf_prev, "p_surf_prev", .false., restart_CSp, & "Previous ocean surface pressure", "Pa") endif - call register_restart_field(CS%ave_ssh, "ave_ssh", .false., CS%restart_CSp, & + call register_restart_field(MS%ave_ssh, "ave_ssh", .false., restart_CSp, & "Time average sea surface height", "meter") ! hML is needed when using the ice shelf module if (use_ice_shelf .and. associated(CS%Hml)) then - call register_restart_field(CS%Hml, "hML", .false., CS%restart_CSp, & + call register_restart_field(CS%Hml, "hML", .false., restart_CSp, & "Mixed layer thickness", "meter") endif @@ -2791,7 +2857,7 @@ end subroutine allocate_surface_state !> This subroutine sets the surface (return) properties of the ocean !! model by setting the appropriate fields in state. Unused fields !! are set to NULL or are unallocated. -subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) +subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, MS, CS) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(surface), intent(inout) :: sfc_state !< ocean surface state @@ -2799,6 +2865,7 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< meridional velocity (m/s) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness (m or kg/m2) real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< time mean surface height (m) + type(MOM_state_type), intent(in) :: MS !< structure describing the MOM state type(MOM_control_struct), intent(inout) :: CS !< control structure ! local @@ -2826,11 +2893,11 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) ! integrals, since the 3-d sums are not negligible in cost. call allocate_surface_state(sfc_state, G, CS%use_temperature, do_integrals=.true.) endif - sfc_state%frazil => CS%tv%frazil - sfc_state%TempxPmE => CS%tv%TempxPmE - sfc_state%internal_heat => CS%tv%internal_heat - sfc_state%T_is_conT = CS%tv%T_is_conT - sfc_state%S_is_absS = CS%tv%S_is_absS + sfc_state%frazil => MS%tv%frazil + sfc_state%TempxPmE => MS%tv%TempxPmE + sfc_state%internal_heat => MS%tv%internal_heat + sfc_state%T_is_conT = MS%tv%T_is_conT + sfc_state%S_is_absS = MS%tv%S_is_absS if (associated(CS%visc%taux_shelf)) sfc_state%taux_shelf => CS%visc%taux_shelf if (associated(CS%visc%tauy_shelf)) sfc_state%tauy_shelf => CS%visc%tauy_shelf @@ -2840,8 +2907,8 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) if (CS%bulkmixedlayer) then if (CS%use_temperature) then ; do j=js,je ; do i=is,ie - sfc_state%SST(i,j) = CS%tv%T(i,j,1) - sfc_state%SSS(i,j) = CS%tv%S(i,j,1) + sfc_state%SST(i,j) = MS%tv%T(i,j,1) + sfc_state%SSS(i,j) = MS%tv%S(i,j,1) enddo ; enddo ; endif do j=js,je ; do I=IscB,IecB sfc_state%u(I,j) = u(I,j,1) @@ -2877,8 +2944,8 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) dh = 0.0 endif if (CS%use_temperature) then - sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * CS%tv%T(i,j,k) - sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * CS%tv%S(i,j,k) + sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * MS%tv%T(i,j,k) + sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * MS%tv%S(i,j,k) else sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * GV%Rlay(k) endif @@ -2962,11 +3029,11 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) endif endif ! end BULKMIXEDLAYER - if (allocated(sfc_state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then + if (allocated(sfc_state%salt_deficit) .and. associated(MS%tv%salt_deficit)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie ! Convert from gSalt to kgSalt - sfc_state%salt_deficit(i,j) = 1000.0 * CS%tv%salt_deficit(i,j) + sfc_state%salt_deficit(i,j) = 1000.0 * MS%tv%salt_deficit(i,j) enddo ; enddo endif @@ -2981,9 +3048,9 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) do j=js,je ; do k=1,nz; do i=is,ie mass = GV%H_to_kg_m2*h(i,j,k) sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass - sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k) + sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*MS%tv%T(i,j,k) sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + & - mass * (1.0e-3*CS%tv%S(i,j,k)) + mass * (1.0e-3*MS%tv%S(i,j,k)) enddo ; enddo ; enddo else if (allocated(sfc_state%ocean_mass)) then @@ -3000,7 +3067,7 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) !$OMP parallel do default(shared) private(mass) do j=js,je ; do k=1,nz ; do i=is,ie mass = GV%H_to_kg_m2*h(i,j,k) - sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k) + sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*MS%tv%T(i,j,k) enddo ; enddo ; enddo endif if (allocated(sfc_state%ocean_salt)) then @@ -3010,7 +3077,7 @@ subroutine calculate_surface_state(sfc_state, u, v, h, ssh, G, GV, CS) do j=js,je ; do k=1,nz ; do i=is,ie mass = GV%H_to_kg_m2*h(i,j,k) sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + & - mass * (1.0e-3*CS%tv%S(i,j,k)) + mass * (1.0e-3*MS%tv%S(i,j,k)) enddo ; enddo ; enddo endif endif @@ -3078,21 +3145,22 @@ end subroutine calculate_surface_state !> End of model -subroutine MOM_end(CS) +subroutine MOM_end(MS, CS) + type(MOM_state_type), pointer :: MS !< structure describing the MOM state type(MOM_control_struct), pointer :: CS !< MOM control structure if (CS%use_ALE_algorithm) then call ALE_end(CS%ALE_CSp) endif - DEALLOC_(CS%u) ; DEALLOC_(CS%v) ; DEALLOC_(CS%h) - DEALLOC_(CS%uh) ; DEALLOC_(CS%vh) + DEALLOC_(MS%u) ; DEALLOC_(MS%v) ; DEALLOC_(MS%h) + DEALLOC_(MS%uh) ; DEALLOC_(MS%vh) if (CS%use_temperature) then - DEALLOC_(CS%T) ; CS%tv%T => NULL() ; DEALLOC_(CS%S) ; CS%tv%S => NULL() + DEALLOC_(MS%T) ; MS%tv%T => NULL() ; DEALLOC_(MS%S) ; MS%tv%S => NULL() endif - if (associated(CS%tv%frazil)) deallocate(CS%tv%frazil) - if (associated(CS%tv%salt_deficit)) deallocate(CS%tv%salt_deficit) + if (associated(MS%tv%frazil)) deallocate(MS%tv%frazil) + if (associated(MS%tv%salt_deficit)) deallocate(MS%tv%salt_deficit) if (associated(CS%Hml)) deallocate(CS%Hml) call tracer_advect_end(CS%tracer_adv_CSp) @@ -3100,11 +3168,11 @@ subroutine MOM_end(CS) call tracer_registry_end(CS%tracer_Reg) call tracer_flow_control_end(CS%tracer_flow_CSp) - if(CS%offline_tracer_mode) then + if (CS%offline_tracer_mode) then call offline_transport_end(CS%offline_CSp) endif - DEALLOC_(CS%uhtr) ; DEALLOC_(CS%vhtr) + DEALLOC_(MS%uhtr) ; DEALLOC_(MS%vhtr) if (CS%split) then call end_dyn_split_RK2(CS%dyn_split_RK2_CSp) elseif (CS%use_RK2) then @@ -3112,12 +3180,13 @@ subroutine MOM_end(CS) else call end_dyn_unsplit(CS%dyn_unsplit_CSp) endif - DEALLOC_(CS%ave_ssh) + DEALLOC_(MS%ave_ssh) if (associated(CS%update_OBC_CSp)) call OBC_register_end(CS%update_OBC_CSp) - call verticalGridEnd(CS%GV) - call MOM_grid_end(CS%G) + call verticalGridEnd(MS%GV) + call MOM_grid_end(MS%G) + deallocate(MS) deallocate(CS) end subroutine MOM_end diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 2440724331..9ff9410f49 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -16,8 +16,7 @@ module MOM_sum_output !* * !* In addition, if the number of velocity truncations since the * !* previous call to write_energy exceeds maxtrunc or the total energy * -!* exceeds a very large threshold, the day is increased to Huge_time * -!* so that the model will gracefully halt itself. * +!* exceeds a very large threshold, a fatal termination is triggered. * !* * !* This file also contains a few miscelaneous initialization * !* calls to FMS-related modules. * @@ -51,8 +50,10 @@ module MOM_sum_output use MOM_io, only : APPEND_FILE, ASCII_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S -use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>), operator(-) -use MOM_time_manager, only : get_calendar_type, NO_CALENDAR +use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>) +use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<) +use MOM_time_manager, only : get_calendar_type, time_type_to_real, NO_CALENDAR use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_stocks use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -114,15 +115,25 @@ module MOM_sum_output net_salt_in_EFP, & ! correspondingly named variables above. net_heat_in_EFP, heat_prev_EFP, salt_prev_EFP, mass_prev_EFP real :: dt ! The baroclinic dynamics time step, in s. + + type(time_type) :: energysavedays !< The interval between writing the energies + !! and other integral quantities of the run. + type(time_type) :: energysavedays_geometric !< The starting interval for computing a geometric + !! progression of time deltas between calls to + !! write_energy. This interval will increase by a factor of 2. + !! after each call to write_energy. + logical :: energysave_geometric !< Logical to control whether calls to write_energy should + !! follow a geometric progression + type(time_type) :: write_energy_time !< The next time to write to the energy file. + type(time_type) :: geometric_end_time !< Time at which to stop the geometric progression + !! of calls to write_energy and revert to the standard + !! energysavedays interval + real :: timeunit ! The length of the units for the time ! axis, in s. logical :: date_stamped_output ! If true, use dates (not times) in messages to stdout. type(time_type) :: Start_time ! The start time of the simulation. ! Start_time is set in MOM_initialization.F90 - type(time_type) :: Huge_time ! A large time, which is used to indicate - ! that an error has been encountered - ! and the run should be terminated with - ! an error code. integer, pointer :: ntrunc ! The number of times the velocity has been ! truncated since the last call to write_energy. real :: max_Energy ! The maximum permitted energy per unit mass; @@ -145,7 +156,7 @@ module MOM_sum_output !> MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module. subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & - Input_start_time, CS) + Input_start_time, CS) 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. @@ -165,6 +176,7 @@ subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & ! (in) Input_start_time - The start time of the simulation. ! (in/out) CS - A pointer that is set to point to the control structure ! for this module + real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS. real :: Rho_0, maxvel ! This include declares and sets the variable "version". #include "version_variable.h" @@ -261,7 +273,26 @@ subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & CS%list_size = 0 endif - CS%Huge_time = set_time(INT(1e9),0) + 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",CS%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, "ENERGYSAVEDAYS_GEOMETRIC",CS%energysavedays_geometric, & + "The starting interval in units of TIMEUNIT for the first call \n"//& + "to save the energies of the run and other globally summed diagnostics. \n"//& + "The interval increases by a factor of 2. after each call to write_energy.",& + default=set_time(seconds=0), timeunit=Time_unit) + + if ((time_type_to_real(CS%energysavedays_geometric) > 0.) .and. & + (CS%energysavedays_geometric < CS%energysavedays)) then + CS%energysave_geometric = .true. + else + CS%energysave_geometric = .false. + endif + CS%Start_time = Input_start_time CS%ntrunc => ntrnc @@ -281,28 +312,29 @@ subroutine MOM_sum_output_end(CS) endif end subroutine MOM_sum_output_end -!> This subroutine calculates and writes the total model energy, the -!! energy and mass of each layer, and other globally integrated -!! physical quantities. -subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid - !! structure. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity, - !! in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H - !! (usually m or kg m-2). - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various - !! thermodynamic variables. - type(time_type), intent(inout) :: day !< The current model time. - integer, intent(in) :: n !< The time step number of the - !! current execution. - type(Sum_output_CS), pointer :: CS !< The control structure returned - !! by a previous call to - !! MOM_sum_output_init. - type(tracer_flow_control_CS), optional, pointer :: tracer_CSp !< tracer constrol structure. - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. +!> This subroutine calculates and writes the total model energy, the energy and +!! mass of each layer, and other globally integrated physical quantities. +subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC, dt_forcing) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & + intent(in) :: u !< The zonal velocity, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: v !< The meridional velocity, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables. + type(time_type), intent(in) :: day !< The current model time. + integer, intent(in) :: n !< The time step number of the + !! current execution. + type(Sum_output_CS), pointer :: CS !< The control structure returned by a + !! previous call to MOM_sum_output_init. + type(tracer_flow_control_CS), & + optional, pointer :: tracer_CSp !< tracer control structure. + type(ocean_OBC_type), & + optional, pointer :: OBC !< Open boundaries control structure. + type(time_type), optional, intent(in) :: dt_forcing !< The forcing time step real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! The height of interfaces, in m. real :: areaTm(SZI_(G),SZJ_(G)) ! A masked version of areaT, in m2. @@ -381,6 +413,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC) character(len=200) :: mesg character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str logical :: date_stamped + type(time_type) :: dt_force real :: Tr_stocks(MAX_FIELDS_) real :: Tr_min(MAX_FIELDS_),Tr_max(MAX_FIELDS_) real :: Tr_min_x(MAX_FIELDS_), Tr_min_y(MAX_FIELDS_), Tr_min_z(MAX_FIELDS_) @@ -398,6 +431,39 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC) ! A description for output of each of the fields. type(vardesc) :: vars(NUM_FIELDS+MAX_FIELDS_) + ! write_energy_time is the next integral multiple of energysavedays. + dt_force = set_time(seconds=2) ; if (present(dt_forcing)) dt_force = dt_forcing + if (CS%previous_calls == 0) then + if (CS%energysave_geometric) then + if (CS%energysavedays_geometric < CS%energysavedays) then + CS%write_energy_time = day + CS%energysavedays_geometric + CS%geometric_end_time = CS%Start_time + CS%energysavedays * & + (1 + (day - CS%Start_time) / CS%energysavedays) + else + CS%write_energy_time = CS%Start_time + CS%energysavedays * & + (1 + (day - CS%Start_time) / CS%energysavedays) + endif + else + CS%write_energy_time = CS%Start_time + CS%energysavedays * & + (1 + (day - CS%Start_time) / CS%energysavedays) + endif + elseif (day + (dt_force/2) <= CS%write_energy_time) then + return ! Do not write this step + else ! Determine the next write time before proceeding + if (CS%energysave_geometric) then + CS%energysavedays_geometric = CS%energysavedays_geometric*2 + if (CS%write_energy_time + CS%energysavedays_geometric >= & + CS%geometric_end_time) then + CS%write_energy_time = CS%geometric_end_time + CS%energysave_geometric = .false. ! stop geometric progression + else + CS%write_energy_time = CS%write_energy_time + CS%energysavedays_geometric + endif + else + CS%write_energy_time = CS%write_energy_time + CS%energysavedays + endif + endif + num_nc_fields = 17 if (.not.CS%use_temperature) num_nc_fields = 11 vars(1) = var_desc("Ntrunc","Nondim","Number of Velocity Truncations",'1','1') @@ -860,15 +926,12 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp, OBC) ! The second (impossible-looking) test looks for a NaN in En_mass. if ((En_mass>CS%max_Energy) .or. & ((En_mass>CS%max_Energy) .and. (En_massCS%maxtrunc) then - day = CS%Huge_time call MOM_error(FATAL, "write_energy : Ocean velocity has been truncated too many times.") endif CS%ntrunc = 0 diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index fd600df72c..0a9c66897d 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -318,7 +318,7 @@ subroutine add_tracer_OBC_values(name, Reg, OBC_inflow, OBC_in_u, OBC_in_v) integer :: m - if (.not. associated(Reg)) call MOM_error(FATAL, "add_tracer_OBC_values :"// & + if (.not. associated(Reg)) call MOM_error(FATAL, "add_tracer_OBC_values :"//& "register_tracer must be called before add_tracer_OBC_values") do m=1,Reg%ntr ; if (Reg%Tr(m)%name == trim(name)) exit ; enddo @@ -374,7 +374,7 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE, diag_ isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - if (.not. associated(Reg)) call MOM_error(FATAL, "register_tracer_diagnostics: "// & + if (.not. associated(Reg)) call MOM_error(FATAL, "register_tracer_diagnostics: "//& "register_tracer must be called before register_tracer_diagnostics") nTr_in = Reg%ntr @@ -486,33 +486,36 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE, diag_ ! Lateral diffusion convergence tendencies if (Tr%diag_form == 1) then - Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion = Tr%conv_scale, x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) + Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & + diag%axesTL, Time, "Lateral or neutral diffusion tracer content tendency for "//trim(shortnm), & + conv_units, conversion=Tr%conv_scale, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer concentration"// & - "tendency for "//trim(shortnm), conv_units, conversion = Tr%conv_scale, & + diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer concentration "//& + "tendency for "//trim(shortnm), conv_units, conversion = Tr%conv_scale, & x_cell_method = 'sum', y_cell_method = 'sum') else - cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& - trim(lowercase(flux_longname))//' content due to parameterized mesoscale diffusion' - Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer concentration tendency for", conv_units, & - conversion = Tr%conv_scale, cmor_field_name = trim(Tr%cmor_tendname)//'pmdiff', & - cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), & + cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//& + ' expressed as '//trim(lowercase(flux_longname))//& + ' content due to parameterized mesoscale diffusion' + Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & + diag%axesTL, Time, "Lateral or neutral diffusion tracer concentration tendency for "//trim(shortnm), & + conv_units, conversion = Tr%conv_scale, cmor_field_name = trim(Tr%cmor_tendname)//'pmdiff', & + cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), & x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& trim(lowercase(flux_longname))//' content due to parameterized mesoscale diffusion' - Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer concentration tendency for", & - conv_units, conversion = Tr%conv_scale, cmor_field_name=trim(Tr%cmor_tendname)//'pmdiff_2d', & - cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), & - x_cell_method = 'sum', y_cell_method = 'sum') + Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer "//& + "concentration tendency for "//trim(shortnm), conv_units, & + conversion=Tr%conv_scale, cmor_field_name=trim(Tr%cmor_tendname)//'pmdiff_2d', & + cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & + x_cell_method='sum', y_cell_method='sum') endif Tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', & - diag%axesTL, Time, "Lateral (neutral) tracer concentration tendency for", units//' s-1') + diag%axesTL, Time, "Lateral (neutral) tracer concentration tendency for "//trim(shortnm), & + trim(units)//' s-1') var_lname = "Net time tendency for "//lowercase(flux_longname) if (len_trim(Tr%cmor_tendname) == 0) then @@ -521,7 +524,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE, diag_ Tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', & diag%axesT1, Time, "Vertical sum of "//trim(lowercase(var_lname)), conv_units) else - cmor_var_lname = "Tendency of "//trim(cmor_longname)//" Expressed as "//trim(flux_longname)//" Content" + cmor_var_lname = "Tendency of "//trim(cmor_longname)//" Expressed as "//& + trim(flux_longname)//" Content" Tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', & diag%axesTL, Time, var_lname, conv_units, & cmor_field_name=trim(Tr%cmor_tendname)//"tend", & @@ -554,14 +558,15 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE, diag_ var_lname = "Vertical remapping tracer concentration tendency for "//trim(Reg%Tr(m)%name) Tr%id_remap_conc= register_diag_field('ocean_model', & trim(Tr%flux_nameroot)//'_tendency_vert_remap', diag%axesTL, Time, var_lname, & - trim(units)//'s-1') + trim(units)//' s-1') var_lname = "Vertical remapping tracer content tendency for "//trim(Reg%Tr(m)%flux_longname) Tr%id_remap_cont = register_diag_field('ocean_model', & trim(Tr%flux_nameroot)//'h_tendency_vert_remap', & diag%axesTL, Time, var_lname, flux_units, v_extensive=.true., conversion = Tr%conv_scale) - var_lname = "Vertical sum of vertical remapping tracer content tendency for "//trim(Reg%Tr(m)%flux_longname) + var_lname = "Vertical sum of vertical remapping tracer content tendency for "//& + trim(Reg%Tr(m)%flux_longname) Tr%id_remap_cont_2d = register_diag_field('ocean_model', & trim(Tr%flux_nameroot)//'h_tendency_vert_remap_2d', & diag%axesT1, Time, var_lname, flux_units, conversion = Tr%conv_scale) @@ -571,8 +576,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE, diag_ if (use_ALE .and. (Reg%ntr 0) unit2 = "("//trim(units)//")2" - Tr%id_tr_vardec = register_diag_field('ocean_model', trim(shortnm)//"_vardec", diag%axesTL, Time, & - "ALE variance decay for "//lowercase(longname), trim(unit2)//" s-1") + Tr%id_tr_vardec = register_diag_field('ocean_model', trim(shortnm)//"_vardec", diag%axesTL, & + Time, "ALE variance decay for "//lowercase(longname), trim(unit2)//" s-1") if (Tr%id_tr_vardec > 0) then ! Set up a new tracer for this tracer squared m2 = Reg%ntr+1