diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 4db2d95330..6095761247 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -47,6 +47,7 @@ program MOM_main use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end 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_tracers 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 @@ -134,6 +135,7 @@ program MOM_main ! representation of time_step. real :: time_step ! The time step of a call to step_MOM in seconds. real :: dt ! The baroclinic dynamics time step, in seconds. + real :: dt_off ! Offline time step in seconds integer :: ntstep ! The number of baroclinic dynamics time steps ! within time_step. @@ -169,6 +171,14 @@ program MOM_main logical :: unit_in_use integer :: initClock, mainClock, termClock + logical :: do_online ! If true, use the model in prognostic mode where + ! the barotropic and baroclinic dynamics, thermodynamics, + ! etc. are stepped forward integrated in time. + ! If false, then all of the above are bypassed with all + ! fields necessary to integrate only the tracer advection + ! and diffusion equation are read in from files stored from + ! a previous integration of the prognostic model + type(MOM_control_struct), pointer :: MOM_CSp => NULL() type(surface_forcing_CS), pointer :: surface_forcing_CSp => NULL() type(sum_output_CS), pointer :: sum_output_CSp => NULL() @@ -250,17 +260,18 @@ program MOM_main else Start_time = set_time(0,days=0) endif - + if (sum(date) >= 0) then ! 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 - call initialize_MOM(Time, param_file, dirs, MOM_CSp, segment_start_time) + ! Note the not before CS%d + call initialize_MOM(Time, param_file, dirs, MOM_CSp, segment_start_time, do_online_out = do_online) 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) + call initialize_MOM(Time, param_file, dirs, MOM_CSp, do_online_out=do_online) endif fluxes%C_p = MOM_CSp%tv%C_p ! Copy the heat capacity for consistency. @@ -298,8 +309,12 @@ program MOM_main call get_param(param_file, mod, "DT_FORCING", time_step, & "The time step for changing forcing, coupling with other \n"//& "components, or potentially writing certain diagnostics. \n"//& - "The default value is given by DT.", units="s", default=dt) - + "The default value is given by DT.", units="s", default=dt) + if (.not. do_online) then + call get_param(param_file, mod, "DT_OFFLINE", time_step, & + "Time step for the offline time step") + dt = time_step + endif ntstep = MAX(1,ceiling(time_step/dt - 0.001)) Time_step_ocean = set_time(int(floor(time_step+0.5))) @@ -350,6 +365,7 @@ program MOM_main "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 log_param(param_file, mod, "ELAPSED TIME AS MASTER", elapsed_time_master) ! Close the param_file. No further parsing of input is possible after this. @@ -399,8 +415,10 @@ program MOM_main call callTree_enter("Main loop, MOM_driver.F90",n) ! Set the forcing for the next steps. - call set_forcing(state, fluxes, Time, Time_step_ocean, grid, & + if (do_online) then + call set_forcing(state, fluxes, Time, Time_step_ocean, grid, & surface_forcing_CSp) + endif if (MOM_CSp%debug) then call MOM_forcing_chksum("After set forcing", fluxes, grid, haloshift=0) endif @@ -423,7 +441,8 @@ program MOM_main ! This call steps the model over a time time_step. Time1 = Master_Time ; Time = Master_Time - call step_MOM(fluxes, state, Time1, time_step, MOM_CSp) + if (do_online) call step_MOM(fluxes, state, Time1, time_step, MOM_CSp) + if (.not. do_online) call step_tracers(fluxes, state, Time1, time_step, MOM_CSp) ! Time = Time + Time_step_ocean ! This is here to enable fractional-second time steps. @@ -450,17 +469,21 @@ program MOM_main surface_forcing_CSp%handles) call disable_averaging(MOM_CSp%diag) - if (fluxes%fluxes_used) then - call enable_averaging(fluxes%dt_buoy_accum, Time, MOM_CSp%diag) - call forcing_diagnostics(fluxes, state, fluxes%dt_buoy_accum, grid, & - MOM_CSp%diag, surface_forcing_CSp%handles) - call accumulate_net_input(fluxes, state, fluxes%dt_buoy_accum, grid, sum_output_CSp) - call disable_averaging(MOM_CSp%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 couping timestep.") + if (do_online) then + if (fluxes%fluxes_used) then + call enable_averaging(fluxes%dt_buoy_accum, Time, MOM_CSp%diag) + call forcing_diagnostics(fluxes, state, fluxes%dt_buoy_accum, grid, & + MOM_CSp%diag, surface_forcing_CSp%handles) + call accumulate_net_input(fluxes, state, fluxes%dt_buoy_accum, grid, sum_output_CSp) + call disable_averaging(MOM_CSp%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%dt_trans == 0.0)) then @@ -502,7 +525,7 @@ program MOM_main if (Restart_control>=0) then if (MOM_CSp%dt_trans > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//& "with a non-zero dt_trans. Additional restart fields are required.") - if (.not.fluxes%fluxes_used) call MOM_error(FATAL, "End of MOM_main reached "//& + if (.not.fluxes%fluxes_used .and. do_online) call MOM_error(FATAL, "End of MOM_main reached "//& "with unused buoyancy fluxes. For conservation, the ocean restart "//& "files can only be created after the buoyancy forcing is applied.") diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a97fdd7dce..fa350742b6 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -97,14 +97,17 @@ module MOM_ALE integer, dimension(:), allocatable :: id_tracer_remap_tendency !< diagnostic id integer, dimension(:), allocatable :: id_Htracer_remap_tendency !< diagnostic id integer, dimension(:), allocatable :: id_Htracer_remap_tendency_2d !< diagnostic id - logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics + logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics + integer :: id_dzRegrid end type ! Publicly available functions public ALE_init public ALE_end -public ALE_main +public ALE_main +public ALE_main_offline +public ALE_offline_tracer_final public ALE_build_grid public ALE_remap_scalar public pressure_gradient_plm @@ -282,6 +285,9 @@ subroutine ALE_register_diags(Time, G, diag, C_p, Reg, CS) CS%id_Htracer_remap_tendency(:) = -1 CS%id_Htracer_remap_tendency_2d(:) = -1 + CS%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,Time, & + 'Change in interface height due to ALE regridding', 'meter') + if(ntr > 0) then do m=1,ntr @@ -375,7 +381,6 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt) type(tracer_registry_type), pointer :: Reg !< Tracer registry structure type(ALE_CS), pointer :: CS !< Regridding parameters and options real, optional, intent(in) :: dt !< Time step between calls to ALE_main() - ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step (m or Pa) @@ -399,6 +404,7 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt) if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)") ! Remap all variables from old grid h onto new grid h_new + call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, -dzRegrid, Reg, & u, v, CS%show_call_tree, dt ) @@ -414,8 +420,114 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt) enddo if (CS%show_call_tree) call callTree_leave("ALE_main()") + + if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag) + + end subroutine ALE_main +!> Takes care of (1) building a new grid and (2) remapping all variables between +!! the old grid and the new grid. The creation of the new grid can be based +!! on z coordinates, target interface densities, sigma coordinates or any +!! arbitrary coordinate system. +subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, dt) + type(ocean_grid_type), intent(in) :: G !< Ocean grid informations + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after last time step (m or Pa) + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure + type(tracer_registry_type), pointer :: Reg !< Tracer registry structure + type(ALE_CS), pointer :: CS !< Regridding parameters and options + real, optional, intent(in) :: dt !< Time step between calls to ALE_main() + ! Local variables + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step (m or Pa) + integer :: nk, i, j, k, isc, iec, jsc, jec + + nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + + if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90") + + if (present(dt)) then + call ALE_update_regrid_weights( dt, CS ) + endif + dzRegrid(:,:,:) = 0.0 + + ! Build new grid. The new grid is stored in h_new. The old grid is h. + ! Both are needed for the subsequent remapping of variables. + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid ) + + call check_grid( G, GV, h, 0. ) + + if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)") + + ! Remap all variables from old grid h onto new grid h_new + + call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, -dzRegrid, Reg, & + debug=CS%show_call_tree, dt=dt ) + + if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)") + + ! Override old grid with new one. The new grid 'h_new' is built in + ! one of the 'build_...' routines above. +!$OMP parallel do default(none) shared(isc,iec,jsc,jec,nk,h,h_new,CS) + do k = 1,nk + do j = jsc-1,jec+1 ; do i = isc-1,iec+1 + h(i,j,k) = h_new(i,j,k) + enddo ; enddo + enddo + + if (CS%show_call_tree) call callTree_leave("ALE_main()") + if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag) + +end subroutine ALE_main_offline + +!> Remaps all tracers from h onto h_target. This is intended to be called when tracers +!! are done offline. In the case where transports don't quite conserve, we still want to +!! make sure that layer thicknesses offline do not drift too far away from the online model +subroutine ALE_offline_tracer_final( G, GV, h, h_target, Reg, CS) + type(ocean_grid_type), intent(in) :: G !< Ocean grid informations + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after last time step (m or Pa) + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_target !< Current 3D grid obtained after last time step (m or Pa) + type(tracer_registry_type), pointer :: Reg !< Tracer registry structure + type(ALE_CS), pointer :: CS !< Regridding parameters and options + ! Local variables + + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions + integer :: nk, i, j, k, isc, iec, jsc, jec + + nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + + if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90") + + ! It does not seem that remap_all_state_vars uses dzRegrid for tracers, only for u, v + dzRegrid(:,:,:) = 0.0 + + call check_grid( G, GV, h, 0. ) + call check_grid( G, GV, h_target, 0. ) + + if (CS%show_call_tree) call callTree_waypoint("Source and target grids checked (ALE_offline_tracer)") + + ! Remap all variables from old grid h onto new grid h_new + + call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_target, -dzRegrid, Reg, & + debug=CS%show_call_tree ) + + if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_offline_tracer)") + + ! Override old grid with new one. The new grid 'h_new' is built in + ! one of the 'build_...' routines above. +!$OMP parallel do default(none) shared(isc,iec,jsc,jec,nk,h,h_target,CS) + do k = 1,nk + do j = jsc-1,jec+1 ; do i = isc-1,iec+1 + h(i,j,k) = h_target(i,j,k) + enddo ; enddo + enddo + + if (CS%show_call_tree) call callTree_leave("ALE_offline_tracer()") + +end subroutine ALE_offline_tracer_final + !> Check grid for negative thicknesses subroutine check_grid( G, GV, h, threshold ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f9d6b35da6..69760b17a3 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -39,7 +39,7 @@ module MOM use MOM_diag_mediator, only : register_scalar_field use MOM_diag_mediator, only : set_axes_info, diag_ctrl, diag_masks_set use MOM_domains, only : MOM_domains_init, clone_MOM_domain -use MOM_domains, only : sum_across_PEs, pass_var +use MOM_domains, only : sum_across_PEs, pass_var, pass_vector use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : start_group_pass, complete_group_pass @@ -63,6 +63,7 @@ module MOM use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags +use MOM_ALE, only : ALE_main_offline, ALE_offline_tracer_final use MOM_continuity, only : continuity, continuity_init, continuity_CS use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS @@ -117,6 +118,7 @@ module MOM use MOM_tracer_registry, only : lock_tracer_registry, tracer_registry_end use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_CS use MOM_tracer_flow_control, only : tracer_flow_control_init, call_tracer_surface_state +use MOM_tracer_flow_control, only : call_tracer_column_fns use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid use MOM_vert_friction, only : vertvisc, vertvisc_remnant use MOM_vert_friction, only : vertvisc_limit_vel, vertvisc_init @@ -124,6 +126,13 @@ module MOM use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units use MOM_wave_speed, only : wave_speed_init, wave_speed_CS +! Offline modules +use MOM_offline_transport, only : offline_transport_CS +use MOM_offline_transport, only : transport_by_files, next_modulo_time +use MOM_offline_transport, only : offline_transport_init, register_diags_offline_transport +use MOM_offline_transport, only : limit_mass_flux_3d, update_h_horizontal_flux, update_h_vertical_flux +use MOM_tracer_diabatic, only : applyTracerBoundaryFluxesInOut + implicit none ; private #include @@ -192,7 +201,8 @@ module MOM !! MOM_regridding module. logical :: do_dynamics !< If false, does not call step_MOM_dyn_*. This is an !! undocumented run-time flag that is fragile. - + logical :: do_online !< If false, step_tracers is called instead of step_MOM. + !! This is intended for running MOM6 in offline tracer mode real :: dt !< (baroclinic) dynamics time step (seconds) real :: dt_therm !< thermodynamics time step (seconds) logical :: thermo_spans_coupling !< If true, thermodynamic and tracer time @@ -343,6 +353,10 @@ module MOM integer :: id_S_preale = -1 integer :: id_e_preale = -1 + ! Diagnostics for tracer horizontal transport + integer :: id_uhtr = -1 + integer :: id_vhtr = -1 + ! The remainder provides pointers to child module control structures. type(MOM_dyn_unsplit_CS), pointer :: dyn_unsplit_CSp => NULL() type(MOM_dyn_unsplit_RK2_CS), pointer :: dyn_unsplit_RK2_CSp => NULL() @@ -368,6 +382,7 @@ module MOM 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() ! These are used for group halo updates. type(group_pass_type) :: pass_tau_ustar_psurf @@ -385,6 +400,7 @@ module MOM public initialize_MOM public finish_MOM_initialization public step_MOM +public step_tracers public MOM_end public calculate_surface_state @@ -469,6 +485,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_predia, eta_preale + real :: tot_wt_ssh, Itot_wt_ssh, I_time_int real :: zos_area_mean, volo, ssh_ga type(time_type) :: Time_local @@ -898,12 +915,19 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) endif ! -------------------------------------------------- end SPLIT + + if (CS%thickness_diffuse .and. .not.CS%thickness_diffuse_first) then call cpu_clock_begin(id_clock_thick_diff) + + if (CS%debug) call hchksum(h,"Pre-thickness_diffuse h", G%HI, haloshift=0) + 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, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) + + if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1) call cpu_clock_end(id_clock_thick_diff) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_h, G%Domain) @@ -977,6 +1001,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call enable_averaging(CS%dt_trans, Time_local, CS%diag) call cpu_clock_begin(id_clock_tracer) + call advect_tracer(h, CS%uhtr, CS%vhtr, CS%OBC, CS%dt_trans, G, GV, & CS%tracer_adv_CSp, CS%tracer_Reg) call tracer_hordiff(h, CS%dt_trans, CS%MEKE, CS%VarMix, G, GV, & @@ -1030,10 +1055,19 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! Regridding/remapping is done here, at the end of the thermodynamics time step ! (that may comprise several dynamical time steps) ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'. + if (CS%id_u_preale > 0) call post_data(CS%id_u_preale, u, CS%diag) + if (CS%id_v_preale > 0) call post_data(CS%id_v_preale, v, CS%diag) + if (CS%id_h_preale > 0) call post_data(CS%id_h_preale, h, CS%diag) + if (CS%id_T_preale > 0) call post_data(CS%id_T_preale, CS%tv%T, CS%diag) + if (CS%id_S_preale > 0) call post_data(CS%id_S_preale, CS%tv%S, CS%diag) + if (CS%id_e_preale > 0) then + call find_eta(h, CS%tv, G%g_Earth, G, GV, eta_preale) + call post_data(CS%id_e_preale, eta_preale, CS%diag) + endif + if ( CS%use_ALE_algorithm ) then ! call pass_vector(u, v, G%Domain) call do_group_pass(CS%pass_T_S_h, G%Domain) - ! update squared quantities if (associated(CS%S_squared)) & CS%S_squared(:,:,:) = CS%tv%S(:,:,:) ** 2 @@ -1155,6 +1189,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%id_Sdiffx_2d > 0) call post_data(CS%id_Sdiffx_2d, CS%S_diffx_2d, CS%diag) if (CS%id_Sdiffy_2d > 0) call post_data(CS%id_Sdiffy_2d, CS%S_diffy_2d, CS%diag) + if (CS%id_uhtr > 0) call post_data(CS%id_uhtr, CS%uhtr, CS%diag) + if (CS%id_vhtr > 0) call post_data(CS%id_vhtr, CS%vhtr, CS%diag) + call post_diags_TS_tendency(G,GV,CS,dtdia) call disable_averaging(CS%diag) @@ -1174,6 +1211,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call cpu_clock_end(id_clock_other) call cpu_clock_begin(id_clock_thermo) + + + ! Reset the accumulated transports to 0. CS%uhtr(:,:,:) = 0.0 CS%vhtr(:,:,:) = 0.0 @@ -1364,16 +1404,431 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) end subroutine step_MOM +!> step_tracers is the main driver for running tracers offline in MOM6. This has been primarily +!! 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_tracers(fluxes, state, Time_start, time_interval, CS) + type(forcing), intent(inout) :: fluxes !< pointers to forcing fields + type(surface), intent(inout) :: 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_control_struct), pointer :: CS !< control structure from initialize_MOM + + ! Local pointers + type(ocean_grid_type), pointer :: G => NULL() ! Pointer to a structure containing + ! metrics and related information + type(verticalGrid_type), pointer :: GV => NULL() ! Pointer to structure containing information + ! about the vertical grid + ! Zonal mass transports + real, dimension(SZIB_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: uhtr, uhtr_sub + ! Zonal diffusive transport + real, dimension(SZIB_(CS%G),SZJ_(CS%G)) :: khdt_x + ! Meridional mass transports + real, dimension(SZI_(CS%G),SZJB_(CS%G),SZK_(CS%G)) :: vhtr, vhtr_sub + ! Meridional diffusive transports + real, dimension(SZI_(CS%G),SZJB_(CS%G)) :: khdt_y + + real :: sum_abs_fluxes, sum_u, sum_v ! Used to keep track of how close to convergence we are + real :: dt_offline, minimum_forcing_depth, evap_CFL_limit ! Shorthand variables from offline CS + + ! Local variables + ! Vertical diffusion related variables + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: & + eatr, & ! Amount of fluid entrained from the layer above within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + ebtr, & ! Amount of fluid entrained from the layer below within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + eatr_sub, & + ebtr_sub + ! Variables used to keep track of layer thicknesses at various points in the code + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: & + h_new, & + h_end, & + h_vol, & + h_pre, & + h_temp + ! Work arrays for temperature and salinity + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: & + temp_old, salt_old, & + zero_3dh ! + integer :: niter, iter + real :: Inum_iter, dt_iter + integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: isv, iev, jsv, jev ! The valid range of the indices. + integer :: IsdB, IedB, JsdB, JedB + logical :: z_first, x_before_y + + ! Fail out if do_online is true + if(CS%do_online) call MOM_error(FATAL,"DO_ONLINE=True when calling step_tracers") + + ! Grid-related pointer assignments + G => CS%G + GV => CS%GV + + ! Initialize some shorthand variables from other structures + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + dt_offline = CS%offline_CSp%dt_offline + evap_CFL_limit = CS%offline_CSp%evap_CFL_limit + minimum_forcing_depth = CS%offline_CSp%minimum_forcing_depth + + niter = CS%offline_CSp%num_off_iter + Inum_iter = 1./real(niter) + dt_iter = dt_offline*Inum_iter + + ! Initialize working arrays + uhtr(:,:,:) = 0.0 + vhtr(:,:,:) = 0.0 + khdt_x(:,:) = 0.0 + khdt_y(:,:) = 0.0 + eatr(:,:,:) = 0.0 + ebtr(:,:,:) = 0.0 + h_pre(:,:,:) = GV%Angstrom + h_new(:,:,:) = GV%Angstrom + h_end(:,:,:) = GV%Angstrom + temp_old(:,:,:) = 0.0 + salt_old(:,:,:) = 0.0 + uhtr_sub(:,:,:) = 0.0 + vhtr_sub(:,:,:) = 0.0 + eatr_sub(:,:,:) = 0.0 + ebtr_sub(:,:,:) = 0.0 + + call cpu_clock_begin(id_clock_tracer) + call enable_averaging(time_interval, Time_start+set_time(int(time_interval)), CS%diag) + + ! Read in all fields that might be used this timestep + call transport_by_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, & + khdt_x, khdt_y, temp_old, salt_old, fluxes, CS%use_ALE_algorithm) + + ! Set the starting layer thicknesses to those from the previous timestep + do k=1,nz ; do j=jsd,jed ; do i=isd,ied + h_pre(i,j,k) = CS%h(i,j,k) + enddo ; enddo; enddo + call pass_var(h_pre,G%Domain) + + x_before_y = (MOD(G%first_direction,2) == 0) + z_first = CS%diabatic_first + + if(CS%use_ALE_algorithm) then + + ! Tracers are transported using the stored mass fluxes. Where possible, operators are Strang-split around + ! the call to + ! 1) Using the layer thicknesses and tracer concentrations from the previous timestep, + ! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to tracer_column_fns. + ! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline + ! 2) Half of the accumulated surface freshwater fluxes are applied + !! START ITERATION + ! 3) Accumulated mass fluxes are used to do horizontal transport. The number of iterations used in + ! advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are stored for later use + ! and resulting layer thicknesses fed into the next step + ! 4) Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for layers which might + ! 'vanish' because of horizontal mass transport to be 'reinflated' + ! 5) Check that transport is done if the remaining mass fluxes equals 0 or if the max number of iterations + ! has been reached + !! END ITERATION + ! 6) Repeat steps 1 and 2 + ! 7) Force a remapping to the stored layer thicknesses that correspond to the snapshot of the online model + ! 8) Reset T/S and h to their stored snapshotted values to prevent model drift + + ! Convert flux rates into explicit mass/height of freshwater flux. Also note, that + ! fluxes are halved because diabatic processes are split before and after advection + + ! Do horizontal diffusion first (but only half of it), remainder will be applied after advection + call tracer_hordiff(h_pre, CS%offline_CSp%dt_offline*0.5, CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv, do_online_flag=CS%do_online, read_khdt_x=khdt_x*0.5, & + read_khdt_y=khdt_y*0.5) + + do j=jsd,jed ; do i=isd,ied + fluxes%netMassOut(i,j) = 0.5*fluxes%netMassOut(i,j) + fluxes%netMassIn(i,j) = 0.5*fluxes%netMassIn(i,j) + enddo ; enddo + + zero_3dh(:,:,:)=0.0 + + ! Copy over the horizontal mass fluxes from the remaining total mass fluxes + do k=1,nz ; do j=jsd,jed ; do i=isdB,iedB + uhtr_sub(i,j,k) = uhtr(i,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do j=jsdB,jedB ; do i=isd,ied + vhtr_sub(i,j,k) = vhtr(i,j,k) + enddo ; enddo ; enddo + + if(CS%debug) then + call uchksum(uhtr_sub,"uhtr_sub before transport",G%HI) + call vchksum(vhtr_sub,"vhtr_sub before transport",G%HI) + call hchksum(h_pre,"h_pre before transport",G%HI) + endif + + ! Note that here, h_new does nto represent any physical, should double check that any individual + ! tracer does not use h_new + call call_tracer_column_fns(h_pre, h_new, eatr*0.5, ebtr*0.5, & + fluxes, CS%offline_CSp%dt_offline*0.5, G, GV, CS%tv, & + CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) + ! Add half of the total freshwater fluxes + call applyTracerBoundaryFluxesInOut(G, GV, zero_3dh, 0.5*dt_offline, fluxes, h_pre, & + evap_CFL_limit, minimum_forcing_depth) + + if(CS%debug) then + call hchksum(h_pre,"h_pre after 1st diabatic",G%HI) + endif + + ! This loop does essentially a flux-limited, nonlinear advection scheme until all mass fluxes + ! are used. ALE is done after the horizontal advection. + do iter=1,CS%offline_CSp%num_off_iter + + do k=1,nz ; do j=jsd,jed ; do i=isd,ied + h_vol(i,j,k) = h_pre(i,j,k)*G%areaT(i,j) + enddo ; enddo ; enddo + + call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & + CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=2, & + uhr_out=uhtr, vhr_out=vhtr, h_out=h_new, x_first_in=x_before_y) + ! Switch the direction every iteration? Maybe not useful + ! x_before_y = .not. x_before_y + + do k=1,nz ; do j=jsd,jed ; do i=isd,ied + h_pre(i,j,k) = h_new(i,j,k)/G%areaT(i,j) + enddo ; enddo ; enddo + + + if(CS%debug) then + call hchksum(h_pre,"h_pre after advect_tracer",G%HI) + endif + + call cpu_clock_begin(id_clock_ALE) + call ALE_main_offline(G, GV, h_pre, CS%tv, & + CS%tracer_Reg, CS%ALE_CSp, CS%offline_CSp%dt_offline) + call cpu_clock_end(id_clock_ALE) + + if(CS%debug) then + call hchksum(h_pre,"h_pre after ALE",G%HI) + endif + + do k=1,nz ; do j=jsd,jed ; do i=isdB,iedB + uhtr_sub(i,j,k) = uhtr(i,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do j=jsdB,jedB ; do i=isd,ied + vhtr_sub(i,j,k) = vhtr(i,j,k) + enddo ; enddo ; enddo + + call pass_vector(uhtr_sub,vhtr_sub,G%Domain) + call pass_var(h_pre, G%Domain) + + if(CS%debug) then + call uchksum(uhtr_sub,"uhtr_sub after adv iteration",G%HI) + call vchksum(vhtr_sub,"vhtr_sub after adv iteration",G%HI) + call hchksum(h_pre,"h_pre after adv iteration",G%HI) + endif + + sum_u = 0.0 + do k=1,nz; do j=js,je ; do i=is-1,ie + sum_u = sum_u + abs(uhtr_sub(i,j,k)) + enddo; enddo; enddo + sum_v = 0.0 + do k=1,nz; do j=js-1,je; do i=is,ie + sum_v = sum_v + abs(vhtr_sub(i,j,k)) + enddo; enddo ; enddo + + call sum_across_PEs(sum_u) + call sum_across_PEs(sum_v) + + if(sum_u+sum_v==0.0) then + if(is_root_pe()) print *, "Converged after iteration", iter + exit +! print *, "Remaining uflux, vflux:", sum(abs(uhtr)), sum(abs(vhtr)) + endif + + enddo + + ! Now do the other half of the vertical mixing and tracer source/sink functions + call call_tracer_column_fns(h_pre, h_new, eatr*0.5, ebtr*0.5, & + fluxes, CS%offline_CSp%dt_offline*0.5, G, GV, CS%tv, & + CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) + call applyTracerBoundaryFluxesInOut(G, GV, zero_3dh, 0.5*dt_offline, fluxes, h_pre, & + evap_CFL_limit, minimum_forcing_depth) + + if(CS%debug) then + call hchksum(h_pre,"h_pre after 2nd diabatic",G%HI) + endif + + ! 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, h_pre, h_end, CS%tracer_Reg, CS%ALE_CSp) + call cpu_clock_end(id_clock_ALE) + + ! Finish with the other half of the tracer horizontal diffusion + call tracer_hordiff(h_pre, CS%offline_CSp%dt_offline*0.5, CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv, do_online_flag=.false., read_khdt_x=khdt_x*0.5, & + read_khdt_y=khdt_y*0.5) + + elseif (.not. CS%use_ALE_algorithm) then + do iter=1,CS%offline_CSp%num_off_iter + + do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + eatr_sub(i,j,k) = eatr(i,j,k) + ebtr_sub(i,j,k) = ebtr(i,j,k) + enddo; enddo ; enddo + + do k = 1, nz ; do j=js-1,je+1 ; do i=is-2,ie+1 + uhtr_sub(I,j,k) = uhtr(I,j,k) + enddo; enddo ; enddo + + do k = 1, nz ; do j=js-2,je+1 ; do i=is-1,ie+1 + vhtr_sub(i,J,k) = vhtr(i,J,k) + enddo; enddo ; enddo + + + ! Calculate 3d mass transports to be used in this iteration + call limit_mass_flux_3d(G, GV, uhtr_sub, vhtr_sub, eatr_sub, ebtr_sub, h_pre, & + CS%offline_CSp%max_off_cfl) + + if (z_first) then + ! First do vertical advection + call update_h_vertical_flux(G, GV, eatr_sub, ebtr_sub, h_pre, h_new) + call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, & + fluxes, dt_iter, G, GV, CS%tv, CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug) + ! We are now done with the vertical mass transports, so now h_new is h_sub + do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + h_pre(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + call pass_var(h_pre,G%Domain) + + ! Second zonal and meridional advection + call update_h_horizontal_flux(G, GV, uhtr_sub, vhtr_sub, h_pre, h_new) + do k = 1, nz ; do i = is-1, ie+1 ; do j=js-1, je+1 + h_vol(i,j,k) = h_pre(i,j,k)*G%areaT(i,j) + enddo; enddo; enddo + call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & + CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) + + ! Done with horizontal so now h_pre should be h_new + do k = 1, nz ; do i=is-1,ie+1 ; do j=js-1,je+1 + h_pre(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + + endif + + if (.not. z_first) then + + ! First zonal and meridional advection + call update_h_horizontal_flux(G, GV, uhtr_sub, vhtr_sub, h_pre, h_new) + do k = 1, nz ; do i = is-1, ie+1 ; do j=js-1, je+1 + h_vol(i,j,k) = h_pre(i,j,k)*G%areaT(i,j) + enddo; enddo; enddo + call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & + CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) + + ! Done with horizontal so now h_pre should be h_new + do k = 1, nz ; do i=is-1,ie+1 ; do j=js-1,je+1 + h_pre(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + + ! Second vertical advection + call update_h_vertical_flux(G, GV, eatr_sub, ebtr_sub, h_pre, h_new) + call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, & + fluxes, dt_iter, G, GV, CS%tv, CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug) + ! We are now done with the vertical mass transports, so now h_new is h_sub + do k = 1, nz ; do i=is-1,ie+1 ; do j=js-1,je+1 + h_pre(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + + + endif + + ! Update remaining transports + do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + eatr(i,j,k) = eatr(i,j,k) - eatr_sub(i,j,k) + ebtr(i,j,k) = ebtr(i,j,k) - ebtr_sub(i,j,k) + enddo; enddo ; enddo + + do k = 1, nz ; do j=js-1,je+1 ; do i=is-2,ie+1 + uhtr(I,j,k) = uhtr(I,j,k) - uhtr_sub(I,j,k) + enddo; enddo ; enddo + + do k = 1, nz ; do j=js-2,je+1 ; do i=is-1,ie+1 + vhtr(i,J,k) = vhtr(i,J,k) - vhtr_sub(i,J,k) + enddo; enddo ; enddo + + call pass_var(eatr,G%Domain) + call pass_var(ebtr,G%Domain) + call pass_var(h_pre,G%Domain) + call pass_vector(uhtr,vhtr,G%Domain) + ! + ! Calculate how close we are to converging by summing the remaining fluxes at each point + sum_abs_fluxes = 0.0 + sum_u = 0.0 + sum_v = 0.0 + do k=1,nz; do j=js,je; do i=is,ie + sum_u = sum_u + abs(uhtr(I-1,j,k))+abs(uhtr(I,j,k)) + sum_v = sum_v + abs(vhtr(i,J-1,k))+abs(vhtr(I,J,k)) + sum_abs_fluxes = sum_abs_fluxes + abs(eatr(i,j,k)) + abs(ebtr(i,j,k)) + abs(uhtr(I-1,j,k)) + & + abs(uhtr(I,j,k)) + abs(vhtr(i,J-1,k)) + abs(vhtr(i,J,k)) + enddo; enddo; enddo + call sum_across_PEs(sum_abs_fluxes) + + print *, "Remaining u-flux, v-flux:", sum_u, sum_v + if (sum_abs_fluxes==0) then + print *, 'Converged after iteration', iter + exit + endif + + ! Switch order of Strang split every iteration + z_first = .not. z_first + x_before_y = .not. x_before_y + + end do + call tracer_hordiff(h_end, CS%offline_CSp%dt_offline*0.5, CS%MEKE, CS%VarMix, G, GV, & + CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv, do_online_flag=.false., read_khdt_x=khdt_x*0.5, & + read_khdt_y=khdt_y*0.5) + endif + + h_temp = h_end-h_pre + + if (CS%offline_CSp%id_hr>0) call post_data(CS%offline_CSp%id_hr, h_temp, CS%diag) + if (CS%offline_CSp%id_uhr>0) call post_data(CS%offline_CSp%id_uhr, uhtr, CS%diag) + if (CS%offline_CSp%id_vhr>0) call post_data(CS%offline_CSp%id_vhr, vhtr, CS%diag) + if (CS%offline_CSp%id_ear>0) call post_data(CS%offline_CSp%id_ear, eatr, CS%diag) + if (CS%offline_CSp%id_ebr>0) call post_data(CS%offline_CSp%id_ebr, ebtr, CS%diag) + + call cpu_clock_end(id_clock_tracer) + + call disable_averaging(CS%diag) + + ! Note here T/S are reset to the stored snap shot to ensure that the offline model + ! densities, used in the neutral diffusion code don't drift too far from the online + ! model + do i = is, ie ; do j = js, je ; do k=1,nz + CS%T(i,j,k) = temp_old(i,j,k) + CS%S(i,j,k) = salt_old(i,j,k) + CS%h(i,j,k) = h_end(i,j,k) + enddo ; enddo; enddo + + call pass_var(CS%h,G%Domain) + call pass_var(CS%T,G%Domain) + call pass_var(CS%S,G%Domain) + + + +end subroutine step_tracers !> This subroutine initializes MOM. -subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) +subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, do_online_out) type(time_type), target, intent(inout) :: Time !< model time, set in this routine 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_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure 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) :: do_online_out !< .false. if tracers are being run offline ! local type(ocean_grid_type), pointer :: G => NULL() ! A pointer to a structure with metrics and related @@ -1486,7 +1941,15 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "If False, skips the dynamics calls that update u & v, as well as\n"//& "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", "DO_ONLINE", CS%do_online, & + "If false, use the model in prognostic mode where\n"//& + "the barotropic and baroclinic dynamics, thermodynamics,\n"//& + "etc. are stepped forward integrated in time.\n"//& + "If true, the all of the above are bypassed with all\n"//& + "fields necessary to integrate only the tracer advection\n"//& + "and diffusion equation are read in from files stored from\n"//& + "a previous integration of the prognostic model\n"//& + "NOTE: This option only used in the ocean_solo_driver.", default=.true.) call get_param(param_file, "MOM", "USE_REGRIDDING", CS%use_ALE_algorithm , & "If True, use the ALE algorithm (regridding/remapping).\n"//& "If False, use the layered isopycnal algorithm.", default=.false. ) @@ -1501,8 +1964,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call get_param(param_file, "MOM", "THICKNESSDIFFUSE", CS%thickness_diffuse, & "If true, interface heights are diffused with a \n"//& "coefficient of KHTH.", default=.false.) - call get_param(param_file, "MOM", "THICKNESSDIFFUSE_FIRST", & - CS%thickness_diffuse_first, & + call get_param(param_file, "MOM", "THICKNESSDIFFUSE_FIRST", & + CS%thickness_diffuse_first, & "If true, do thickness diffusion before dynamics.\n"//& "This is only used if THICKNESSDIFFUSE is true.", & default=.false.) @@ -2050,6 +2513,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif + ! If need a diagnostic field, then would have been allocated in register_diags. if (CS%use_temperature) then call add_tracer_diagnostics("T", CS%tracer_Reg, CS%T_adx, CS%T_ady, & @@ -2068,6 +2532,14 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) cmor_long_name ="Sea Water Salinity") endif + ! If running in offline tracer mode, initialize the necessary control structure and + ! parameters + if(present(do_online_out)) do_online_out=CS%do_online + + if(.not. CS%do_online) then + call offline_transport_init(param_file, CS%offline_CSp, CS%diabatic_CSp%diabatic_aux_CSp, G, GV) + call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp) + endif ! This subroutine initializes any tracer packages. new_sim = ((dirs%input_filename(1:1) == 'n') .and. & @@ -2376,7 +2848,7 @@ subroutine register_diags(Time, G, GV, CS, ADp) CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, & 'Meridional velocity before remapping', 'meter second-1') CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, & - 'Layer Thickness before remapping', thickness_units, v_cell_method='sum') + 'Layer Thickness before remapping', thickness_units) CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, & 'Temperature before remapping', 'degC') CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, & @@ -2391,7 +2863,13 @@ subroutine register_diags(Time, G, GV, CS, ADp) CS%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, Time, & 'Salinity', 'PPT') endif - + + ! Diagnostics related to tracer transport + CS%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, Time, & + 'Accumulated zonal thickness fluxes to advect tracers', 'kg') + CS%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, Time, & + 'Accumulated meridional thickness fluxes to advect tracers', 'kg') + end subroutine register_diags diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 11ecb7a3db..f7a1aac006 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -175,6 +175,7 @@ module MOM_forcing_type integer :: id_lprec = -1, id_fprec = -1 integer :: id_lrunoff = -1, id_frunoff = -1 integer :: id_net_massout = -1, id_net_massin = -1 + integer :: id_massout_flux = -1, id_massin_flux = -1 integer :: id_seaice_melt = -1 ! global area integrated mass flux diagnostic handles @@ -1055,7 +1056,12 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles) handles%id_net_massin = register_diag_field('ocean_model', 'net_massin', diag%axesT1, Time, & 'Net mass entering ocean due to precip, runoff, ice melt', 'kilogram meter-2 second-1') + handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, & + 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', & + 'kilogram meter-2') + handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, & + 'Net mass mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kilogram meter-2') !========================================================================= ! area integrated surface mass transport @@ -1796,6 +1802,8 @@ subroutine forcing_diagnostics(fluxes, state, dt, G, diag, handles) call post_data(handles%id_total_net_massout, total_transport, diag) endif endif + + if(handles%id_massout_flux > 0) call post_data(handles%id_massout_flux,fluxes%netMassOut,diag) if(handles%id_net_massin > 0 .or. handles%id_total_net_massin > 0) then sum(:,:) = 0.0 @@ -1812,6 +1820,8 @@ subroutine forcing_diagnostics(fluxes, state, dt, G, diag, handles) call post_data(handles%id_total_net_massin, total_transport, diag) endif endif + + if(handles%id_massin_flux > 0) call post_data(handles%id_massin_flux,fluxes%netMassIn,diag) if ((handles%id_evap > 0) .and. ASSOCIATED(fluxes%evap)) & call post_data(handles%id_evap, fluxes%evap, diag) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index d93b13ca53..d7379f7a0e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1043,6 +1043,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, & endif ! Change in state due to forcing + dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) ) dTemp = fractionOfForcing*netHeat(i) ! ### The 0.9999 here should become a run-time parameter? @@ -1071,6 +1072,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, & ! Update state by the appropriate increment. hOld = h2d(i,k) ! Keep original thickness in hand h2d(i,k) = h2d(i,k) + dThickness ! New thickness + if (h2d(i,k) > 0.) then if (calculate_energetics) then ! Calculate the energy required to mix the newly added water over diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 521beaf7ad..0ffc550558 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -62,6 +62,7 @@ module MOM_diabatic_driver use MOM_wave_speed, only : wave_speeds use time_manager_mod, only : increment_time ! for testing itides (BDM) + implicit none ; private #include @@ -73,7 +74,7 @@ module MOM_diabatic_driver public adiabatic_driver_init !> Control structure for this module -type, public :: diabatic_CS ; private +type, public :: diabatic_CS ; logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with !! nkml sublayers (and additional buffer layers). logical :: use_energetic_PBL !< If true, use the implicit energetics planetary @@ -742,7 +743,6 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & h, tv, CS%aggregate_FW_forcing, cTKE, dSV_dT, dSV_dS) - if (CS%debug) then call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0) call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0) @@ -1136,9 +1136,9 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) endif ; endif enddo ; enddo do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo + enddo - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied ! so hold should be h_orig @@ -1170,6 +1170,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent eatr(i,j,k) = ea(i,j,k) + add_ent enddo ; enddo ; enddo + if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, dt, G, GV, tv, & @@ -1195,6 +1196,8 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) endif ! (CS%mix_boundary_tracers) + + call cpu_clock_end(id_clock_tracers) @@ -1404,8 +1407,9 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) - if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag) if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag) if (CS%id_wd > 0) call post_data(CS%id_wd, CDp%diapyc_vel, CS%diag) diff --git a/src/tracer/MOM_offline_control.F90 b/src/tracer/MOM_offline_control.F90 new file mode 100644 index 0000000000..2d971e2b1e --- /dev/null +++ b/src/tracer/MOM_offline_control.F90 @@ -0,0 +1,566 @@ +!> Contains routines related to offline transport of tracers +module MOM_offline_transport +! This file is part of MOM6. See LICENSE.md for the license. + + use data_override_mod, only : data_override_init, data_override + use MOM_time_manager, only : time_type + use MOM_domains, only : pass_var, pass_vector, To_All + use MOM_error_handler, only : callTree_enter, callTree_leave, MOM_error, FATAL, WARNING, is_root_pe + use MOM_grid, only : ocean_grid_type + use MOM_verticalGrid, only : verticalGrid_type + use MOM_io, only : read_data + use MOM_file_parser, only : get_param, log_version, param_file_type + use MOM_diag_mediator, only : diag_ctrl, register_diag_field + use mpp_domains_mod, only : CENTER, CORNER, NORTH, EAST + use MOM_variables, only : vertvisc_type + use MOM_forcing_type, only : forcing + use MOM_shortwave_abs, only : optics_type + use MOM_diag_mediator, only : post_data + use MOM_forcing_type, only : forcing + use MOM_diabatic_aux, only : diabatic_aux_CS + + implicit none + +#include + + type, public :: offline_transport_CS + + !> Variables related to reading in fields from online run + integer :: start_index ! Timelevel to start + integer :: numtime ! How many timelevels in the input fields + integer :: & ! Index of each of the variables to be read in + ridx_sum = -1, & ! Separate indices for each variabile if they are + ridx_snap = -1 ! setoff from each other in time + character(len=200) :: offlinedir ! Directory where offline fields are stored + character(len=200) :: & ! ! Names of input files + snap_file, & + sum_file + logical :: fields_are_offset ! True if the time-averaged fields and snapshot fields are + ! offset by one time level + + !> Variables controlling some of the numerical considerations of offline transport + integer :: num_off_iter + real :: dt_offline ! Timestep used for offline tracers + real :: max_off_cfl=0.5 ! Hardcoded for now, only used in non-ALE mode + real :: evap_CFL_limit, minimum_forcing_depth + + !> Diagnostic manager IDs for use in the online model of additional fields necessary + !> for offline tracer modeling + integer :: & + id_uhtr_preadv = -1, & + id_vhtr_preadv = -1, & + !> Diagnostic manager IDs for some fields that may be of interest when doing offline transport + id_uhr = -1, & + id_vhr = -1, & + id_ear = -1, & + id_ebr = -1, & + id_hr = -1 + + end type offline_transport_CS + +#include "MOM_memory.h" +#include "version_variable.h" + public offline_transport_init + public transport_by_files + public register_diags_offline_transport + public update_h_horizontal_flux + public update_h_vertical_flux + public limit_mass_flux_3d + +contains + !> Controls the reading in 3d mass fluxes, diffusive fluxes, and other fields stored + !! in a previous integration of the online model + subroutine transport_by_files(G, GV, CS, h_end, eatr, ebtr, uhtr, vhtr, khdt_x, khdt_y, & + temp, salt, fluxes, do_ale_in) + + type(ocean_grid_type), intent(inout) :: G + type(verticalGrid_type), intent(inout) :: GV + type(offline_transport_CS), intent(inout) :: CS + logical, optional :: do_ale_in + + !! Mandatory variables + ! Fields at U-points + ! 3D + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uhtr + ! 2D + real, dimension(SZIB_(G),SZJ_(G)) :: khdt_x + ! Fields at V-points + ! 3D + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vhtr + ! 2D + real, dimension(SZI_(G),SZJB_(G)) :: khdt_y + ! Fields at T-point + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & + h_end, & + eatr, ebtr, & + temp, salt + type(forcing) :: fluxes + logical :: do_ale + integer :: i, j, k, is, ie, js, je, nz + + do_ale = .false.; + if (present(do_ale_in) ) do_ale = do_ale_in + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + + call callTree_enter("transport_by_files, MOM_offline_control.F90") + + + !! Time-summed fields + ! U-grid + call read_data(CS%sum_file, 'uhtr_sum', uhtr,domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=EAST) + call read_data(CS%sum_file, 'khdt_x_sum', khdt_x,domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=EAST) + ! V-grid + call read_data(CS%sum_file, 'vhtr_sum', vhtr, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=NORTH) + call read_data(CS%sum_file, 'khdt_y_sum', khdt_y, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=NORTH) + ! T-grid + call read_data(CS%sum_file, 'ea_sum', eatr, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=CENTER) + call read_data(CS%sum_file, 'eb_sum', ebtr, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=CENTER) + + !! Time-averaged fields + call read_data(CS%snap_file, 'temp', temp, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=CENTER) + call read_data(CS%snap_file, 'salt', salt, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_sum,position=CENTER) + + !! Read snapshot fields (end of time interval timestamp) + call read_data(CS%snap_file, 'h_end', h_end, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_snap,position=CENTER) + + ! Apply masks at T, U, and V points + do k=1,nz ; do j=js,je ; do i=is,ie + if(G%mask2dT(i,j)<1.0) then + h_end(i,j,k) = GV%Angstrom + eatr(i,j,k) = 0.0 + ebtr(i,j,k) = 0.0 + endif + enddo; enddo ; enddo + + do k=1,nz ; do j=js-1,je ; do i=is,ie + if(G%mask2dCv(i,j)<1.0) then + khdt_y(i,j) = 0.0 + vhtr(i,j,k) = 0.0 + endif + enddo; enddo ; enddo + + do k=1,nz ; do j=js,je ; do i=is-1,ie + if(G%mask2dCu(i,j)<1.0) then + khdt_x(i,j) = 0.0 + uhtr(i,j,k) = 0.0 + endif + enddo; enddo ; enddo + + ! This block makes sure that the fluxes control structure, which may not be used in the solo_driver, + ! contains netMassIn and netMassOut which is necessary for the applyTracerBoundaryFluxesInOut routine + if (do_ale) then + if (.not. ASSOCIATED(fluxes%netMassOut)) then + ALLOCATE(fluxes%netMassOut(G%isd:G%ied,G%jsd:G%jed)) + fluxes%netMassOut(:,:) = 0.0 + endif + if (.not. ASSOCIATED(fluxes%netMassIn)) then + ALLOCATE(fluxes%netMassIn(G%isd:G%ied,G%jsd:G%jed)) + fluxes%netMassIn(:,:) = 0.0 + endif + + call read_data(CS%sum_file,'massout_flux_sum',fluxes%netMassOut, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_snap,position=center) + call read_data(CS%sum_file,'massin_flux_sum', fluxes%netMassIn, domain=G%Domain%mpp_domain, & + timelevel=CS%ridx_snap,position=center) + endif + + !! Make sure all halos have been updated + ! Vector fields + call pass_vector(uhtr, vhtr, G%Domain) + call pass_vector(khdt_x, khdt_y, G%Domain) + + ! Scalar fields + call pass_var(h_end, G%Domain) + call pass_var(eatr, G%Domain) + call pass_var(ebtr, G%Domain) + call pass_var(temp, G%Domain) + call pass_var(salt, G%Domain) + + if (do_ale) then + call pass_var(fluxes%netMassOut,G%Domain) + call pass_var(fluxes%netMassIn,G%Domain) + endif + + ! Update the read indices + CS%ridx_snap = next_modulo_time(CS%ridx_snap,CS%numtime) + CS%ridx_sum = next_modulo_time(CS%ridx_sum,CS%numtime) + + call callTree_leave("transport_by_file") + + end subroutine transport_by_files + + !> Initialize additional diagnostics required for offline tracer transport + subroutine register_diags_offline_transport(Time, diag, CS) + + type(offline_transport_CS), pointer :: CS !< control structure for MOM + type(time_type), intent(in) :: Time !< current model time + type(diag_ctrl) :: diag + + + ! U-cell fields + CS%id_uhr = register_diag_field('ocean_model', 'uhr', diag%axesCuL, Time, & + 'Zonal thickness fluxes remaining at end of timestep', 'kg') + + ! V-cell fields + CS%id_vhr = register_diag_field('ocean_model', 'vhr', diag%axesCvL, Time, & + 'Meridional thickness fluxes remaining at end of timestep', 'kg') + + ! T-cell fields + CS%id_hr = register_diag_field('ocean_model', 'hdiff', diag%axesTL, Time, & + 'Difference between the stored and calculated layer thickness', 'm') + CS%id_ear = register_diag_field('ocean_model', 'ear', diag%axesTL, Time, & + 'Remaining thickness entrained from above', 'm') + CS%id_ebr = register_diag_field('ocean_model', 'ebr', diag%axesTL, Time, & + 'Remaining thickness entrained from below', 'm') + + end subroutine register_diags_offline_transport + + !> Initializes the control structure for offline transport and reads in some of the + ! run time parameters from MOM_input + subroutine offline_transport_init(param_file, CS, diabatic_aux_CSp, G, GV) + + type(param_file_type), intent(in) :: param_file + type(offline_transport_CS), pointer, intent(inout) :: CS + type(diabatic_aux_CS), pointer, intent(in) :: diabatic_aux_CSp + type(ocean_grid_type), intent(in) :: G + type(verticalGrid_type), intent(in) :: GV + + character(len=40) :: mod = "offline_transport" + + integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: IsdB, IedB, JsdB, JedB + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + call callTree_enter("offline_transport_init, MOM_offline_control.F90") + + if (associated(CS)) then + call MOM_error(WARNING, "offline_transport_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + call log_version(param_file,mod,version, & + "This module allows for tracers to be run offline") + + ! Parse MOM_input for offline control + call get_param(param_file, mod, "OFFLINEDIR", CS%offlinedir, & + "Input directory where the offline fields can be found", default=" ") + call get_param(param_file, mod, "OFF_SUM_FILE", CS%sum_file, & + "Filename where the accumulated fields can be found", default = " ") + call get_param(param_file, mod, "OFF_SNAP_FILE", CS%snap_file, & + "Filename where snapshot fields can be found",default=" ") + call get_param(param_file, mod, "START_INDEX", CS%start_index, & + "Which time index to start from", default=1) + call get_param(param_file, mod, "NUMTIME", CS%numtime, & + "Number of timelevels in offline input files", default=0) + call get_param(param_file, mod, "FIELDS_ARE_OFFSET", CS%fields_are_offset, & + "True if the time-averaged fields and snapshot fields are offset by one time level", & + default=.false.) + call get_param(param_file, mod, "NUM_OFF_ITER", CS%num_off_iter, & + "Number of iterations to subdivide the offline tracer advection and diffusion" ) + call get_param(param_file, mod, "DT_OFFLINE", CS%dt_offline, & + "Length of the offline timestep") + + ! Concatenate offline directory and file names + CS%snap_file = trim(CS%offlinedir)//trim(CS%snap_file) + CS%sum_file = trim(CS%offlinedir)//trim(CS%sum_file) + + ! Set the starting read index for time-averaged and snapshotted fields + CS%ridx_sum = CS%start_index + if(CS%fields_are_offset) CS%ridx_snap = next_modulo_time(CS%start_index,CS%numtime) + if(.not. CS%fields_are_offset) CS%ridx_snap = CS%start_index + + ! Copy over parameters from other control structures + if(associated(diabatic_aux_CSp)) then + CS%evap_CFL_limit = diabatic_aux_CSp%evap_CFL_limit + CS%minimum_forcing_depth = diabatic_aux_CSp%minimum_forcing_depth + endif + + call callTree_leave("offline_transport_init") + + end subroutine offline_transport_init + + !> Calculates the next timelevel to read from the input fields. This allows the 'looping' + !! of the fields + function next_modulo_time(inidx, numtime) + ! Returns the next time interval to be read + integer :: numtime ! Number of time levels in input fields + integer :: inidx ! The current time index + + integer :: read_index ! The index in the input files that corresponds + ! to the current timestep + + integer :: next_modulo_time + + read_index = mod(inidx+1,numtime) + if (read_index < 0) read_index = inidx-read_index + if (read_index == 0) read_index = numtime + + next_modulo_time = read_index + + end function next_modulo_time + + !> This updates thickness based on the convergence of horizontal mass fluxes + !! NOTE: Only used in non-ALE mode + subroutine update_h_horizontal_flux(G, GV, uhtr, vhtr, h_pre, h_new) + type(ocean_grid_type), pointer :: G + type(verticalGrid_type), pointer :: GV + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhtr + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhtr + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: h_pre + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h_new + + ! Local variables + integer :: i, j, k, m, is, ie, js, je, nz + ! Set index-related variables for fields on T-grid + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + do k = 1, nz + do i=is-1,ie+1 ; do j=js-1,je+1 + + h_new(i,j,k) = max(0.0, G%areaT(i,j)*h_pre(i,j,k) + & + ((uhtr(I-1,j,k) - uhtr(I,j,k)) + (vhtr(i,J-1,k) - vhtr(i,J,k)))) + + ! In the case that the layer is now dramatically thinner than it was previously, + ! add a bit of mass to avoid truncation errors. This will lead to + ! non-conservation of tracers + h_new(i,j,k) = h_new(i,j,k) + & + max(GV%Angstrom, 1.0e-13*h_new(i,j,k) - G%areaT(i,j)*h_pre(i,j,k)) + + ! Convert back to thickness + h_new(i,j,k) = h_new(i,j,k)/G%areaT(i,j) + + enddo ; enddo + enddo + end subroutine update_h_horizontal_flux + + !> Updates layer thicknesses due to vertical mass transports + !! NOTE: Only used in non-ALE configuration + subroutine update_h_vertical_flux(G, GV, ea, eb, h_pre, h_new) + type(ocean_grid_type), pointer :: G + type(verticalGrid_type), pointer :: GV + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: ea + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: eb + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: h_pre + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h_new + + ! Local variables + integer :: i, j, k, m, is, ie, js, je, nz + ! Set index-related variables for fields on T-grid + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + ! Update h_new with convergence of vertical mass transports + do j=js-1,je+1 + do i=is-1,ie+1 + + ! Top layer + h_new(i,j,1) = max(0.0, h_pre(i,j,1) + (eb(i,j,1) - ea(i,j,2) + ea(i,j,1) )) + h_new(i,j,1) = h_new(i,j,1) + & + max(0.0, 1.0e-13*h_new(i,j,1) - h_pre(i,j,1)) + + ! Bottom layer +! h_new(i,j,nz) = h_pre(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1)+eb(i,j,nz)) + h_new(i,j,nz) = max(0.0, h_pre(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1)+eb(i,j,nz))) + h_new(i,j,nz) = h_new(i,j,nz) + & + max(0.0, 1.0e-13*h_new(i,j,nz) - h_pre(i,j,nz)) + + enddo + + ! Interior layers + do k=2,nz-1 ; do i=is-1,ie+1 + + h_new(i,j,k) = max(0.0, h_pre(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + & + (eb(i,j,k) - ea(i,j,k+1)))) + h_new(i,j,k) = h_new(i,j,k) + & + max(0.0, 1.0e-13*h_new(i,j,k) - h_pre(i,j,k)) + + enddo ; enddo + + enddo + + end subroutine update_h_vertical_flux + + !> This routine limits the mass fluxes so that the a layer cannot be completely depleted. + !! NOTE: Only used in non-ALE mode + subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre, max_off_cfl) + type(ocean_grid_type), pointer :: G + type(verticalGrid_type), pointer :: GV + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: ea + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: eb + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: h_pre + real, intent(in) :: max_off_cfl + + ! Local variables + integer :: i, j, k, m, is, ie, js, je, nz + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: top_flux, bottom_flux + real :: pos_flux, hvol, h_neglect, scale_factor + + + ! In this subroutine, fluxes out of the box are scaled away if they deplete + ! the layer, note that we define the positive direction as flux out of the box. + ! Hence, uh(I-1) is multipled by negative one, but uh(I) is not + + ! Set index-related variables for fields on T-grid + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + ! Calculate top and bottom fluxes from ea and eb. Note the explicit negative signs + ! to enforce the positive out convention + k = 1 + do j=js-1,je+1 ; do i=is-1,ie+1 + top_flux(i,j,k) = -ea(i,j,k) + bottom_flux(i,j,k) = -(eb(i,j,k)-ea(i,j,k+1)) + enddo ; enddo + + do k=2, nz-1 ; do j=js-1,je+1 ; do i=is-1,ie+1 + top_flux(i,j,k) = -(ea(i,j,k)-eb(i,j,k-1)) + bottom_flux(i,j,k) = -(eb(i,j,k)-ea(i,j,k+1)) + enddo ; enddo ; enddo + + k=nz + do j=js-1,je+1 ; do i=is-1,ie+1 + top_flux(i,j,k) = -(ea(i,j,k)-eb(i,j,k-1)) + bottom_flux(i,j,k) = -eb(i,j,k) + enddo ; enddo + + + ! Calculate sum of positive fluxes (negatives applied to enforce convention) + ! in a given cell and scale it back if it would deplete a layer + do k = 1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + + hvol = h_pre(i,j,k)*G%areaT(i,j) + pos_flux = max(0.0,-uh(I-1,j,k)) + max(0.0, -vh(i,J-1,k)) + & + max(0.0, uh(I,j,k)) + max(0.0, vh(i,J,k)) + & + max(0.0, top_flux(i,j,k)*G%areaT(i,j)) + max(0.0, bottom_flux(i,j,k)*G%areaT(i,j)) + + if (pos_flux>hvol .and. pos_flux>0.0) then + scale_factor = ( hvol )/pos_flux*max_off_cfl + else ! Don't scale + scale_factor = 1.0 + endif + + ! Scale horizontal fluxes + if (-uh(I-1,j,k)>0) uh(I-1,j,k) = uh(I-1,j,k)*scale_factor + if (uh(I,j,k)>0) uh(I,j,k) = uh(I,j,k)*scale_factor + if (-vh(i,J-1,k)>0) vh(i,J-1,k) = vh(i,J-1,k)*scale_factor + if (vh(i,J,k)>0) vh(i,J,k) = vh(i,J,k)*scale_factor + + if (k>1 .and. k0.0) then + ea(i,j,k) = ea(i,j,k)*scale_factor + eb(i,j,k-1) = eb(i,j,k-1)*scale_factor + endif + if(bottom_flux(i,j,k)>0.0) then + eb(i,j,k) = eb(i,j,k)*scale_factor + ea(i,j,k+1) = ea(i,j,k+1)*scale_factor + endif + ! Scale top layer + elseif (k==1) then + if(top_flux(i,j,k)>0.0) ea(i,j,k) = ea(i,j,k)*scale_factor + if(bottom_flux(i,j,k)>0.0) then + eb(i,j,k) = eb(i,j,k)*scale_factor + ea(i,j,k+1) = ea(i,j,k+1)*scale_factor + endif + ! Scale bottom layer + elseif (k==nz) then + if(top_flux(i,j,k)>0.0) then + ea(i,j,k) = ea(i,j,k)*scale_factor + eb(i,j,k-1) = eb(i,j,k-1)*scale_factor + endif + if (bottom_flux(i,j,k)>0.0) eb(i,j,k)=eb(i,j,k)*scale_factor + endif + enddo ; enddo ; enddo + + end subroutine limit_mass_flux_3d + +!> \namespace mom_offline_transport +!! \section offline_overview Offline Tracer Transport in MOM6 +!! 'Offline tracer modeling' uses physical fields (e.g. mass transports and layer thicknesses) saved +!! from a previous integration of the physical model to transport passive tracers. These fields are +!! accumulated or averaged over a period of time (in this test case, 1 day) and used to integrate +!! portions of the MOM6 code base that handle the 3d advection and diffusion of passive tracers. +!! +!! The distribution of tracers in the ocean modeled offline should not be expected to match an online +!! simulation. Accumulating transports over more than one online model timestep implicitly assumes +!! homogeneity over that time period and essentially aliases over processes that occur with higher +!! frequency. For example, consider the case of a surface boundary layer with a strong diurnal cycle. +!! An offline simulation with a 1 day timestep, captures the net transport into or out of that layer, +!! but not the exact cycling. This effective aliasing may also complicate online model configurations +!! which strongly-eddying regions. In this case, the offline model timestep must be limited to some +!! fraction of the eddy correlation timescale. Lastly, the nonlinear advection scheme which applies +!! limited mass-transports over a sequence of iterations means that tracers are not transported along +!! exactly the same path as they are in the online model. +!! +!! This capability has currently targeted the Baltic_ALE_z test case, though some work has also been +!! done with the OM4 1/2 degree configuration. Work is ongoing to develop recommendations and best +!! practices for investigators seeking to use MOM6 for offline tracer modeling. +!! +!! \section offline_technical Implementation of offline routine in MOM6 +!! +!! The subroutine step_tracers that coordinates this can be found in MOM.F90 and is only called +!! using the solo ocean driver. This is to avoid issues with coupling to other climate components +!! that may be relying on fluxes from the ocean to be coupled more often than the offline time step. +!! Other routines related to offline tracer modeling can be found in tracers/MOM_offline_control.F90 +!! +!! As can also be seen in the comments for the step_tracers subroutine, an offline time step +!! comprises the following steps: +!! -# Using the layer thicknesses and tracer concentrations from the previous timestep, +!! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to +!! tracer_column_fns. +!! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline +!! -# Half of the accumulated surface freshwater fluxes are applied +!! START ITERATION +!! -# Accumulated mass fluxes are used to do horizontal transport. The number of iterations +!! used in advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are +!! stored for later use and resulting layer thicknesses fed into the next step +!! -# Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for +!! layers which might 'vanish' because of horizontal mass transport to be 'reinflated' +!! and essentially allows for the vertical transport of tracers +!! -# Check that transport is done if the remaining mass fluxes equals 0 or if the max +!! number of iterations has been reached +!! END ITERATION +!! -# Repeat steps 1 and 2 +!! -# Force a remapping to the stored layer thicknesses that correspond to the snapshot of +!! the online model at the end of an accumulation interval +!! -# Reset T/S and h to their stored snapshotted values to prevent model drift +!! +!! \section offline_evaluation Evaluating the utility of an offline tracer model +!! How well an offline tracer model can be used as an alternative to integrating tracers online +!! with the prognostic model must be evaluated for each application. This efficacy may be related +!! to the native coordinate of the online model, to the length of the offline timestep, and to the +!! behavior of the tracer itself. +!! +!! A framework for formally regression testing the offline capability still needs to be developed. +!! However, as a simple way of testing whether the offline model is nominally behaving as expected, +!! the total inventory of the advection test tracers (tr1, tr2, etc.) should be conserved between +!! time steps except for the last 4 decimal places. As a general guideline, an offline timestep of +!! 5 days or less. +!! +!! \section offline_parameters Runtime parameters for offline tracers +!! - OFFLINEDIR: Input directory where the offline fields can be found +!! - OFF_SUM_FILE: Filename where the accumulated fields can be found (e.g. horizontal mass transports) +!! - OFF_SNAP_FILE: Filename where snapshot fields can be found (e.g. end of timestep layer thickness) +!! - START_INDEX: Which timelevel of the input files to read first +!! - NUMTIME: How many timelevels to read before 'looping' back to 1 +!! - FIELDS_ARE_OFFSET: True if the time-averaged fields and snapshot fields are offset by one +!! time level, probably not needed +!! -NUM_OFF_ITER: Maximum number of iterations to do for the nonlinear advection scheme + +end module MOM_offline_transport + diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index b8668bd89b..b7b66f8b60 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -9,7 +9,7 @@ module MOM_tracer_advect use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type use MOM_domains, only : sum_across_PEs, max_across_PEs -use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type +use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type, pass_var use MOM_checksums, only : hchksum use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type @@ -18,7 +18,6 @@ module MOM_tracer_advect use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum use MOM_verticalGrid, only : verticalGrid_type - implicit none ; private #include @@ -45,7 +44,7 @@ module MOM_tracer_advect !> This routine time steps the tracer concentration using a !! monotonic, conservative, weakly diffusive scheme. -subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) +subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg, h_prev_opt, max_iter_in, x_first_in, uhr_out, vhr_out, h_out) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_end !< layer thickness after advection (m or kg m-2) @@ -54,8 +53,13 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used real, intent(in) :: dt !< time increment (seconds) type(tracer_advect_CS), pointer :: CS !< control structure for module - type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry - + type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional :: h_prev_opt !< layer thickness before advection (m or kg m-2) + integer, optional :: max_iter_in + logical, optional :: x_first_in + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: uhr_out !< accumulated volume/mass flux through zonal face (m3 or kg) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(out) :: vhr_out !< accumulated volume/mass flux through merid face (m3 or kg) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional :: h_out !< layer thickness before advection (m or kg m-2) type(tracer_type) :: Tr(MAX_FIELDS_) ! The array of registered tracers real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -104,6 +108,8 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) max_iter = 2*INT(CEILING(dt/CS%dt)) + 1 + if(present(max_iter_in)) max_iter = max_iter_in + if(present(x_first_in)) x_first = x_first_in call cpu_clock_begin(id_clock_pass) call create_group_pass(CS%pass_uhr_vhr_t_hprev, uhr, vhr, G%Domain) call create_group_pass(CS%pass_uhr_vhr_t_hprev, hprev, G%Domain) @@ -119,6 +125,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) ! This initializes the halos of uhr and vhr because pass_vector might do ! calculations on them, even though they are never used. !$OMP do + do k = 1, nz do j = jsd, jed; do i = IsdB, IedB; uhr(i,j,k) = 0.0; enddo ; enddo do j = jsdB, jedB; do i = Isd, Ied; vhr(i,j,k) = 0.0; enddo ; enddo @@ -127,20 +134,27 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) ! Put the remaining (total) thickness fluxes into uhr and vhr. do j=js,je ; do I=is-1,ie ; uhr(I,j,k) = uhtr(I,j,k) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; vhr(i,J,k) = vhtr(i,J,k) ; enddo ; enddo -! This loop reconstructs the thickness field the last time that the -! tracers were updated, probably just after the diabatic forcing. A useful -! diagnostic could be to compare this reconstruction with that older value. - do i=is,ie ; do j=js,je - hprev(i,j,k) = max(0.0, G%areaT(i,j)*h_end(i,j,k) + & - ((uhr(I,j,k) - uhr(I-1,j,k)) + (vhr(i,J,k) - vhr(i,J-1,k)))) -! In the case that the layer is now dramatically thinner than it was previously, -! add a bit of mass to avoid truncation errors. This will lead to -! non-conservation of tracers - hprev(i,j,k) = hprev(i,j,k) + & - max(0.0, 1.0e-13*hprev(i,j,k) - G%areaT(i,j)*h_end(i,j,k)) - enddo ; enddo + if (.not. present(h_prev_opt)) then + ! This loop reconstructs the thickness field the last time that the + ! tracers were updated, probably just after the diabatic forcing. A useful + ! diagnostic could be to compare this reconstruction with that older value. + do i=is,ie ; do j=js,je + hprev(i,j,k) = max(0.0, G%areaT(i,j)*h_end(i,j,k) + & + ((uhr(I,j,k) - uhr(I-1,j,k)) + (vhr(i,J,k) - vhr(i,J-1,k)))) + ! In the case that the layer is now dramatically thinner than it was previously, + ! add a bit of mass to avoid truncation errors. This will lead to + ! non-conservation of tracers + hprev(i,j,k) = hprev(i,j,k) + & + max(0.0, 1.0e-13*hprev(i,j,k) - G%areaT(i,j)*h_end(i,j,k)) + enddo ; enddo + else + do i=is,ie ; do j=js,je + hprev(i,j,k) = h_prev_opt(i,j,k); + enddo ; enddo + endif enddo + !$OMP do do j=jsd,jed ; do I=isd,ied-1 uh_neglect(I,j) = GV%H_subroundoff*MIN(G%areaT(i,j),G%areaT(i+1,j)) @@ -265,7 +279,9 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) endif ; enddo ! End of k-loop ! If the advection just isn't finishing after max_iter, move on. - if (itt >= max_iter) exit + if (itt >= max_iter) then + exit + endif ! Exit if there are no layers that need more iterations. if (isv > is-stencil) then @@ -274,11 +290,18 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg) call sum_across_PEs(domore_k(:), nz) call cpu_clock_end(id_clock_sync) do k=1,nz ; do_any = do_any + domore_k(k) ; enddo - if (do_any == 0) exit + if (do_any == 0) then + exit + endif + endif enddo ! Iterations loop + if(present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:) + if(present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:) + if(present(h_out)) h_out(:,:,:) = hprev(:,:,:) + call cpu_clock_end(id_clock_advect) end subroutine advect_tracer @@ -363,7 +386,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & endif ! usePLMslope ! Calculate the i-direction fluxes of each tracer, using as much - ! the minimum of the remaining mass flux (uhr) and the half the mass + ! the minimum of the remaining mass flux (uhr) and the half the mass ! in the cell plus whatever part of its half of the mass flux that ! the flux through the other side does not require. do I=is-1,ie @@ -396,6 +419,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & CFL(I) = uhh(I)/(hprev(i,j,k)+h_neglect) ! CFL is positive endif enddo + + if (usePPM) then do m=1,ntr ; do I=is-1,ie if (uhh(I) >= 0.0) then @@ -628,6 +653,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! the flux through the other side does not require. do J=js-1,je ; if (domore_v(J,k)) then domore_v(J,k) = .false. + do i=is,ie if (vhr(i,J,k) == 0.0) then vhh(i,J) = 0.0 @@ -658,6 +684,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & CFL(i) = vhh(i,J) / (hprev(i,j,k)+h_neglect) ! CFL is positive endif enddo + if (usePPM) then do m=1,ntr ; do i=is,ie if (vhh(i,J) >= 0.0) then diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index ee2b6daafa..de6a3595d5 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -248,8 +248,8 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, & is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - ! Only apply forcing if fluxes%sw is associated. - if (.not.ASSOCIATED(fluxes%sw)) return +! ! Only apply forcing if fluxes%sw is associated. +! if (.not.ASSOCIATED(fluxes%sw)) return in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0 if(present(in_flux_optional)) then diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index a07d05e61e..30a1ddebf7 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -43,7 +43,6 @@ module MOM_tracer_flow_control use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type - #include ! Add references to other user-provide tracer modules here. diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 91c4fdf1f7..643404e670 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -9,7 +9,8 @@ module MOM_tracer_hor_diff use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type use MOM_domains, only : sum_across_PEs, max_across_PEs use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type -use MOM_checksums, only : hchksum +use MOM_domains, only : pass_vector +use MOM_checksums, only : hchksum, uchksum, vchksum use MOM_EOS, only : calculate_density use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -61,6 +62,8 @@ module MOM_tracer_hor_diff integer :: id_KhTr_v = -1 integer :: id_KhTr_h = -1 integer :: id_CFL = -1 + integer :: id_khdt_x = -1 + integer :: id_khdt_y = -1 type(group_pass_type) :: pass_t !For group halo pass, used in both !tracer_hordiff and tracer_epipycnal_ML_diff @@ -81,7 +84,7 @@ module MOM_tracer_hor_diff !! using the diffusivity in CS%KhTr, or using space-dependent diffusivity. !! Multiple iterations are used (if necessary) so that there is no limit !! on the acceptable time increment. -subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) +subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y) type(ocean_grid_type), intent(inout) :: G !< Grid type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2) real, intent(in) :: dt !< time step (seconds) @@ -96,6 +99,12 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) !! NULL ptrs, and these may (probably will) point to !! some of the same arrays as Tr does. tv is required !! for epipycnal mixing between mixed layer and the interior. + ! Optional inputs for offline tracer transport + logical, optional :: do_online_flag + real, dimension(SZIB_(G),SZJ_(G)), optional, intent(in) :: read_khdt_x + real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: read_khdt_y + + real, dimension(SZI_(G),SZJ_(G)) :: & Ihdxdy, & ! The inverse of the volume or mass of fluid in a layer in a ! grid cell, in m-3 or kg-1. @@ -103,6 +112,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) CFL, & ! A diffusive CFL number for each cell, nondim. dTr ! The change in a tracer's concentration, in units of ! concentration. + real, dimension(SZIB_(G),SZJ_(G)) :: & khdt_x, & ! The value of Khtr*dt times the open face width divided by ! the distance between adjacent tracer points, in m2. @@ -117,7 +127,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) Kh_v ! Tracer mixing coefficient at u-points, in m2 s-1. real :: max_CFL ! The global maximum of the diffusive CFL number. - logical :: use_VarMix, Resoln_scaled + logical :: use_VarMix, Resoln_scaled, do_online integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts real :: I_numitts ! The inverse of the number of iterations, num_itts. real :: scale ! The fraction of khdt_x or khdt_y that is applied in this @@ -131,6 +141,10 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) real :: normalize ! normalization used for diagnostic Kh_h; diffusivity averaged to h-points. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + do_online = .true. + if (present(do_online_flag)) do_online = do_online_flag + if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "register_tracer must be called before tracer_hordiff.") if (LOC(Reg)==0) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & @@ -170,95 +184,105 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) call cpu_clock_end(id_clock_pass) if (CS%show_call_tree) call callTree_waypoint("Calculating diffusivity (tracer_hordiff)") - if (use_VarMix) then -!$OMP parallel default(none) shared(is,ie,js,je,CS,VarMix,MEKE,Resoln_scaled, & -!$OMP Kh_u,Kh_v,khdt_x,dt,G,khdt_y) & -!$OMP private(Kh_loc,Rd_dx) -!$OMP do - do j=js,je ; do I=is-1,ie - Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) - if (associated(MEKE%Kh)) & - Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) - if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) - if (Resoln_scaled) & - Kh_Loc = Kh_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) - Kh_u(I,j) = max(Kh_loc, CS%KhTr_min) - if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity - Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points - Kh_loc=Kh_u(I,j)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) - if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max - Kh_u(I,j) = max(Kh_loc, CS%KhTr_min) ! Re-apply min - endif - enddo ; enddo -!$OMP do - do J=js-1,je ; do i=is,ie - Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) - if (associated(MEKE%Kh)) & - Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) - if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) - if (Resoln_scaled) & - Kh_Loc = Kh_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) - Kh_v(i,J) = max(Kh_loc, CS%KhTr_min) - if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity - Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points - Kh_loc=Kh_v(I,j)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) - if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max - Kh_v(i,J) = max(Kh_loc, CS%KhTr_min) ! Re-apply min - endif - enddo ; enddo -!$OMP do - do j=js,je ; do I=is-1,ie - khdt_x(I,j) = dt*(Kh_u(I,j)*(G%dy_Cu(I,j)*G%IdxCu(I,j))) - enddo ; enddo -!$OMP do - do J=js-1,je ; do i=is,ie - khdt_y(i,J) = dt*(Kh_v(i,J)*(G%dx_Cv(i,J)*G%IdyCv(i,J))) - enddo ; enddo -!$OMP end parallel - elseif (Resoln_scaled) then -!$OMP parallel default(none) shared(is,ie,js,je,VarMix,Kh_u,Kh_v,khdt_x,khdt_y,CS,dt,G) & -!$OMP private(Res_fn) -!$OMP do - do j=js,je ; do I=is-1,ie - Res_fn = 0.5 * (VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) - Kh_u(I,j) = max(CS%KhTr * Res_fn, CS%KhTr_min) - khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) * Res_fn - enddo ; enddo -!$OMP do - do J=js-1,je ; do i=is,ie - Res_fn = 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) - Kh_v(i,J) = max(CS%KhTr * Res_fn, CS%KhTr_min) - khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) * Res_fn - enddo ; enddo -!$OMP end parallel - else -!$OMP parallel default(none) shared(is,ie,js,je,Kh_u,Kh_v,khdt_x,khdt_y,CS,G,dt) - if (CS%id_KhTr_u > 0) then -!$OMP do - do j=js,je ; do I=is-1,ie - Kh_u(I,j) = CS%KhTr - khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) - enddo ; enddo - else -!$OMP do - do j=js,je ; do I=is-1,ie - khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) - enddo ; enddo - endif - if (CS%id_KhTr_v > 0) then -!$OMP do - do J=js-1,je ; do i=is,ie - Kh_v(i,J) = CS%KhTr - khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) - enddo ; enddo - else -!$OMP do - do J=js-1,je ; do i=is,ie - khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) - enddo ; enddo - endif -!$OMP end parallel - endif + + if (do_online) then + if (use_VarMix) then + !$OMP parallel default(none) shared(is,ie,js,je,CS,VarMix,MEKE,Resoln_scaled, & + !$OMP Kh_u,Kh_v,khdt_x,dt,G,khdt_y) & + !$OMP private(Kh_loc,Rd_dx) + !$OMP do + do j=js,je ; do I=is-1,ie + Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) + if (associated(MEKE%Kh)) & + Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) + if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) + if (Resoln_scaled) & + Kh_Loc = Kh_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) + Kh_u(I,j) = max(Kh_loc, CS%KhTr_min) + if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity + Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points + Kh_loc=Kh_u(I,j)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) + if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max + Kh_u(I,j) = max(Kh_loc, CS%KhTr_min) ! Re-apply min + endif + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + Kh_loc = CS%KhTr + CS%KhTr_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) + if (associated(MEKE%Kh)) & + Kh_Loc = Kh_Loc + MEKE%KhTr_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) + if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) + if (Resoln_scaled) & + Kh_Loc = Kh_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) + Kh_v(i,J) = max(Kh_loc, CS%KhTr_min) + if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity + Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points + Kh_loc=Kh_v(I,j)*max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeff*Rd_dx ) + if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max + Kh_v(i,J) = max(Kh_loc, CS%KhTr_min) ! Re-apply min + endif + enddo ; enddo + + !$OMP do + do j=js,je ; do I=is-1,ie + khdt_x(I,j) = dt*(Kh_u(I,j)*(G%dy_Cu(I,j)*G%IdxCu(I,j))) + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + khdt_y(i,J) = dt*(Kh_v(i,J)*(G%dx_Cv(i,J)*G%IdyCv(i,J))) + enddo ; enddo + !$OMP end parallel + elseif (Resoln_scaled) then + !$OMP parallel default(none) shared(is,ie,js,je,VarMix,Kh_u,Kh_v,khdt_x,khdt_y,CS,dt,G) & + !$OMP private(Res_fn) + !$OMP do + do j=js,je ; do I=is-1,ie + Res_fn = 0.5 * (VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) + Kh_u(I,j) = max(CS%KhTr * Res_fn, CS%KhTr_min) + khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) * Res_fn + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + Res_fn = 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) + Kh_v(i,J) = max(CS%KhTr * Res_fn, CS%KhTr_min) + khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) * Res_fn + enddo ; enddo + !$OMP end parallel + else + !$OMP parallel default(none) shared(is,ie,js,je,Kh_u,Kh_v,khdt_x,khdt_y,CS,G,dt) + if (CS%id_KhTr_u > 0) then + !$OMP do + do j=js,je ; do I=is-1,ie + Kh_u(I,j) = CS%KhTr + khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) + enddo ; enddo + else + !$OMP do + do j=js,je ; do I=is-1,ie + khdt_x(I,j) = dt*(CS%KhTr*(G%dy_Cu(I,j)*G%IdxCu(I,j))) + enddo ; enddo + endif + if (CS%id_KhTr_v > 0) then + !$OMP do + do J=js-1,je ; do i=is,ie + Kh_v(i,J) = CS%KhTr + khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) + enddo ; enddo + else + !$OMP do + do J=js-1,je ; do i=is,ie + khdt_y(i,J) = dt*(CS%KhTr*(G%dx_Cv(i,J)*G%IdyCv(i,J))) + enddo ; enddo + endif + !$OMP end parallel + endif ! VarMix + else ! .not. do_online + khdt_x = read_khdt_x + khdt_y = read_khdt_y + call pass_vector(khdt_x,khdt_y,G%Domain) + endif ! do_online + + if (CS%check_diffusive_CFL) then if (CS%show_call_tree) call callTree_waypoint("Checking diffusive CFL (tracer_hordiff)") @@ -448,6 +472,16 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv) endif + if( CS%debug ) then + call uchksum(khdt_x,"After tracer diffusion khdt_x", G%HI, haloshift=2) + call vchksum(khdt_y,"After tracer diffusion khdt_y", G%HI, haloshift=2) + call uchksum(Coef_x,"After tracer diffusion Coef_x", G%HI, haloshift=2) + call vchksum(Coef_y,"After tracer diffusion Coef_y", G%HI, haloshift=2) + endif + + if (CS%id_khdt_x > 0) call post_data(CS%id_khdt_x, khdt_x, CS%diag) + if (CS%id_khdt_y > 0) call post_data(CS%id_khdt_y, khdt_y, CS%diag) + if (CS%show_call_tree) call callTree_leave("tracer_hordiff()") end subroutine tracer_hordiff @@ -1395,11 +1429,17 @@ subroutine tracer_hor_diff_init(Time, G, param_file, diag, CS, CSnd) cmor_field_name='diftrelo', cmor_units='m2 sec-1', & cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', & cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity') + + CS%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, Time, & + 'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'meter2') + CS%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, Time, & + 'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'meter2') if (CS%check_diffusive_CFL) then CS%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, Time,& 'Grid CFL number for lateral/neutral tracer diffusion', 'dimensionless') endif + end subroutine tracer_hor_diff_init subroutine tracer_hor_diff_end(CS) diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index 46c7264757..95d80933d9 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/src/tracer/advection_test_tracer.F90 @@ -84,13 +84,14 @@ module advection_test_tracer public advection_test_tracer_column_physics, advection_test_stock ! ntr is the number of tracers in this module. -integer, parameter :: ntr = 11 +integer, parameter :: NTR = 11 type p3d real, dimension(:,:,:), pointer :: p => NULL() end type p3d type, public :: advection_test_tracer_CS ; private + integer :: ntr = NTR ! Number of tracers in this module logical :: coupled_tracers = .false. ! These tracers are not offered to the ! coupler. character(len=200) :: tracer_IC_file ! The full path to the IC file, or " " @@ -116,7 +117,6 @@ module advection_test_tracer integer, dimension(NTR) :: ind_tr ! Indices returned by aof_set_coupler_flux ! if it is used and the surface tracer concentrations are to be ! provided to the coupler. - integer :: ntr = NTR type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the ! timing of diagnostic output. @@ -284,7 +284,7 @@ subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS h_neglect = GV%H_subroundoff CS%diag => diag - + CS%ntr = NTR if (.not.restart) then do m=1,NTR do k=1,nz ; do j=js,je ; do i=is,ie