diff --git a/.DoxygenLayout.xml b/.DoxygenLayout.xml index ebe1bea767..80047a4105 100644 --- a/.DoxygenLayout.xml +++ b/.DoxygenLayout.xml @@ -1,31 +1,26 @@ - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + diff --git a/.doxygen b/.doxygen index fceb7fea9b..8b1e2aafab 100644 --- a/.doxygen +++ b/.doxygen @@ -872,7 +872,7 @@ EXCLUDE_SYMBOLS = # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = config_src src # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and @@ -892,7 +892,7 @@ EXAMPLE_RECURSIVE = NO # that contain images that are to be included in the documentation (see the # \image command). -IMAGE_PATH = +IMAGE_PATH = .images config_src src # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program @@ -2028,7 +2028,7 @@ SEARCH_INCLUDES = YES # preprocessor. # This tag requires that the tag SEARCH_INCLUDES is set to YES. -INCLUDE_PATH = +INCLUDE_PATH = src/framework # You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard # patterns (like *.h and *.hpp) to filter out the header-files in the diff --git a/.images/Arakawa_C_grid.png b/.images/Arakawa_C_grid.png new file mode 100644 index 0000000000..a765bf9338 Binary files /dev/null and b/.images/Arakawa_C_grid.png differ diff --git a/.images/Grid_metrics.png b/.images/Grid_metrics.png new file mode 100644 index 0000000000..22b1f0d55f Binary files /dev/null and b/.images/Grid_metrics.png differ diff --git a/.images/Horizontal_NE_indexing_nonsym.png b/.images/Horizontal_NE_indexing_nonsym.png new file mode 100644 index 0000000000..6316a64e7a Binary files /dev/null and b/.images/Horizontal_NE_indexing_nonsym.png differ diff --git a/.images/Horizontal_NE_indexing_sym.png b/.images/Horizontal_NE_indexing_sym.png new file mode 100644 index 0000000000..c2aaa742ac Binary files /dev/null and b/.images/Horizontal_NE_indexing_sym.png differ diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 8be4c72bfe..ecf05f9956 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -194,6 +194,9 @@ module MOM_surface_forcing real, pointer, dimension(:,:) :: fprec =>NULL() ! mass flux of frozen precip (kg/m2/s) real, pointer, dimension(:,:) :: runoff =>NULL() ! mass flux of liquid runoff (kg/m2/s) real, pointer, dimension(:,:) :: calving =>NULL() ! mass flux of frozen runoff (kg/m2/s) + real, pointer, dimension(:,:) :: ustar_berg =>NULL() ! frictional velocity beneath icebergs (m/s) + real, pointer, dimension(:,:) :: area_berg =>NULL() ! area covered by icebergs(m2/m2) + real, pointer, dimension(:,:) :: mass_berg =>NULL() ! mass of icebergs(kg/m2) real, pointer, dimension(:,:) :: runoff_hflx =>NULL() ! heat content of liquid runoff (W/m2) real, pointer, dimension(:,:) :: calving_hflx =>NULL() ! heat content of frozen runoff (W/m2) real, pointer, dimension(:,:) :: p =>NULL() ! pressure of overlying ice and atmosphere @@ -469,6 +472,20 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, if (ASSOCIATED(IOB%calving)) & fluxes%frunoff(i,j) = IOB%calving(i-i0,j-j0) * G%mask2dT(i,j) + if (((ASSOCIATED(IOB%ustar_berg) .and. (.not. ASSOCIATED(fluxes%ustar_berg))) & + .or. (ASSOCIATED(IOB%area_berg) .and. (.not. ASSOCIATED(fluxes%area_berg)))) & + .or. (ASSOCIATED(IOB%mass_berg) .and. (.not. ASSOCIATED(fluxes%mass_berg)))) & + call allocate_forcing_type(G, fluxes, iceberg=.true.) + + if (ASSOCIATED(IOB%ustar_berg)) & + fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%area_berg)) & + fluxes%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%mass_berg)) & + fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) + if (ASSOCIATED(IOB%runoff_hflx)) & fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) @@ -1197,6 +1214,12 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) write(outunit,100) 'iobt%runoff ', mpp_chksum( iobt%runoff ) write(outunit,100) 'iobt%calving ', mpp_chksum( iobt%calving ) write(outunit,100) 'iobt%p ', mpp_chksum( iobt%p ) + if (ASSOCIATED(iobt%ustar_berg)) & + write(outunit,100) 'iobt%ustar_berg ', mpp_chksum( iobt%ustar_berg ) + if (ASSOCIATED(iobt%area_berg)) & + write(outunit,100) 'iobt%area_berg ', mpp_chksum( iobt%area_berg ) + if (ASSOCIATED(iobt%mass_berg)) & + write(outunit,100) 'iobt%mass_berg ', mpp_chksum( iobt%mass_berg ) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) do n = 1, iobt%fluxes%num_bcs !{ diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 3543b5d6c0..6cd379ea9c 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -67,8 +67,10 @@ module ocean_model_mod use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux +use MOM_forcing_type, only : allocate_forcing_type use fms_mod, only : stdout use mpp_mod, only : mpp_chksum +use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE #include @@ -153,6 +155,9 @@ module ocean_model_mod integer :: nstep = 0 ! The number of calls to update_ocean. logical :: use_ice_shelf ! If true, the ice shelf model is enabled. + logical :: icebergs_apply_rigid_boundary ! If true, the icebergs can change ocean bd condition. + real :: kv_iceberg ! The viscosity of the icebergs in m2/s (for ice rigidity) + real :: density_iceberg ! A typical density of icebergs in kg/m3 (for ice rigidity) type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() logical :: restore_salinity ! If true, the coupled MOM driver adds a term to ! restore salinity to a specified value. @@ -160,7 +165,7 @@ module ocean_model_mod ! restore sst to a specified value. real :: press_to_z ! A conversion factor between pressure and ocean ! depth in m, usually 1/(rho_0*g), in m Pa-1. - real :: C_p ! The heat capacity of seawater, in J K-1 kg-1. + real :: C_p ! The heat capacity of seawater, in J K-1 kg-1. type(directories) :: dirs ! A structure containing several relevant directory paths. type(forcing) :: fluxes ! A structure containing pointers to @@ -280,6 +285,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in) call get_param(param_file, mod, "ICE_SHELF", OS%use_ice_shelf, & "If true, enables the ice shelf model.", default=.false.) + call get_param(param_file, mod, "ICEBERGS_APPLY_RIGID_BOUNDARY", OS%icebergs_apply_rigid_boundary, & + "If true, allows icebergs to change boundary condition felt by ocean", default=.false.) + + if (OS%icebergs_apply_rigid_boundary) then + call get_param(param_file, mod, "KV_ICEBERG", OS%kv_iceberg, & + "The viscosity of the icebergs", units="m2 s-1",default=1.0e10) + call get_param(param_file, mod, "DENSITY_ICEBERGS", OS%density_iceberg, & + "A typical density of icebergs.", units="kg m-3", default=917.0) + endif + OS%press_to_z = 1.0/(Rho0*G_Earth) call surface_forcing_init(Time_in, OS%grid, param_file, OS%MOM_CSp%diag, & @@ -289,6 +304,11 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in) call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & OS%MOM_CSp%diag, OS%fluxes) endif + if (OS%icebergs_apply_rigid_boundary) then + !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.) + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + if (.not. OS%use_ice_shelf) call allocate_forcing_type(OS%grid, OS%fluxes, ustar=.true., shelf=.true.) + endif call MOM_sum_output_init(OS%grid, param_file, OS%dirs%output_directory, & OS%MOM_CSp%ntrunc, Time_init, OS%sum_output_CSp) @@ -394,12 +414,14 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes #endif - ! Add ice shelf fluxes - + ! Add ice shelf fluxes if (OS%use_ice_shelf) then call shelf_calc_flux(OS%State, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp) endif - + if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + call add_berg_flux_to_shelf(OS%grid, OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg) + endif ! Indicate that there are new unused fluxes. OS%fluxes%fluxes_used = .false. OS%fluxes%dt_buoy_accum = time_step @@ -407,10 +429,13 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%flux_tmp%C_p = OS%fluxes%C_p call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) - if (OS%use_ice_shelf) then call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) endif + if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg) + endif call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight) #ifdef _USE_GENERIC_TRACER @@ -477,6 +502,66 @@ end subroutine update_ocean_model ! the any restart file name as a prefix. ! ! + +subroutine add_berg_flux_to_shelf(G, fluxes, use_ice_shelf, density_ice, kv_ice) + type(ocean_grid_type), intent(inout) :: G + type(forcing), intent(inout) :: fluxes + logical, intent(in) :: use_ice_shelf + real, intent(in) :: kv_ice ! The viscosity of ice, in m2 s-1. + real, intent(in) :: density_ice ! A typical density of ice, in kg m-3. +! Arguments: +! (in) fluxes - A structure of surface fluxes that may be used. +! (in) G - The ocean's grid structure. + integer :: i, j, is, ie, js, je, isd, ied, jsd, jed + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed + !This routine adds iceberg data to the ice shelf data (if ice shelf is used) + !which can then be used to change the top of ocean boundary condition used in + !the ocean model. This routine is taken from the add_shelf_flux subroutine + !within the ice shelf model. + + if (.not. (((associated(fluxes%frac_shelf_h) .and. associated(fluxes%frac_shelf_u)) & + .and.(associated(fluxes%frac_shelf_v) .and. associated(fluxes%ustar_shelf)))& + .and.(associated(fluxes%rigidity_ice_u) .and. associated(fluxes%rigidity_ice_v)))) return + + if (.not. ((associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg)) .and. associated(fluxes%mass_berg) ) ) return + + if (.not. use_ice_shelf) then + fluxes%frac_shelf_h(:,:)=0. + fluxes%frac_shelf_u(:,:)=0. + fluxes%frac_shelf_v(:,:)=0. + fluxes%ustar_shelf(:,:)=0. + fluxes%rigidity_ice_u(:,:)=0. + fluxes%rigidity_ice_v(:,:)=0. + endif + + do j=jsd,jed ; do i=isd,ied + if (G%areaT(i,j) > 0.0) & + fluxes%frac_shelf_h(i,j) = fluxes%frac_shelf_h(i,j) + fluxes%area_berg(i,j) + fluxes%ustar_shelf(i,j) = fluxes%ustar_shelf(i,j) + fluxes%ustar_berg(i,j) + enddo ; enddo + !do I=isd,ied-1 ; do j=isd,jed + do j=jsd,jed ; do i=isd,ied-1 ! ### changed stride order; i->ied-1? + fluxes%frac_shelf_u(I,j) = 0.0 + if ((G%areaT(i,j) + G%areaT(i+1,j) > 0.0)) & ! .and. (G%dxdy_u(I,j) > 0.0)) & + fluxes%frac_shelf_u(I,j) = fluxes%frac_shelf_u(I,j) + (((fluxes%area_berg(i,j)*G%areaT(i,j)) + (fluxes%area_berg(i+1,j)*G%areaT(i+1,j))) / & + (G%areaT(i,j) + G%areaT(i+1,j))) + fluxes%rigidity_ice_u(I,j) = fluxes%rigidity_ice_u(I,j) +((kv_ice / density_ice) * & + min(fluxes%mass_berg(i,j), fluxes%mass_berg(i+1,j))) + enddo ; enddo + do j=jsd,jed-1 ; do i=isd,ied ! ### change stride order; j->jed-1? + !do i=isd,ied ; do J=isd,jed-1 + fluxes%frac_shelf_v(i,J) = 0.0 + if ((G%areaT(i,j) + G%areaT(i,j+1) > 0.0)) & ! .and. (G%dxdy_v(i,J) > 0.0)) & + fluxes%frac_shelf_v(i,J) = fluxes%frac_shelf_v(i,J) + (((fluxes%area_berg(i,j)*G%areaT(i,j)) + (fluxes%area_berg(i,j+1)*G%areaT(i,j+1))) / & + (G%areaT(i,j) + G%areaT(i,j+1) )) + fluxes%rigidity_ice_v(i,J) = fluxes%rigidity_ice_v(i,J) +((kv_ice / density_ice) * & + max(fluxes%mass_berg(i,j), fluxes%mass_berg(i,j+1))) + enddo ; enddo + call pass_vector(fluxes%frac_shelf_u, fluxes%frac_shelf_v, G%domain, TO_ALL, CGRID_NE) + +end subroutine add_berg_flux_to_shelf + subroutine ocean_model_restart(OS, timestamp) type(ocean_state_type), pointer :: OS character(len=*), intent(in), optional :: timestamp diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 0e9fca8edf..3b39e753a6 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -2870,7 +2870,7 @@ end subroutine regridding_memory_deallocation !! Most calculations in this module start with the coordinate at the bottom !! of the column set to -depth, and use a increasing value of coordinate with !! decreasing k. This is consistent with the rest of MOM6 that uses position, -!! f$z\f$ which is a negative quantity for most of the ocean. +!! \f$z\f$ which is a negative quantity for most of the ocean. !! !! A change in grid is define through a change in position of the interfaces: !! \f[ diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3062f638dd..0298e53925 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -544,8 +544,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call create_group_pass(CS%pass_tau_ustar_psurf, fluxes%ustar(:,:), G%Domain) if (ASSOCIATED(fluxes%p_surf)) & call create_group_pass(CS%pass_tau_ustar_psurf, fluxes%p_surf(:,:), G%Domain) - if (CS%thickness_diffuse .OR. CS%mixedlayer_restrat) & + if ((CS%thickness_diffuse .and. (.not.CS%thickness_diffuse_first .or. CS%dt_trans == 0) ) .OR. CS%mixedlayer_restrat) & call create_group_pass(CS%pass_h, h, G%Domain) + if (CS%diabatic_first) then if (associated(CS%visc%Ray_u) .and. associated(CS%visc%Ray_v)) & call create_group_pass(CS%pass_ray, CS%visc%Ray_u, CS%visc%Ray_v, G%Domain, & @@ -4079,9 +4080,9 @@ end subroutine MOM_end !! * FRAZIL_HEAT_TENDENCY generally has 3d structure, since MOM6 frazil calculation checks the !! full ocean column. !! -!! * FRAZIL_HEAT_TENDENCY[k=@sum] = HFSIFRAZIL = column integrated frazil heating. +!! * FRAZIL_HEAT_TENDENCY[k=\@sum] = HFSIFRAZIL = column integrated frazil heating. !! -!! * HFDS = FRAZIL_HEAT_TENDENCY[k=@sum] + BOUNDARY_FORCING_HEAT_TENDENCY[k=@sum] +!! * HFDS = FRAZIL_HEAT_TENDENCY[k=\@sum] + BOUNDARY_FORCING_HEAT_TENDENCY[k=\@sum] !! !! Here is an example 2d heat budget (depth summed) diagnostic for MOM. !! @@ -4107,7 +4108,7 @@ end subroutine MOM_end !! * BOUNDARY_FORCING_SALT_TENDENCY generally has 3d structure, with k > 1 contributions from !! the case when layers are tiny, in which case MOM6 partitions tendencies into k > 1 layers. !! -!! * SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=@sum] +!! * SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=\@sum] !! !! Here is an example 2d salt budget (depth summed) diagnostic for MOM. !! diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index d53d1c9728..a175eeed2f 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -1,66 +1,9 @@ +!> Accelerations due to the Coriolis force and momentum advection module MOM_CoriolisAdv -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** - -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - June 2002 * -!* * -!* This file contains the subroutine that calculates the time * -!* derivatives of the velocities due to Coriolis acceleration and * -!* momentum advection. This subroutine uses either a vorticity * -!* advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or * -!* Sadourny's (JAS 1975) energy conserving scheme. Both have been * -!* modified to use general orthogonal coordinates as described in * -!* Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second * -!* order accurate, and allow for vanishingly small layer thicknesses. * -!* The Arakawa and Hsu scheme globally conserves both total energy * -!* and potential enstrophy in the limit of nondivergent flow. * -!* Sadourny's energy conserving scheme conserves energy if the flow * -!* is nondivergent or centered difference thickness fluxes are used. * -!* * -!* Two sets of boundary conditions have been coded in the * -!* definition of relative vorticity. These are written as: * -!* NOSLIP defined (in spherical coordinates): * -!* relvort = dv/dx (east & west), with v = 0. * -!* relvort = -sec(Q) * d(u cos(Q))/dy (north & south), with u = 0. * -!* NOSLIP not defined (free slip): * -!* relvort = 0 (all boundaries) * -!* with Q temporarily defined as latitude. The free slip boundary * -!* condition is much more natural on a C-grid. * -!* * -!* Macros written all in capital letters are defined in MOM_memory.h. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q, CoriolisBu * -!* j+1 > o > o > At ^: v, CAv, vh * -!* j x ^ x ^ x At >: u, CAu, uh, a, b, c, d * -!* j > o > o > At o: h, KE * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** + +! This file is part of MOM6. See LICENSE.md for the license. + +!> \author Robert Hallberg, April 1994 - June 2002 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 @@ -78,60 +21,54 @@ module MOM_CoriolisAdv #include +!> Control structure for mom_coriolisadv type, public :: CoriolisAdv_CS ; private - integer :: Coriolis_Scheme ! CORIOLIS_SCHEME selects the discretization for - ! the Coriolis terms. Valid values are: - ! SADOURNY75_ENERGY - Sadourny, 1975 - ! ARAKAWA_HSU90 - Arakawa & Hsu, 1990 - ! Energy & non-div. Enstrophy - ! ROBUST_ENSTRO - Pseudo-enstrophy scheme - ! SADOURNY75_ENSTRO - Sadourny, JAS 1975 - ! Enstrophy - ! ARAKAWA_LAMB81 - Arakawa & Lamb, MWR 1981 - ! Energy & Enstrophy - ! ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb - ! with Arakawa & Hsu and - ! Sadourny energy. - ! The default, SADOURNY75_ENERGY, is the safest - ! choice when the deformation radius is poorly - ! resolved. - integer :: KE_Scheme ! KE_SCHEME selects the discretization for - ! the kinetic energy. Valid values are: - ! KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV - integer :: PV_Adv_Scheme ! PV_ADV_SCHEME selects the discretization for - ! PV advection. Valid values are: - ! PV_ADV_CENTERED - centered (aka Sadourny, 75) - ! PV_ADV_UPWIND1 - upwind, first order - real :: F_eff_max_blend ! The factor by which the maximum effective Coriolis - ! acceleration from any point can be increased when - ! blending different discretizations with the - ! ARAKAWA_LAMB_BLEND Coriolis scheme. This must be - ! greater than 2.0, and is 4.0 by default. - real :: wt_lin_blend ! A weighting value beyond which the blending between - ! Sadourny and Arakawa & Hsu goes linearly to 0. - ! This must be between 1 and 1e-15, often 1/8. - logical :: no_slip ! If true, no slip boundary conditions are used. - ! Otherwise free slip boundary conditions are assumed. - ! The implementation of the free slip boundary - ! conditions on a C-grid is much cleaner than the - ! no slip boundary conditions. The use of free slip - ! b.c.s is strongly encouraged. The no slip b.c.s - ! are not implemented with the biharmonic viscosity. - logical :: bound_Coriolis ! If true, the Coriolis terms at u points are - ! bounded by the four estimates of (f+rv)v from the - ! four neighboring v points, and similarly at v - ! points. This option would have no effect on the - ! SADOURNY75_ENERGY scheme if it were possible to - ! use centered difference thickness fluxes. - logical :: Coriolis_En_Dis ! If CORIOLIS_EN_DIS is defined, two estimates of - ! the thickness fluxes are used to estimate the - ! Coriolis term, and the one that dissipates energy - ! relative to the other one is used. This is only - ! available at present if Coriolis scheme is - ! SADOURNY75_ENERGY. - type(time_type), pointer :: Time ! A pointer to the ocean model's clock. - type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the - ! timing of diagnostic output. + integer :: Coriolis_Scheme !< Selects the discretization for the Coriolis terms. + !! Valid values are: + !! - SADOURNY75_ENERGY - Sadourny, 1975 + !! - ARAKAWA_HSU90 - Arakawa & Hsu, 1990, Energy & non-div. Enstrophy + !! - ROBUST_ENSTRO - Pseudo-enstrophy scheme + !! - SADOURNY75_ENSTRO - Sadourny, JAS 1975, Enstrophy + !! - ARAKAWA_LAMB81 - Arakawa & Lamb, MWR 1981, Energy & Enstrophy + !! - ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with Arakawa & Hsu and Sadourny energy. + !! The default, SADOURNY75_ENERGY, is the safest choice then the + !! deformation radius is poorly resolved. + integer :: KE_Scheme !< KE_SCHEME selects the discretization for + !! the kinetic energy. Valid values are: + !! KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV + integer :: PV_Adv_Scheme !< PV_ADV_SCHEME selects the discretization for PV advection + !! Valid values are: + !! - PV_ADV_CENTERED - centered (aka Sadourny, 75) + !! - PV_ADV_UPWIND1 - upwind, first order + real :: F_eff_max_blend !< The factor by which the maximum effective Coriolis + !! acceleration from any point can be increased when + !! blending different discretizations with the + !! ARAKAWA_LAMB_BLEND Coriolis scheme. This must be + !! greater than 2.0, and is 4.0 by default. + real :: wt_lin_blend !< A weighting value beyond which the blending between + !! Sadourny and Arakawa & Hsu goes linearly to 0. + !! This must be between 1 and 1e-15, often 1/8. + logical :: no_slip !< If true, no slip boundary conditions are used. + !! Otherwise free slip boundary conditions are assumed. + !! The implementation of the free slip boundary + !! conditions on a C-grid is much cleaner than the + !! no slip boundary conditions. The use of free slip + !! b.c.s is strongly encouraged. The no slip b.c.s + !! are not implemented with the biharmonic viscosity. + logical :: bound_Coriolis !< If true, the Coriolis terms at u points are + !! bounded by the four estimates of (f+rv)v from the + !! four neighboring v points, and similarly at v + !! points. This option would have no effect on the + !! SADOURNY75_ENERGY scheme if it were possible to + !! use centered difference thickness fluxes. + logical :: Coriolis_En_Dis !< If CORIOLIS_EN_DIS is defined, two estimates of + !! the thickness fluxes are used to estimate the + !! Coriolis term, and the one that dissipates energy + !! relative to the other one is used. This is only + !! available at present if Coriolis scheme is + !! SADOURNY75_ENERGY. + type(time_type), pointer :: Time !< A pointer to the ocean model's clock. + type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output. integer :: id_rv = -1, id_PV = -1, id_gKEu = -1, id_gKEv = -1 integer :: id_rvxu = -1, id_rvxv = -1 end type CoriolisAdv_CS @@ -164,46 +101,23 @@ module MOM_CoriolisAdv contains +!> Calculates the Coriolis and momentum advection contributions to the acceleration. subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: CAu - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: CAv - type(ocean_OBC_type), pointer :: OBC - type(accel_diag_ptrs), intent(inout) :: AD - type(CoriolisAdv_CS), pointer :: CS -! This subroutine calculates the Coriolis and momentum advection -! contributions to the acceleration. -! -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) h - Layer thickness, in m or kg m-2. -! (in) uh - Volume flux through zonal faces = u*h*dy, m3 s-1 or kg s-1. -! (in) vh - Volume flux through meridional faces = v*h*dx, in units of -! m3 s-1 or kg s-1. -! (out) CAu - Zonal acceleration due to Coriolis and momentum -! advection terms, in m s-2. -! (out) CAv - Meridional acceleration due to Coriolis and -! momentum advection terms, in m s-2. -! (in) AD - A structure pointing to the various accelerations in -! the momentum equations. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! CoriolisAdv_init. -! -! To work, the following fields must be set outside of the usual -! is to ie range before this subroutine is called: -! v[is-1,ie+1,ie+2], u[is-1,ie+1], vh[ie+1], uh[is-1], and -! h[is-1,ie+1,ie+2]. -! In the y-direction, the following fields must be set: -! v[js-1,je+1], u[js-1,je+1,je+2], vh[js-1], uh[je+1], and -! h[js-1,je+1,je+2]. + type(ocean_grid_type), intent(in) :: G !< Ocen grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity (m/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg/m2) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh !< Zonal transport u*h*dy (m3/s or kg/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh !< Meridional transport v*h*dx (m3/s or kg/s) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: CAu !< Zonal acceleration due to Coriolis + !! and momentum advection, in m/s2. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: CAv !< Meridional acceleration due to Coriolis + !! and momentum advection, in m/s2. + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics + type(CoriolisAdv_CS), pointer :: CS !< Control structure for MOM_CoriolisAdv + ! Local variables real, dimension(SZIB_(G),SZJB_(G)) :: & q, & ! Layer potential vorticity, in m-1 s-1. @@ -283,6 +197,14 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) real :: QUHeff,QVHeff ! More temporary variables, in m3 s-2 or kg s-2. integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz +! To work, the following fields must be set outside of the usual +! is to ie range before this subroutine is called: +! v[is-1,ie+1,ie+2], u[is-1,ie+1], vh[ie+1], uh[is-1], and +! h[is-1,ie+1,ie+2]. +! In the y-direction, the following fields must be set: +! v[js-1,je+1], u[js-1,je+1,je+2], vh[js-1], uh[je+1], and +! h[js-1,je+1,je+2]. + if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_CoriolisAdv: Module must be initialized before it is used.") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -311,10 +233,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) !$OMP do do k=1,nz -! Here the second order accurate layer potential vorticities, q, -! are calculated. hq is second order accurate in space. Relative -! vorticity is second order accurate everywhere with free slip b.c.s, -! but only first order accurate at boundaries with no slip b.c.s. + ! Here the second order accurate layer potential vorticities, q, + ! are calculated. hq is second order accurate in space. Relative + ! vorticity is second order accurate everywhere with free slip b.c.s, + ! but only first order accurate at boundaries with no slip b.c.s. do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 if (CS%no_slip ) then relative_vorticity = (2.0-G%mask2dBu(I,J)) * & @@ -361,9 +283,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) q2(I,J) = relative_vorticity * Ih enddo ; enddo -! a, b, c, and d are combinations of neighboring potential -! vorticities which form the Arakawa and Hsu vorticity advection -! scheme. All are defined at u grid points. + ! a, b, c, and d are combinations of neighboring potential + ! vorticities which form the Arakawa and Hsu vorticity advection + ! scheme. All are defined at u grid points. if (CS%Coriolis_Scheme == ARAKAWA_HSU90) then do j=Jsq,Jeq+1 @@ -480,12 +402,12 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) enddo ; enddo endif -! Calculate KE and the gradient of KE + ! Calculate KE and the gradient of KE call gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) -! Calculate the tendencies of zonal velocity due to the Coriolis -! force and momentum advection. On a Cartesian grid, this is -! CAu = q * vh - d(KE)/dx. + ! Calculate the tendencies of zonal velocity due to the Coriolis + ! force and momentum advection. On a Cartesian grid, this is + ! CAu = q * vh - d(KE)/dx. if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then if (CS%Coriolis_En_Dis) then ! Energy dissipating biased scheme, Hallberg 200x @@ -592,9 +514,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) enddo ; enddo -! Calculate the tendencies of meridional velocity due to the Coriolis -! force and momentum advection. On a Cartesian grid, this is -! CAv = - q * uh - d(KE)/dy. + ! Calculate the tendencies of meridional velocity due to the Coriolis + ! force and momentum advection. On a Cartesian grid, this is + ! CAv = - q * uh - d(KE)/dy. if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then if (CS%Coriolis_En_Dis) then ! Energy dissipating biased scheme, Hallberg 200x @@ -698,7 +620,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) enddo ; enddo if (ASSOCIATED(AD%rv_x_u) .or. ASSOCIATED(AD%rv_x_v)) then -! Calculate the Coriolis-like acceleration due to relative vorticity. + ! Calculate the Coriolis-like acceleration due to relative vorticity. if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then if (ASSOCIATED(AD%rv_x_u)) then do J=Jsq,Jeq ; do i=is,ie @@ -741,8 +663,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) enddo ! k-loop. !$OMP end parallel -! Here the various Coriolis-related derived quantities are offered -! for averaging. + ! Here the various Coriolis-related derived quantities are offered for averaging. if (query_averaging_enabled(CS%diag)) then if (CS%id_rv > 0) call post_data(CS%id_rv, RV, CS%diag) if (CS%id_PV > 0) call post_data(CS%id_PV, PV, CS%diag) @@ -754,38 +675,24 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) end subroutine CorAdCalc -! ========================================================================================= +!> Calculates the acceleration due to the gradient of kinetic energy. subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZIB_(G),SZJ_(G) ,SZK_(G)), intent(in) :: u - real, dimension(SZI_(G) ,SZJB_(G),SZK_(G)), intent(in) :: v - real, dimension(SZI_(G) ,SZJ_(G) ,SZK_(G)), intent(in) :: h - real, dimension(SZIB_(G),SZJ_(G) ,SZK_(G)), intent(in) :: uh - real, dimension(SZI_(G),SZJB_(G) ,SZK_(G)), intent(in) :: vh - real, dimension(SZI_(G) ,SZJ_(G) ), intent(out) :: KE - real, dimension(SZIB_(G),SZJ_(G) ), intent(out) :: KEx - real, dimension(SZI_(G) ,SZJB_(G)), intent(out) :: KEy - integer, intent(in) :: k - type(ocean_OBC_type), pointer :: OBC - type(CoriolisAdv_CS), pointer :: CS -! This subroutine calculates the acceleration due to the gradient of kinetic energy. -! -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) h - Layer thickness, in m. -! (in) uh - Volume flux through zonal faces = u*h*dy, m3 s-1. -! (in) vh - Volume flux through meridional faces = v*h*dx, in m3 s-1. -! (out) KE - Kinetic energy used by s/r, in m2 s-2. -! (out) KEx - Zonal acceleration due to Coriolis and momentum -! advection terms, in m s-2. -! (out) KEy - Meridional acceleration due to Coriolis and -! momentum advection terms, in m s-2. -! (in) k - layer number -! (in) G - The ocean's grid structure. -! (in) CS - The control structure returned by a previous call to -! CoriolisAdv_init. - + type(ocean_grid_type), intent(in) :: G !< Ocen grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity (m/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg/m2) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh !< Zonal transport u*h*dy (m3/s or kg/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh !< Meridional transport v*h*dx (m3/s or kg/s) + real, dimension(SZI_(G) ,SZJ_(G) ), intent(out) :: KE !< Kinetic energy (m2/s2) + real, dimension(SZIB_(G),SZJ_(G) ), intent(out) :: KEx !< Zonal acceleration due to kinetic + !! energy gradient (m/s2) + real, dimension(SZI_(G) ,SZJB_(G)), intent(out) :: KEy !< Meridional acceleration due to kinetic + !! energy gradient (m/s2) + integer, intent(in) :: k !< Layer number to calculate for + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(CoriolisAdv_CS), pointer :: CS !< Control structure for MOM_CoriolisAdv + ! Local variables real :: um, up, vm, vp ! Temporary variables with units of m s-1. real :: um2, up2, vm2, vp2 ! Temporary variables with units of m2 s-2. real :: um2a, up2a, vm2a, vp2a ! Temporary variables with units of m4 s-2. @@ -795,11 +702,11 @@ subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB -! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term). + ! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term). if (CS%KE_Scheme.eq.KE_ARAKAWA) then -! The following calculation of Kinetic energy includes the metric terms -! identified in Arakawa & Lamb 1982 as important for KE conservation. It -! also includes the possibility of partially-blocked tracer cell faces. + ! The following calculation of Kinetic energy includes the metric terms + ! identified in Arakawa & Lamb 1982 as important for KE conservation. It + ! also includes the possibility of partially-blocked tracer cell faces. do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE(i,j) = ( ( G%areaCu( I ,j)*(u( I ,j,k)*u( I ,j,k)) & +G%areaCu(I-1,j)*(u(I-1,j,k)*u(I-1,j,k)) ) & @@ -808,8 +715,8 @@ subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) )*0.25*G%IareaT(i,j) enddo ; enddo elseif (CS%KE_Scheme.eq.KE_SIMPLE_GUDONOV) then -! The following discretization of KE is based on the one-dimensinal Gudonov -! scheme which does not take into account any geometric factors + ! The following discretization of KE is based on the one-dimensinal Gudonov + ! scheme which does not take into account any geometric factors do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 up = 0.5*( u(I-1,j,k) + ABS( u(I-1,j,k) ) ) ; up2 = up*up um = 0.5*( u( I ,j,k) - ABS( u( I ,j,k) ) ) ; um2 = um*um @@ -818,8 +725,8 @@ subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) KE(i,j) = ( max(up2,um2) + max(vp2,vm2) ) *0.5 enddo ; enddo elseif (CS%KE_Scheme.eq.KE_GUDONOV) then -! The following discretization of KE is based on the one-dimensinal Gudonov -! scheme but has been adapted to take horizontal grid factors into account + ! The following discretization of KE is based on the one-dimensinal Gudonov + ! scheme but has been adapted to take horizontal grid factors into account do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 up = 0.5*( u(I-1,j,k) + ABS( u(I-1,j,k) ) ) ; up2a = up*up*G%areaCu(I-1,j) um = 0.5*( u( I ,j,k) - ABS( u( I ,j,k) ) ) ; um2a = um*um*G%areaCu( I ,j) @@ -829,12 +736,12 @@ subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) enddo ; enddo endif -! Term - d(KE)/dx. + ! Term - d(KE)/dx. do j=js,je ; do I=Isq,Ieq KEx(I,j) = (KE(i+1,j) - KE(i,j)) * G%IdxCu(i,j) enddo ; enddo -! Term - d(KE)/dy. + ! Term - d(KE)/dy. do J=Jsq,Jeq ; do i=is,ie KEy(i,J) = (KE(i,j+1) - KE(i,j)) * G%IdyCv(i,J) enddo ; enddo @@ -850,31 +757,22 @@ subroutine gradKE(u, v, h, uh, vh, KE, KEx, KEy, k, OBC, G, CS) end subroutine gradKE -! ========================================================================================= - +!> Initializes the control structure for coriolisadv_cs subroutine CoriolisAdv_init(Time, G, param_file, diag, AD, CS) - type(time_type), target, intent(in) :: Time - type(ocean_grid_type), intent(in) :: G - type(param_file_type), intent(in) :: param_file - type(diag_ctrl), target, intent(inout) :: diag - type(accel_diag_ptrs), target, intent(inout) :: AD - type(CoriolisAdv_CS), pointer :: CS -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (inout) AD - A structure pointing to the various accelerations in -! the momentum equations. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module - + type(time_type), target, intent(in) :: Time !< Current model time + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Runtime parameter handles + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(accel_diag_ptrs), target, intent(inout) :: AD !< Strorage for acceleration diagnostics + type(CoriolisAdv_CS), pointer :: CS !< Control structure fro MOM_CoriolisAdv + ! Local variables ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_CoriolisAdv" ! This module's name. character(len=20) :: tmpstr character(len=400) :: mesg integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -1029,9 +927,51 @@ subroutine CoriolisAdv_init(Time, G, param_file, diag, AD, CS) end subroutine CoriolisAdv_init +!> Destructor for coriolisadv_cs subroutine CoriolisAdv_end(CS) - type(CoriolisAdv_CS), pointer :: CS + type(CoriolisAdv_CS), pointer :: CS !< Control structure fro MOM_CoriolisAdv deallocate(CS) end subroutine CoriolisAdv_end +!> \namespace mom_coriolisadv +!! +!! This file contains the subroutine that calculates the time +!! derivatives of the velocities due to Coriolis acceleration and +!! momentum advection. This subroutine uses either a vorticity +!! advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or +!! Sadourny's (JAS 1975) energy conserving scheme. Both have been +!! modified to use general orthogonal coordinates as described in +!! Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second +!! order accurate, and allow for vanishingly small layer thicknesses. +!! The Arakawa and Hsu scheme globally conserves both total energy +!! and potential enstrophy in the limit of nondivergent flow. +!! Sadourny's energy conserving scheme conserves energy if the flow +!! is nondivergent or centered difference thickness fluxes are used. +!! +!! Two sets of boundary conditions have been coded in the +!! definition of relative vorticity. These are written as: +!! NOSLIP defined (in spherical coordinates): +!! relvort = dv/dx (east & west), with v = 0. +!! relvort = -sec(Q) * d(u cos(Q))/dy (north & south), with u = 0. +!! +!! NOSLIP not defined (free slip): +!! relvort = 0 (all boundaries) +!! +!! with Q temporarily defined as latitude. The free slip boundary +!! condition is much more natural on a C-grid. +!! +!! A small fragment of the grid is shown below: +!! \verbatim +!! +!! j+1 x ^ x ^ x At x: q, CoriolisBu +!! j+1 > o > o > At ^: v, CAv, vh +!! j x ^ x ^ x At >: u, CAu, uh, a, b, c, d +!! j > o > o > At o: h, KE +!! j-1 x ^ x ^ x +!! i-1 i i+1 At x & ^: +!! i i+1 At > & o: +!! \endverbatim +!! +!! The boundaries always run through q grid points (x). + end module MOM_CoriolisAdv diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 912dde87a7..0f068da864 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -1,44 +1,7 @@ +!> Solve the layer continuity equation. module MOM_continuity -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** - -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg and Alistair Adcroft, September 2006. * -!* * -!* This file contains the driver routine which selects which * -!* continuity solver will be called, based on run-time input. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v, vh * -!* j x ^ x ^ x At >: u, uh * -!* j > o > o > At o: h, hin * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** + +! This file is part of MOM6. See LICENSE.md for the license. use MOM_continuity_PPM, only : continuity_PPM, continuity_PPM_init use MOM_continuity_PPM, only : continuity_PPM_end, continuity_PPM_CS @@ -57,85 +20,66 @@ module MOM_continuity public continuity, continuity_init, continuity_end -integer :: id_clock_pass, id_clock_vertvisc - +!> Control structure for mom_continuity type, public :: continuity_CS ; private - integer :: continuity_scheme ! CONTINUITY_SCHEME selects the discretization - ! for the continuity solver. Valid values are: - ! PPM - A directionally split peicewise - ! parabolic reconstruction solver. - ! The default, PPM, seems most appropriate for - ! use with our current time-splitting strategies. - type(continuity_PPM_CS), pointer :: PPM_CSp => NULL() + integer :: continuity_scheme !< Selects the discretization for the continuity solver. + !! Valid values are: + !! - PPM - A directionally split piecewise parabolic reconstruction solver. + !! The default, PPM, seems most appropriate for use with our current + !! time-splitting strategies. + type(continuity_PPM_CS), pointer :: PPM_CSp => NULL() !< Control structure for mom_continuity_ppm end type continuity_CS -integer, parameter :: PPM_SCHEME = 1 -character(len=20), parameter :: PPM_STRING = "PPM" +integer, parameter :: PPM_SCHEME = 1 !< Enumerated constant to select PPM +character(len=20), parameter :: PPM_STRING = "PPM" !< String to select PPM contains +!> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme, +!! based on Lin (1994). subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) - type(ocean_grid_type), intent(inout) :: G - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh - real, intent(in) :: dt - type(continuity_CS), pointer :: CS - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt - type(ocean_OBC_type), pointer, optional :: OBC - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout), optional :: u_cor_aux - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout), optional :: v_cor_aux - type(BT_cont_type), pointer, optional :: BT_cont -! This subroutine time steps the layer thicknesses, using a monotonically -! limit, directionally split PPM scheme, based on Lin (1994). - -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) hin - Initial layer thickness, in m. -! (out) h - Final layer thickness, in m. -! (out) uh - Volume flux through zonal faces = u*h*dy, m3 s-1. -! (out) vh - Volume flux through meridional faces = v*h*dx, -! in m3 s-1. -! (in) dt - Time increment in s. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! continuity_init. -! (in, opt) uhbt - The summed volume flux through zonal faces, m3 s-1. -! (in, opt) vhbt - The summed volume flux through meridional faces, m3 s-1. -! (in, opt) OBC - This open boundary condition type specifies whether, where, -! and what open boundary conditions are used. -! (in, opt) visc_rem_u - Both the fraction of the momentum originally in a -! (in, opt) visc_rem_v - layer that remains after a time-step of viscosity, -! and the fraction of a time-step's worth of a -! barotropic acceleration that a layer experiences -! after viscosity is applied, in the zonal (_u) and -! meridional (_v) directions. Nondimensional between -! 0 (at the bottom) and 1 (far above the bottom). -! (out, opt) u_cor - The zonal velocities that give uhbt as the depth- -! integrated transport, in m s-1. -! (out, opt) v_cor - The meridional velocities that give vhbt as the -! depth-integrated transport, in m s-1. -! (in, opt) uhbt_aux - A second set of summed volume fluxes through zonal -! (in, opt) vhbt_aux - and meridional faces, both in m3 s-1. -! (out, opt) u_cor_aux - The zonal and meridional velocities that give uhbt_aux -! (out, opt) v_cor_aux - and vhbt_aux as the depth-integrated transports, -! both in m s-1. -! (out, opt) BT_cont - A structure with elements that describe the effective -! open face areas as a function of barotropic flow. + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< Zonal velocity, in m/s. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< Meridional velocity, in m/s. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin !< Initial layer thickness, in m or kg/m2. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Final layer thickness, in m or kg/m2. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Volume flux through zonal faces = + !! u*h*dy, in m3/s. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh !< Volume flux through meridional faces = + !! v*h*dx, in m3/s. + real, intent(in) :: dt !< Time increment, in s. + type(continuity_CS), pointer :: CS !< Control structure for mom_continuity. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt !< The vertically summed volume + !! flux through zonal faces, in m3/s. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt !< The vertically summed volume + !! flux through meridional faces, in m3/s. + type(ocean_OBC_type), pointer, optional :: OBC !< Open boundaries control structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u !< Both the fraction of + !! zonal momentum that remains after a time-step of viscosity, and the fraction of a time-step's + !! worth of a barotropic acceleration that a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v !< Both the fraction of + !! meridional momentum that remains after a time-step of viscosity, and the fraction of a time-step's + !! worth of a barotropic acceleration that a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor !< The zonal velocities that + !! give uhbt as the depth-integrated transport, in m/s. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor !< The meridional velocities that + !! give vhbt as the depth-integrated transport, in m/s. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux !< A second summed zonal + !! volume flux in m3/s. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux !< A second summed meridional + !! volume flux in m3/s. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout), optional :: u_cor_aux !< The zonal velocities + !! that give uhbt_aux as the depth-integrated transport, in m/s. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout), optional :: v_cor_aux !< The meridional velocities + !! that give vhbt_aux as the depth-integrated transport, in m/s. + type(BT_cont_type), pointer, optional :: BT_cont !< A structure with elements + !! that describe the effective open face areas as a function of barotropic flow. + if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, & "MOM_continuity: Either both visc_rem_u and visc_rem_v or neither"// & " one must be present in call to continuity.") @@ -163,21 +107,14 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & end subroutine continuity +!> Initializes continuity_cs subroutine continuity_init(Time, G, GV, param_file, diag, CS) - type(time_type), target, intent(in) :: Time - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - type(param_file_type), intent(in) :: param_file - type(diag_ctrl), target, intent(inout) :: diag - type(continuity_CS), pointer :: CS -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + type(time_type), target, intent(in) :: Time !< Current model time. + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Parameter file handles. + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(continuity_CS), pointer :: CS !< Control structure for mom_continuity. ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_continuity" ! This module's name. @@ -215,8 +152,9 @@ subroutine continuity_init(Time, G, GV, param_file, diag, CS) end subroutine continuity_init +!> Destructor for continuity_cs. subroutine continuity_end(CS) - type(continuity_CS), pointer :: CS + type(continuity_CS), pointer :: CS !< Control structure for mom_continuity. if (CS%continuity_scheme == PPM_SCHEME) then call continuity_PPM_end(CS%PPM_CSp) diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index a0097250cb..46820a7889 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -1,23 +1,7 @@ +!> Solve the layer continuity equation using the PPM method for layer fluxes. module MOM_continuity_PPM -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** + +! This file is part of MOM6. See LICENSE.md for the license. use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : time_type, diag_ctrl @@ -37,56 +21,60 @@ module MOM_continuity_PPM integer :: id_clock_update, id_clock_correct +!> Control structure for mom_continuity_ppm type, public :: continuity_PPM_CS ; private - type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the - ! timing of diagnostic output. - logical :: upwind_1st ! If true, use a first-order upwind scheme. - logical :: monotonic ! If true, use the Colella & Woodward monotonic - ! limiter; otherwise use a simple positive - ! definite limiter. - logical :: simple_2nd ! If true, use a simple second order (arithmetic - ! mean) interpolation of the edge values instead - ! of the higher order interpolation. - real :: tol_eta ! The tolerance for free-surface height - ! discrepancies between the barotropic solution and - ! the sum of the layer thicknesses, in m. - real :: tol_vel ! The tolerance for barotropic velocity - ! discrepancies between the barotropic solution and - ! the sum of the layer thicknesses, in m s-1. - real :: tol_eta_aux ! The tolerance for free-surface height - ! discrepancies between the barotropic solution and - ! the sum of the layer thicknesses when calculating - ! the auxiliary corrected velocities, in m. - real :: CFL_limit_adjust ! The maximum CFL of the adjusted velocities, ND. - logical :: aggress_adjust ! If true, allow the adjusted velocities to have a - ! relative CFL change up to 0.5. False by default. - logical :: vol_CFL ! If true, use the ratio of the open face lengths - ! to the tracer cell areas when estimating CFL - ! numbers. Without aggress_adjust, the default is - ! false; it is always true with. - logical :: better_iter ! If true, stop corrective iterations using a - ! velocity-based criterion and only stop if the - ! iteration is better than all predecessors. - logical :: use_visc_rem_max ! If true, use more appropriate limiting bounds - ! for corrections in strongly viscous columns. - logical :: marginal_faces ! If true, use the marginal face areas from the - ! continuity solver for use as the weights in the - ! barotropic solver. Otherwise use the transport - ! averaged areas. + type(diag_ctrl), pointer :: diag !< Diagnostics control structure. + logical :: upwind_1st !< If true, use a first-order upwind scheme. + logical :: monotonic !< If true, use the Colella & Woodward monotonic + !! limiter; otherwise use a simple positive + !! definite limiter. + logical :: simple_2nd !< If true, use a simple second order (arithmetic + !! mean) interpolation of the edge values instead + !! of the higher order interpolation. + real :: tol_eta !< The tolerance for free-surface height + !! discrepancies between the barotropic solution and + !! the sum of the layer thicknesses, in m. + real :: tol_vel !< The tolerance for barotropic velocity + !! discrepancies between the barotropic solution and + !! the sum of the layer thicknesses, in m s-1. + real :: tol_eta_aux !< The tolerance for free-surface height + !! discrepancies between the barotropic solution and + !! the sum of the layer thicknesses when calculating + !! the auxiliary corrected velocities, in m. + real :: CFL_limit_adjust !< The maximum CFL of the adjusted velocities, ND. + logical :: aggress_adjust !< If true, allow the adjusted velocities to have a + !! relative CFL change up to 0.5. False by default. + logical :: vol_CFL !< If true, use the ratio of the open face lengths + !! to the tracer cell areas when estimating CFL + !! numbers. Without aggress_adjust, the default is + !! false; it is always true with. + logical :: better_iter !< If true, stop corrective iterations using a + !! velocity-based criterion and only stop if the + !! iteration is better than all predecessors. + logical :: use_visc_rem_max !< If true, use more appropriate limiting bounds + !! for corrections in strongly viscous columns. + logical :: marginal_faces !< If true, use the marginal face areas from the + !! continuity solver for use as the weights in the + !! barotropic solver. Otherwise use the transport + !! averaged areas. end type continuity_PPM_CS +!> A container for loop bounds type :: loop_bounds_type ; private + !>@{ + !! Loop bounds integer :: ish, ieh, jsh, jeh + !>@} end type loop_bounds_type contains -!> This subroutine time steps the layer thicknesses, using a monotonically -!! limit, directionally split PPM scheme, based on Lin (1994). In the following -!! documentation, H is used for the units of thickness (usually m or kg m-2.) +!> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, +!! based on Lin (1994). subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) + ! In the following documentation, H is used for the units of thickness (usually m or kg m-2.) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(continuity_PPM_CS), pointer :: CS !< Module's control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< Zonal velocity, in m s-1. @@ -99,42 +87,37 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, !! v*h*dx, H m2 s-1. real, intent(in) :: dt !< Time increment in s. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt !< - !! The summed volume flux through zonal faces, H m2 s-1. - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt !< - !! The summed volume flux through meridional faces, H m2 s-1. - type(ocean_OBC_type), pointer, optional :: OBC !< - !! This open boundary condition type specifies whether, where, - !! and what open boundary conditions are used. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u !< - !! Both the fraction of the momentum originally in a layer that remains after a time-step - !! of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that - !! a layer experiences after viscosity is applied, in the zonal (_u) - !! direction. Nondimensional between 0 (at the bottom) and 1 (far above the bottom). - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v !< - !! Same as above for meridional direction. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor !< - !! The zonal velocities that give uhbt as the depth- - !! integrated transport, in m s-1. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor !< - !! The meridional velocities that give vhbt as the - !! depth-integrated transport, in m s-1. - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux !< - !! A second set of summed volume fluxes through zonal faces, in H m2 s-1. - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux !< - !! A second set of summed volume fluxes through meridional faces, in H m2 s-1. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux !< - !! The zonal velocities that give uhbt_aux as the - !! depth-integrated transports, in m s-1. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux !< - !! The meridional velocities that give vhbt_aux as the - !! depth-integrated transports, in m s-1. - type(BT_cont_type), pointer, optional :: BT_cont !< - !! A structure with elements that describe the effective - !! open face areas as a function of barotropic flow. - - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - h_input ! Left and right face thicknesses, in H. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt + !< The summed volume flux through zonal faces, H m2 s-1. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt + !< The summed volume flux through meridional faces, H m2 s-1. + type(ocean_OBC_type), pointer, optional :: OBC !< Open boundaries control structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u + !< The fraction of zonal momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that + !! a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v + !< The fraction of meridional momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that + !! a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor + !< The zonal velocities that give uhbt as the depth-integrated transport, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor + !< The meridional velocities that give vhbt as the depth-integrated transport, in m s-1. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux + !< A second set of summed volume fluxes through zonal faces, in H m2 s-1. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux + !< A second set of summed volume fluxes through meridional faces, in H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux + !< The zonal velocities that give uhbt_aux as the depth-integrated transports, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux + !< The meridional velocities that give vhbt_aux as the depth-integrated transports, in m s-1. + type(BT_cont_type), pointer, optional :: BT_cont !< A structure with + !! elements that describe the effective open face areas as a function of barotropic flow. + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_input ! Left and right face thicknesses, in H. real :: h_min type(loop_bounds_type) :: LB integer :: is, ie, js, je, nz, stencil @@ -321,8 +304,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, end subroutine continuity_PPM -!> This subroutine calculates the mass or volume fluxes through the zonal -!! faces, and other related quantities. +!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities. subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & visc_rem_u, u_cor, uhbt_aux, u_cor_aux, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. @@ -335,35 +317,28 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & real, intent(in) :: dt !< Time increment in s. type(continuity_PPM_CS), pointer :: CS !< This module's control structure. type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. - type(ocean_OBC_type), pointer, optional :: OBC !< - !! This open boundary condition type specifies whether, where, - !! and what open boundary conditions are used. + type(ocean_OBC_type), pointer, optional :: OBC !< Open boundaries control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u !< - !! Both the fraction of the momentum originally in a - !! layer that remains after a time-step of viscosity, - !! and the fraction of a time-step's worth of a - !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between - !! 0 (at the bottom) and 1 (far above the bottom). - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt !< - !! The summed volume flux through zonal faces, H m2 s-1. - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux !< - !! A second set of summed volume fluxes through zonal - !! faces, in H m2 s-1. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor !< - !! The zonal velocitiess (u with a barotropic correction) + !< The fraction of zonal momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that + !! a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt + !< The summed volume flux through zonal faces, H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux + !< A second set of summed volume fluxes through zonal faces, in H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor + !< The zonal velocitiess (u with a barotropic correction) !! that give uhbt as the depth-integrated transport, m s-1. - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux !< - !! The zonal velocities (u with a barotropic correction) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux + !< The zonal velocities (u with a barotropic correction) !! that give uhbt_aux as the depth-integrated transports, in m s-1. type(BT_cont_type), pointer, optional :: BT_cont !< - !! A structure with elements that describe the effective + !< A structure with elements that describe the effective !! open face areas as a function of barotropic flow. - - real, dimension(SZIB_(G),SZK_(G)) :: & - duhdu ! Partial derivative of uh with u, in H m. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - h_L, h_R ! Left and right face thicknesses, in H. + ! Local variables + real, dimension(SZIB_(G),SZK_(G)) :: duhdu ! Partial derivative of uh with u, in H m. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R ! Left and right face thicknesses, in H. real, dimension(SZIB_(G)) :: & du, & ! Corrective barotropic change in the velocity, in m s-1. du_min_CFL, & ! Min/max limits on du correction @@ -611,16 +586,16 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & end subroutine zonal_mass_flux -!> This subroutines evaluates the zonal mass or volume fluxes in a layer. +!> Evaluates the zonal mass or volume fluxes in a layer. subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, j, & ish, ieh, do_I, vol_CFL) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity, in m s-1. real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the - !! momentum originally in a layer that remains after a time-step - !! of viscosity, and the fraction of a time-step's worth of a barotropic - !! acceleration that a layer experiences after viscosity is applied. - !! Nondimensional between 0 (at the bottom) and 1 (far above the bottom). + !! momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic + !! acceleration that a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZI_(G)), intent(in) :: h !< Layer thickness, in H. real, dimension(SZI_(G)), intent(in) :: h_L !< Left thickness, in H. real, dimension(SZI_(G)), intent(in) :: h_R !< Right thickness, in H. @@ -634,8 +609,8 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, j, & integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I !< Which i values to work on. logical, intent(in) :: vol_CFL !< If true, rescale the - !! ratio of face areas to the cell areas when estimating the CFL number. - + !! ratio of face areas to the cell areas when estimating the CFL number. + ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing, ND. real :: curv_3 ! A measure of the thickness curvature over a grid length, ! with the same units as h_in. @@ -667,8 +642,7 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, j, & end subroutine zonal_flux_layer -!> This subroutines sets the effective interface thickness at each zonal -!! velocity point. +!> Sets the effective interface thickness at each zonal velocity point. subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, & marginal, visc_rem_u) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. @@ -694,9 +668,9 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, & !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). - + ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing, ND. real :: curv_3 ! A measure of the thickness curvature over a grid length, ! with the same units as h_in. @@ -745,7 +719,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, & end subroutine zonal_face_thickness -!> This subroutine returns the barotropic velocity adjustment that gives the +!> Returns the barotropic velocity adjustment that gives the !! desired barotropic (layer-summed) transport. subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & du, du_max_CFL, du_min_CFL, dt, G, CS, visc_rem, & @@ -763,7 +737,7 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZIB_(G)), intent(in), optional :: uhbt !< !! The summed volume flux through zonal faces, H m2 s-1. @@ -783,13 +757,13 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & integer, intent(in) :: ish !< Start of index range. integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I_in !< - !! A logical flag indiciating which I values to work on. + !! A logical flag indicating which I values to work on. logical, intent(in), optional :: full_precision !< !! A flag indicating how carefully to iterate. The !! default is .true. (more accurate). real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout), optional :: uh_3d !< !! Volume flux through zonal faces = u*h*dy, H m2 s-1. - + ! Local variables real, dimension(SZIB_(G),SZK_(G)) :: & uh_aux, & ! An auxiliary zonal volume flux, in H m s-1. duhdu ! Partial derivative of uh with u, in H m. @@ -908,9 +882,8 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & end subroutine zonal_flux_adjust -!> This subroutine sets of a structure that describes the zonal barotropic -!! volume or mass fluxes as a function of barotropic flow to agree closely with -!! the sum of the layer's transports. +!> Sets a structure that describes the zonal barotropic volume or mass fluxes as a +!! function of barotropic flow to agree closely with the sum of the layer's transports. subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, & du_max_CFL, du_min_CFL, dt, G, CS, visc_rem, & visc_rem_max, j, ish, ieh, do_I) @@ -940,15 +913,15 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZIB_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem. integer, intent(in) :: j !< Spatial index. integer, intent(in) :: ish !< Start of index range. integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I !< - !! A logical flag indiciating which I values to work on. - + !! A logical flag indicating which I values to work on. + ! Local variables real, dimension(SZIB_(G)) :: & du0, & ! The barotropic velocity increment that gives 0 transport, m s-1. duL, duR, & ! The barotropic velocity increments that give the westerly @@ -1080,8 +1053,7 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, end subroutine set_zonal_BT_cont -!> This subroutine calculates the mass or volume fluxes through the meridional -!! faces, and other related quantities. +!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities. subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & visc_rem_v, v_cor, vhbt_aux, v_cor_aux, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. @@ -1117,7 +1089,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & !! that give vhbt_aux as the depth-integrated transports, in m s-1. type(BT_cont_type), pointer, optional :: BT_cont !< !! A structure with elements that describe the effective - + ! Local variables real, dimension(SZI_(G),SZK_(G)) :: & dvhdv ! Partial derivative of vh with v, in m2. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -1366,16 +1338,16 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & end subroutine meridional_mass_flux -!> This subroutines evaluates the meridional mass or volume fluxes in a layer. +!> Evaluates the meridional mass or volume fluxes in a layer. subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, J, & ish, ieh, do_I, vol_CFL) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity, in m s-1. real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the - !! momentum originally in a layer that remains after a time-step - !! of viscosity, and the fraction of a time-step's worth of a barotropic - !! acceleration that a layer experiences after viscosity is applied. - !! Nondimensional between 0 (at the bottom) and 1 (far above the bottom). + !! momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic + !! acceleration that a layer experiences after viscosity is applied. + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Layer thickness used to !! calculate fluxes, in H. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_L !< Left thickness in the @@ -1392,8 +1364,8 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, J, & integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on. logical, intent(in) :: vol_CFL !< If true, rescale the - !! ratio of face areas to the cell areas when estimating the CFL number. - + !! ratio of face areas to the cell areas when estimating the CFL number. + ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing, ND. real :: curv_3 ! A measure of the thickness curvature over a grid length, ! with the same units as h_in. @@ -1426,8 +1398,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, J, & end subroutine merid_flux_layer -!> This subroutines sets the effective interface thickness at each meridional -!! velocity point. +!> Sets the effective interface thickness at each meridional velocity point. subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, LB, vol_CFL, & marginal, visc_rem_v) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. @@ -1453,9 +1424,9 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, LB, vol_CFL, & !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). - + ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing, ND. real :: curv_3 ! A measure of the thickness curvature over a grid length, ! with the same units as h_in. @@ -1506,8 +1477,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, LB, vol_CFL, & end subroutine merid_face_thickness -!> This subroutine returns the barotropic velocity adjustment that gives the -!! desired barotropic (layer-summed) transport. +!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport. subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, & dv, dv_max_CFL, dv_min_CFL, dt, G, CS, visc_rem, & j, ish, ieh, do_I_in, full_precision, vh_3d) @@ -1524,7 +1494,7 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZI_(G)), intent(in), optional :: vhbt !< !! The summed volume flux through meridional faces, H m2 s-1. @@ -1544,15 +1514,15 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 integer, intent(in) :: ish !< Start of index range. integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), intent(in) :: do_I_in !< - !! A logical flag indiciating which I values to work on. + !! A logical flag indicating which I values to work on. logical, intent(in), optional :: full_precision !< !! full_precision - A flag indicating how carefully to iterate. The !! default is .true. (more accurate). real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout), optional :: vh_3d !< !! Volume flux through meridional faces = v*h*dx, H m2 s-1. - + ! Local variables real, dimension(SZI_(G),SZK_(G)) :: & - vh_aux, & ! An auxiliary meridonal volume flux, in H m s-1. + vh_aux, & ! An auxiliary meridional volume flux, in H m s-1. dvhdv ! Partial derivative of vh with v, in H m. real, dimension(SZI_(G)) :: & vh_err, & ! Difference between vhbt and the summed vh, in H m2 s-1. @@ -1669,9 +1639,8 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 end subroutine meridional_flux_adjust -!> This subroutine sets of a structure that describes the meridional -!! barotropic volume or mass fluxes as a function of barotropic flow to agree -!! closely with the sum of the layer's transports. +!> Sets of a structure that describes the meridional barotropic volume or mass fluxes as a +!! function of barotropic flow to agree closely with the sum of the layer's transports. subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, & dv_max_CFL, dv_min_CFL, dt, G, CS, visc_rem, & visc_rem_max, j, ish, ieh, do_I) @@ -1701,15 +1670,15 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, !! layer that remains after a time-step of viscosity, !! and the fraction of a time-step's worth of a !! barotropic acceleration that a layer experiences - !! after viscosity is applied. Nondimensional between + !! after viscosity is applied. Non-dimensional between !! 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZI_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem. integer, intent(in) :: j !< Spatial index. integer, intent(in) :: ish !< Start of index range. integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), intent(in) :: do_I !< - !! A logical flag indiciating which I values to work on. - + !! A logical flag indicating which I values to work on. + ! Local variables real, dimension(SZI_(G)) :: & dv0, & ! The barotropic velocity increment that gives 0 transport, m s-1. dvL, dvR, & ! The barotropic velocity increments that give the southerly @@ -1838,7 +1807,7 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, end subroutine set_merid_BT_cont -!> This subroutine calculates left/right edge values for PPM reconstruction. +!> Calculates left/right edge values for PPM reconstruction. subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd) type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness, in H. @@ -1856,7 +1825,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. -! Local variables with useful mnemonic names. + ! Local variables with useful mnemonic names. real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes. real, parameter :: oneSixth = 1./6. real :: h_ip1, h_im1 @@ -1929,7 +1898,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ return end subroutine PPM_reconstruction_x -!> This subroutine calculates left/right edge valus for PPM reconstruction. +!> Calculates left/right edge values for PPM reconstruction. subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd) type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness, in H. @@ -1947,7 +1916,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. -! Local variables with useful mnemonic names. + ! Local variables with useful mnemonic names. real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes. real, parameter :: oneSixth = 1./6. real :: h_jp1, h_jm1 @@ -2077,7 +2046,7 @@ subroutine PPM_limit_CW84(h_in, h_L, h_R, G, iis, iie, jis, jie) integer, intent(in) :: jis !< Start of j index range. integer, intent(in) :: jie !< End of j index range. -! Local variables + ! Local variables real :: h_i, RLdiff, RLdiff2, RLmean, FunFac character(len=256) :: mesg integer :: i,j @@ -2101,7 +2070,7 @@ subroutine PPM_limit_CW84(h_in, h_L, h_R, G, iis, iie, jis, jie) return end subroutine PPM_limit_CW84 -!> Return the maxiumum ratio of a/b or maxrat. +!> Return the maximum ratio of a/b or maxrat. function ratio_max(a, b, maxrat) result(ratio) real, intent(in) :: a !< Numerator real, intent(in) :: b !< Denominator @@ -2115,7 +2084,7 @@ function ratio_max(a, b, maxrat) result(ratio) endif end function ratio_max -!> This include declares and sets the variable "version". +!> Initializes continuity_ppm_cs subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS) type(time_type), target, intent(in) :: Time !< Time increment in s. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -2125,7 +2094,7 @@ subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS) type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to !! regulate diagnostic output. type(continuity_PPM_CS), pointer :: CS !< Module's control structure. - +!> This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_continuity_PPM" ! This module's name. @@ -2210,35 +2179,16 @@ subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS) end subroutine continuity_PPM_init -!> Deconstructor for CS object. +!> Destructor for continuity_ppm_cs subroutine continuity_PPM_end(CS) type(continuity_PPM_CS), pointer :: CS !< Module's control structure. deallocate(CS) end subroutine continuity_PPM_end -!> \class MOM_continuity_PPM -!!*******+*********+*********+*********+*********+*********+*********+** -!! * -!! By Robert Hallberg and Alistair Adcroft, September 2006 - . * -!! * -!! This program contains the subroutine that advects layer * -!! thickness. The scheme here uses a Piecewise-Parabolic method with * -!! a positive definite limiter. * -!! * -!! Macros written all in capital letters are defined in MOM_memory.h. * -!! * -!! A small fragment of the grid is shown below: * -!! * -!! j+1 x ^ x ^ x At x: q * -!! j+1 > o > o > At ^: v, vh * -!! j x ^ x ^ x At >: u, uh * -!! j > o > o > At o: h, hin * -!! j-1 x ^ x ^ x * -!! i-1 i i+1 At x & ^: * -!! i i+1 At > & o: * -!! * -!! The boundaries always run through q grid points (x). * -!! * -!!*******+*********+*********+*********+*********+*********+*********+** +!> \namespace mom_continuity_ppm +!! +!! This module contains the subroutines that advect layer +!! thickness. The scheme here uses a Piecewise-Parabolic method with +!! a positive definite limiter. end module MOM_continuity_PPM diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 3f7969f0b4..f7a1aac006 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -109,7 +109,7 @@ module MOM_forcing_type p_surf_full => NULL(), & !< Pressure at the top ocean interface (Pa). !! if there is sea-ice, then p_surf_flux is at ice-ocean interface p_surf => NULL(), & !< Pressure at the top ocean interface (Pa) as used - !! to drive the ocean model. If p_surf is limited, + !! to drive the ocean model. If p_surf is limited, !! p_surf may be smaller than p_surf_full, !! otherwise they are the same. p_surf_SSH => NULL() !< Pressure at the top ocean interface that is used @@ -122,6 +122,12 @@ module MOM_forcing_type TKE_tidal => NULL(), & !< tidal energy source driving mixing in bottom boundary layer (W/m^2) ustar_tidal => NULL() !< tidal contribution to bottom ustar (m/s) + ! iceberg related inputs + real, pointer, dimension(:,:) :: & + ustar_berg => NULL(),& !< iceberg contribution to top ustar (m/s) + area_berg => NULL(),& !< area of ocean surface covered by icebergs (m2/m2) + mass_berg => NULL() !< mass of icebergs (kg/m2) + ! land ice-shelf related inputs real, pointer, dimension(:,:) :: & ustar_shelf => NULL(), & !< friction velocity under ice-shelves (m/s) @@ -196,7 +202,7 @@ module MOM_forcing_type integer :: id_heat_content_cond = -1, id_heat_content_surfwater= -1 integer :: id_heat_content_vprec = -1, id_heat_content_massout = -1 integer :: id_heat_added = -1, id_heat_content_massin = -1 - integer :: id_hfrainds = -1, id_hfrunoffds = -1 + integer :: id_hfrainds = -1, id_hfrunoffds = -1 ! global area integrated heat flux diagnostic handles @@ -245,6 +251,15 @@ module MOM_forcing_type ! clock id handle integer :: id_clock_forcing + ! iceberg id handle + integer :: id_ustar_berg + integer :: id_area_berg + integer :: id_mass_berg + + !Iceberg + Ice shelf + integer :: id_ustar_ice_cover + integer :: id_frac_ice_cover + end type forcing_diags contains @@ -629,7 +644,7 @@ subroutine extractFluxes2d(G, GV, fluxes, optics, nsw, dt, !! Units of net_heat are (K * H). real, dimension(SZI_(G),SZJ_(G)), intent(out) :: net_salt !< surface salt flux into the ocean accumulated !! over a time step (ppt * H) - real, dimension(:,:,:), intent(out) :: pen_SW_bnd !! penetrating shortwave flux, split into bands. + real, dimension(:,:,:), intent(out) :: pen_SW_bnd !< penetrating shortwave flux, split into bands. !! Units (deg K * H) & array size nsw x SZI_(G), !! where nsw=number of SW bands in pen_SW_bnd. !! This heat flux is not in net_heat. @@ -864,14 +879,13 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, haloshift) end subroutine MOM_forcing_chksum -!> Write out values of the fluxes arrays at the i,j location +!> Write out values of the fluxes arrays at the i,j location. This is a debugging tool. subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) - - type(forcing), intent(in) :: fluxes !< fluxes type - type(ocean_grid_type), intent(in) :: G !< grid type - character(len=*), intent(in) :: mesg !< message - integer, intent(in) :: i, j !< horizontal indices - + type(forcing), intent(in) :: fluxes !< Fluxes type + type(ocean_grid_type), intent(in) :: G !< Grid type + character(len=*), intent(in) :: mesg !< Message + integer, intent(in) :: i !< i-index + integer, intent(in) :: j !< j-index write(0,'(2a)') 'MOM_forcing_type, forcing_SinglePointPrint: Called from ',mesg write(0,'(a,2es15.3)') 'MOM_forcing_type, forcing_SinglePointPrint: lon,lat = ',G%geoLonT(i,j),G%geoLatT(i,j) @@ -910,9 +924,10 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) call locMsg(fluxes%heat_content_cond,'heat_content_massout') contains + !> Format and write a message depending on associated state of array subroutine locMsg(array,aname) - real, dimension(:,:), pointer :: array - character(len=*) :: aname + real, dimension(:,:), pointer :: array !< Array to write element from + character(len=*) :: aname !< Name of array if (associated(array)) then write(0,'(3a,es15.3)') 'MOM_forcing_type, forcing_SinglePointPrint: ',trim(aname),' = ',array(i,j) @@ -952,6 +967,21 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles) handles%id_ustar = register_diag_field('ocean_model', 'ustar', diag%axesT1, Time, & 'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', 'meter second-1') + handles%id_ustar_berg = register_diag_field('ocean_model', 'ustar_berg', diag%axesT1, Time, & + 'Friction velocity below iceberg ', 'meter second-1') + + handles%id_area_berg = register_diag_field('ocean_model', 'area_berg', diag%axesT1, Time, & + 'Area of grid cell covered by iceberg ', 'm2/m2') + + handles%id_mass_berg = register_diag_field('ocean_model', 'mass_berg', diag%axesT1, Time, & + 'Mass of icebergs ', 'kg/m2') + + handles%id_ustar_ice_cover = register_diag_field('ocean_model', 'ustar_ice_cover', diag%axesT1, Time, & + 'Friction velocity below iceberg and ice shelf together', 'meter second-1') + + handles%id_frac_ice_cover = register_diag_field('ocean_model', 'frac_ice_cover', diag%axesT1, Time, & + 'Area of grid cell below iceberg and ice shelf together ', 'm2/m2') + handles%id_psurf = register_diag_field('ocean_model', 'p_surf', diag%axesT1, Time, & 'Pressure at ice-ocean or atmosphere-ocean interface', 'Pascal', cmor_field_name='pso',& cmor_long_name='Sea Water Pressure at Sea Water Surface', cmor_units='Pa', & @@ -1146,7 +1176,7 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles) handles%id_hfrunoffds = register_diag_field('ocean_model', 'hfrunoffds', & diag%axesT1, Time, 'Heat content (relative to 0C) of liquid+solid runoff into ocean', 'W m-2',& standard_name='temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water') - + handles%id_heat_content_lprec = register_diag_field('ocean_model', 'heat_content_lprec', & diag%axesT1,Time,'Heat content (relative to 0degC) of liquid precip entering ocean', & 'W/m^2') @@ -1628,7 +1658,7 @@ subroutine forcing_accumulate(flux_tmp, fluxes, dt, G, wt2) fluxes%ustar_shelf(i,j) = flux_tmp%ustar_shelf(i,j) enddo ; enddo endif - + if (associated(fluxes%iceshelf_melt) .and. associated(flux_tmp%iceshelf_melt)) then do i=isd,ied ; do j=jsd,jed fluxes%iceshelf_melt(i,j) = flux_tmp%iceshelf_melt(i,j) @@ -1690,6 +1720,16 @@ subroutine mech_forcing_diags(fluxes, dt, G, diag, handles) call post_data(handles%id_tauy, fluxes%tauy, diag) if ((handles%id_ustar > 0) .and. ASSOCIATED(fluxes%ustar)) & call post_data(handles%id_ustar, fluxes%ustar, diag) + if ((handles%id_ustar_berg > 0) .and. ASSOCIATED(fluxes%ustar_berg)) & + call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) + if ((handles%id_area_berg > 0) .and. ASSOCIATED(fluxes%area_berg)) & + call post_data(handles%id_area_berg, fluxes%area_berg, diag) + if ((handles%id_mass_berg > 0) .and. ASSOCIATED(fluxes%mass_berg)) & + call post_data(handles%id_mass_berg, fluxes%mass_berg, diag) + if ((handles%id_frac_ice_cover > 0) .and. ASSOCIATED(fluxes%frac_shelf_h)) & + call post_data(handles%id_frac_ice_cover, fluxes%frac_shelf_h, diag) + if ((handles%id_ustar_ice_cover > 0) .and. ASSOCIATED(fluxes%ustar_shelf)) & + call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) endif @@ -1987,31 +2027,31 @@ subroutine forcing_diagnostics(fluxes, state, dt, G, diag, handles) endif ! for OMIP, hfrunoffds = heat content of liquid plus frozen runoff - if (handles%id_hfrunoffds > 0) then - sum(:,:) = 0.0 - if(ASSOCIATED(fluxes%heat_content_lrunoff)) then + if (handles%id_hfrunoffds > 0) then + sum(:,:) = 0.0 + if(ASSOCIATED(fluxes%heat_content_lrunoff)) then sum(:,:) = sum(:,:) + fluxes%heat_content_lrunoff(:,:) - endif - if(ASSOCIATED(fluxes%heat_content_frunoff)) then + endif + if(ASSOCIATED(fluxes%heat_content_frunoff)) then sum(:,:) = sum(:,:) + fluxes%heat_content_frunoff(:,:) - endif + endif call post_data(handles%id_hfrunoffds, sum, diag) - endif + endif - ! for OMIP, hfrainds = heat content of lprec + fprec + cond - if (handles%id_hfrainds > 0) then - sum(:,:) = 0.0 - if(ASSOCIATED(fluxes%heat_content_lprec)) then + ! for OMIP, hfrainds = heat content of lprec + fprec + cond + if (handles%id_hfrainds > 0) then + sum(:,:) = 0.0 + if(ASSOCIATED(fluxes%heat_content_lprec)) then sum(:,:) = sum(:,:) + fluxes%heat_content_lprec(:,:) - endif - if(ASSOCIATED(fluxes%heat_content_fprec)) then + endif + if(ASSOCIATED(fluxes%heat_content_fprec)) then sum(:,:) = sum(:,:) + fluxes%heat_content_fprec(:,:) - endif - if(ASSOCIATED(fluxes%heat_content_cond)) then + endif + if(ASSOCIATED(fluxes%heat_content_cond)) then sum(:,:) = sum(:,:) + fluxes%heat_content_cond(:,:) - endif + endif call post_data(handles%id_hfrainds, sum, diag) - endif + endif if ((handles%id_LwLatSens > 0) .and. ASSOCIATED(fluxes%lw) .and. & ASSOCIATED(fluxes%latent) .and. ASSOCIATED(fluxes%sens)) then @@ -2152,7 +2192,7 @@ subroutine forcing_diagnostics(fluxes, state, dt, G, diag, handles) call post_data(handles%id_netFWGlobalScl, fluxes%netFWGlobalScl, diag) - ! remamining boundary terms ================================================== + ! remaining boundary terms ================================================== if ((handles%id_psurf > 0) .and. ASSOCIATED(fluxes%p_surf)) & call post_data(handles%id_psurf, fluxes%p_surf, diag) @@ -2171,7 +2211,7 @@ end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type -subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, press) +subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, press, iceberg) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< Forcing fields structure logical, optional, intent(in) :: stress !< If present and true, allocate taux, tauy @@ -2180,6 +2220,7 @@ subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, p logical, optional, intent(in) :: heat !< If present and true, allocate heat fluxes logical, optional, intent(in) :: shelf !< If present and true, allocate fluxes for ice-shelf logical, optional, intent(in) :: press !< If present and true, allocate p_surf + logical, optional, intent(in) :: iceberg !< If present and true, allocate fluxes for icebergs ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -2232,11 +2273,19 @@ subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, p call myAlloc(fluxes%p_surf,isd,ied,jsd,jed, press) + !These fields should only on allocated when iceberg area is being passed through the coupler. + call myAlloc(fluxes%ustar_berg,isd,ied,jsd,jed, iceberg) + call myAlloc(fluxes%area_berg,isd,ied,jsd,jed, iceberg) + call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) contains + !> Allocates and zeroes-out array. subroutine myAlloc(array, is, ie, js, je, flag) - real, dimension(:,:), pointer :: array - integer, intent(in) :: is, ie, js, je !< Bounds + real, dimension(:,:), pointer :: array !< Array to be allocated + integer, intent(in) :: is !< Start i-index + integer, intent(in) :: ie !< End i-index + integer, intent(in) :: js !< Start j-index + integer, intent(in) :: je !< End j-index logical, optional, intent(in) :: flag !< Flag to indicate to allocate if (present(flag)) then @@ -2299,6 +2348,9 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%rigidity_ice_u)) deallocate(fluxes%rigidity_ice_u) if (associated(fluxes%rigidity_ice_v)) deallocate(fluxes%rigidity_ice_v) if (associated(fluxes%tr_fluxes)) deallocate(fluxes%tr_fluxes) + if (associated(fluxes%ustar_berg)) deallocate(fluxes%ustar_berg) + if (associated(fluxes%area_berg)) deallocate(fluxes%area_berg) + if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) end subroutine deallocate_forcing_type diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index cae09ffcbd..57293a65fa 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -1,24 +1,7 @@ +!> Provides the ocean grid type module MOM_grid -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** +! This file is part of MOM6. See LICENSE.md for the license. use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent @@ -32,124 +15,162 @@ module MOM_grid public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction public isPointInCell, hor_index_type +!> Ocean grid type. See mom_grid for details. type, public :: ocean_grid_type - type(MOM_domain_type), pointer :: Domain => NULL() - type(MOM_domain_type), pointer :: Domain_aux => NULL() ! A non-symmetric auxiliary domain type. - type(hor_index_type) :: HI - integer :: isc, iec, jsc, jec ! The range of the computational domain indices - integer :: isd, ied, jsd, jed ! and data domain indices at tracer cell centers. - integer :: isg, ieg, jsg, jeg ! The range of the global domain tracer cell indices. - integer :: IscB, IecB, JscB, JecB ! The range of the computational domain indices - integer :: IsdB, IedB, JsdB, JedB ! and data domain indices at tracer cell vertices. - integer :: IsgB, IegB, JsgB, JegB ! The range of the global domain vertex indices. - integer :: isd_global ! The values of isd and jsd in the global - integer :: jsd_global ! (decomposition invariant) index space. - integer :: idg_offset ! The offset between the corresponding global - integer :: jdg_offset ! and local array indices; add to local to get global. - integer :: ke ! The number of layers in the vertical. - logical :: symmetric ! True if symmetric memory is used. - logical :: nonblocking_updates ! If true, non-blocking halo updates are - ! allowed. The default is .false. (for now). - integer :: first_direction ! An integer that indicates which direction is - ! to be updated first in directionally split - ! parts of the calculation. This can be altered - ! during the course of the run via calls to - ! set_first_direction. + type(MOM_domain_type), pointer :: Domain => NULL() !< Ocean model domain + type(MOM_domain_type), pointer :: Domain_aux => NULL() !< A non-symmetric auxiliary domain type. + type(hor_index_type) :: HI !< Horizontal index ranges + + integer :: isc !< The start i-index of cell centers within the computational domain + integer :: iec !< The end i-index of cell centers within the computational domain + integer :: jsc !< The start j-index of cell centers within the computational domain + integer :: jec !< The end j-index of cell centers within the computational domain + + integer :: isd !< The start i-index of cell centers within the data domain + integer :: ied !< The end i-index of cell centers within the data domain + integer :: jsd !< The start j-index of cell centers within the data domain + integer :: jed !< The end j-index of cell centers within the data domain + + integer :: isg !< The start i-index of cell centers within the global domain + integer :: ieg !< The end i-index of cell centers within the global domain + integer :: jsg !< The start j-index of cell centers within the global domain + integer :: jeg !< The end j-index of cell centers within the global domain + + integer :: IscB !< The start i-index of cell vertices within the computational domain + integer :: IecB !< The end i-index of cell vertices within the computational domain + integer :: JscB !< The start j-index of cell vertices within the computational domain + integer :: JecB !< The end j-index of cell vertices within the computational domain + + integer :: IsdB !< The start i-index of cell vertices within the data domain + integer :: IedB !< The end i-index of cell vertices within the data domain + integer :: JsdB !< The start j-index of cell vertices within the data domain + integer :: JedB !< The end j-index of cell vertices within the data domain + + integer :: IsgB !< The start i-index of cell vertices within the global domain + integer :: IegB !< The end i-index of cell vertices within the global domain + integer :: JsgB !< The start j-index of cell vertices within the global domain + integer :: JegB !< The end j-index of cell vertices within the global domain + + integer :: isd_global !< The value of isd in the global index space (decompoistion invariant). + integer :: jsd_global !< The value of isd in the global index space (decompoistion invariant). + integer :: idg_offset !< The offset between the corresponding global and local i-indices. + integer :: jdg_offset !< The offset between the corresponding global and local j-indices. + integer :: ke !< The number of layers in the vertical. + logical :: symmetric !< True if symmetric memory is used. + logical :: nonblocking_updates !< If true, non-blocking halo updates are + !! allowed. The default is .false. (for now). + integer :: first_direction !< An integer that indicates which direction is + !! to be updated first in directionally split + !! parts of the calculation. This can be altered + !! during the course of the run via calls to + !! set_first_direction. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - mask2dT, & ! 0 for land points and 1 for ocean points on the h-grid. Nd. - geoLatT, & ! The geographic latitude at q points in degrees of latitude or m. - geoLonT, & ! The geographic longitude at q points in degrees of longitude or m. - dxT, IdxT, & ! dxT is delta x at h points, in m, and IdxT is 1/dxT in m-1. - dyT, IdyT, & ! dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1. - areaT, & ! areaT is the area of an h-cell, in m2. - IareaT, & ! IareaT = 1/areaT, in m-2. - sin_rot, & ! The sine and cosine of the angular rotation between the local - cos_rot ! model grid's northward and the true northward directions. + mask2dT, & !< 0 for land points and 1 for ocean points on the h-grid. Nd. + geoLatT, & !< The geographic latitude at q points in degrees of latitude or m. + geoLonT, & !< The geographic longitude at q points in degrees of longitude or m. + dxT, & !< dxT is delta x at h points, in m. + IdxT, & !< 1/dxT in m-1. + dyT, & !< dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1. + IdyT, & !< dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1. + areaT, & !< The area of an h-cell, in m2. + IareaT, & !< 1/areaT, in m-2. + sin_rot, & !< The sine of the angular rotation between the local model grid's northward + !! and the true northward directions. + cos_rot !< The cosine of the angular rotation between the local model grid's northward + !! and the true northward directions. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & - mask2dCu, & ! 0 for boundary points and 1 for ocean points on the u grid. Nondim. - geoLatCu, & ! The geographic latitude at u points in degrees of latitude or m. - geoLonCu, & ! The geographic longitude at u points in degrees of longitude or m. - dxCu, IdxCu, & ! dxCu is delta x at u points, in m, and IdxCu is 1/dxCu in m-1. - dyCu, IdyCu, & ! dyCu is delta y at u points, in m, and IdyCu is 1/dyCu in m-1. - dy_Cu, & ! The unblocked lengths of the u-faces of the h-cell in m. - dy_Cu_obc, & ! The unblocked lengths of the u-faces of the h-cell in m for OBC. - IareaCu, & ! The masked inverse areas of u-grid cells in m2. - areaCu ! The areas of the u-grid cells in m2. + mask2dCu, & !< 0 for boundary points and 1 for ocean points on the u grid. Nondim. + geoLatCu, & !< The geographic latitude at u points in degrees of latitude or m. + geoLonCu, & !< The geographic longitude at u points in degrees of longitude or m. + dxCu, & !< dxCu is delta x at u points, in m. + IdxCu, & !< 1/dxCu in m-1. + dyCu, & !< dyCu is delta y at u points, in m. + IdyCu, & !< 1/dyCu in m-1. + dy_Cu, & !< The unblocked lengths of the u-faces of the h-cell in m. + dy_Cu_obc, & !< The unblocked lengths of the u-faces of the h-cell in m for OBC. + IareaCu, & !< The masked inverse areas of u-grid cells in m2. + areaCu !< The areas of the u-grid cells in m2. + !> \todo dy_Cu_obc is not used? real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & - mask2dCv, & ! 0 for boundary points and 1 for ocean points on the v grid. Nondim. - geoLatCv, & ! The geographic latitude at v points in degrees of latitude or m. - geoLonCv, & ! The geographic longitude at v points in degrees of longitude or m. - dxCv, IdxCv, & ! dxCv is delta x at v points, in m, and IdxCv is 1/dxCv in m-1. - dyCv, IdyCv, & ! dyCv is delta y at v points, in m, and IdyCv is 1/dyCv in m-1. - dx_Cv, & ! The unblocked lengths of the v-faces of the h-cell in m. - dx_Cv_obc, & ! The unblocked lengths of the v-faces of the h-cell in m for OBC. - IareaCv, & ! The masked inverse areas of v-grid cells in m2. - areaCv ! The areas of the v-grid cells in m2. + mask2dCv, & !< 0 for boundary points and 1 for ocean points on the v grid. Nondim. + geoLatCv, & !< The geographic latitude at v points in degrees of latitude or m. + geoLonCv, & !< The geographic longitude at v points in degrees of longitude or m. + dxCv, & !< dxCv is delta x at v points, in m. + IdxCv, & !< 1/dxCv in m-1. + dyCv, & !< dyCv is delta y at v points, in m. + IdyCv, & !< 1/dyCv in m-1. + dx_Cv, & !< The unblocked lengths of the v-faces of the h-cell in m. + dx_Cv_obc, & !< The unblocked lengths of the v-faces of the h-cell in m for OBC. + IareaCv, & !< The masked inverse areas of v-grid cells in m2. + areaCv !< The areas of the v-grid cells in m2. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & - mask2dBu, & ! 0 for boundary points and 1 for ocean points on the q grid. Nondim. - geoLatBu, & ! The geographic latitude at q points in degrees of latitude or m. - geoLonBu, & ! The geographic longitude at q points in degrees of longitude or m. - dxBu, IdxBu, & ! dxBu is delta x at q points, in m, and IdxBu is 1/dxBu in m-1. - dyBu, IdyBu, & ! dyBu is delta y at q points, in m, and IdyBu is 1/dyBu in m-1. - areaBu, & ! areaBu is the area of a q-cell, in m2 - IareaBu ! IareaBu = 1/areaBu in m-2. + mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid. Nondim. + geoLatBu, & !< The geographic latitude at q points in degrees of latitude or m. + geoLonBu, & !< The geographic longitude at q points in degrees of longitude or m. + dxBu, & !< dxBu is delta x at q points, in m. + IdxBu, & !< 1/dxBu in m-1. + dyBu, & !< dyBu is delta y at q points, in m. + IdyBu, & !< 1/dyBu in m-1. + areaBu, & !< areaBu is the area of a q-cell, in m2 + IareaBu !< IareaBu = 1/areaBu in m-2. real, pointer, dimension(:) :: & - gridLatT => NULL(), gridLatB => NULL() ! The latitude of T or B points for - ! the purpose of labeling the output axes. - ! On many grids these are the same as geoLatT & geoLatBu. + gridLatT => NULL(), & !< The latitude of T points for the purpose of labeling the output axes. + !! On many grids this is the same as geoLatT. + gridLatB => NULL() !< The latitude of B points for the purpose of labeling the output axes. + !! On many grids this is the same as geoLatBu. real, pointer, dimension(:) :: & - gridLonT => NULL(), gridLonB => NULL() ! The longitude of T or B points for - ! the purpose of labeling the output axes. - ! On many grids these are the same as geoLonT & geoLonBu. + gridLonT => NULL(), & !< The longitude of T points for the purpose of labeling the output axes. + !! On many grids this is the same as geoLonT. + gridLonB => NULL() !< The longitude of B points for the purpose of labeling the output axes. + !! On many grids this is the same as geoLonBu. character(len=40) :: & - x_axis_units, & ! The units that are used in labeling the coordinate - y_axis_units ! axes. Except on a Cartesian grid, these are usually - ! some variant of "degrees". + x_axis_units, & !< The units that are used in labeling the x coordinate axes. + y_axis_units !< The units that are used in labeling the y coordinate axes. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - bathyT ! Ocean bottom depth at tracer points, in m. + bathyT !< Ocean bottom depth at tracer points, in m. - logical :: bathymetry_at_vel ! If true, there are separate values for the - ! basin depths at velocity points. Otherwise the effects of - ! of topography are entirely determined from thickness points. + logical :: bathymetry_at_vel !< If true, there are separate values for the + !! basin depths at velocity points. Otherwise the effects of + !! of topography are entirely determined from thickness points. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & - Dblock_u, & ! Topographic depths at u-points at which the flow is blocked - Dopen_u ! (Dblock_u) and open at width dy_Cu (Dopen_u), both in m. + Dblock_u, & !< Topographic depths at u-points at which the flow is blocked, in m. + Dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu, in m. real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & - Dblock_v, & ! Topographic depths at v-points at which the flow is blocked - Dopen_v ! (Dblock_v) and open at width dx_Cv (Dopen_v), both in m. + Dblock_v, & !< Topographic depths at v-points at which the flow is blocked, in m. + Dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv, in m. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & - CoriolisBu ! The Coriolis parameter at corner points, in s-1. + CoriolisBu !< The Coriolis parameter at corner points, in s-1. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - dF_dx, dF_dy ! Derivatives of f (Coriolis parameter) at h-points, in s-1 m-1. - real :: g_Earth ! The gravitational acceleration in m s-2. + dF_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points, in s-1 m-1. + dF_dy !< Derivative d/dy f (Coriolis parameter) at h-points, in s-1 m-1. + real :: g_Earth !< The gravitational acceleration in m s-2. ! These variables are global sums that are useful for 1-d diagnostics - real :: areaT_global ! Global sum of h-cell area in m2 - real :: IareaT_global ! Global sum of inverse h-cell area (1/areaT_global) - ! in m2 + real :: areaT_global !< Global sum of h-cell area in m2 + real :: IareaT_global !< Global sum of inverse h-cell area (1/areaT_global) in m2. + ! These variables are for block structures. - integer :: nblocks + integer :: nblocks type(hor_index_type), pointer :: Block(:) => NULL() ! store indices for each block ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) - real :: south_lat ! The latitude (or y-coordinate) of the first v-line - real :: west_lon ! The longitude (or x-coordinate) of the first u-line - real :: len_lat = 0. ! The latitudinal (or y-coord) extent of physical domain - real :: len_lon = 0. ! The longitudinal (or x-coord) extent of physical domain - real :: Rad_Earth = 6.378e6 ! The radius of the planet in meters. - real :: max_depth ! The maximum depth of the ocean in meters. + real :: south_lat !< The latitude (or y-coordinate) of the first v-line + real :: west_lon !< The longitude (or x-coordinate) of the first u-line + real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain + real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth = 6.378e6 !< The radius of the planet in meters. + real :: max_depth !< The maximum depth of the ocean in meters. end type ocean_grid_type contains -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> MOM_grid_init initializes the ocean grid array sizes and grid memory. subroutine MOM_grid_init(G, param_file, HI, global_indexing, bathymetry_at_vel) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type @@ -382,10 +403,12 @@ end function Adcroft_reciprocal !> Returns true if the coordinates (x,y) are within the h-cell (i,j) logical function isPointInCell(G, i, j, x, y) - type(ocean_grid_type), intent(in) :: G !< Grid type - integer, intent(in) :: i, j !< i,j indices of cell to test - real, intent(in) :: x, y !< x,y coordinates of point -! This is a crude calculation that assume a geographic coordinate system + type(ocean_grid_type), intent(in) :: G !< Grid type + integer, intent(in) :: i !< i index of cell to test + integer, intent(in) :: j !< j index of cell to test + real, intent(in) :: x !< x coordinate of point + real, intent(in) :: y !< y coordinate of point + ! Local variables real :: xNE, xNW, xSE, xSW, yNE, yNW, ySE, ySW real :: p0, p1, p2, p3, l0, l1, l2, l3 isPointInCell = .false. @@ -393,6 +416,7 @@ logical function isPointInCell(G, i, j, x, y) xNW = G%geoLonBu(i-1,j ) ; yNW = G%geoLatBu(i-1,j ) xSE = G%geoLonBu(i ,j-1) ; ySE = G%geoLatBu(i ,j-1) xSW = G%geoLonBu(i-1,j-1) ; ySW = G%geoLatBu(i-1,j-1) + ! This is a crude calculation that assume a geographic coordinate system if (xmax(xNE,xNW,xSE,xSW) .or. & ymax(yNE,yNW,ySE,ySW) ) then return ! Avoid the more complicated calculation @@ -419,7 +443,6 @@ subroutine set_first_direction(G, y_first) G%first_direction = y_first end subroutine set_first_direction -!--------------------------------------------------------------------- !> Allocate memory used by the ocean_grid_type and related structures. subroutine allocate_metrics(G) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type @@ -493,7 +516,6 @@ subroutine allocate_metrics(G) end subroutine allocate_metrics -!--------------------------------------------------------------------- !> Release memory used by the ocean_grid_type and related structures. subroutine MOM_grid_end(G) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type @@ -537,4 +559,26 @@ subroutine MOM_grid_end(G) end subroutine MOM_grid_end +!> \namespace mom_grid +!! +!! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid. +!! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction. +!! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between two corners of a T-cell. +!! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points. +!! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point. +!! +!! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered location on an T-cell and around a q-point. +!! +!! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively. +!! +!! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name. +!! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`. +!! +!! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in `geoLatT`, `geoLonT` for T-points. +!! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively. +!! +!! Each location also has a 2D mask indicating whether the entire column is land or ocean. +!! `mask2dT` is 1. is the column is wet or 0. is the T-cell is land. +!! `mask2dCu` is 1. if both neighboring column are ocean, and 0. if either is land. + end module MOM_grid diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9557b82509..76f6805d54 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -720,12 +720,12 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & h_new, h_old, G) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_new - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< New u values on open boundaries + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old !< Original unadjusted u + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< New v values on open boundaries + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old !< Original unadjusted v + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_new !< New h values on open boundaries + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old !< Original h values ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: grad real :: dhdt, dhdx, dhdy, gamma_u, gamma_h, gamma_v diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index b84ec40432..264500066d 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -1,46 +1,7 @@ +!> Routines for calculating baroclinic wave speeds module MOM_wave_speed -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** -! -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg and Ben Mater, May 2008 - December 2015 * -!* * -!* The subroutines in this module calculates the first baroclinic * -!* mode internal wave speed, or a set of N wave speeds. * -!* * -!* Macros written all in capital letters are defined in MOM_memory.h. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v, vh, vav * -!* j x ^ x ^ x At >: u, uh, uav * -!* j > o > o > At o: h * -!* j-1 x ^ x ^ x * -!* i-1 i i+1 At x & ^: * -!* i i+1 At > & o: * -!* * -!* The boundaries always run through q grid points (x). * -!* * -!********+*********+*********+*********+*********+*********+*********+** + +! This file is part of MOM6. See LICENSE.md for the license. 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 @@ -57,64 +18,24 @@ module MOM_wave_speed public wave_speed, wave_speeds, wave_speed_init +!> Control structure for MOM_wave_speed type, public :: wave_speed_CS ; private - type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the - ! timing of diagnostic output. + type(diag_ctrl), pointer :: diag ! Diagnostics control structure end type wave_speed_CS contains +!> Calculates the wave speed of the first baroclinic mode. subroutine wave_speed(h, tv, G, GV, cg1, CS, full_halos) - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h - type(thermo_var_ptrs), intent(in) :: tv - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 - type(wave_speed_CS), optional, pointer :: CS - logical, optional, intent(in) :: full_halos -! This subroutine determines the first mode internal wave speed. -! Arguments: h - Layer thickness, in m or kg m-2. -! (in) tv - A structure containing the thermobaric variables. -! (out) cg1 - The first mode internal gravity wave speed, in m s-1. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! wave_speed_init. -! (in,opt) full_halos - If true, do the calculation over the entire -! computational domain. - -! This subroutine solves for the first baroclinic mode wave speed. (It could -! solve for all the wave speeds, but the iterative approach taken here means -! that this is not particularly efficient.) - -! If e(k) is the perturbation interface height, this means solving for the -! smallest eigenvalue (lam = 1/c^2) of the system -! -! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0 -! -! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving -! -! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0 -! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0 -! -! Here -! Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1)) -! -! Alternately, these same eigenvalues can be found from the second smallest -! eigenvalue of the Montgomery potential (M(k)) calculation: -! -! -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0 -! -! with rigid lid and flat bottom boundary conditions -! -! (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0 -! -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0 -! -! Note that the barotropic mode has been eliminated from the rigid lid -! interface height equations, hence the matrix is one row smaller. Without -! the rigid lid, the top boundary condition is simpler to implement with -! the M equations. - + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg/m2) + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed (m/s) + type(wave_speed_CS), optional, pointer :: CS !< Control structure for MOM_wave_speed + logical, optional, intent(in) :: full_halos !< If true, do the calculation + !! over the entire computational domain. + ! Local variables real, dimension(SZK_(G)+1) :: & dRho_dT, dRho_dS, & pres, T_int, S_int, & @@ -388,57 +309,18 @@ subroutine wave_speed(h, tv, G, GV, cg1, CS, full_halos) end subroutine wave_speed +!> Calculates the wave speeds for the first few barolinic modes. subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) - type(ocean_grid_type), intent(in) :: G - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h - type(thermo_var_ptrs), intent(in) :: tv - integer, intent(in) :: nmodes - real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn - type(wave_speed_CS), optional, pointer :: CS - logical, optional, intent(in) :: full_halos -! This subroutine determines the first mode internal wave speed. -! Arguments: h - Layer thickness, in m or kg m-2. -! (in) tv - A structure containing the thermobaric variables. -! (out) cn - The internal gravity wave mode speeds, in m s-1. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! wave_speed_init. -! (in,opt) full_halos - If true, do the calculation over the entire -! computational domain. -! -! This subroutine solves for the baroclinic mode wave speed for the first -! few modes. (It is an iterative approach so is not particularly efficient.) - -! If e(k) is the perturbation interface height, this means solving for the -! eigenvalues (lam = 1/c^2) of the system -! -! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0 -! -! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving -! -! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0 -! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0 -! -! Here -! Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1)) -! -! Alternately, these same eigenvalues can be found from the second smallest -! eigenvalue of the Montgomery potential (M(k)) calculation: -! -! -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0 -! -! with rigid lid and flat bottom boundary conditions -! -! (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0 -! -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0 -! -! Note that the barotropic mode has been eliminated from the rigid lid -! interface height equations, hence the matrix is one row smaller. Without -! the rigid lid, the top boundary condition is simpler to implement with -! the M equations. - + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg/m2) + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + integer, intent(in) :: nmodes !< Number of modes + real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn !< Waves speeds (m/s) + type(wave_speed_CS), optional, pointer :: CS !< Control structure for MOM_wave_speed + logical, optional, intent(in) :: full_halos !< If true, do the calculation + !! over the entire computational domain. + ! Local variables real, dimension(SZK_(G)+1) :: & dRho_dT, dRho_dS, & pres, T_int, S_int, & @@ -969,31 +851,21 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) end subroutine wave_speeds +!> Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c where lam is constant across rows. subroutine tridiag_det(a,b,c,nrows,lam,det_out,ddet_out) - real, dimension(:), intent(in) :: a - real, dimension(:), intent(in) :: b - real, dimension(:), intent(in) :: c - integer, intent(in) :: nrows - real, intent(in) :: lam - real, intent(out):: det_out - real, intent(out):: ddet_out - ! Arguments: - ! (in) a - lower diagonal with first entry equal to zero - ! (in) b - middle diagonal - ! (in) c - upper diagonal with last entry equal to zero - ! (in) nrows - number of rows - ! (in) lam - value to subtract from middle diagonal - ! (out) det_out - determinate evaluated for given lam value - ! (out) ddet_out - derivative of determinate with respect to lam - ! evaluated for given lam value - real, dimension(nrows) :: & - det ! value of recursion function - real, dimension(nrows) :: & - ddet ! value of recursion function for derivative - real, parameter:: & - rescale = 1024.0**4 ! max value of determinant allowed before rescaling - real :: I_rescale ! inverse of rescale - integer :: n ! row (layer interface) index + real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry = 0) + real, dimension(:), intent(in) :: b !< Leading diagonal of matrix (excluding lam) + real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry = 0) + integer, intent(in) :: nrows !< Size of matrix + real, intent(in) :: lam !< Value subtracted from b + real, intent(out):: det_out !< Determinant + real, intent(out):: ddet_out !< Derivative of determinant w.r.t. lam + ! Local variables + real, dimension(nrows) :: det ! value of recursion function + real, dimension(nrows) :: ddet ! value of recursion function for derivative + real, parameter:: rescale = 1024.0**4 ! max value of determinant allowed before rescaling + real :: I_rescale ! inverse of rescale + integer :: n ! row (layer interface) index if (size(b) .ne. nrows) call MOM_error(WARNING, "Diagonal b must be same length as nrows.") if (size(a) .ne. nrows) call MOM_error(WARNING, "Diagonal a must be same length as nrows.") @@ -1017,19 +889,13 @@ subroutine tridiag_det(a,b,c,nrows,lam,det_out,ddet_out) end subroutine tridiag_det +!> Initialize control structure for MOM_wave_speed subroutine wave_speed_init(Time, G, param_file, diag, CS) - type(time_type), intent(in) :: Time - type(ocean_grid_type), intent(in) :: G - type(param_file_type), intent(in) :: param_file - type(diag_ctrl), target, intent(inout) :: diag - type(wave_speed_CS), pointer :: CS -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + type(time_type), intent(in) :: Time !< Current model time + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file handles + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(wave_speed_CS), pointer :: CS !< Control structure for MOM_wave_speed ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_wave_speed" ! This module's name. @@ -1047,4 +913,48 @@ subroutine wave_speed_init(Time, G, param_file, diag, CS) end subroutine wave_speed_init +!> \namespace mom_wave_speed +!! +!! Subroutine wave_speed() solves for the first baroclinic mode wave speed. (It could +!! solve for all the wave speeds, but the iterative approach taken here means +!! that this is not particularly efficient.) +!! +!! If `e(k)` is the perturbation interface height, this means solving for the +!! smallest eigenvalue (`lam` = 1/c^2) of the system +!! +!! \verbatim +!! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0 +!! \endverbatim +!! +!! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving +!! +!! \verbatim +!! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0 +!! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0 +!! \endverbatim +!! +!! Here +!! \verbatim +!! Igl(k) = 1.0/(gprime(k)*h(k)) ; Igu(k) = 1.0/(gprime(k)*h(k-1)) +!! \endverbatim +!! +!! Alternately, these same eigenvalues can be found from the second smallest +!! eigenvalue of the Montgomery potential (M(k)) calculation: +!! +!! \verbatim +!! -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0 +!! \endverbatim +!! +!! with rigid lid and flat bottom boundary conditions +!! +!! \verbatim +!! (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0 +!! -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0 +!! \endverbatim +!! +!! Note that the barotropic mode has been eliminated from the rigid lid +!! interface height equations, hence the matrix is one row smaller. Without +!! the rigid lid, the top boundary condition is simpler to implement with +!! the M equations. + end module MOM_wave_speed diff --git a/src/framework/MOM_hor_index.F90 b/src/framework/MOM_hor_index.F90 index db22d8fc84..4326693957 100644 --- a/src/framework/MOM_hor_index.F90 +++ b/src/framework/MOM_hor_index.F90 @@ -1,24 +1,7 @@ +!> Defines the horizontal index type (hor_index_type) used for providing index ranges module MOM_hor_index -!*********************************************************************** -!* GNU General Public License * -!* This file is a part of MOM. * -!* * -!* MOM is free software; you can redistribute it and/or modify it and * -!* are expected to follow the terms of the GNU General Public License * -!* as published by the Free Software Foundation; either version 2 of * -!* the License, or (at your option) any later version. * -!* * -!* MOM is distributed in the hope that it will be useful, but WITHOUT * -!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * -!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * -!* License for more details. * -!* * -!* For the full text of the GNU General Public License, * -!* write to: Free Software Foundation, Inc., * -!* 675 Mass Ave, Cambridge, MA 02139, USA. * -!* or see: http://www.gnu.org/licenses/gpl.html * -!*********************************************************************** +! This file is part of MOM6. See LICENSE.md for the license. use MOM_domains, only : MOM_domain_type, get_domain_extent use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL @@ -28,23 +11,48 @@ module MOM_hor_index public :: hor_index_init, assignment(=) +!> Container for horizontal index ranges for data, computational and global domains type, public :: hor_index_type - integer :: isc, iec, jsc, jec ! The range of the computational indices and - integer :: isd, ied, jsd, jed ! data indices at tracer cell centers. - integer :: isg, ieg, jsg, jeg ! The range of the global domain tracer cell indices. - integer :: IscB, IecB, JscB, JecB ! The range of the computational indices and - integer :: IsdB, IedB, JsdB, JedB ! data indices at tracer cell vertices. - integer :: IsgB, IegB, JsgB, JegB ! The range of the global domain vertex indices. - integer :: idg_offset ! The offset between the corresponding global - integer :: jdg_offset ! and local array indices. - logical :: symmetric ! True if symmetric memory is used. + integer :: isc !< The start i-index of cell centers within the computational domain + integer :: iec !< The end i-index of cell centers within the computational domain + integer :: jsc !< The start j-index of cell centers within the computational domain + integer :: jec !< The end j-index of cell centers within the computational domain + + integer :: isd !< The start i-index of cell centers within the data domain + integer :: ied !< The end i-index of cell centers within the data domain + integer :: jsd !< The start j-index of cell centers within the data domain + integer :: jed !< The end j-index of cell centers within the data domain + + integer :: isg !< The start i-index of cell centers within the global domain + integer :: ieg !< The end i-index of cell centers within the global domain + integer :: jsg !< The start j-index of cell centers within the global domain + integer :: jeg !< The end j-index of cell centers within the global domain + + integer :: IscB !< The start i-index of cell vertices within the computational domain + integer :: IecB !< The end i-index of cell vertices within the computational domain + integer :: JscB !< The start j-index of cell vertices within the computational domain + integer :: JecB !< The end j-index of cell vertices within the computational domain + + integer :: IsdB !< The start i-index of cell vertices within the data domain + integer :: IedB !< The end i-index of cell vertices within the data domain + integer :: JsdB !< The start j-index of cell vertices within the data domain + integer :: JedB !< The end j-index of cell vertices within the data domain + + integer :: IsgB !< The start i-index of cell vertices within the global domain + integer :: IegB !< The end i-index of cell vertices within the global domain + integer :: JsgB !< The start j-index of cell vertices within the global domain + integer :: JegB !< The end j-index of cell vertices within the global domain + + integer :: idg_offset !< The offset between the corresponding global and local i-indices. + integer :: jdg_offset !< The offset between the corresponding global and local j-indices. + logical :: symmetric !< True if symmetric memory is used. end type hor_index_type interface assignment(=); module procedure HIT_assign ; end interface contains -!> hor_index_init sets various index values in a hor_index_type. +!> Sets various index values in a hor_index_type. subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset) type(MOM_domain_type), intent(in) :: Domain !< The MOM domain from which to extract information. type(hor_index_type), intent(inout) :: HI !< A horizontal index type to populate with data @@ -82,10 +90,10 @@ subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset) end subroutine hor_index_init !> HIT_assign copies one hor_index_type into another. It is accessed via an -!! assignment (=) operator. +!! assignment (=) operator. subroutine HIT_assign(HI1, HI2) - type(hor_index_type), intent(out) :: HI1 - type(hor_index_type), intent(in) :: HI2 + type(hor_index_type), intent(out) :: HI1 !< Horizontal index type to copy to + type(hor_index_type), intent(in) :: HI2 !< Horizontal index type to copy from ! This subroutine copies all components of the horizontal array index type ! variable on the RHS (HI2) to the variable on the LHS (HI1). @@ -102,4 +110,20 @@ subroutine HIT_assign(HI1, HI2) end subroutine HIT_assign +!> \namespace mom_hor_index +!! +!! The hor_index_type provides the decalarations and loop ranges for almost all data with horizontal extent. +!! +!! Declarations and loop ranges should always be coded with the symmetric memory model in mind. +!! The non-symmetric memory mode will then also work, albeit with a different (less efficient) communication pattern. +!! +!! Using the hor_index_type HI: +!! - declaration of h-point data is of the form `h(HI%%isd:HI%%ied,HI%%jsd:HI%%jed)`; +!! - declaration of q-point data is of the form `q(HI%%IsdB:HI%%IedB,HI%%JsdB:HI%%JedB)`; +!! - declaration of u-point data is of the form `u(HI%%IsdB:HI%%IedB,HI%%jsd:HI%%jed)`; +!! - declaration of v-point data is of the form `v(HI%%isd:HI%%ied,HI%%JsdB:HI%%JedB)`. +!! +!! For more detail explanation of horizontal indexing see \ref Horizontal_indexing. + + end module MOM_hor_index diff --git a/src/framework/_Horizontal_indexing.dox b/src/framework/_Horizontal_indexing.dox new file mode 100644 index 0000000000..9eea1e9698 --- /dev/null +++ b/src/framework/_Horizontal_indexing.dox @@ -0,0 +1,96 @@ +/*! \page Horizontal_indexing Horizontal indexing and memory + +MOM6 is written in Fortran90 and uses the `i,j,k` order of indexing. +`i` corresponds to the fastest index (stride-1 in memory) and thus should be the inner-most loop variable. +We often refer to the i-direction as the x- or zonal direction, and similarly to the j-direction as y- or meridional direction. +The model can use curvilinear grids/coordinates in the horizontal and so these labels have loose meanings but convenient. + +\section Staggering Loops and staggered variables + +Many variables are staggered horizontally with respect to each other. +The dynamics and tracer equations are discretized on an Arakawa C grid. +Staggered variables must still have integer indices and we use a north-east convention centered on the h-points. +These means a variable with indices `i,j` will be either collocated, to the east, to the north, or to the north-east of the h-point with the same indices. + +\image html Arakawa_C_grid.png MOM6 uses an Arakawa C grid staggering of variables with a North-East indexing convention. "Cells" refer to the control volumes around tracer- or h-point located variables unless labelled otherwise. + +\subsection Soft_convention Soft convention for loop variables + +To ease reading the code we use a "soft" convection (soft because there is no syntax checking) where an upper-case index variable can be interpreted as the lower-case index variable plus \f$\frac{1}{2}\f$. + +For example, when a loop is over h-points collocated variables +- the do-loop statements will be for lower-case `i,j` variables +- references to h-point variables will be `h(i,j)`, `D(i+1,j)`, etc. +- references to u-point variables will be `u(I,j)` (meaning \f$u_{i+\frac{1}{2},j}\f$), `u(I-1,j)` (meaning \f$u_{i-\frac{1}{2},j}\f$), etc. +- references to v-point variables will be `v(i,J)` (meaning \f$v_{i,j+\frac{1}{2}}\f$), `u(I-1,j)` (meaning \f$u_{i,j-\frac{1}{2}}\f$), etc. +- references to q-point variables will be `q(I,J)` (meaning \f$q_{i+\frac{1}{2},j+\frac{1}{2}}\f$), etc. + +In contrast, when a loop is over u-points collocated variables +- the do-loop statements will be for upper-case `I` and lower-case `j` variables +- the expression \f$ u_{i+\frac{1}{2},j} ( h_{i,j} + h_{i+1,j} ) \f$ is `u(I,j) * ( h(i,j) + h(i+1,j)`. + + +\section Memory Declaration of variables + +\image html Horizontal_NE_indexing_nonsym.png Non-symmetric mode: All arrays are declared with the same shape `(isd:ied,jsd:jed)`. + +\image html Horizontal_NE_indexing_sym.png Symmetric mode: Arrays have different shapes depending on their staggering location on the Arakawa C grid. + +A field is described by 2D or 3D arrays which are distributed across parallel processors. +Each processor only sees a small window of the global field. +The processor "owns" the computational domain (red in above figure) but arrays are extended horizontally with halos which are intermittently updated with the values from neighboring processors. +The halo regions (blue in above figure) may not always be up-to-date. +Data in halo regions (blue in above figure) will be overwritten my mpp_updates. + +MOM6 has two memory models, "symmetric" and "non-symmetric". +In non-symmetric mode all arrays are given the same shape. +The consequence of this is that there are fewer staggered variables to the south-west of the computational domain. +An operator applied at h-point locations involving u- or v- point data can not have as wide a stencil on the south-west side of the processor domain as it can on the north-east side. + +In symmetric mode, declarations are dependent on the variables staggered location on the Arakawa C grid. +This allows loops to be symmetric and stencils to be applied more uniformly. + +In the code, declarations are consistent with the symmetric memory model. +The non-symmetric mode is implemented by setting the start values of the staggered data domain to be the same as the non-staggered start value. + +The horizontal index type (mom_hor_index::hor_index_type) provides the data domain start values. +The values are also copied into the mom_grid::ocean_grid_type for convenience although we might deprecate this convenience in the future. + +Declarations of h-point data take the form: +- `real, dimension(HI%%isd:HI%%ied, HI%%jsd:HI%%jed) :: D !< Depth at h-points (m)` +- `real, dimension(HI%%isd:HI%%ied, HI%%jsd:HI%%jed, GV%%ke) :: h !< Layer thickness (H units)` + +Declarations of u-point data take the form: +- `real, dimension(HI%%IsdB:HI%%IedB, HI%%jsd:HI%%jed) :: Du !< Depth at u-points (m)` +- `real, dimension(HI%%IsdB:HI%%IedB, HI%%jsd:HI%%jed, GV%%ke) :: h !< Zonal flow (m/s)` + +Declarations of v-point data take the form: +- `real, dimension(HI%%isd:HI%%ied, HI%%JsdB:HI%%JedB) :: Dv !< Depth at v-points (m)` +- `real, dimension(HI%%isd:HI%%ied, HI%%JsdB:HI%%JedB, GV%%ke) :: h !< Zonal flow (m/s)` + +Declarations of q-point data take the form: +- `real, dimension(HI%%IsdB:HI%%IedB, HI%%JsdB:HI%%JedB) :: Dq !< Depth at q-points (m)` +- `real, dimension(HI%%IsdB:HI%%IedB, HI%%JsdB:HI%%JedB, GV%%ke) :: vort !< Vertical componentof vorticity (s-1)` + +The file MOM_memory_macros.h provides the macros `SZI_`, `SZJ_`, `SZIB_` and `SZJB_` that help make the above more concise: +- `real, dimension(SZI_(HI), SZJ_(HI)) :: D !< Depth at h-points (m)` +- `real, dimension(SZIB_(HI), SZJ_(HI)) :: Du !< Depth at u-points (m)` +- `real, dimension(SZI_(HI), SZJB_(HI)) :: Dv !< Depth at v-points (m)` +- `real, dimension(SZIB_(HI), SZJB_(HI)) :: Dq !< Depth at q-points (m)` + +\section Global_index Calculating a global index + +For the most part MOM6 code should be independent of an equivalent absolute global index. +There are exceptions and when the global index of a cell `i,j` is needed is can be calculated as follows: + + `i_global = i + HI%%idg_offset` + +Before the mom_hor_index::hor_index_type was introduced, this conversion was done use variables in the mom_grid::ocean_grid_type: + + `i_global = (i-G%%isd) + G%%isd_global` + +which is no longer preferred. + +Note that a global index only makes sense for a rectangular global domain. If the domain is a Mosaic of connected tiles (e.g. size tiles of a cube) the global indices (i,j) become meaningless. + +*/ diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index b4685674af..135952c86b 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -34,7 +34,7 @@ module MOM_diapyc_energy_req use MOM_grid, only : ocean_grid_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_specific_vol_derivs +use MOM_EOS, only : calculate_specific_vol_derivs, calculate_density implicit none ; private @@ -50,9 +50,9 @@ module MOM_diapyc_energy_req ! profile in place of any that might be passed ! in as an argument. type(diag_ctrl), pointer :: diag ! Structure used to regulate timing of diagnostic output - integer :: id_ERt=-1, id_ERb=-1, id_ERc=-1, id_ERh=-1, id_Kddt=-1 + integer :: id_ERt=-1, id_ERb=-1, id_ERc=-1, id_ERh=-1, id_Kddt=-1, id_Kd=-1 integer :: id_CHCt=-1, id_CHCb=-1, id_CHCc=-1, id_CHCh=-1 - integer :: id_T0=-1, id_Tf=-1, id_S0=-1, id_Sf=-1 + integer :: id_T0=-1, id_Tf=-1, id_S0=-1, id_Sf=-1, id_N2_0=-1, id_N2_f=-1 integer :: id_h=-1, id_zInt=-1 end type diapyc_energy_req_CS @@ -211,6 +211,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & real, dimension(GV%ke+1) :: & pres, & ! Interface pressures in Pa. z_Int, & ! Interface heights relative to the surface, in m. + N2, & ! An estimate of the buoyancy frequency in s-2. Kddt_h_a, & ! Kddt_h_b, & ! Kddt_h, & ! The diapycnal diffusivity times a timestep divided by the @@ -238,6 +239,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & ! converted to turbulent kinetic energy if the velocity in ! the layer below an interface were homogenized with all of ! the water above the interface, in J m-2 = kg s-2. + real :: rho_here ! The in-situ density, in kg m-3. real :: PE_change, ColHt_cor real :: htot real :: dT_k, dT_km1, dS_k, dS_km1 ! Temporary arrays @@ -930,6 +932,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & if (CS%id_ERc>0) call post_data_1d_k(CS%id_ERc, PE_chg_k(:,3), CS%diag) if (CS%id_ERh>0) call post_data_1d_k(CS%id_ERh, PE_chg_k(:,4), CS%diag) if (CS%id_Kddt>0) call post_data_1d_k(CS%id_Kddt, GV%H_to_m*Kddt_h, CS%diag) + if (CS%id_Kd>0) call post_data_1d_k(CS%id_Kd, Kd, CS%diag) if (CS%id_h>0) call post_data_1d_k(CS%id_h, GV%H_to_m*h_tr, CS%diag) if (CS%id_zInt>0) call post_data_1d_k(CS%id_zInt, Z_int, CS%diag) if (CS%id_CHCt>0) call post_data_1d_k(CS%id_CHCt, ColHt_cor_k(:,1), CS%diag) @@ -940,6 +943,28 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & if (CS%id_Tf>0) call post_data_1d_k(CS%id_Tf, Tf, CS%diag) if (CS%id_S0>0) call post_data_1d_k(CS%id_S0, S0, CS%diag) if (CS%id_Sf>0) call post_data_1d_k(CS%id_Sf, Sf, CS%diag) + if (CS%id_N2_0>0) then + N2(1) = 0.0 ; N2(nz+1) = 0.0 + do K=2,nz + call calculate_density(0.5*(T0(k-1) + T0(k)), 0.5*(S0(k-1) + S0(k)), & + pres(K), rho_here, tv%eqn_of_state) + N2(K) = (GV%g_Earth * rho_here / (0.5*GV%H_to_m*(h_tr(k-1) + h_tr(k)))) * & + ( 0.5*(dSV_dT(k-1) + dSV_dT(k)) * (T0(k-1) - T0(k)) + & + 0.5*(dSV_dS(k-1) + dSV_dS(k)) * (S0(k-1) - S0(k)) ) + enddo + call post_data_1d_k(CS%id_N2_0, N2, CS%diag) + endif + if (CS%id_N2_f>0) then + N2(1) = 0.0 ; N2(nz+1) = 0.0 + do K=2,nz + call calculate_density(0.5*(Tf(k-1) + Tf(k)), 0.5*(Sf(k-1) + Sf(k)), & + pres(K), rho_here, tv%eqn_of_state) + N2(K) = (GV%g_Earth * rho_here / (0.5*GV%H_to_m*(h_tr(k-1) + h_tr(k)))) * & + ( 0.5*(dSV_dT(k-1) + dSV_dT(k)) * (Tf(k-1) - Tf(k)) + & + 0.5*(dSV_dS(k-1) + dSV_dS(k)) * (Sf(k-1) - Sf(k)) ) + enddo + call post_data_1d_k(CS%id_N2_f, N2, CS%diag) + endif endif end subroutine diapyc_energy_req_calc @@ -1038,7 +1063,7 @@ subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & ! The expression for the change in potential energy used here is derived ! from the expression for the final estimates of the changes in temperature ! and salinities, and then extensively manipulated to get it into its most - ! succint form. The derivation is not necessarily obvious, but it demonstably + ! succint form. The derivation is not necessarily obvious, but it demonstrably ! works by comparison with separate calculations of the energy changes after ! the tridiagonal solver for the final changes in temperature and salinity are ! applied. @@ -1092,11 +1117,14 @@ subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & end subroutine find_PE_chg +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interfaces's diapycnal diffusivity times a timestep +!! using the original form used in the first version of ePBL. subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & - dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, & - dT_to_dPEa, dS_to_dPEa, pres, dT_to_dColHt_k, & - dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, & - PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0) + dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, & + dT_to_dPEa, dS_to_dPEa, pres, dT_to_dColHt_k, & + dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, & + PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0) real, intent(in) :: Kddt_h !< The diffusivity at an interface times the time step and !! divided by the average of the thicknesses around the !! interface, in units of H (m or kg-2). @@ -1167,64 +1195,7 @@ subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & ! The comments describing these arguments are for a downward mixing pass, but ! this routine can also be used for an upward pass with the sense of direction ! reversed. -! -! Arguments: Kddt_h - The diffusivity at an interface times the time step and -! divided by the average of the thicknesses around the -! interface, in units of H (m or kg-2). -! (in) h_k - The thickness of the layer below the interface, in H. -! (in) b_den_1 - The first first term in the denominator of the pivot -! for the tridiagonal solver, given by h_k plus a term that -! is a fraction (determined from the tridiagonal solver) of -! Kddt_h for the interface above, in H. -! (in) dTe_term - A diffusivity-independent term related to the -! temperature change in the layer below the interface, in K H. -! (in) dSe_term - A diffusivity-independent term related to the -! salinity change in the layer below the interface, in ppt H. -! (in) dT_km1_t2 - A diffusivity-independent term related to the -! temperature change in the layer above the interface, in K. -! (in) dS_km1_t2 - A diffusivity-independent term related to the -! salinity change in the layer above the interface, in ppt. -! (in) dT_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating -! a layer's temperature change to the change in column -! potential energy, in J m-2 K-1. -! (in) dS_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating -! a layer's salinity change to the change in column -! potential energy, in J m-2 ppt-1. -! (in) dT_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating -! a layer's temperature change to the change in column -! potential energy, including all implicit diffusive changes -! in the temperatures of all the layers above, in J m-2 K-1. -! (in) dS_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating -! a layer's salinity change to the change in column -! potential energy, including all implicit diffusive changes -! in the salinities of all the layers above, in J m-2 ppt-1. -! (in) pres - The hydrostatic interface pressure, which is used to relate -! the changes in column thickness to the energy that is radiated -! as gravity waves and unavailable to drive mixing, in Pa. -! (in) dT_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dT) relating -! a layer's temperature change to the change in column -! height, in m K-1. -! (in) dS_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dS) relating -! a layer's salinity change to the change in column -! height, in m ppt-1. -! (in) dT_to_dColHta - A factor (mass_lay*dSColHtc_vol/dT) relating -! a layer's temperature change to the change in column -! height, including all implicit diffusive changes -! in the temperatures of all the layers above, in m K-1. -! (in) dS_to_dColHta - A factor (mass_lay*dSColHtc_vol/dS) relating -! a layer's salinity change to the change in column -! height, including all implicit diffusive changes -! in the salinities of all the layers above, in m ppt-1. -! (out,opt) PE_chg - The change in column potential energy from applying -! Kddt_h at the present interface, in J m-2. -! (out,opt) dPE_chg_dKd - The partial derivative of PE_chg with Kddt_h, -! in units of J m-2 H-1. -! (out,opt) PE_max - The maximum change in column potential energy that could -! be realizedd by applying a huge value of Kddt_h at the -! present interface, in J m-2. -! (out,opt) dPEc_dKc0 - The partial derivative of PE_chg with Kddt_h in the -! limit where Kddt_h = 0, in J m-2 H-1. - ! b_den_1 - The first term in the denominator of b1, in m or kg m-2. + real :: b1 ! b1 is used by the tridiagonal solver, in H-1. real :: b1Kd ! Temporary array (nondim.) real :: ColHt_chg ! The change in column thickness in m. @@ -1345,6 +1316,8 @@ subroutine diapyc_energy_req_init(Time, G, param_file, diag, CS) "Diffusivity Energy Requirements, halves", "J m-2") CS%id_Kddt = register_diag_field('ocean_model', 'EnReqTest_Kddt', diag%axesZi, Time, & "Implicit diffusive coupling coefficient", "m") + CS%id_Kd = register_diag_field('ocean_model', 'EnReqTest_Kd', diag%axesZi, Time, & + "Diffusivity in test", "m2 s-1") CS%id_h = register_diag_field('ocean_model', 'EnReqTest_h_lay', diag%axesZL, Time, & "Test column layer thicknesses", "m") CS%id_zInt = register_diag_field('ocean_model', 'EnReqTest_z_int', diag%axesZi, Time, & @@ -1365,6 +1338,10 @@ subroutine diapyc_energy_req_init(Time, G, param_file, diag, CS) "Salinity before mixing", "g kg-1") CS%id_Sf = register_diag_field('ocean_model', 'EnReqTest_Sf', diag%axesZL, Time, & "Salinity after mixing", "g kg-1") + CS%id_N2_0 = register_diag_field('ocean_model', 'EnReqTest_N2_0', diag%axesZi, Time, & + "Squared buoyancy frequency before mixing", "second-2") + CS%id_N2_f = register_diag_field('ocean_model', 'EnReqTest_N2_f', diag%axesZi, Time, & + "Squared buoyancy frequency after mixing", "second-2") end subroutine diapyc_energy_req_init diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index be1da1087b..4a374d5262 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -114,9 +114,19 @@ module MOM_energetic_PBL ! the inhibition of the diffusive length scale by ! rotation. Making this larger decreases the ! diffusivity in the planetary boundary layer. + real :: transLay_scale ! A scale for the mixing length in the transition layer + ! at the edge of the boundary layer as a fraction of the + ! boundary layer thickness. The default is 0, but a + ! value of 0.1 might be better justified by observations. + real :: MLD_tol ! A tolerance for determining the boundary layer + ! thickness when Use_MLD_iteration is true, in m. + real :: min_mix_len ! The minimum mixing length scale that will be + ! used by ePBL, in m. The default (0) does not + ! set a minimum. type(time_type), pointer :: Time ! A pointer to the ocean model's clock. logical :: TKE_diagnostics = .false. - LOGICAL :: Use_MLD_ITERATION=.false.!False to use old ePBL method. + logical :: orig_PE_calc = .true. + logical :: Use_MLD_iteration=.false. ! False to use old ePBL method. logical :: Mixing_Diagnostics = .false. ! Will be true when outputing mixing ! length and velocity scale type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the @@ -125,7 +135,7 @@ module MOM_energetic_PBL ! These are terms in the mixed layer TKE budget, all in J m-2 = kg s-2. real, allocatable, dimension(:,:) :: & ML_depth, & ! The mixed layer depth in m. - ML_depth2, & ! The mixed layer depth in m. + ML_depth2, & ! The mixed layer depth in m. diag_TKE_wind, & ! The wind source of TKE. diag_TKE_MKE, & ! The resolved KE source of TKE. diag_TKE_conv, & ! The convective source of TKE. @@ -224,20 +234,21 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & real, dimension(SZI_(G),SZK_(GV)+1) :: & Kd, & ! The diapycnal diffusivity, in m2 s-1. pres, & ! Interface pressures in Pa. - hb_hs ! The distance from the bottom over the thickness of the water, - ! times a conversion factor from H to m, in m H-1. + hb_hs ! The distance from the bottom over the thickness of the + ! water column, nondim. real, dimension(SZI_(G)) :: & mech_TKE, & ! The mechanically generated turbulent kinetic energy ! available for mixing over a time step, in J m-2 = kg s-2. - conv_PErel, & ! The potential energy that has been convectively released + conv_PErel, & ! The potential energy that has been convectively released ! during this timestep, in J m-2 = kg s-2. A portion nstar_FC ! of conv_PErel is available to drive mixing. htot, & ! The total depth of the layers above an interface, in H. uhtot, & ! The depth integrated zonal and meridional velocities in the vhtot, & ! layers above, in H m s-1. + mech_TKE_top, & ! The value of mech_TKE at the top of the column, in J m-2. + conv_PErel_top, & ! The value of conv_PErel at the top of the column, in J m-2. Idecay_len_TKE, & ! The inverse of a turbulence decay length scale, in H-1. - h_bot, & ! The distance from the bottom, in H. h_sum, & ! The total thickness of the water column, in H. absf ! The absolute value of f, in s-1. @@ -259,9 +270,23 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & Te, Se, & ! Estimated final values of T and S in the column, in K and ppt. c1, & ! c1 is used by the tridiagonal solver, ND. dTe, dSe ! Running (1-way) estimates of temperature and salinity change. + real, dimension(SZK_(GV)) :: & + Th_a, & ! An effective temperature times a thickness in the layer above, + ! including implicit mixing effects with other yet higher layers, in K H. + Sh_a, & ! An effective salinity times a thickness in the layer above, + ! including implicit mixing effects with other yet higher layers, in K H. + Th_b, & ! An effective temperature times a thickness in the layer below, + ! including implicit mixing effects with other yet lower layers, in K H. + Sh_b ! An effective salinity times a thickness in the layer below, + ! including implicit mixing effects with other yet lower layers, in K H. real, dimension(SZI_(G)) :: & - b_den_1 ! The first term in the denominator of b1, in H. + hp_a ! An effective pivot thickness of the layer including the effects + ! of coupling with layers above, in H. This is the first term + ! in the denominator of b1 in a downward-oriented tridiagonal solver. real, dimension(SZK_(GV)+1) :: & + MixLen_shape, & ! A nondimensional shape factor for the mixing length that + ! gives it an appropriate assymptotic value at the bottom of + ! the boundary layer. Kddt_h ! The diapycnal diffusivity times a timestep divided by the ! average thicknesses around a layer, in H (m or kg m-2). real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver, in H-1. @@ -279,7 +304,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & real :: dt_h ! The timestep divided by the averages of the thicknesses around ! a layer, times a thickness conversion factor, in H s m-2. - real :: I_hs ! The inverse of h_sum. + real :: h_bot ! The distance from the bottom, in H. + real :: h_rsum ! The running sum of h from the top, in H. + real :: I_hs ! The inverse of h_sum, in H-1. + real :: I_mld ! The inverse of the current value of MLD, in H-1. real :: h_tt ! The distance from the surface or up to the next interface ! that did not exhibit turbulent mixing from this scheme plus ! a surface mixing roughness length given by h_tt_min, in H. @@ -291,6 +319,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! used convert TKE back into ustar^3. real :: U_star ! The surface friction velocity, in m s-1. real :: vstar ! An in-situ turbulent velocity, in m s-1. + real :: hbs_here ! The local minimum of hb_hs and MixLen_shape, times a + ! conversion factor from H to M, in m H-1. real :: nstar_FC ! The fraction of conv_PErel that can be converted to mixing, nondim. real :: TKE_reduc ! The fraction by which TKE and other energy fields are ! reduced to support mixing, nondim. between 0 and 1. @@ -342,11 +372,14 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & logical :: write_diags ! If true, write out diagnostics with this step. logical :: reset_diags ! If true, zero out the accumulated diagnostics. ! detrainment, in units of m. + ! Local column copies of energy change diagnostics, all in J m-2. + real :: dTKE_conv, dTKE_forcing, dTKE_mixing + real :: dTKE_MKE, dTKE_mech_decay, dTKE_conv_decay !---------------------------------------------------------------------- !/BGR added Aug24,2016 for adding iteration to get boundary layer depth ! - needed to compute new mixing length. - real :: MLD_GUESS, MLD_FOUND ! Mixing Layer depth guessed/found for iteration - real :: MAX_MLD, MIN_MLD ! Iteration bounds which are adjusted at each step + real :: MLD_guess, MLD_found ! Mixing Layer depth guessed/found for iteration, in m. + real :: max_MLD, min_MLD ! Iteration bounds, in m, which are adjusted at each step ! - These are initialized based on surface/bottom ! 1. The iteration guesses a value (possibly from ! prev step or neighbor). @@ -358,44 +391,43 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! is very quick (e.g. if MLD doesn't change ! over timestep). Otherwise it takes 5-10 ! passes, but has a high convergence rate. - ! Other iteration may be tried, but this + ! Other iteration may be tried, but this ! method seems to rarely fail and the added ! cost is likely not significant. Additionally, - ! when it fails it does so in a reasonable - ! manner giving a usable guess. When it + ! when it fails it does so in a reasonable + ! manner giving a usable guess. When it ! does fail, it is due to convection within ! the boundary. Likely, a new method e.g. ! surface_disconnect, can improve this. logical :: FIRST_OBL ! Flag for computing "found" Mixing layer depth logical :: OBL_CONVERGED ! Flag for convergence of MLD - INTEGER :: OBL_IT ! Iteration counter - REAL :: MinMixLen=1.0 ! Minimum value for mixing length (must be non-zero - ! for iteration to converge when OBL shoals). - ! At present, this is the value used for convection - ! in the ocean interior. - INTEGER :: MAX_OBL_IT=20 ! Set maximum number of iterations. Probably + integer :: OBL_IT ! Iteration counter +!### These need to be made into run-time parameters. + integer :: MAX_OBL_IT=20 ! Set maximum number of iterations. Probably ! best as an input parameter, but then may want ! to use allocatable arrays if storing ! guess/found (as diagnostic); skipping for now. - ! In reality, the maximum number of guesses - ! needed is set by: + ! In reality, the maximum number of guesses + ! needed is set by: ! DEPTH/2^M < DZ ! where M is the number of guesses ! e.g. M=12 for DEPTH=4000m and DZ=1m real, dimension(SZK_(GV)+1) :: Vstar_Used, & ! 1D arrays used to store Mixing_Length_Used ! Vstar and Mixing_Length + !/BGR - remaining variables are related to tracking iteration statistics. logical :: OBL_IT_STATS=.false. ! Flag for computing OBL iteration statistics REAL :: ITguess(20), ITresult(20),ITmax(20),ITmin(20) ! Flag for storing guess/result ! should have dim=MAX_OBL_IT integer, save :: MAXIT=0 ! Stores maximum number of iterations - integer, save :: MINIT=1e8 ! Stokes minimum number of iterations + integer, save :: MINIT=1e8 ! Stores minimum number of iterations integer, save :: SUMIT=0 ! Stores total iterations (summed over all) integer, save :: NUMIT=0 ! Stores number of times iterated !e.g. Average iterations = SUMIT/NUMIT integer, save :: CONVERGED! integer, save :: NOTCONVERGED! !-End BGR iteration parameters----------------------------------------- + logical :: debug=.false. ! Change this hard-coded value for debugging. ! The following arrays are used only for debugging purposes. real :: dPE_debug, mixing_debug @@ -421,7 +453,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & h_neglect = GV%H_subroundoff - if(.not.CS%Use_MLD_Iteration) MAX_OBL_IT=1 + if(.not.CS%Use_MLD_Iteration) MAX_OBL_IT=1 C1_3 = 1.0 / 3.0 dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag IdtdR0 = 1.0 / (dt__diag * GV%Rho0) @@ -448,9 +480,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & CS%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced_forcing(i,j) = 0.0 enddo ; enddo endif - IF (CS%Mixing_Diagnostics) then - CS%Mixing_Length(:,:,:)=0.0 - CS%Velocity_Scale(:,:,:)=0.0 + if (CS%Mixing_Diagnostics) then + CS%Mixing_Length(:,:,:) = 0.0 + CS%Velocity_Scale(:,:,:) = 0.0 endif !!OMP end parallel endif @@ -465,17 +497,19 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !!OMP U_Star,absf,mech_TKE,conv_PErel,nstar_k, & !!OMP h_sum,I_hs,h_bot,hb_hs,T0,S0,num_itts, & !!OMP pres,dMass,dPres,dT_to_dPE,dS_to_dPE, & -!!OMP dT_to_dColHt,dS_to_dColHt,Kddt_h, & -!!OMP b_den_1,dT_to_dPE_a,dT_to_dColHt_a, & -!!OMP dS_to_dColHt_a,htot,uhtot,vhtot, & +!!OMP dT_to_dColHt,dS_to_dColHt,Kddt_h,hp_a, & +!!OMP Th_a,Sh_a,Th_b,Sh_b,dT_to_dPE_a,htot, & +!!OMP dT_to_dColHt_a,dS_to_dColHt_a,uhtot,vhtot, & !!OMP Idecay_len_TKE,exp_kh,nstar_FC,tot_TKE, & !!OMP TKE_reduc,dTe_t2,dSe_t2,dTe,dSe,dt_h, & !!OMP Convectively_stable,sfc_disconnect,b1, & !!OMP c1,dT_km1_t2,dS_km1_t2,dTe_term, & -!!OMP dSe_term,MKE2_Hharm,vstar,h_tt, & +!!OMP dSe_term,MKE2_Hharm,vstar,h_tt,h_rsum, & !!OMP Kd_guess0,Kddt_h_g0,dPEc_dKd_Kd0, & !!OMP PE_chg_max,dPEa_dKd_g0,PE_chg_g0, & !!OMP MKE_src,dPE_conv,Kddt_h_max,Kddt_h_min, & +!!OMP dTKE_conv, dTKE_forcing, dTKE_mixing, & +!!OMP dTKE_MKE,dTKE_mech_decay,dTKE_conv_decay,& !!OMP TKE_left_max,TKE_left_min,Kddt_h_guess, & !!OMP TKE_left_itt,dPEa_dKd_itt,PE_chg_itt, & !!OMP MKE_src_itt,Kddt_h_itt,dPEc_dKd,PE_chg, & @@ -490,10 +524,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & T(i,k) = tv%T(i,j,k) ; S(i,k) = tv%S(i,j,k) Kd(i,K) = 0.0 enddo ; enddo - do i=is,ie - CS%ML_depth(i,j) = h(i,1)*GV%H_to_m ; - !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m ; - sfc_connected(i) = .true. ; + do i=is,ie + CS%ML_depth(i,j) = h(i,1)*GV%H_to_m + !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m + sfc_connected(i) = .true. enddo if (debug) then @@ -507,39 +541,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! and ustar and wstar available to drive mixing at the first interior ! interface. do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then - !/The following lines are for the iteration over MLD - !{ - OBL_CONVERGED=.false.!Initialize convergence state to 'false' - MAX_MLD = 0.0;!MAX_MLD will initialized as ocean bottom depth - do k=1,nz ; - MAX_MLD = MAX_MLD + h(i,k)*GV%H_to_m; - enddo - MIN_MLD = 0.0 !MIN_MLD will initialize as 0. - !/BGR: Add MLD_GUESS based on stored previous value. - ! note that this is different from ML_Depth already - ! computed by EPBL, need to figure out why. - if (CS%ML_Depth2(i,j).gt.1.) then - !If prev value is present use for guess. - MLD_GUESS=CS%ML_Depth2(i,j) - else - !Otherwise guess middle of water column. - MLD_GUESS = 0.5 * (MIN_MLD+MAX_MLD) - endif - !/BGR: May add user-input bounds for max/min MLD - DO OBL_IT=1,MAX_OBL_IT ! Iterate up to MAX_OBL_IT times - IF (.not. OBL_CONVERGED) THEN !Cycle through EPBL if not converged - CS%ML_depth(i,j) = h(i,1)*GV%H_to_m ; - !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m ; - sfc_connected(i) = .true. ; - ! Store in 1D arrays cleared out each iteration. Only write in - ! 3D arrays after convergence. - VSTAR_USED(:)=0.0 - MIXING_LENGTH_USED(:)=0.0 - IF (.not.CS%Use_MLD_Iteration) OBL_CONVERGED=.true. - !/This ends the new code for the iteration. - !} U_Star = fluxes%ustar(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & @@ -581,40 +584,15 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! endif ; enddo ! do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then - if (debug) then - mech_TKE_k(i,1) = mech_TKE(i) ; conv_PErel_k(i,1) = conv_PErel(i) - nstar_k(:) = 0.0 ; nstar_k(1) = CS%nstar - endif - - if (.not.CS%Use_MLD_Iteration) then - h_sum(i) = H_neglect - do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo - I_hs = 0.0 ; if (h_sum(i) > 0.0) I_hs = 1.0 / h_sum(i) - h_bot(i) = 0.0 ; hb_hs(i,nz+1) = 0.0 - do k=nz,1,-1 - h_bot(i) = h_bot(i) + h(i,k) - hb_hs(i,K) = GV%H_to_m * (h_bot(i)*I_hs) - enddo - else - !New method where mixing length is reduced based on MLD - I_hs=1.0/MLD_GUESS - h_sum(i) = 0.0 - h_bot(i) = 0.0 ; hb_hs(i,1:nz+1) = 0.0 - do k=2,nz - h_sum(i)=h_sum(i)+h(i,k) - h_bot(i) = max(0.0,MLD_GUESS - h_sum(i)) - hb_hs(i,K) = GV%H_to_m * (h_bot(i)*I_hs)**2!Notice the square - ! makes KPP-like. - enddo - endif -! endif ; enddo + h_sum(i) = H_neglect ; do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo + I_hs = 0.0 ; if (h_sum(i) > 0.0) I_hs = 1.0 / h_sum(i) - ! Note the outer i-loop and inner k-loop loop order!!! -! do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then - do k=1,nz ; T0(k) = T(i,k) ; S0(k) = S(i,k) ; enddo - - num_itts(:) = -1 + h_bot = 0.0 ; hb_hs(i,nz+1) = 0.0 + do k=nz,1,-1 + h_bot = h_bot + h(i,k) + hb_hs(i,K) = h_bot * I_hs + enddo pres(i,1) = 0.0 do k=1,nz @@ -628,517 +606,645 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & pres(i,K+1) = pres(i,K) + dPres enddo - Kd(i,1) = 0.0 ; Kddt_h(1) = 0.0 - b_den_1(i) = h(i,1) - dT_to_dPE_a(i,1) = dT_to_dPE(i,1) - dS_to_dPE_a(i,1) = dS_to_dPE(i,1) - dT_to_dColHt_a(i,1) = dT_to_dColHt(i,1) - dS_to_dColHt_a(i,1) = dS_to_dColHt(i,1) - - htot(i) = h(i,1) - uhtot(i) = u(i,1)*h(i,1) - vhtot(i) = v(i,1)*h(i,1) - - do K=2,nz - ! Apply dissipation to the TKE, here applied as an exponential decay - ! due to 3-d turbulent energy being lost to inefficient rotational modes. - - ! There should be several different "flavors" of TKE that decay at - ! different rates. The following form is often used for mechanical - ! stirring from the surface, perhaps due to breaking surface gravity - ! waves and wind-driven turbulence. - Idecay_len_TKE(i) = (CS%TKE_decay * absf(i) / U_Star) * GV%H_to_m - exp_kh = 1.0 - if (Idecay_len_TKE(i) > 0.0) exp_kh = exp(-h(i,k-1)*Idecay_len_TKE(i)) - if (CS%TKE_diagnostics) CS%diag_TKE_mech_decay(i,j) = & - CS%diag_TKE_mech_decay(i,j) + (exp_kh-1.0) * mech_TKE(i) * IdtdR0 - mech_TKE(i) = mech_TKE(i) * exp_kh - - - ! Accumulate any convectively released potential energy to contribute - ! to wstar and to drive penetrating convection. - if (TKE_forced(i,j,k) > 0.0) then - conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,k) - if (CS%TKE_diagnostics) then - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + CS%nstar*TKE_forced(i,j,k) * IdtdR0 - endif - endif +! endif ; enddo - if (debug) then - mech_TKE_k(i,K) = mech_TKE(i) ; conv_PErel_k(i,K) = conv_PErel(i) - endif + ! Note the outer i-loop and inner k-loop loop order!!! +! do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + do k=1,nz ; T0(k) = T(i,k) ; S0(k) = S(i,k) ; enddo - ! Determine the total energy - nstar_FC = CS%nstar - if (CS%nstar * conv_PErel(i) > 0.0) then - ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based - ! on a curve fit from the data of Wang (GRL, 2003). - ! Note: Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot(i))**3 / conv_PErel(i)) - nstar_FC = CS%nstar * conv_PErel(i) / (conv_PErel(i) + 0.2 * & - sqrt(0.5 * dt * GV%Rho0 * (absf(i)*(htot(i)*GV%H_to_m))**3 * conv_PErel(i))) + ! Store the initial mechanical TKE and convectively released PE to + ! enable multiple iterations. + mech_TKE_top(i) = mech_TKE(i) ; conv_PErel_top(i) = conv_PErel(i) + + !/The following lines are for the iteration over MLD + !{ + ! max_MLD will initialized as ocean bottom depth + max_MLD = 0.0 ; do k=1,nz ; max_MLD = max_MLD + h(i,k)*GV%H_to_m ; enddo + min_MLD = 0.0 !min_MLD will initialize as 0. + !/BGR: May add user-input bounds for max/min MLD + + !/BGR: Add MLD_guess based on stored previous value. + ! note that this is different from ML_Depth already + ! computed by EPBL, need to figure out why. + if (CS%ML_Depth2(i,j) > 1.) then + !If prev value is present use for guess. + MLD_guess=CS%ML_Depth2(i,j) + else + !Otherwise guess middle of water column. + MLD_guess = 0.5 * (min_MLD+max_MLD) + endif + + ! Iterate up to MAX_OBL_IT times to determine a converged EPBL depth. + OBL_CONVERGED = .false. + do OBL_IT=1,MAX_OBL_IT ; if (.not. OBL_CONVERGED) then + CS%ML_depth(i,j) = h(i,1)*GV%H_to_m + !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m + + sfc_connected(i) = .true. + mech_TKE(i) = mech_TKE_top(i) ; conv_PErel(i) = conv_PErel_top(i) + + + if (CS%TKE_diagnostics) then + dTKE_conv = 0.0 ; dTKE_forcing = 0.0 ; dTKE_mixing = 0.0 + dTKE_MKE = 0.0 ; dTKE_mech_decay = 0.0 ; dTKE_conv_decay = 0.0 endif - if (debug) nstar_k(K) = nstar_FC - - tot_TKE = mech_TKE(i) + nstar_FC * conv_PErel(i) - - ! For each interior interface, first discard the TKE to account for - ! mixing of shortwave radiation through the next denser cell. - if (TKE_forced(i,j,k) < 0.0) then - if (TKE_forced(i,j,k) + tot_TKE < 0.0) then - ! The shortwave requirements deplete all the energy in this layer. - if (CS%TKE_diagnostics) then - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) + tot_TKE * IdtdR0 - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) - tot_TKE * IdtdR0 - ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + & - ! (TKE_forced(i,j,k) + tot_TKE) * IdtdR0 - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + & - (CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 - endif - tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 + + ! Store in 1D arrays cleared out each iteration. Only write in + ! 3D arrays after convergence. + do k=1,nz + Vstar_Used(k) = 0.0 ; Mixing_Length_Used(k) = 0.0 + enddo + if (.not.CS%Use_MLD_Iteration) OBL_CONVERGED = .true. + + if ((.not.CS%Use_MLD_Iteration) .or. & + (CS%transLay_scale >= 1.0) .or. (CS%transLay_scale < 0.0) ) then + do K=1,nz+1 ; MixLen_shape(K) = 1.0 ; enddo + elseif (MLD_guess <= 0.0) then + if (CS%transLay_scale > 0.0) then + do K=1,nz+1 ; MixLen_shape(K) = CS%transLay_scale ; enddo else - ! Reduce the mechanical and convective TKE proportionately. - TKE_reduc = (tot_TKE + TKE_forced(i,j,k)) / tot_TKE - if (CS%TKE_diagnostics) then - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) - TKE_forced(i,j,k) * IdtdR0 - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + TKE_forced(i,j,k) * IdtdR0 - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + & - (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 - endif - tot_TKE = TKE_reduc*tot_TKE ! = tot_TKE + TKE_forced(i,j,k) - mech_TKE(i) = TKE_reduc*mech_TKE(i) - conv_PErel(i) = TKE_reduc*conv_PErel(i) + do K=1,nz+1 ; MixLen_shape(K) = 1.0 ; enddo endif + else + ! Reduce the mixing length based on MLD, with a quadratic + ! expression that follows KPP. + I_MLD = 1.0 / MLD_guess ; h_rsum = 0.0 + MixLen_shape(1) = 1.0 + do K=2,nz+1 + h_rsum = h_rsum + h(i,k-1) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**2 + enddo endif - ! Precalculate some temporary expressions that are independent of Kddt_h(K). - if (K==2) then - dTe_t2 = 0.0 ; dSe_t2 = 0.0 - else - dTe_t2 = Kddt_h(K-1) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) - dSe_t2 = Kddt_h(K-1) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + + Kd(i,1) = 0.0 ; Kddt_h(1) = 0.0 + hp_a(i) = h(i,1) + dT_to_dPE_a(i,1) = dT_to_dPE(i,1) ; dT_to_dColHt_a(i,1) = dT_to_dColHt(i,1) + dS_to_dPE_a(i,1) = dS_to_dPE(i,1) ; dS_to_dColHt_a(i,1) = dS_to_dColHt(i,1) + + htot(i) = h(i,1) ; uhtot(i) = u(i,1)*h(i,1) ; vhtot(i) = v(i,1)*h(i,1) + + if (debug) then + mech_TKE_k(i,1) = mech_TKE(i) ; conv_PErel_k(i,1) = conv_PErel(i) + nstar_k(:) = 0.0 ; nstar_k(1) = CS%nstar ; num_itts(:) = -1 endif - dt_h = (GV%m_to_H**2*dt) / max(0.5*(h(i,k-1)+h(i,k)), 1e-15*h_sum(i)) - - ! This tests whether the layers above and below this interface are in - ! a convetively stable configuration, without considering any effects of - ! mixing at higher interfaces. It is an approximation to the more - ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of - ! mixing across interface K-1. The dT_to_dColHt here are effectively - ! mass-weigted estimates of dSV_dT. - Convectively_stable = ( 0.0 <= & - ( (dT_to_dColHt(i,k) + dT_to_dColHt(i,k-1) ) * (T0(k-1)-T0(k)) + & - (dS_to_dColHt(i,k) + dS_to_dColHt(i,k-1) ) * (S0(k-1)-S0(k)) ) ) - - if ((mech_TKE(i) + conv_PErel(i)) <= 0.0 .and. Convectively_stable) then - ! Energy is already exhausted, so set Kd = 0 and cycle or exit? - tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 - Kd(i,K) = 0.0 ; Kddt_h(K) = 0.0 - sfc_disconnect = .true. - ! if (.not.debug) exit - - ! The estimated properties for layer k-1 can be calculated, using - ! greatly simplified expressions when Kddt_h = 0. This enables the - ! tridiagonal solver for the whole column to be completed for debugging - ! purposes, and also allows for something akin to convective adjustment - ! in unstable interior regions? - b1 = 1.0 / (b_den_1(i)) - c1(K) = 0.0 - dTe(k-1) = b1 * ( dTe_t2 ) - dSe(k-1) = b1 * ( dSe_t2 ) - - b_den_1(i) = h(i,k) - dT_to_dPE_a(i,k) = dT_to_dPE(i,k) - dS_to_dPE_a(i,k) = dS_to_dPE(i,k) - dT_to_dColHt_a(i,k) = dT_to_dColHt(i,k) - dS_to_dColHt_a(i,k) = dS_to_dColHt(i,k) - - else ! tot_TKE > 0.0 or this is a potentially convectively unstable profile. - sfc_disconnect = .false. - - ! Precalculate some more temporary expressions that are independent of - ! Kddt_h(K). - if (K==2) then - dT_km1_t2 = (T0(k)-T0(k-1)) - dS_km1_t2 = (S0(k)-S0(k-1)) - else - dT_km1_t2 = (T0(k)-T0(k-1)) - & - (Kddt_h(K-1) / b_den_1(i)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) - dS_km1_t2 = (S0(k)-S0(k-1)) - & - (Kddt_h(K-1) / b_den_1(i)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + + do K=2,nz + ! Apply dissipation to the TKE, here applied as an exponential decay + ! due to 3-d turbulent energy being lost to inefficient rotational modes. + + ! There should be several different "flavors" of TKE that decay at + ! different rates. The following form is often used for mechanical + ! stirring from the surface, perhaps due to breaking surface gravity + ! waves and wind-driven turbulence. + Idecay_len_TKE(i) = (CS%TKE_decay * absf(i) / U_Star) * GV%H_to_m + exp_kh = 1.0 + if (Idecay_len_TKE(i) > 0.0) exp_kh = exp(-h(i,k-1)*Idecay_len_TKE(i)) + if (CS%TKE_diagnostics) & + dTKE_mech_decay = dTKE_mech_decay + (exp_kh-1.0) * mech_TKE(i) * IdtdR0 + mech_TKE(i) = mech_TKE(i) * exp_kh + + + ! Accumulate any convectively released potential energy to contribute + ! to wstar and to drive penetrating convection. + if (TKE_forced(i,j,k) > 0.0) then + conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,k) + if (CS%TKE_diagnostics) & + dTKE_forcing = dTKE_forcing + CS%nstar*TKE_forced(i,j,k) * IdtdR0 endif - dTe_term = dTe_t2 + b_den_1(i) * (T0(k-1)-T0(k)) - dSe_term = dSe_t2 + b_den_1(i) * (S0(k-1)-S0(k)) - - ! Using Pr=1 and the diffusivity at the bottom interface (once it is - ! known), determine how much resolved mean kinetic energy (MKE) will be - ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of - ! this to the mTKE budget available for mixing in the next layer. - - if ((CS%MKE_to_TKE_effic > 0.0) .and. (htot(i)*h(i,k) > 0.0)) then - ! This is the energy that would be available from homogenizing the - ! velocities between layer k and the layers above. - dMKE_max = (GV%H_to_kg_m2 * CS%MKE_to_TKE_effic) * 0.5 * & - (h(i,k) / ((htot(i) + h(i,k))*htot(i))) * & - ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2) - ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be - ! extracted by mixing with a finite viscosity. - MKE2_Hharm = (htot(i) + h(i,k) + 2.0*h_neglect) / & - ((htot(i)+h_neglect) * (h(i,k)+h_neglect)) - else - dMKE_max = 0.0 ; MKE2_Hharm = 0.0 + + if (debug) then + mech_TKE_k(i,K) = mech_TKE(i) ; conv_PErel_k(i,K) = conv_PErel(i) endif - ! At this point, Kddt_h(K) will be unknown because its value may depend - ! on how much energy is available. mech_TKE might be negative due to - ! contributions from TKE_forced. - h_tt = htot(i) + h_tt_min - TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*conv_PErel(i) - if (TKE_here > 0.0) then - vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 - Mixing_Length_Used(k) = MAX(MinMixLen,((h_tt*hb_hs(i,K))*vstar) / & - ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar)) - !Note setting Kd_guess0 to Mixing_Length_Used(K) here will - ! change the answers. Therefore, skipping that. - if (.not.CS%Use_MLD_Iteration) then - Kd_guess0 = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & - ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar) - else - Kd_guess0 = vstar * vonKar * Mixing_Length_Used(k) - endif - else - vstar = 0.0 ; Kd_guess0 = 0.0 + ! Determine the total energy + nstar_FC = CS%nstar + if (CS%nstar * conv_PErel(i) > 0.0) then + ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based + ! on a curve fit from the data of Wang (GRL, 2003). + ! Note: Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot(i))**3 / conv_PErel(i)) + nstar_FC = CS%nstar * conv_PErel(i) / (conv_PErel(i) + 0.2 * & + sqrt(0.5 * dt * GV%Rho0 * (absf(i)*(htot(i)*GV%H_to_m))**3 * conv_PErel(i))) endif - Vstar_Used(k) = vstar!Track vstar - Kddt_h_g0 = Kd_guess0*dt_h - - call find_PE_chg(Kddt_h_g0, h(i,k), b_den_1(i), dTe_term, dSe_term, & - dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & - dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & - pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & - dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & - PE_chg=PE_chg_g0, dPEc_dKd=dPEa_dKd_g0, dPE_max=PE_chg_max, & - dPEc_dKd_0=dPEc_dKd_Kd0 ) - - MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) - - if ((PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0))) then - ! This column is convectively unstable. - if (PE_chg_max <= 0.0) then - ! Does MKE_src need to be included in the calculation of vstar here? - TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*(conv_PErel(i)-PE_chg_max) - if (TKE_here > 0.0) then - vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 - Mixing_Length_Used(k) = max(MinMixLen,((h_tt*hb_hs(i,K))*vstar) / & - ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar)) - if (.not.CS%Use_MLD_Iteration) then - ! Note again (as prev) that using Mixing_Length_Used here - ! instead of redoing the computation will change answers... - Kd(i,k) = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & - ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar) - else - Kd(i,k) = vstar * vonKar * Mixing_Length_Used(k) - endif - else - vstar = 0.0 ; Kd(i,k) = 0.0 - endif - Vstar_Used(k) = vstar - - call find_PE_chg(Kd(i,k)*dt_h, h(i,k), b_den_1(i), dTe_term, dSe_term, & - dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & - dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & - pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & - dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & - PE_chg=dPE_conv) - ! Should this be iterated to convergence for Kd? - if (dPE_conv > 0.0) then - Kd(i,k) = Kd_guess0 ; dPE_conv = PE_chg_g0 - else - MKE_src = dMKE_max*(1.0 - exp(-(Kd(i,k)*dt_h) * MKE2_Hharm)) + if (debug) nstar_k(K) = nstar_FC + + tot_TKE = mech_TKE(i) + nstar_FC * conv_PErel(i) + + ! For each interior interface, first discard the TKE to account for + ! mixing of shortwave radiation through the next denser cell. + if (TKE_forced(i,j,k) < 0.0) then + if (TKE_forced(i,j,k) + tot_TKE < 0.0) then + ! The shortwave requirements deplete all the energy in this layer. + if (CS%TKE_diagnostics) then + dTKE_mixing = dTKE_mixing + tot_TKE * IdtdR0 + dTKE_forcing = dTKE_forcing - tot_TKE * IdtdR0 + ! dTKE_unbalanced_forcing = dTKE_unbalanced_forcing + & + ! (TKE_forced(i,j,k) + tot_TKE) * IdtdR0 + dTKE_conv_decay = dTKE_conv_decay + & + (CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 endif + tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 else - ! The energy change does not vary monotonically with Kddt_h. Find the maximum? - Kd(i,k) = Kd_guess0 ; dPE_conv = PE_chg_g0 - endif - conv_PErel(i) = conv_PErel(i) - dPE_conv - mech_TKE(i) = mech_TKE(i) + MKE_src - if (CS%TKE_diagnostics) then - CS%diag_TKE_conv(i,j) = CS%diag_TKE_conv(i,j) - CS%nstar*dPE_conv * IdtdR0 - CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + MKE_src * IdtdR0 - endif - if (sfc_connected(i)) then - CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) - !CS%ML_depth2(i,j) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) + ! Reduce the mechanical and convective TKE proportionately. + TKE_reduc = (tot_TKE + TKE_forced(i,j,k)) / tot_TKE + if (CS%TKE_diagnostics) then + dTKE_mixing = dTKE_mixing - TKE_forced(i,j,k) * IdtdR0 + dTKE_forcing = dTKE_forcing + TKE_forced(i,j,k) * IdtdR0 + dTKE_conv_decay = dTKE_conv_decay + & + (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 + endif + tot_TKE = TKE_reduc*tot_TKE ! = tot_TKE + TKE_forced(i,j,k) + mech_TKE(i) = TKE_reduc*mech_TKE(i) + conv_PErel(i) = TKE_reduc*conv_PErel(i) endif + endif - Kddt_h(K) = Kd(i,k)*dt_h - elseif (tot_TKE + (MKE_src - PE_chg_g0) >= 0.0) then - ! There is energy to support the suggested mixing. Keep that estimate. - Kd(i,k) = Kd_guess0 - Kddt_h(K) = Kddt_h_g0 - - ! Reduce the mechanical and convective TKE proportionately. - tot_TKE = tot_TKE + MKE_src - TKE_reduc = 0.0 ! tot_TKE could be 0 if Convectively_stable is false. - if (tot_TKE > 0.0) TKE_reduc = (tot_TKE - PE_chg_g0) / tot_TKE - - if (CS%TKE_diagnostics) then - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) - PE_chg_g0 * IdtdR0 - CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + MKE_src * IdtdR0 - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + & - (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 - endif - tot_TKE = TKE_reduc*tot_TKE - mech_TKE(i) = TKE_reduc*(mech_TKE(i) + MKE_src) - conv_PErel(i) = TKE_reduc*conv_PErel(i) - if (sfc_connected(i)) then - CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) - !CS%ML_depth2(i,J) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) + ! Precalculate some temporary expressions that are independent of Kddt_h(K). + if (CS%orig_PE_calc) then + if (K==2) then + dTe_t2 = 0.0 ; dSe_t2 = 0.0 + else + dTe_t2 = Kddt_h(K-1) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + dSe_t2 = Kddt_h(K-1) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) endif - elseif (tot_TKE == 0.0) then - ! This can arise if nstar_FC = 0. - Kd(i,k) = 0.0 ; Kddt_h(K) = 0.0 - tot_TKE = 0.0 ; conv_PErel(i) = 0.0 ; mech_TKE(i) = 0.0 + endif + dt_h = (GV%m_to_H**2*dt) / max(0.5*(h(i,k-1)+h(i,k)), 1e-15*h_sum(i)) + + ! This tests whether the layers above and below this interface are in + ! a convetively stable configuration, without considering any effects of + ! mixing at higher interfaces. It is an approximation to the more + ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of + ! mixing across interface K-1. The dT_to_dColHt here are effectively + ! mass-weigted estimates of dSV_dT. + Convectively_stable = ( 0.0 <= & + ( (dT_to_dColHt(i,k) + dT_to_dColHt(i,k-1) ) * (T0(k-1)-T0(k)) + & + (dS_to_dColHt(i,k) + dS_to_dColHt(i,k-1) ) * (S0(k-1)-S0(k)) ) ) + + if ((mech_TKE(i) + conv_PErel(i)) <= 0.0 .and. Convectively_stable) then + ! Energy is already exhausted, so set Kd = 0 and cycle or exit? + tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 + Kd(i,K) = 0.0 ; Kddt_h(K) = 0.0 sfc_disconnect = .true. - else - ! There is not enough energy to support the mixing, so reduce the - ! diffusivity to what can be supported. - Kddt_h_max = Kddt_h_g0 ; Kddt_h_min = 0.0 - TKE_left_max = tot_TKE + (MKE_src - PE_chg_g0) ; TKE_left_min = tot_TKE - - ! As a starting guess, take the minimum of a false position estimate - ! and a Newton's method estimate starting from Kddt_h = 0.0. - Kddt_h_guess = tot_TKE * Kddt_h_max / max( PE_chg_g0 - MKE_src, & - Kddt_h_max * (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) - ! The above expression is mathematically the same as the following - ! except it is not susceptible to division by zero when - ! dPEc_dKd_Kd0 = dMKE_max = 0 . - ! Kddt_h_guess = tot_TKE * min( Kddt_h_max / (PE_chg_g0 - MKE_src), & - ! 1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) - if (debug) then - TKE_left_itt(:) = 0.0 ; dPEa_dKd_itt(:) = 0.0 ; PE_chg_itt(:) = 0.0 - MKE_src_itt(:) = 0.0 ; Kddt_h_itt(:) = 0.0 + ! if (.not.debug) exit + + ! The estimated properties for layer k-1 can be calculated, using + ! greatly simplified expressions when Kddt_h = 0. This enables the + ! tridiagonal solver for the whole column to be completed for debugging + ! purposes, and also allows for something akin to convective adjustment + ! in unstable interior regions? + b1 = 1.0 / hp_a(i) + c1(K) = 0.0 + if (CS%orig_PE_calc) then + dTe(k-1) = b1 * ( dTe_t2 ) + dSe(k-1) = b1 * ( dSe_t2 ) endif - do itt=1,max_itt - call find_PE_chg(Kddt_h_guess, h(i,k), b_den_1(i), dTe_term, dSe_term, & - dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & - dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & - pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & - dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & - PE_chg=PE_chg, dPEc_dKd=dPEc_dKd ) - MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) - dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) - - TKE_left = tot_TKE + (MKE_src - PE_chg) - if (debug) then - Kddt_h_itt(itt) = Kddt_h_guess ; MKE_src_itt(itt) = MKE_src - PE_chg_itt(itt) = PE_chg ; TKE_left_itt(itt) = TKE_left - dPEa_dKd_itt(itt) = dPEc_dKd + + hp_a(i) = h(i,k) + dT_to_dPE_a(i,k) = dT_to_dPE(i,k) + dS_to_dPE_a(i,k) = dS_to_dPE(i,k) + dT_to_dColHt_a(i,k) = dT_to_dColHt(i,k) + dS_to_dColHt_a(i,k) = dS_to_dColHt(i,k) + + else ! tot_TKE > 0.0 or this is a potentially convectively unstable profile. + sfc_disconnect = .false. + + ! Precalculate some more temporary expressions that are independent of + ! Kddt_h(K). + if (CS%orig_PE_calc) then + if (K==2) then + dT_km1_t2 = (T0(k)-T0(k-1)) + dS_km1_t2 = (S0(k)-S0(k-1)) + else + dT_km1_t2 = (T0(k)-T0(k-1)) - & + (Kddt_h(K-1) / hp_a(i)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + dS_km1_t2 = (S0(k)-S0(k-1)) - & + (Kddt_h(K-1) / hp_a(i)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) endif - ! Store the new bounding values, bearing in mind that min and max - ! here refer to Kddt_h and dTKE_left/dKddt_h < 0: - if (TKE_left >= 0.0) then - Kddt_h_min = Kddt_h_guess ; TKE_left_min = TKE_left + dTe_term = dTe_t2 + hp_a(i) * (T0(k-1)-T0(k)) + dSe_term = dSe_t2 + hp_a(i) * (S0(k-1)-S0(k)) + else + if (K<=2) then + Th_a(k-1) = h(i,k-1) * T0(k-1) ; Sh_a(k-1) = h(i,k-1) * S0(k-1) else - Kddt_h_max = Kddt_h_guess ; TKE_left_max = TKE_left + Th_a(k-1) = h(i,k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2) + Sh_a(k-1) = h(i,k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2) endif + Th_b(k) = h(i,k) * T0(k) ; Sh_b(k) = h(i,k) * S0(k) + endif + + ! Using Pr=1 and the diffusivity at the bottom interface (once it is + ! known), determine how much resolved mean kinetic energy (MKE) will be + ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of + ! this to the mTKE budget available for mixing in the next layer. + + if ((CS%MKE_to_TKE_effic > 0.0) .and. (htot(i)*h(i,k) > 0.0)) then + ! This is the energy that would be available from homogenizing the + ! velocities between layer k and the layers above. + dMKE_max = (GV%H_to_kg_m2 * CS%MKE_to_TKE_effic) * 0.5 * & + (h(i,k) / ((htot(i) + h(i,k))*htot(i))) * & + ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2) + ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be + ! extracted by mixing with a finite viscosity. + MKE2_Hharm = (htot(i) + h(i,k) + 2.0*h_neglect) / & + ((htot(i)+h_neglect) * (h(i,k)+h_neglect)) + else + dMKE_max = 0.0 ; MKE2_Hharm = 0.0 + endif - ! Try to use Newton's method, but if it would go outside the bracketed - ! values use the false-position method instead. - use_Newt = .true. - if (dPEc_dKd - dMKE_src_dK <= 0.0) then - use_Newt = .false. + ! At this point, Kddt_h(K) will be unknown because its value may depend + ! on how much energy is available. mech_TKE might be negative due to + ! contributions from TKE_forced. + h_tt = htot(i) + h_tt_min + TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*conv_PErel(i) + if (TKE_here > 0.0) then + vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 + hbs_here = GV%H_to_m * min(hb_hs(i,K), MixLen_shape(K)) + Mixing_Length_Used(k) = MAX(CS%min_mix_len,((h_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)) + !Note setting Kd_guess0 to Mixing_Length_Used(K) here will + ! change the answers. Therefore, skipping that. + if (.not.CS%Use_MLD_Iteration) then + Kd_guess0 = vstar * vonKar * ((h_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar) else - dKddt_h_Newt = TKE_left / (dPEc_dKd - dMKE_src_dK) - Kddt_h_Newt = Kddt_h_guess + dKddt_h_Newt - if ((Kddt_h_Newt > Kddt_h_max) .or. (Kddt_h_Newt < Kddt_h_min)) & - use_Newt = .false. + Kd_guess0 = vstar * vonKar * Mixing_Length_Used(k) endif + else + vstar = 0.0 ; Kd_guess0 = 0.0 + endif + Vstar_Used(k) = vstar ! Track vstar + Kddt_h_g0 = Kd_guess0*dt_h + + if (CS%orig_PE_calc) then + call find_PE_chg_orig(Kddt_h_g0, h(i,k), hp_a(i), dTe_term, dSe_term, & + dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & + pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + PE_chg=PE_chg_g0, dPEc_dKd=dPEa_dKd_g0, dPE_max=PE_chg_max, & + dPEc_dKd_0=dPEc_dKd_Kd0 ) + else + call find_PE_chg(0.0, Kddt_h_g0, hp_a(i), h(i,k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), dT_to_dPE(i,k), dS_to_dPE(i,k), & + pres(i,K), dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + PE_chg=PE_chg_g0, dPEc_dKd=dPEa_dKd_g0, dPE_max=PE_chg_max, & + dPEc_dKd_0=dPEc_dKd_Kd0 ) + endif - if (use_Newt) then - Kddt_h_next = Kddt_h_guess + dKddt_h_Newt - dKddt_h = dKddt_h_Newt + MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) + + if ((PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0))) then + ! This column is convectively unstable. + if (PE_chg_max <= 0.0) then + ! Does MKE_src need to be included in the calculation of vstar here? + TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*(conv_PErel(i)-PE_chg_max) + if (TKE_here > 0.0) then + vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 + hbs_here = GV%H_to_m * min(hb_hs(i,K), MixLen_shape(K)) + Mixing_Length_Used(k) = max(CS%min_mix_len,((h_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)) + if (.not.CS%Use_MLD_Iteration) then + ! Note again (as prev) that using Mixing_Length_Used here + ! instead of redoing the computation will change answers... + Kd(i,k) = vstar * vonKar * ((h_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar) + else + Kd(i,k) = vstar * vonKar * Mixing_Length_Used(k) + endif + else + vstar = 0.0 ; Kd(i,k) = 0.0 + endif + Vstar_Used(k) = vstar + + if (CS%orig_PE_calc) then + call find_PE_chg_orig(Kd(i,k)*dt_h, h(i,k), hp_a(i), dTe_term, dSe_term, & + dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & + pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + PE_chg=dPE_conv) + else + call find_PE_chg(0.0, Kd(i,k)*dt_h, hp_a(i), h(i,k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), dT_to_dPE(i,k), dS_to_dPE(i,k), & + pres(i,K), dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + PE_chg=dPE_conv) + endif + ! Should this be iterated to convergence for Kd? + if (dPE_conv > 0.0) then + Kd(i,k) = Kd_guess0 ; dPE_conv = PE_chg_g0 + else + MKE_src = dMKE_max*(1.0 - exp(-(Kd(i,k)*dt_h) * MKE2_Hharm)) + endif else - Kddt_h_next = (TKE_left_max * Kddt_h_min - Kddt_h_max * TKE_left_min) / & - (TKE_left_max - TKE_left_min) - dKddt_h = Kddt_h_next - Kddt_h_guess + ! The energy change does not vary monotonically with Kddt_h. Find the maximum? + Kd(i,k) = Kd_guess0 ; dPE_conv = PE_chg_g0 + endif + conv_PErel(i) = conv_PErel(i) - dPE_conv + mech_TKE(i) = mech_TKE(i) + MKE_src + if (CS%TKE_diagnostics) then + dTKE_conv = dTKE_conv - CS%nstar*dPE_conv * IdtdR0 + dTKE_MKE = dTKE_MKE + MKE_src * IdtdR0 + endif + if (sfc_connected(i)) then + CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) + !CS%ML_depth2(i,j) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) endif - if ((abs(dKddt_h) < 1e-9*Kddt_h_guess) .or. (itt==max_itt)) then - ! Use the old value so that the energy calculation does not need to be repeated. - if (debug) num_itts(K) = itt - exit - else - Kddt_h_guess = Kddt_h_next + Kddt_h(K) = Kd(i,k)*dt_h + elseif (tot_TKE + (MKE_src - PE_chg_g0) >= 0.0) then + ! There is energy to support the suggested mixing. Keep that estimate. + Kd(i,k) = Kd_guess0 + Kddt_h(K) = Kddt_h_g0 + + ! Reduce the mechanical and convective TKE proportionately. + tot_TKE = tot_TKE + MKE_src + TKE_reduc = 0.0 ! tot_TKE could be 0 if Convectively_stable is false. + if (tot_TKE > 0.0) TKE_reduc = (tot_TKE - PE_chg_g0) / tot_TKE + + if (CS%TKE_diagnostics) then + dTKE_mixing = dTKE_mixing - PE_chg_g0 * IdtdR0 + dTKE_MKE = dTKE_MKE + MKE_src * IdtdR0 + dTKE_conv_decay = dTKE_conv_decay + & + (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 + endif + tot_TKE = TKE_reduc*tot_TKE + mech_TKE(i) = TKE_reduc*(mech_TKE(i) + MKE_src) + conv_PErel(i) = TKE_reduc*conv_PErel(i) + if (sfc_connected(i)) then + CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) + !CS%ML_depth2(i,J) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) endif - enddo - Kd(i,K) = Kddt_h_guess / dt_h ; Kddt_h(K) = Kd(i,K)*dt_h - - ! All TKE should have been consumed. - if (CS%TKE_diagnostics) then - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) - (tot_TKE + MKE_src) * IdtdR0 - CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + MKE_src * IdtdR0 - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + & - (CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 + elseif (tot_TKE == 0.0) then + ! This can arise if nstar_FC = 0. + Kd(i,k) = 0.0 ; Kddt_h(K) = 0.0 + tot_TKE = 0.0 ; conv_PErel(i) = 0.0 ; mech_TKE(i) = 0.0 + sfc_disconnect = .true. + else + ! There is not enough energy to support the mixing, so reduce the + ! diffusivity to what can be supported. + Kddt_h_max = Kddt_h_g0 ; Kddt_h_min = 0.0 + TKE_left_max = tot_TKE + (MKE_src - PE_chg_g0) ; TKE_left_min = tot_TKE + + ! As a starting guess, take the minimum of a false position estimate + ! and a Newton's method estimate starting from Kddt_h = 0.0. + Kddt_h_guess = tot_TKE * Kddt_h_max / max( PE_chg_g0 - MKE_src, & + Kddt_h_max * (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + ! The above expression is mathematically the same as the following + ! except it is not susceptible to division by zero when + ! dPEc_dKd_Kd0 = dMKE_max = 0 . + ! Kddt_h_guess = tot_TKE * min( Kddt_h_max / (PE_chg_g0 - MKE_src), & + ! 1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + if (debug) then + TKE_left_itt(:) = 0.0 ; dPEa_dKd_itt(:) = 0.0 ; PE_chg_itt(:) = 0.0 + MKE_src_itt(:) = 0.0 ; Kddt_h_itt(:) = 0.0 + endif + do itt=1,max_itt + if (CS%orig_PE_calc) then + call find_PE_chg_orig(Kddt_h_guess, h(i,k), hp_a(i), dTe_term, dSe_term, & + dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), & + pres(i,K), dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + PE_chg=PE_chg, dPEc_dKd=dPEc_dKd ) + else + call find_PE_chg(0.0, Kddt_h_guess, hp_a(i), h(i,k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(i,k-1), dS_to_dPE_a(i,k-1), dT_to_dPE(i,k), dS_to_dPE(i,k), & + pres(i,K), dT_to_dColHt_a(i,k-1), dS_to_dColHt_a(i,k-1), & + dT_to_dColHt(i,k), dS_to_dColHt(i,k), & + PE_chg=dPE_conv) + endif + MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) + dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) + + TKE_left = tot_TKE + (MKE_src - PE_chg) + if (debug) then + Kddt_h_itt(itt) = Kddt_h_guess ; MKE_src_itt(itt) = MKE_src + PE_chg_itt(itt) = PE_chg ; TKE_left_itt(itt) = TKE_left + dPEa_dKd_itt(itt) = dPEc_dKd + endif + ! Store the new bounding values, bearing in mind that min and max + ! here refer to Kddt_h and dTKE_left/dKddt_h < 0: + if (TKE_left >= 0.0) then + Kddt_h_min = Kddt_h_guess ; TKE_left_min = TKE_left + else + Kddt_h_max = Kddt_h_guess ; TKE_left_max = TKE_left + endif + + ! Try to use Newton's method, but if it would go outside the bracketed + ! values use the false-position method instead. + use_Newt = .true. + if (dPEc_dKd - dMKE_src_dK <= 0.0) then + use_Newt = .false. + else + dKddt_h_Newt = TKE_left / (dPEc_dKd - dMKE_src_dK) + Kddt_h_Newt = Kddt_h_guess + dKddt_h_Newt + if ((Kddt_h_Newt > Kddt_h_max) .or. (Kddt_h_Newt < Kddt_h_min)) & + use_Newt = .false. + endif + + if (use_Newt) then + Kddt_h_next = Kddt_h_guess + dKddt_h_Newt + dKddt_h = dKddt_h_Newt + else + Kddt_h_next = (TKE_left_max * Kddt_h_min - Kddt_h_max * TKE_left_min) / & + (TKE_left_max - TKE_left_min) + dKddt_h = Kddt_h_next - Kddt_h_guess + endif + + if ((abs(dKddt_h) < 1e-9*Kddt_h_guess) .or. (itt==max_itt)) then + ! Use the old value so that the energy calculation does not need to be repeated. + if (debug) num_itts(K) = itt + exit + else + Kddt_h_guess = Kddt_h_next + endif + enddo + Kd(i,K) = Kddt_h_guess / dt_h ; Kddt_h(K) = Kd(i,K)*dt_h + + ! All TKE should have been consumed. + if (CS%TKE_diagnostics) then + dTKE_mixing = dTKE_mixing - (tot_TKE + MKE_src) * IdtdR0 + dTKE_MKE = dTKE_MKE + MKE_src * IdtdR0 + dTKE_conv_decay = dTKE_conv_decay + & + (CS%nstar-nstar_FC) * conv_PErel(i) * IdtdR0 + endif + + if (sfc_connected(i)) CS%ML_depth(i,J) = CS%ML_depth(i,J) + & + (PE_chg / PE_chg_g0) * GV%H_to_m * h(i,k) + tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 + sfc_disconnect = .true. endif - if (sfc_connected(i)) CS%ML_depth(i,J) = CS%ML_depth(i,J) + & - (PE_chg / PE_chg_g0) * GV%H_to_m * h(i,k) - tot_TKE = 0.0 ; mech_TKE(i) = 0.0 ; conv_PErel(i) = 0.0 - sfc_disconnect = .true. + Kddt_h(K) = Kd(i,K)*dt_h + ! At this point, the final value of Kddt_h(K) is known, so the + ! estimated properties for layer k-1 can be calculated. + b1 = 1.0 / (hp_a(i) + Kddt_h(K)) + c1(K) = Kddt_h(K) * b1 + if (CS%orig_PE_calc) then + dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) + dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) + endif + + hp_a(i) = h(i,k) + (hp_a(i) * b1) * Kddt_h(K) + dT_to_dPE_a(i,k) = dT_to_dPE(i,k) + c1(K)*dT_to_dPE_a(i,k-1) + dS_to_dPE_a(i,k) = dS_to_dPE(i,k) + c1(K)*dS_to_dPE_a(i,k-1) + dT_to_dColHt_a(i,k) = dT_to_dColHt(i,k) + c1(K)*dT_to_dColHt_a(i,k-1) + dS_to_dColHt_a(i,k) = dS_to_dColHt(i,k) + c1(K)*dS_to_dColHt_a(i,k-1) + + endif ! tot_TKT > 0.0 branch. Kddt_h(K) has been set. + + ! Store integrated velocities and thicknesses for MKE conversion calculations. + if (sfc_disconnect) then + ! There is no turbulence at this interface, so zero out the running sums. + uhtot(i) = u(i,k)*h(i,k) + vhtot(i) = v(i,k)*h(i,k) + htot(i) = h(i,k) + sfc_connected(i) = .false. + else + uhtot(i) = uhtot(i) + u(i,k)*h(i,k) + vhtot(i) = vhtot(i) + v(i,k)*h(i,k) + htot(i) = htot(i) + h(i,k) endif - Kddt_h(K) = Kd(i,K)*dt_h - ! At this point, the final value of Kddt_h(K) is known, so the - ! estimated properties for layer k-1 can be calculated. - b1 = 1.0 / (b_den_1(i) + Kddt_h(K)) - c1(K) = Kddt_h(K) * b1 - dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) - dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) - - b_den_1(i) = h(i,k) + (b_den_1(i) * b1) * Kddt_h(K) - dT_to_dPE_a(i,k) = dT_to_dPE(i,k) + c1(K)*dT_to_dPE_a(i,k-1) - dS_to_dPE_a(i,k) = dS_to_dPE(i,k) + c1(K)*dS_to_dPE_a(i,k-1) - dT_to_dColHt_a(i,k) = dT_to_dColHt(i,k) + c1(K)*dT_to_dColHt_a(i,k-1) - dS_to_dColHt_a(i,k) = dS_to_dColHt(i,k) + c1(K)*dS_to_dColHt_a(i,k-1) - - endif ! tot_TKT > 0.0 branch. Kddt_h(K) has been set. - - ! Store integrated velocities and thicknesses for MKE conversion calculations. - if (sfc_disconnect) then - ! There is no turbulence at this interface, so zero out the running sums. - uhtot(i) = u(i,k)*h(i,k) - vhtot(i) = v(i,k)*h(i,k) - htot(i) = h(i,k) - sfc_connected(i) = .false. - else - uhtot(i) = uhtot(i) + u(i,k)*h(i,k) - vhtot(i) = vhtot(i) + v(i,k)*h(i,k) - htot(i) = htot(i) + h(i,k) - endif + if (debug) then + if (k==2) then + Te(1) = b1*(h(i,1)*T0(1)) + Se(1) = b1*(h(i,1)*S0(1)) + else + Te(k-1) = b1 * (h(i,k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2)) + Se(k-1) = b1 * (h(i,k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2)) + endif + endif + enddo + Kd(i,nz+1) = 0.0 if (debug) then - if (k==2) then - Te(1) = b1*(h(i,1)*T0(1)) - Se(1) = b1*(h(i,1)*S0(1)) - else - Te(k-1) = b1 * (h(i,k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2)) - Se(k-1) = b1 * (h(i,k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2)) - endif + ! Complete the tridiagonal solve for Te. + b1 = 1.0 / hp_a(i) + Te(nz) = b1 * (h(i,nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) + Se(nz) = b1 * (h(i,nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) + do k=nz-1,1,-1 + Te(k) = Te(k) + c1(K+1)*Te(k+1) + Se(k) = Se(k) + c1(K+1)*Se(k+1) + enddo endif - enddo - Kd(i,nz+1) = 0.0 - - if (debug) then - ! Complete the tridiagonal solve for Te. - b1 = 1.0 / (b_den_1(i)) - Te(nz) = b1 * (h(i,nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) - Se(nz) = b1 * (h(i,nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) - do k=nz-1,1,-1 - Te(k) = Te(k) + c1(K+1)*Te(k+1) - Se(k) = Se(k) + c1(K+1)*Se(k+1) + if (present(dT_expected)) then + do k=1,nz ; dT_expected(i,j,k) = Te(k) - T0(k) ; enddo + endif + if (present(dS_expected)) then + do k=1,nz ; dS_expected(i,j,k) = Se(k) - S0(k) ; enddo + endif + if (debug) then + dPE_debug = 0.0 + do k=1,nz + dPE_debug = dPE_debug + (dT_to_dPE(i,k) * (Te(k) - T0(k)) + & + dS_to_dPE(i,k) * (Se(k) - S0(k))) + enddo + mixing_debug = dPE_debug * IdtdR0 + endif + k = nz ! This is here to allow a breakpoint to be set. + + !/BGR: The following lines are used for the iteration + ITmax(obl_it) = max_MLD ! Track max } + ITmin(obl_it) = min_MLD ! Track min } For debug purpose + ITguess(obl_it) = MLD_guess ! Track guess } + MLD_FOUND=0.0 ; FIRST_OBL=.true. + ! MLD_FOUND=CS%ML_depth2(i,J) + do k=2,nz + if (FIRST_OBL) then !Breaks when OBL found + if (Vstar_Used(k) > 1.e-10 .and. k < nz) then + !If within OBL, keep integrating to find OBL + !/ Add thickness of level above to OBL depth. + MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m + !1. Check if guess was too shallow + !Adding -1 m as cushion, will help avoid + ! non-convergence flag when nearly converged. + elseif (MLD_FOUND-CS%MLD_tol > MLD_guess) then + !elseif (MLD_FOUND > MLD_guess) then + !/ Guess was too shallow, set new minimum guess + min_MLD = MLD_guess + FIRST_OBL = .false. !Break OBL loop + !2. Check if Guess minus found MLD + ! is less than thickness of level (=converged) + ! - We could try to add a more precise + ! value for found MLD, but seems difficult to + ! to contrain beyond within a level. + elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then + ! Converged. Exit iteration. + FIRST_OBL = .false.!Break OBL loop + OBL_CONVERGED = .true.!Break convergence loop + ! Testing Output, comment for use. + !print*,'Converged--------' + !print*,MLD_FOUND,MLD_guess + !/ + if (OBL_IT_STATS) then !Compute iteration statistics + MAXIT = max(MAXIT,obl_it) + MINIT = min(MINIT,obl_it) + SUMIT = SUMIT+obl_it + NUMIT = NUMIT+1 + print*,MAXIT,MINIT,SUMIT/NUMIT + endif + !BGR this can be where MLD is stored for next time... + CS%ML_Depth2(i,j) = MLD_guess + !/ + !2. If not, guess was too deep + else + !Guess was too deep, set new maximum guess + max_MLD = MLD_guess !We know this guess was too deep + FIRST_OBL = .false.!Break OBL loop + endif + endif enddo + ! For next pass, guess average of minimum and maximum values. + MLD_guess = min_MLD*0.5 + max_MLD*0.5 + ITresult(obl_it) = MLD_FOUND + endif ; enddo ! Iteration loop for converged boundary layer thickness. + + if (.not.OBL_CONVERGED) then + !/Temp output, warn that EPBL didn't converge + !/Print guess/found for every iteration step + !print*,'EPBL MLD DID NOT CONVERGE' + NOTCONVERGED=NOTCONVERGED+1 + !do obl_it=1,max_obl_it + ! print*,ITmin(obl_it),ITmax(obl_it) + ! print*,obl_it,ITguess(obl_it),ITresult(obl_it) + !enddo + !Activate to print out some output when not converged + !{ + !print*,'Min/Max: ',ITmin(50),ITmax(50) + !print*,'Guess/result: ',ITguess(50),ITresult(50) + !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,& + ! real(NOTCONVERGED)/real(CONVERGED) + !} + !stop !Kill if not converged during testing. + else + CONVERGED=CONVERGED+1 endif - if (present(dT_expected)) then - do k=1,nz ; dT_expected(i,j,k) = Te(k) - T0(k) ; enddo - endif - if (present(dS_expected)) then - do k=1,nz ; dS_expected(i,j,k) = Se(k) - S0(k) ; enddo + + if (CS%TKE_diagnostics) then + CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + dTKE_MKE + CS%diag_TKE_conv(i,j) = CS%diag_TKE_conv(i,j) + dTKE_conv + CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + dTKE_forcing + CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) + dTKE_mixing + CS%diag_TKE_mech_decay(i,j) = CS%diag_TKE_mech_decay(i,j) + dTKE_mech_decay + CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + dTKE_conv_decay + ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + dTKE_unbalanced endif - if (debug) then - dPE_debug = 0.0 + if (CS%Mixing_Diagnostics) then + !Write to 3-D for outputing Mixing length and + ! velocity scale. do k=1,nz - dPE_debug = dPE_debug + (dT_to_dPE(i,k) * (Te(k) - T0(k)) + & - dS_to_dPE(i,k) * (Se(k) - S0(k))) + CS%Mixing_Length(i,j,k) = Mixing_Length_Used(k) + CS%Velocity_Scale(i,j,k) = Vstar_Used(k) enddo - mixing_debug = dPE_debug * IdtdR0 endif - k = nz ! This is here to allow a breakpoint to be set. - - !/BGR: The following lines are used for the iteration - ITmax(obl_it)=max_MLD ! Track max } - ITmin(obl_it)=min_MLD ! Track min } For debug purpose - ITguess(obl_it) = MLD_GUESS ! Track guess } - MLD_FOUND=0.0;FIRST_OBL=.true.; - ! MLD_FOUND=CS%ML_depth2(i,J) - do k=2,nz - IF (FIRST_OBL) then !Breaks when OBL found - if (Vstar_Used(k).gt.1.e-10 .and. k.lt.nz) then - !If within OBL, keep integrating to find OBL - !/ Add thickness of level above to OBL depth. - MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m - !1. Check if guess was too shallow - !Adding -1 m as cushion, will help avoid - ! non-convergence flag when nearly converged. - elseif (MLD_FOUND-1.0.GT.MLD_GUESS) then - !elseif (MLD_FOUND.GT.MLD_GUESS) then - !/ Guess was too shallow, set new minimum guess - MIN_MLD=MLD_GUESS - FIRST_OBL=.false. !Break OBL loop - !2. Check if Guess minus found MLD - ! is less than thickness of level (=converged) - ! - We could try to add a more precise - ! value for found MLD, but seems difficult to - ! to contrain beyond within a level. - elseif ((MLD_GUESS-MLD_FOUND).lt.max(1.,h(i,k-1)*GV%H_to_m)) then - ! Converged. Exit iteration. - FIRST_OBL=.false.!Break OBL loop - OBL_CONVERGED=.true.!Break convergence loop - ! Testing Output, comment for use. - !print*,'Converged--------' - !print*,MLD_FOUND,MLD_GUESS - !/ - IF (OBL_IT_STATS) then !Compute iteration statistics - MAXIT=max(MAXIT,obl_it) - MINIT=min(MINIT,obl_it) - SUMIT=SUMIT+obl_it - NUMIT=NUMIT+1 - print*,MAXIT,MINIT,SUMIT/NUMIT - ENDIF - !BGR this can be where MLD is stored for next time... - CS%ML_Depth2(i,j)=MLD_GUESS - !/ - !2. If not, guess was too deep - else - !Guess was too deep, set new maximum guess - MAX_MLD=MLD_GUESS !We know this guess was too deep - FIRST_OBL=.false.!Break OBL loop - endif - endif - enddo - ! For next pass, guess average of minimum and maximum values. - MLD_GUESS = MIN_MLD*0.5 + MAX_MLD*0.5 - ITresult(obl_it)=MLD_FOUND - ENDIF!Check for convergence loop - ENDDO!Iteration loop - if (.not.OBL_CONVERGED) then - !/Temp output, warn that EPBL didn't converge - !/Print guess/found for every iteration step - !print*,'EPBL MLD DID NOT CONVERGE' - NOTCONVERGED=NOTCONVERGED+1 - !do obl_it=1,max_obl_it - ! print*,ITmin(obl_it),ITmax(obl_it) - ! print*,obl_it,ITguess(obl_it),ITresult(obl_it) - !enddo - !Activate to print out some output when not converged - !{ - !print*,'Min/Max: ',ITmin(50),ITmax(50) - !print*,'Guess/result: ',ITguess(50),ITresult(50) - !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,& - ! real(NOTCONVERGED)/real(CONVERGED) - !} - !stop !Kill if not converged during testing. - else - CONVERGED=CONVERGED+1 - endif - if (cs%Mixing_Diagnostics) then - !Write to 3-D for outputing Mixing length and - ! velocity scale. - do k=1,nz - cs%Mixing_Length(i,j,k)=Mixing_Length_Used(k) - cs%Velocity_Scale(i,j,k)=Vstar_Used(k) - enddo - endif - + else ! For masked points, Kd_int must still be set (to 0) because it has intent(out). do K=1,nz+1 @@ -1192,16 +1298,222 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & end subroutine energetic_PBL -subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interfaces's diapycnal diffusivity times a timestep. +subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & + dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, & + pres, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & + PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, ColHt_cor) + real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface, in units of H (m or kg-2). + real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface, in units of H (m or kg-2). + real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers, in K H. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers, in K H. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers, in K H. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers, in K H. + real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers above, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers above, in J m-2 ppt-1. + real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers below, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers below, in J m-2 ppt-1. + real, intent(in) :: pres !< The hydrostatic interface pressure, which is used to relate + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing, in Pa. + real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above, in m K-1. + real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above, in m ppt-1. + real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below, in m K-1. + real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below, in m ppt-1. + + real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying + !! Kddt_h at the present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h, + !! in units of J m-2 H-1. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realizedd by applying a huge value of Kddt_h at the + !! present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the + !! limit where Kddt_h = 0, in J m-2 H-1. + real, optional, intent(out) :: ColHt_cor !< The correction to PE_chg that is made due to a net + !! change in the column height, in J m-2. + + real :: hps ! The sum of the two effective pivot thicknesses, in H. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term, in H2. + real :: dT_c ! The core term in the expressions for the temperature changes, in K H2. + real :: dS_c ! The core term in the expressions for the salinity changes, in psu H2. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes, J m-3. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes, J m-3. + real :: ColHt_chg ! The change in the column height, in m. + real :: y1 ! A local temporary term, in units of H-3 or H-4 in various contexts. + + ! The expression for the change in potential energy used here is derived + ! from the expression for the final estimates of the changes in temperature + ! and salinities, and then extensively manipulated to get it into its most + ! succint form. The derivation is not necessarily obvious, but it demonstrably + ! works by comparison with separate calculations of the energy changes after + ! the tridiagonal solver for the final changes in temperature and salinity are + ! applied. + + hps = hp_a + hp_b + bdt1 = hp_a * hp_b + Kddt_h0 * hps + dT_c = hp_a * Th_b - hp_b * Th_a + dS_c = hp_a * Sh_b - hp_b * Sh_a + PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & + hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) + ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & + hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) + + if (present(PE_chg)) then + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKddt_h. + y1 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + PE_chg = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) PE_chg = PE_chg - pres * ColHt_chg + if (present(ColHt_cor)) ColHt_cor = -pres * min(ColHt_chg, 0.0) + else if (present(ColHt_cor)) then + y1 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + ColHt_cor = -pres * min(ColHt_core * y1, 0.0) + endif + + if (present(dPEc_dKd)) then + ! Find the derivative of the potential energy change with dKddt_h. + y1 = 1.0 / (bdt1 + dKddt_h * hps)**2 + dPEc_dKd = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres * ColHt_chg + endif + + if (present(dPE_max)) then + ! This expression is the limit of PE_chg for infinite dKddt_h. + y1 = 1.0 / (bdt1 * hps) + dPE_max = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPE_max = dPE_max - pres * ColHt_chg + endif + + if (present(dPEc_dKd_0)) then + ! This expression is the limit of dPEc_dKd for dKddt_h = 0. + y1 = 1.0 / bdt1**2 + dPEc_dKd_0 = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres * ColHt_chg + endif + +end subroutine find_PE_chg + +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interfaces's diapycnal diffusivity times a timestep +!! using the original form used in the first version of ePBL. +subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, & dT_to_dPEa, dS_to_dPEa, pres, dT_to_dColHt_k, & dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, & PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0) - real, intent(in) :: Kddt_h, h_k, b_den_1, dTe_term, dSe_term - real, intent(in) :: dT_km1_t2, dS_km1_t2, pres - real, intent(in) :: dT_to_dPE_k, dS_to_dPE_k, dT_to_dPEa, dS_to_dPEa - real, intent(in) :: dT_to_dColHt_k, dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta - real, optional, intent(out) :: PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0 + real, intent(in) :: Kddt_h !< The diffusivity at an interface times the time step and + !! divided by the average of the thicknesses around the + !! interface, in units of H (m or kg-2). + real, intent(in) :: h_k !< The thickness of the layer below the interface, in H. + real, intent(in) :: b_den_1 !< The first term in the denominator of the pivot + !! for the tridiagonal solver, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: dTe_term !< A diffusivity-independent term related to the + !! temperature change in the layer below the interface, in K H. + real, intent(in) :: dSe_term !< A diffusivity-independent term related to the + !! salinity change in the layer below the interface, in ppt H. + real, intent(in) :: dT_km1_t2 !< A diffusivity-independent term related to the + !! temperature change in the layer above the interface, in K. + real, intent(in) :: dS_km1_t2 !< A diffusivity-independent term related to the + !! salinity change in the layer above the interface, in ppt. + real, intent(in) :: pres !< The hydrostatic interface pressure, which is used to relate + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing, in Pa. + real, intent(in) :: dT_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers below, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers below, in J m-2 ppt-1. + real, intent(in) :: dT_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers above, in J m-2 K-1. + real, intent(in) :: dS_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers above, in J m-2 ppt-1. + real, intent(in) :: dT_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below, in m K-1. + real, intent(in) :: dS_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below, in m ppt-1. + real, intent(in) :: dT_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above, in m K-1. + real, intent(in) :: dS_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above, in m ppt-1. + + real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying + !! Kddt_h at the present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h, + !! in units of J m-2 H-1. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realizedd by applying a huge value of Kddt_h at the + !! present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the + !! limit where Kddt_h = 0, in J m-2 H-1. ! This subroutine determines the total potential energy change due to mixing ! at an interface, including all of the implicit effects of the prescribed @@ -1211,64 +1523,7 @@ subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & ! The comments describing these arguments are for a downward mixing pass, but ! this routine can also be used for an upward pass with the sense of direction ! reversed. -! -! Arguments: Kddt_h - The diffusivity at an interface times the time step and -! divided by the average of the thicknesses around the -! interface, in units of H (m or kg-2). -! (in) h_k - The thickness of the layer below the interface, in H. -! (in) b_den_1 - The first first term in the denominator of the pivot -! for the tridiagonal solver, given by h_k plus a term that -! is a fraction (determined from the tridiagonal solver) of -! Kddt_h for the interface above, in H. -! (in) dTe_term - A diffusivity-independent term related to the -! temperature change in the layer below the interface, in K H. -! (in) dSe_term - A diffusivity-independent term related to the -! salinity change in the layer below the interface, in ppt H. -! (in) dT_km1_t2 - A diffusivity-independent term related to the -! temperature change in the layer above the interface, in K. -! (in) dS_km1_t2 - A diffusivity-independent term related to the -! salinity change in the layer above the interface, in ppt. -! (in) dT_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating -! a layer's temperature change to the change in column -! potential energy, in J m-2 K-1. -! (in) dS_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating -! a layer's salinity change to the change in column -! potential energy, in J m-2 ppt-1. -! (in) dT_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating -! a layer's temperature change to the change in column -! potential energy, including all implicit diffusive changes -! in the temperatures of all the layers above, in J m-2 K-1. -! (in) dS_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating -! a layer's salinity change to the change in column -! potential energy, including all implicit diffusive changes -! in the salinities of all the layers above, in J m-2 ppt-1. -! (in) pres - The hydrostatic interface pressure, which is used to relate -! the changes in column thickness to the energy that is radiated -! as gravity waves and unavailable to drive mixing, in Pa. -! (in) dT_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dT) relating -! a layer's temperature change to the change in column -! height, in m K-1. -! (in) dS_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dS) relating -! a layer's salinity change to the change in column -! height, in m ppt-1. -! (in) dT_to_dColHta - A factor (mass_lay*dSColHtc_vol/dT) relating -! a layer's temperature change to the change in column -! height, including all implicit diffusive changes -! in the temperatures of all the layers above, in m K-1. -! (in) dS_to_dColHta - A factor (mass_lay*dSColHtc_vol/dS) relating -! a layer's salinity change to the change in column -! height, including all implicit diffusive changes -! in the salinities of all the layers above, in m ppt-1. -! (out,opt) PE_chg - The change in column potential energy from applying -! Kddt_h at the present interface, in J m-2. -! (out,opt) dPE_chg_dKd - The partial derivative of PE_chg with Kddt_h, -! in units of J m-2 H-1. -! (out,opt) PE_max - The maximum change in column potential energy that could -! be realizedd by applying a huge value of Kddt_h at the -! present interface, in J m-2. -! (out,opt) dPEc_dKc0 - The partial derivative of PE_chg with Kddt_h in the -! limit where Kddt_h = 0, in J m-2 H-1. - ! b_den_1 - The first term in the denominator of b1, in m or kg m-2. + real :: b1 ! b1 is used by the tridiagonal solver, in H-1. real :: b1Kd ! Temporary array (nondim.) real :: ColHt_chg ! The change in column thickness in m. @@ -1342,7 +1597,7 @@ subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & if (dColHt_dKd < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres*dColHt_dKd endif -end subroutine find_PE_chg +end subroutine find_PE_chg_orig !> Copies the ePBL active mixed layer depth into MLD subroutine energetic_PBL_get_MLD(CS, MLD, G) @@ -1430,29 +1685,47 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", & units="nondim", default=omega_frac_dflt) call get_param(param_file, mod, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, & - "A ratio relating the efficiency with which convectively\n"//& - "released energy is converted to a turbulent velocity,\n"//& - "relative to mechanically forced TKE. Making this larger\n"//& - "increases the BL diffusivity", & - "units=nondim", default=1.0) + "A ratio relating the efficiency with which convectively \n"//& + "released energy is converted to a turbulent velocity, \n"//& + "relative to mechanically forced TKE. Making this larger \n"//& + "increases the BL diffusivity", units="nondim", default=1.0) call get_param(param_file, mod, "VSTAR_SCALE_FACTOR", CS%vstar_scale_fac, & "An overall nondimensional scaling factor for v*. \n"//& "Making this larger decreases the PBL diffusivity.", & - "units=nondim", default=1.0) + units="nondim", default=1.0) call get_param(param_file, mod, "EKMAN_SCALE_COEF", CS%Ekman_scale_coef, & - "A nondimensional scaling factor controlling the inhibition\n"//& - "of the diffusive length scale by rotation. Making this larger\n"//& - "decreases the PBL diffusivity.", & - "units=nondim", default=1.0) + "A nondimensional scaling factor controlling the inhibition \n"//& + "of the diffusive length scale by rotation. Making this larger \n"//& + "decreases the PBL diffusivity.", units="nondim", default=1.0) call get_param(param_file, mod, "USE_MLD_ITERATION", CS%USE_MLD_ITERATION, & - "A logical that determines whether or not to use the\n"//& - "MLD to set the EPBL length scale.", default=.false.) - - + "A logical that specifies whether or not to use the \n"//& + "distance to the bottom of the actively turblent boundary \n"//& + "layer to help set the EPBL length scale.", default=.false.) + call get_param(param_file, mod, "EPBL_MLD_TOLERANCE", CS%MLD_tol, & + "The tolerance for the iteratively determined mixed \n"//& + "layer depth. This is only used with USE_MLD_ITERATION.", & + units="meter", default=1.0) + call get_param(param_file, mod, "EPBL_MIN_MIX_LEN", CS%min_mix_len, & + "The minimum mixing length scale that will be used \n"//& + "by ePBL. The default (0) does not set a minimum.", & + units="meter", default=0.0) + call get_param(param_file, mod, "EPBL_ORIGINAL_PE_CALC", CS%orig_PE_calc, & + "If true, the ePBL code uses the original form of the \n"//& + "potential energy change code. Otherwise, the newer \n"//& + "version that can work with successive increments to the \n"//& + "diffusivity in upward or downward passes is used.", default=.true.) + call get_param(param_file, mod, "EPBL_TRANSITION_SCALE", CS%transLay_scale, & + "A scale for the mixing length in the transition layer \n"//& + "at the edge of the boundary layer as a fraction of the \n"//& + "boundary layer thickness. The default is 0, but a \n"//& + "value of 0.1 might be better justified by observations.", & + units="nondim", default=0.0) ! This gives a minimum decay scale that is typically much less than Angstrom. CS%ustar_min = 2e-4*CS%omega*(GV%Angstrom_z + GV%H_to_m*GV%H_subroundoff) - ! NOTE from AJA: The above parameter is not logged? + call log_param(param_file, mod, "EPBL_USTAR_MIN", CS%ustar_min, & + "The (tiny) minimum friction velocity used within the \n"//& + "ePBL code, derived from OMEGA and ANGSTROM.", units="meter second-1") CS%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, & Time, 'Surface mixed layer depth', 'meter') @@ -1495,17 +1768,16 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) CS%TKE_diagnostics = .true. endif - if (max(CS%id_Mixing_Length,CS%id_Velocity_Scale)>0) then - call safe_alloc_alloc(CS%Velocity_Scale,isd,ied,jsd,jed,SZK_(GV)+1) - call safe_alloc_alloc(CS%Mixing_Length,isd,ied,jsd,jed,SZK_(GV)+1) - CS%Velocity_Scale(:,:,:)=0.0 - CS%Mixing_Length(:,:,:)=0.0 - CS%mixing_diagnostics=.true. + if ((CS%id_Mixing_Length>0) .or. (CS%id_Velocity_Scale>0)) then + call safe_alloc_alloc(CS%Velocity_Scale,isd,ied,jsd,jed,GV%ke+1) + call safe_alloc_alloc(CS%Mixing_Length,isd,ied,jsd,jed,GV%ke+1) + CS%Velocity_Scale(:,:,:) = 0.0 + CS%Mixing_Length(:,:,:) = 0.0 + CS%mixing_diagnostics = .true. endif call safe_alloc_alloc(CS%ML_depth, isd, ied, jsd, jed) call safe_alloc_alloc(CS%ML_depth2, isd, ied, jsd, jed) - end subroutine energetic_PBL_init subroutine energetic_PBL_end(CS) @@ -1514,7 +1786,7 @@ subroutine energetic_PBL_end(CS) if (.not.associated(CS)) return if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) - if (allocated(CS%ML_depth2)) deallocate(CS%ML_depth2) + if (allocated(CS%ML_depth2)) deallocate(CS%ML_depth2) if (allocated(CS%diag_TKE_wind)) deallocate(CS%diag_TKE_wind) if (allocated(CS%diag_TKE_MKE)) deallocate(CS%diag_TKE_MKE) if (allocated(CS%diag_TKE_conv)) deallocate(CS%diag_TKE_conv) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index ffa9536c86..bac6c0eab8 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -636,9 +636,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G !Calculate tendencies (i.e., field changes at dt) from the sources / sinks ! - call generic_tracer_source(tv%T,tv%S,sosga,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,& + call generic_tracer_source(tv%T,tv%S,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,& G%areaT,get_diag_time_end(CS%diag),& - optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band) + optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, sosga=sosga) ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes ! usually in ALE mode @@ -910,10 +910,9 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS) call generic_tracer_coupler_set(state%tr_fields,& ST=state%SST,& SS=state%SSS,& - sosga=sosga, & rho=rho0,& !nnz: required for MOM5 and previous versions. ilb=G%isd, jlb=G%jsd,& - tau=1,model_time=get_diag_time_end(CS%diag)) + tau=1,sosga=sosga,model_time=get_diag_time_end(CS%diag)) !Output diagnostics via diag_manager for all tracers in this module ! if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&