From 0add693a8970bf206fa6879915b56c3170ae33da Mon Sep 17 00:00:00 2001 From: Philip Pegion Date: Tue, 5 May 2020 09:41:38 -0500 Subject: [PATCH] additions for stochastic physics and ePBL perts --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 2 + src/core/MOM.F90 | 39 +++++++++--- src/core/MOM_forcing_type.F90 | 19 ++++-- src/diagnostics/MOM_diagnostics.F90 | 22 +++++-- src/framework/MOM_domains.F90 | 59 ++++--------------- .../vertical/MOM_energetic_PBL.F90 | 2 +- 6 files changed, 79 insertions(+), 64 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index c704214930..10add0f8d0 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -315,6 +315,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%p_surf ,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) + print*,'allocate fluxes%t_rp' + call safe_alloc_ptr(fluxes%t_rp,isd,ied,jsd,jed) if (CS%use_limited_P_SSH) then fluxes%p_surf_SSH => fluxes%p_surf else diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c36c0545e1..ae535dcbeb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -30,6 +30,7 @@ module MOM use MOM_domains, only : To_All, Omit_corners, 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, Omit_Corners +use MOM_domains, only : root_PE,PE_here,Get_PElist,num_PEs use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -154,8 +155,7 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline -use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf -use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end +use stochastic_physics, only : init_stochastic_physics_ocn,run_stochastic_physics_ocn implicit none ; private @@ -248,6 +248,8 @@ module MOM logical :: offline_tracer_mode = .false. !< If true, step_offline() is called instead of step_MOM(). !! This is intended for running MOM6 in offline tracer mode + logical :: do_stochy = .false. + !< If true, call stochastic physics pattern generator type(time_type), pointer :: Time !< pointer to the ocean clock real :: dt !< (baroclinic) dynamics time step [T ~> s] @@ -826,6 +828,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS enddo ; enddo endif + print*,'calling run_stochastic_physics_ocn',CS%do_stochy + if (CS%do_stochy) call run_stochastic_physics_ocn(forces%t_rp) call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, & dt_therm_here, bbl_time_int, CS, & @@ -978,7 +982,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%time_in_thermo_cycle > 0.0) then call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag) call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, & - sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc) + sfc_state, CS%tv, ssh, fluxes%t_rp, CS%ave_ssh_ibc) endif call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) @@ -1811,10 +1815,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: nkml, nkbl, verbosity, write_geom integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. - real :: conv2watt ! A conversion factor from temperature fluxes to heat - ! fluxes [J m-2 H-1 degC-1 ~> J m-3 degC-1 or J kg-1 degC-1] - real :: conv2salt ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1] - real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors + integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean + integer :: num_procs +! model + integer :: me ! my pe + integer :: master ! root pe + real :: conv2watt, conv2salt character(len=48) :: flux_units, S_flux_units type(vardesc) :: vd_T, vd_S ! Structures describing temperature and salinity variables. @@ -2490,6 +2496,25 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call callTree_waypoint("returned from ALE_init() (initialize_MOM)") endif + ! Shift from using the temporary dynamic grid type to using the final + ! (potentially static) ocean-specific grid type. + ! The next line would be needed if G%Domain had not already been init'd above: + ! call clone_MOM_domain(dG%Domain, G%Domain) + call MOM_grid_init(G, param_file, US, HI, bathymetry_at_vel=bathy_at_vel) + call copy_dyngrid_to_MOM_grid(dG, G, US) + call destroy_dyn_horgrid(dG) + + num_procs=num_PEs() + allocate(pelist(num_procs)) + call Get_PElist(pelist) + me=PE_here() + master=root_PE() + + !call init_stochastic_physics_ocn(CS%dt_therm,G,me,master,pelist,CS%do_stochy) + print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) + call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,CS%do_stochy) + print*,'back from init_stochastic_physics_ocn',CS%do_stochy + ! Set a few remaining fields that are specific to the ocean grid type. if (CS%rotate_index) then call set_first_direction(G, modulo(first_direction + turns, 2)) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 3248c09fa4..e11b6a39ce 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -145,6 +145,8 @@ module MOM_forcing_type !< Pressure at the top ocean interface [R L2 T-2 ~> Pa] that is used in corrections to the sea surface !! height field that is passed back to the calling routines. !! p_surf_SSH may point to p_surf or to p_surf_full. + real, pointer, dimension(:,:) :: t_rp => NULL() + !< random pattern at t-points logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere !! and various types of ice needs to be accumulated, and the !! surface pressure explicitly reset to zero at the driver level @@ -248,10 +250,10 @@ module MOM_forcing_type !! nondimensional from 0 to 1 [nondim]. This is only associated if ice shelves are enabled, !! and is exactly 0 away from shelves or on land. real, pointer, dimension(:,:) :: & - rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at - !! u-points [L4 Z-1 T-1 ~> m3 s-1] - rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at - !! v-points [L4 Z-1 T-1 ~> m3 s-1] + rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at u-points [m3 s-1] + rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at v-points [m3 s-1] + real, pointer, dimension(:,:) :: t_rp => NULL() + !< random pattern at t-points real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes !! have been averaged [s]. logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided. @@ -2126,6 +2128,12 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) do_pres = .true. ; if (present(skip_pres)) do_pres = .not.skip_pres + if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then + do j=js,je ; do i=is,ie + fluxes%t_rp(i,j) = forces%t_rp(i,j) + enddo ; enddo + endif + if (associated(forces%ustar) .and. associated(fluxes%ustar)) then do j=js,je ; do i=is,ie fluxes%ustar(i,j) = forces%ustar(i,j) @@ -3088,6 +3096,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & call myAlloc(forces%p_surf,isd,ied,jsd,jed, press) call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press) call myAlloc(forces%net_mass_src,isd,ied,jsd,jed, press) + call myAlloc(forces%t_rp,isd,ied,jsd,jed, press) call myAlloc(forces%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(forces%rigidity_ice_v,isd,ied,JsdB,JedB, shelf) @@ -3252,6 +3261,7 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%netMassIn)) deallocate(fluxes%netMassIn) if (associated(fluxes%salt_flux)) deallocate(fluxes%salt_flux) if (associated(fluxes%p_surf_full)) deallocate(fluxes%p_surf_full) + if (associated(fluxes%t_rp)) deallocate(fluxes%t_rp) if (associated(fluxes%p_surf)) deallocate(fluxes%p_surf) if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) @@ -3280,6 +3290,7 @@ subroutine deallocate_mech_forcing(forces) if (associated(forces%ustar)) deallocate(forces%ustar) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) + if (associated(forces%t_rp)) deallocate(forces%t_rp) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) if (associated(forces%rigidity_ice_u)) deallocate(forces%rigidity_ice_u) if (associated(forces%rigidity_ice_v)) deallocate(forces%rigidity_ice_v) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 8d667503d7..fca2ed869f 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -139,7 +139,9 @@ module MOM_diagnostics integer :: id_salt_deficit = -1 integer :: id_Heat_PmE = -1 integer :: id_intern_heat = -1 - !>@} +! stochastic pattern + integer :: id_t_rp = -1 + !!@} end type surface_diag_IDs @@ -1277,7 +1279,7 @@ end subroutine post_surface_dyn_diags !> This routine posts diagnostics of various ocean surface and integrated !! quantities at the time the ocean state is reported back to the caller subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, & - ssh, ssh_ibc) + ssh, t_rp, ssh_ibc) type(surface_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs. type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1286,8 +1288,10 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv real, intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s]. type(surface), intent(in) :: sfc_state !< structure describing the ocean surface state type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< Time mean surface height without corrections - !! for ice displacement [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: t_rp!< random pattern for stochastic proceeses real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections !! for ice displacement and the inverse barometer [Z ~> m] @@ -1414,6 +1418,11 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv call post_data(IDs%id_sss_sq, work_2d, diag, mask=G%mask2dT) endif + if (IDs%id_t_rp > 0) then + !call post_data(IDs%id_t_rp, t_rp, diag, mask=G%mask2dT) + call post_data(IDs%id_t_rp, t_rp, diag) + endif + call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag)) end subroutine post_surface_thermo_diags @@ -1901,8 +1910,9 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) 'Heat flux into ocean from mass flux into ocean', & 'W m-2', conversion=US%QRZ_T_to_W_m2) IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time,& - 'Heat flux into ocean from geothermal or other internal sources', & - 'W m-2', conversion=US%QRZ_T_to_W_m2) + 'Heat flux into ocean from geothermal or other internal sources', 'W m-2') + IDs%id_t_rp = register_diag_field('ocean_model', 'random_pattern', diag%axesT1, Time, & + 'random pattern for stochastics', 'None') end subroutine register_surface_diags diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 0cdcc455fc..acbc1ccaea 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -3,56 +3,23 @@ module MOM_domains ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_coms_infra, only : MOM_infra_init, MOM_infra_end -use MOM_coms_infra, only : PE_here, root_PE, num_PEs, broadcast -use MOM_coms_infra, only : sum_across_PEs, min_across_PEs, max_across_PEs -use MOM_domain_infra, only : MOM_domain_type, domain2D, domain1D, group_pass_type -use MOM_domain_infra, only : create_MOM_domain, clone_MOM_domain, deallocate_MOM_domain -use MOM_domain_infra, only : get_domain_extent, get_domain_components, same_domain -use MOM_domain_infra, only : compute_block_extent, get_global_shape -use MOM_domain_infra, only : pass_var, pass_vector, fill_symmetric_edges, global_field_sum -use MOM_domain_infra, only : pass_var_start, pass_var_complete -use MOM_domain_infra, only : pass_vector_start, pass_vector_complete -use MOM_domain_infra, only : create_group_pass, do_group_pass -use MOM_domain_infra, only : start_group_pass, complete_group_pass -use MOM_domain_infra, only : rescale_comp_data, global_field, redistribute_array, broadcast_domain -use MOM_domain_infra, only : MOM_thread_affinity_set, set_MOM_thread_affinity -use MOM_domain_infra, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM -use MOM_domain_infra, only : CORNER, CENTER, NORTH_FACE, EAST_FACE -use MOM_domain_infra, only : To_East, To_West, To_North, To_South, To_All, Omit_Corners -use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_io_infra, only : file_exists +use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist +use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs +use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end +use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL, is_root_pe +use MOM_file_parser, only : get_param, log_param, log_version +use MOM_file_parser, only : param_file_type use MOM_string_functions, only : slasher implicit none ; private -public :: MOM_infra_init, MOM_infra_end -! Domain types and creation and destruction routines -public :: MOM_domain_type, domain2D, domain1D -public :: MOM_domains_init, create_MOM_domain, clone_MOM_domain, deallocate_MOM_domain -public :: MOM_thread_affinity_set, set_MOM_thread_affinity -! Domain query routines -public :: get_domain_extent, get_domain_components, get_global_shape, same_domain -public :: PE_here, root_PE, num_PEs -! Blocks are not actively used in MOM6, so this routine could be deprecated. -public :: compute_block_extent -! Single call communication routines -public :: pass_var, pass_vector, fill_symmetric_edges, broadcast -! Non-blocking communication routines -public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete -! Multi-variable group communication routines and type -public :: create_group_pass, do_group_pass, group_pass_type, start_group_pass, complete_group_pass -! Global reduction routines -public :: sum_across_PEs, min_across_PEs, max_across_PEs -public :: global_field, redistribute_array, broadcast_domain -! Simple index-convention-invariant array manipulation routine -public :: rescale_comp_data -!> These encoding constants are used to indicate the staggering of scalars and vectors -public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR -!> These encoding constants are used to indicate the discretization position of a variable -public :: CORNER, CENTER, NORTH_FACE, EAST_FACE -!> These encoding constants indicate communication patterns. In practice they can be added. +public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 +public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain +public :: pass_var, pass_vector, PE_here, root_PE, num_PEs, Get_PElist +public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast +public :: pass_vector_start, pass_vector_complete +public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs +public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER, CENTER public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners ! These are no longer used by MOM6 because the reproducing sum works so well, but they are ! still referenced by some of the non-GFDL couplers. diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 99dd38135d..c9ae6e43ed 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -401,7 +401,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS do K=1,nz+1 ; Kd(K) = 0.0 ; enddo ! Make local copies of surface forcing and process them. - u_star = fluxes%ustar(i,j) + u_star = fluxes%ustar(i,j)*(fluxes%t_rp(i,j)) u_star_Mean = fluxes%ustar_gustless(i,j) B_flux = buoy_flux(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then