From df0eaf0c608d40c8ef45adc91991490bfc21111e Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Wed, 19 Jul 2023 21:15:57 -0400 Subject: [PATCH 01/14] Fms2 io read3d slice (#399) * Restore functionality for reading slices from 3d volumes in MOM_io - The recent MOM_io modifications in support of FMS2_io accidentally removed support for reading on-grid data (same horizontal grid as model) k-slices. This is needed in some configurations in the model state initialization. * Add FMS1 interfaces * Additional patches to enable reading ongrid state initialization data - read local 3d volume rather than attempting to slice ongrid data vertically. - Related bugfixes in MOM_io --- config_src/infra/FMS1/MOM_coms_infra.F90 | 24 ++- config_src/infra/FMS1/MOM_io_infra.F90 | 41 +++++- config_src/infra/FMS2/MOM_coms_infra.F90 | 23 ++- config_src/infra/FMS2/MOM_io_infra.F90 | 70 ++++++++- src/framework/MOM_horizontal_regridding.F90 | 20 ++- src/framework/MOM_io.F90 | 154 +++++++++++++++++++- 6 files changed, 322 insertions(+), 10 deletions(-) diff --git a/config_src/infra/FMS1/MOM_coms_infra.F90 b/config_src/infra/FMS1/MOM_coms_infra.F90 index 939161875e..13f8006184 100644 --- a/config_src/infra/FMS1/MOM_coms_infra.F90 +++ b/config_src/infra/FMS1/MOM_coms_infra.F90 @@ -25,7 +25,7 @@ module MOM_coms_infra !> Communicate an array, string or scalar from one PE to others interface broadcast module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D - module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D + module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D end interface broadcast !> Compute a checksum for a field distributed over a PE list. If no PE list is @@ -260,6 +260,28 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking) end subroutine broadcast_real2D + +!> Communicate a 3-D array of reals from one PE to others +subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking) + real, dimension(:,:,:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The total number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real3D + ! field_chksum wrappers !> Compute a checksum for a field distributed over a PE list. If no PE list is diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index c0ccfcbcc8..e37e5db3cb 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -57,7 +57,7 @@ module MOM_io_infra !> Read a data field from a file interface read_field module procedure read_field_4d - module procedure read_field_3d + module procedure read_field_3d, read_field_3d_region module procedure read_field_2d, read_field_2d_region module procedure read_field_1d, read_field_1d_int module procedure read_field_0d, read_field_0d_int @@ -696,6 +696,45 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, & endif ; endif end subroutine read_field_3d +!> This routine uses the fms_io subroutine read_data to read a region from a distributed or +!! global 3-D data field named "fieldname" from file "filename". +subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, & + no_domain, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data + !! should be read + integer, dimension(:), intent(in) :: start !< The starting index to read in each of 4 + !! dimensions. For this 3-d read, the + !! 4th values are always 1. + integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4 + !! dimensions. For this 3-d read, the + !! 4th values are always 1. + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + logical, optional, intent(in) :: no_domain !< If present and true, this variable does not + !! use domain decomposion. + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + if (present(MOM_Domain)) then + call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, & + no_domain=no_domain) + else + call read_data(filename, fieldname, data, start, nread, no_domain=no_domain) + endif + + if (present(scale)) then ; if (scale /= 1.0) then + if (present(MOM_Domain)) then + call rescale_comp_data(MOM_Domain, data, scale) + else + ! Dangerously rescale the whole array + data(:,:,:) = scale*data(:,:,:) + endif + endif ; endif +end subroutine read_field_3d_region + + !> This routine uses the fms_io subroutine read_data to read a distributed !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. diff --git a/config_src/infra/FMS2/MOM_coms_infra.F90 b/config_src/infra/FMS2/MOM_coms_infra.F90 index 939161875e..cf9a724734 100644 --- a/config_src/infra/FMS2/MOM_coms_infra.F90 +++ b/config_src/infra/FMS2/MOM_coms_infra.F90 @@ -25,7 +25,7 @@ module MOM_coms_infra !> Communicate an array, string or scalar from one PE to others interface broadcast module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D - module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D + module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D end interface broadcast !> Compute a checksum for a field distributed over a PE list. If no PE list is @@ -260,6 +260,27 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking) end subroutine broadcast_real2D +!> Communicate a 3-D array of reals from one PE to others +subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking) + real, dimension(:,:,:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The total number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real3D + ! field_chksum wrappers !> Compute a checksum for a field distributed over a PE list. If no PE list is diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 99d0ac3345..a43b4e9344 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -66,7 +66,7 @@ module MOM_io_infra !> Read a data field from a file interface read_field module procedure read_field_4d - module procedure read_field_3d + module procedure read_field_3d, read_field_3d_region module procedure read_field_2d, read_field_2d_region module procedure read_field_1d, read_field_1d_int module procedure read_field_0d, read_field_0d_int @@ -1030,6 +1030,74 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, & end subroutine read_field_3d +!> This routine uses the fms_io subroutine read_data to read a region from a distributed or +!! global 3-D data field named "fieldname" from file "filename". +subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, & + no_domain, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data + !! should be read + integer, dimension(:), intent(in) :: start !< The starting index to read in each of 3 + !! dimensions. For this 3-d read, the + !! 4th value is always 1. + integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4 + !! dimensions. For this 3-d read, the + !! 4th values are always 1. + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + logical, optional, intent(in) :: no_domain !< If present and true, this variable does not + !! use domain decomposion. + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + ! Local variables + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file + type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object + character(len=96) :: var_to_read ! Name of variable to read from the netcdf file + logical :: success ! True if the file was successfully opened + + if (present(MOM_Domain)) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain) + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file and prepare to read it. + call prepare_to_read_var(fileobj_DD, fieldname, "read_field_2d_region: ", & + filename, var_to_read) + + ! Read the data. + call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3)) + + ! Close the file-set. + if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD) + else + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "read_field_2d_region: ", filename, var_to_read) + + ! Read the data. + call fms2_read_data(fileobj, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3)) + + ! Close the file-set. + if (check_if_open(fileobj)) call fms2_close_file(fileobj) + endif + + if (present(scale)) then ; if (scale /= 1.0) then + if (present(MOM_Domain)) then + call rescale_comp_data(MOM_Domain, data, scale) + else + ! Dangerously rescale the whole array + data(:,:,:) = scale*data(:,:,:) + endif + endif ; endif + +end subroutine read_field_3d_region + !> This routine uses the fms_io subroutine read_data to read a distributed !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 883653d715..34d0b73cb9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -309,6 +309,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its !! native horizontal grid, with units that change !! as the input data is interpreted [a] then [A ~> a] + real, dimension(:,:,:), allocatable :: tr_in_full !< A 3-d array for holding input data on the + !! model horizontal grid, with units that change + !! as the input data is interpreted [a] then [A ~> a] real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles !! with units that change as the input data is !! interpreted [a] then [A ~> a] @@ -448,6 +451,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr if (is_ongrid) then allocate(tr_in(is:ie,js:je), source=0.0) + allocate(tr_in_full(is:ie,js:je,kd), source=0.0) allocate(mask_in(is:ie,js:je), source=0.0) else call horizontal_interp_init() @@ -470,14 +474,19 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr ! Loop through each data level and interpolate to model grid. ! After interpolating, fill in points which will be needed to define the layers. + + if (is_ongrid) then + start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = 1 + count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = kd ; start(4) = 1 ; count(4) = 1 + call MOM_read_data(trim(filename), trim(varnam), tr_in_full, start, count, G%Domain) + endif + do k=1,kd mask_in(:,:) = 0.0 tr_out(:,:) = 0.0 if (is_ongrid) then - start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = k - count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = 1 ; start(4) = 1 ; count(4) = 1 - call MOM_read_data(trim(filename), trim(varnam), tr_in, start, count, G%Domain) + tr_in(is:ie,js:je) = tr_in_full(is:ie,js:je,k) do j=js,je do i=is,ie if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then @@ -594,7 +603,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr enddo ! kd - deallocate(lon_in, lat_in) + if (allocated(lat_inp)) deallocate(lat_inp) + deallocate(tr_in) + if (allocated(tr_inp)) deallocate(tr_inp) + if (allocated(tr_in_full)) deallocate(tr_in_full) end subroutine horiz_interp_and_extrap_tracer_record diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index bebce6f502..220a7d6bcf 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -100,6 +100,7 @@ module MOM_io module procedure MOM_read_data_2d module procedure MOM_read_data_2d_region module procedure MOM_read_data_3d + module procedure MOM_read_data_3d_region module procedure MOM_read_data_4d end interface MOM_read_data @@ -137,7 +138,7 @@ module MOM_io interface read_variable module procedure read_variable_0d, read_variable_0d_int module procedure read_variable_1d, read_variable_1d_int - module procedure read_variable_2d + module procedure read_variable_2d, read_variable_3d end interface read_variable !> Read a global or variable attribute from a named netCDF file using netCDF calls @@ -1161,7 +1162,7 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) allocate(field_nread(field_ndims)) field_nread(:2) = field_shape(:2) field_nread(3:) = 1 - if (present(nread)) field_shape(:2) = nread(:2) + if (present(nread)) field_nread(:2) = nread(:2) rc = nf90_get_var(ncid, varid, var, field_start, field_nread) @@ -1182,6 +1183,119 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) call broadcast(var, size(var), blocking=.true.) end subroutine read_variable_2d + +subroutine read_variable_3d(filename, varname, var, start, nread, ncid_in) + character(len=*), intent(in) :: filename !< Name of file to be read + character(len=*), intent(in) :: varname !< Name of variable to be read + real, intent(out) :: var(:,:,:) !< Output array of variable [arbitrary] + integer, optional, intent(in) :: start(:) !< Starting index on each axis. + integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis + integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file. + !! If absent, the file is opened and closed within this routine. + + integer :: ncid, varid + integer :: field_ndims, dim_len + integer, allocatable :: field_dimids(:), field_shape(:) + integer, allocatable :: field_start(:), field_nread(:) + integer :: i, rc + character(len=*), parameter :: hdr = "read_variable_3d: " + + ! Validate shape of start and nread + if (present(start)) then + if (size(start) < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " start must have at least two dimensions.") + endif + + if (present(nread)) then + if (size(nread) < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " nread must have at least two dimensions.") + + if (any(nread(3:) > 1)) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " nread may only read a single level in higher dimensions.") + endif + + ! Since start and nread may be reshaped, we cannot rely on netCDF to ensure + ! that their lengths are equivalent, and must do it here. + if (present(start) .and. present(nread)) then + if (size(start) /= size(nread)) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " start and nread must have the same length.") + endif + + ! Open and read `varname` from `filename` + if (is_root_pe()) then + if (present(ncid_in)) then + ncid = ncid_in + else + call open_file_to_Read(filename, ncid) + endif + + call get_varid(varname, ncid, filename, varid, match_case=.false.) + if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//& + " in "//trim(filename)) + + ! Query for the dimensionality of the input field + rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + ": Difficulties reading "//trim(varname)//" from "//trim(filename)) + + ! Confirm that field is at least 2d + if (field_ndims < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) // " " // & + trim(varname) // " from " // trim(filename) // " is not a 2D field.") + + ! If start and nread are present, then reshape them to match field dims + if (present(start) .or. present(nread)) then + allocate(field_shape(field_ndims)) + allocate(field_dimids(field_ndims)) + + rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + ": Difficulties reading "//trim(varname)//" from "//trim(filename)) + + do i = 1, field_ndims + rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len) + if (rc /= NF90_NOERR) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // ": Difficulties reading dimensions from " // trim(filename)) + field_shape(i) = dim_len + enddo + + ! Reshape start(:) and nreads(:) in case ranks differ + allocate(field_start(field_ndims)) + field_start(:) = 1 + if (present(start)) then + dim_len = min(size(start), size(field_start)) + field_start(:dim_len) = start(:dim_len) + endif + + allocate(field_nread(field_ndims)) + field_nread(:3) = field_shape(:3) + !field_nread(3:) = 1 + if (present(nread)) field_nread(:3) = nread(:3) + + rc = nf90_get_var(ncid, varid, var, field_start, field_nread) + + deallocate(field_start) + deallocate(field_nread) + deallocate(field_shape) + deallocate(field_dimids) + else + rc = nf90_get_var(ncid, varid, var) + endif + + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + " Difficulties reading "//trim(varname)//" from "//trim(filename)) + + if (.not.present(ncid_in)) call close_file_to_read(ncid, filename) + endif + + call broadcast(var, size(var), blocking=.true.) +end subroutine read_variable_3d + !> Read a character-string global or variable attribute subroutine read_attribute_str(filename, attname, att_val, varname, found, all_read, ncid_in) character(len=*), intent(in) :: filename !< Name of the file to read @@ -2198,6 +2312,42 @@ subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & endif end subroutine MOM_read_data_3d +!> Read a 3d region array from file using infrastructure I/O. +subroutine MOM_read_data_3d_region(filename, fieldname, data, start, nread, MOM_domain, & + no_domain, scale, turns) + character(len=*), intent(in) :: filename !< Input filename + character(len=*), intent(in) :: fieldname !< Field variable name + real, dimension(:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a] + integer, dimension(:), intent(in) :: start !< Starting index for each axis. + integer, dimension(:), intent(in) :: nread !< Number of values to read along each axis. + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition + logical, optional, intent(in) :: no_domain !< If true, field does not use + !! domain decomposion. + real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by + !! before it is returned to convert from the units in the file + !! to the internal units for this variable [A a-1 ~> 1] + integer, optional, intent(in) :: turns !< Number of quarter turns from + !! input to model grid + + integer :: qturns ! Number of quarter turns + real, allocatable :: data_in(:,:,:) ! Field array on the input grid in arbitrary units [A ~> a] + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call read_field(filename, fieldname, data, start, nread, & + MOM_Domain=MOM_Domain, no_domain=no_domain, scale=scale & + ) + else + call allocate_rotated_array(data, [1,1,1], -qturns, data_in) + call read_field(filename, fieldname, data_in, start, nread, & + MOM_Domain=MOM_Domain%domain_in, no_domain=no_domain, scale=scale & + ) + call rotate_array(data_in, qturns, data) + deallocate(data_in) + endif +end subroutine MOM_read_data_3d_region !> Read a 4d array from file using infrastructure I/O. subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & From 25feaf2ce017c087a06980ee8106f17f1ea605f2 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 18 Jul 2023 09:44:17 -0400 Subject: [PATCH 02/14] Fix logic for setting KV_ML_INVZ2 from KVML - We were reading KV_ML_INVZ2 without logging, then checking for KVML and finally logging based on a combination of the two. This had the side affect that we get warnings about not using KVML even if KVML was not present. - The fix checks for KVML first, and then changes the default so that when KVML=1e-4 is replaced by KV_ML_INVZ2=1e-4 we end up with no warnings and KVML can be obsoleted safely. Note: this commit alone does not remove all warnings from the MOM6-examples suite because we still need to fix the MOM_input that still use KVML - KVML needs to be unscaled since it is the default for KV_ML_INVZ2 - tc3 used KVML and has been corrected. --- .testing/tc3/MOM_input | 8 +++--- .../vertical/MOM_vert_friction.F90 | 25 +++++++------------ 2 files changed, 13 insertions(+), 20 deletions(-) diff --git a/.testing/tc3/MOM_input b/.testing/tc3/MOM_input index a034960d1e..6963feee98 100644 --- a/.testing/tc3/MOM_input +++ b/.testing/tc3/MOM_input @@ -283,10 +283,10 @@ HMIX_FIXED = 20.0 ! [m] KV = 1.0E-04 ! [m2 s-1] ! The background kinematic viscosity in the interior. ! The molecular value, ~1e-6 m2 s-1, may be used. -KVML = 0.01 ! [m2 s-1] default = 1.0E-04 - ! The kinematic viscosity in the mixed layer. A typical - ! value is ~1e-2 m2 s-1. KVML is not used if - ! BULKMIXEDLAYER is true. The default is set by KV. +KV_ML_INVZ2 = 0.01 ! [m2 s-1] default = 0.0 + ! An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, with + ! the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the distance + ! from the surface, to allow for finite wind stresses to be transmitted through. HBBL = 10.0 ! [m] ! The thickness of a bottom boundary layer with a ! viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 0ab283a90e..496012c3d9 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -2156,6 +2156,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & # include "version_variable.h" character(len=40) :: mdl = "MOM_vert_friction" ! This module's name. character(len=40) :: thickness_units + real :: Kv_mks ! KVML in MKS if (associated(CS)) then call MOM_error(WARNING, "vertvisc_init called with an associated "// & @@ -2339,30 +2340,22 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & CS%Kvml_invZ2 = 0.0 if (GV%nkml < 1) then - call get_param(param_file, mdl, "KV_ML_INVZ2", CS%Kvml_invZ2, & - "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//& - "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//& - "distance from the surface, to allow for finite wind stresses to be "//& - "transmitted through infinitesimally thin surface layers. This is an "//& - "older option for numerical convenience without a strong physical basis, "//& - "and its use is now discouraged.", & - units="m2 s-1", default=-1.0, scale=US%m2_s_to_Z2_T, do_not_log=.true.) - if (CS%Kvml_invZ2 < 0.0) then - call get_param(param_file, mdl, "KVML", CS%Kvml_invZ2, & + call get_param(param_file, mdl, "KVML", Kv_mks, & "The scale for an extra kinematic viscosity in the mixed layer", & - units="m2 s-1", default=0.0, scale=US%m2_s_to_Z2_T, do_not_log=.true.) - if (CS%Kvml_invZ2 >= 0.0) & - call MOM_error(WARNING, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.") + units="m2 s-1", default=-1.0, do_not_log=.true.) + if (Kv_mks >= 0.0) then + call MOM_error(WARNING, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.") + else + Kv_mks = 0.0 endif - if (CS%Kvml_invZ2 < 0.0) CS%Kvml_invZ2 = 0.0 - call log_param(param_file, mdl, "KV_ML_INVZ2", CS%Kvml_invZ2, & + call get_param(param_file, mdl, "KV_ML_INVZ2", CS%Kvml_invZ2, & "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//& "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//& "distance from the surface, to allow for finite wind stresses to be "//& "transmitted through infinitesimally thin surface layers. This is an "//& "older option for numerical convenience without a strong physical basis, "//& "and its use is now discouraged.", & - units="m2 s-1", default=0.0, unscale=US%Z2_T_to_m2_s) + units="m2 s-1", default=Kv_mks, scale=US%m2_s_to_Z2_T) endif if (.not.CS%bottomdraglaw) then From a5129ca3772d461b8bcb12b1aa23b2b942a734bd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 30 Jun 2023 10:07:18 -0400 Subject: [PATCH 03/14] +Use RHO_PGF_REF for the pressure gradient forces Use the new runtime parameter RHO_PGF_REF instead of RHO_0 to set the reference density that is subtracted off from the other densities when calculating the finite volume pressure gradient forces. Although the answers are mathematically equivalent for any value of this parameter, a judicious choice can reduce the impacts of roundoff errors by about 2 orders of magnitude. By default, RHO_PGF_REF is set to RHO_0, and all answers are bitwise identical. However, there is a new runtime parameter that appears in many of the MOM_parameter_doc.all files. --- src/core/MOM_PressureForce_FV.F90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 14c9b2e6dc..281623ae84 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -519,7 +519,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS h_neglect = GV%H_subroundoff - dz_neglect = GV%H_subroundoff * GV%H_to_Z + dz_neglect = GV%dZ_subroundoff I_Rho0 = 1.0 / GV%Rho0 G_Rho0 = GV%g_Earth / GV%Rho0 rho_ref = CS%Rho0 @@ -827,12 +827,11 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS mdl = "MOM_PressureForce_FV" call log_version(param_file, mdl, version, "") - call get_param(param_file, mdl, "RHO_0", CS%Rho0, & - "The mean ocean density used with BOUSSINESQ true to "//& - "calculate accelerations and the mass for conservation "//& - "properties, or with BOUSSINSEQ false to convert some "//& - "parameters from vertical units of m to kg m-2.", & - units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) + call get_param(param_file, mdl, "RHO_PGF_REF", CS%Rho0, & + "The reference density that is subtracted off when calculating pressure "//& + "gradient forces. Its inverse is subtracted off of specific volumes when "//& + "in non-Boussinesq mode. The default is RHO_0.", & + units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, & From bb71c346ac4dc1d9c87f2d4e265ecfc1bed90403 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 20 Jul 2023 09:18:21 -0400 Subject: [PATCH 04/14] Converted warning about depth_list to a note The message that a file is being created was issued as a WARNING when we all agree it should really be a NOTE. Depth_list.nc is read if it is present to avoid recomputing a sorted list, but the absence of the file is not an error and does not warrant a warning. Changes: - Changed WARNING to NOTE. - Removed MOM_mesg from imports since it wasn't being used. --- src/diagnostics/MOM_sum_output.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index fd957d0a44..aae32a7862 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -8,7 +8,7 @@ module MOM_sum_output use MOM_coms, only : sum_across_PEs, PE_here, root_PE, num_PEs, max_across_PEs, field_chksum use MOM_coms, only : reproducing_sum, reproducing_sum_EFP, EFP_to_real, real_to_EFP use MOM_coms, only : EFP_type, operator(+), operator(-), assignment(=), EFP_sum_across_PEs -use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, MOM_mesg +use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type @@ -1077,7 +1077,7 @@ subroutine depth_list_setup(G, GV, US, DL, CS) valid_DL_read = .true. ! Otherwise there would have been a fatal error. endif else - if (is_root_pe()) call MOM_error(WARNING, "depth_list_setup: "// & + if (is_root_pe()) call MOM_error(NOTE, "depth_list_setup: "// & trim(CS%depth_list_file)//" does not exist. Creating a new file.") valid_DL_read = .false. endif From b9c7c8689707b0a4e6c9d1d95fa9c57fcd5202f6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 19 Jul 2023 09:32:28 -0400 Subject: [PATCH 05/14] Correct diagnostic coordinate interpolation scheme The interpolation scheme for state-dependent diagnostic coordinates was incorrectly registering as the same parameter as the main model. This meant it was never possible to change the interpolation scheme from the default (which was not the same as the main model). Fix registers the generated parameter name which was always computed but not used. A typical example of the generated parameter is "DIAG_COORD_INTERP_SCHEME_RHO2". --- src/ALE/MOM_regridding.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index c5f5807f66..4cc60d16b2 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -263,7 +263,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m param_name = create_coord_param(param_prefix, "INTERP_SCHEME", param_suffix) string2 = 'PPM_H4' ! Default for diagnostics endif - call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, & + call get_param(param_file, mdl, param_name, string, & "This sets the interpolation scheme to use to "//& "determine the new grid. These parameters are "//& "only relevant when REGRIDDING_COORDINATE_MODE is "//& From ccd3ded7b6279f433115e16f183cf29a1c38063e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 20 Jul 2023 07:23:09 -0400 Subject: [PATCH 06/14] +(*)Fix wave_speed_init mono_N2_depth bug Fixed a bug in which wave_speed_init was effectively discarding any values of mono_N2_depth passed to it via the optional argument mono_N2_depth, but also changed the default value of RESOLN_N2_FILTER_DEPTH, which was previously being discarded, to disable the monotonization and replicate the previous results. There were also clarifying additions made to the description how to disable RESOLN_N2_FILTER_DEPTH. This will change some entries in MOM_parameter_doc files, and it will change solutions in cases that set RESOLN_N2_FILTER_DEPTH to a non-default value and have parameter settings that use the resolution function to scale their horizontal mixing. There are, however, no known active simulations where the answers are expected to change. --- src/diagnostics/MOM_wave_speed.F90 | 4 +++- src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index bb1b381c15..5757e25cd1 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -39,6 +39,7 @@ module MOM_wave_speed !! wave_speed() which can be overridden by optional arguments. real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of !! calculating the equivalent barotropic wave speed [Z ~> m]. + !! If this parameter is negative, this limiting does not occur. !! This parameter controls the default behavior of wave_speed() which !! can be overridden by optional arguments. real :: min_speed2 = 0. !< The minimum mode 1 internal wave speed squared [L2 T-2 ~> m2 s-2] @@ -1465,7 +1466,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de call log_version(mdl, version) call wave_speed_set_param(CS, use_ebt_mode=use_ebt_mode, mono_N2_column_fraction=mono_N2_column_fraction, & - better_speed_est=better_speed_est, min_speed=min_speed, wave_speed_tol=wave_speed_tol, & + mono_N2_depth=mono_N2_depth, better_speed_est=better_speed_est, & + min_speed=min_speed, wave_speed_tol=wave_speed_tol, & remap_answers_2018=remap_answers_2018, remap_answer_date=remap_answer_date, & c1_thresh=c1_thresh) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7d71a62e25..8f0aa02b12 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1239,8 +1239,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) in_use = .true. call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, & "The depth below which N2 is monotonized to avoid stratification "//& - "artifacts from altering the equivalent barotropic mode structure.",& - units="m", default=2000., scale=US%m_to_Z) + "artifacts from altering the equivalent barotropic mode structure. "//& + "This monotonzization is disabled if this parameter is negative.", & + units="m", default=-1.0, scale=US%m_to_Z) allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif From 6102be250c7ad59ffb3afe8816058a0d29f2cae7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 30 Jun 2023 10:35:09 -0400 Subject: [PATCH 07/14] *+Revise non-Boussinesq gprime expressions Revised the calculation of gprime and the coordinate densities (GV%Rlay) in fully non-Boussinesq mode to use the arithmetic mean of adjacent coordinate densities in the denominator of the expression for g_prime in place of RHO_0. Also use LIGHTEST_DENSITY in place of RHO_0 to specify the top-level coordinate density in certain coordinate modes. Also made corresponding changes to the fully non-Boussinesq APE calculation when CALCULATE_APE is true, and eliminated an incorrect calculation of the layer volumes in non-Boussinesq mode using the Boussinesq reference density that was never actually being used when CALCULATE_APE is false. This commit will change answers in some fully non-Boussinesq calculations, and an existing runtime parameter is used and logged in some new cases, changing the MOM_parameter_doc file in those cases. --- src/diagnostics/MOM_sum_output.F90 | 38 +++++--- .../MOM_coord_initialization.F90 | 94 +++++++++++++++---- 2 files changed, 101 insertions(+), 31 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index aae32a7862..fb95b79a91 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -510,24 +510,18 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci do k=1,nz ; vol_lay(k) = (US%m_to_L**2*GV%H_to_Z/GV%H_to_kg_m2)*mass_lay(k) ; enddo else tmp1(:,:,:) = 0.0 - if (CS%do_APE_calc) then - do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j) - enddo ; enddo ; enddo - mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) + do k=1,nz ; do j=js,je ; do i=is,ie + tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j) + enddo ; enddo ; enddo + mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) + if (CS%do_APE_calc) then call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) do k=1,nz ; do j=js,je ; do i=is,ie tmp1(i,j,k) = US%Z_to_m*US%L_to_m**2*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay) do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2 * vol_lay(k) ; enddo - else - do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j) - enddo ; enddo ; enddo - mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) - do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2*US%kg_m3_to_R * (mass_lay(k) / GV%Rho0) ; enddo endif endif ! Boussinesq @@ -643,7 +637,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci if (GV%Boussinesq) then do j=js,je ; do i=is,ie hbelow = 0.0 - do k=nz,1,-1 + do K=nz,1,-1 hbelow = hbelow + h(i,j,k) * GV%H_to_Z hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref)) hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref) @@ -652,14 +646,28 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci (hint * hint - hbot * hbot) enddo enddo ; enddo - else + elseif (GV%semi_Boussinesq) then do j=js,je ; do i=is,ie - do k=nz,1,-1 + do K=nz,1,-1 hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * & - (hint * hint - hbot * hbot) + (hint * hint - hbot * hbot) + enddo + enddo ; enddo + else + do j=js,je ; do i=is,ie + do K=nz,2,-1 + hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. + hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) + PE_pt(i,j,K) = (0.25 * PE_scale_factor * areaTm(i,j) * & + ((GV%Rlay(k)+GV%Rlay(k-1))*GV%g_prime(K))) * & + (hint * hint - hbot * hbot) enddo + hint = Z_0APE(1) + eta(i,j,1) ! eta and H_0 have opposite signs. + hbot = max(Z_0APE(1) - (G%bathyT(i,j) + G%Z_ref), 0.0) + PE_pt(i,j,1) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * & + (hint * hint - hbot * hbot) enddo ; enddo endif diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index 8af8cd3bc6..37c719209b 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -126,6 +126,7 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file) ! Local variables real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2]. real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. + real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3]. character(len=40) :: mdl = "set_coord_from_gprime" ! This subroutine's name. integer :: k, nz nz = GV%ke @@ -138,11 +139,20 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file) call get_param(param_file, mdl, "GINT", g_int, & "The reduced gravity across internal interfaces.", & units="m s-2", fail_if_missing=.true., scale=US%m_s_to_L_T**2*US%Z_to_m) + call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & + "The reference potential density used for layer 1.", & + units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) g_prime(1) = g_fs do k=2,nz ; g_prime(k) = g_int ; enddo - Rlay(1) = GV%Rho0 - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + Rlay(1) = Rlay_Ref + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + else + do k=2,nz + Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) + enddo + endif call callTree_leave(trim(mdl)//'()') @@ -184,9 +194,15 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file) enddo ! These statements set the interface reduced gravities. ! g_prime(1) = g_fs - do k=2,nz - g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) - enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz + g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) + enddo + else + do k=2,nz + g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) + enddo + endif call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_layer_density @@ -237,7 +253,13 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state) ! These statements set the layer densities. ! - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + else + do k=2,nz + Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) + enddo + endif call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_ref @@ -294,7 +316,15 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s g_prime(1) = g_fs do k=1,nz ; Pref(k) = P_Ref ; enddo call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, (/1,nz/) ) - do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz + g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) + enddo + else + do k=2,nz + g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) + enddo + endif call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_profile @@ -387,7 +417,15 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta do k=k_light-1,1,-1 Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2) enddo - do k=2,nz ; g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) ; enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz + g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) + enddo + else + do k=2,nz + g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) + enddo + endif call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_range @@ -429,7 +467,15 @@ subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file) call MOM_read_data(filename, coord_var, Rlay, scale=US%kg_m3_to_R) g_prime(1) = g_fs - do k=2,nz ; g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) ; enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz + g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) + enddo + else + do k=2,nz + g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) + enddo + endif do k=1,nz ; if (g_prime(k) <= 0.0) then call MOM_error(FATAL, "MOM_initialization set_coord_from_file: "//& "Zero or negative g_primes read from variable "//"Layer"//" in file "//& @@ -479,9 +525,15 @@ subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file) enddo ! These statements set the interface reduced gravities. g_prime(1) = g_fs - do k=2,nz - g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) - enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz + g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) + enddo + else + do k=2,nz + g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) + enddo + endif call callTree_leave(trim(mdl)//'()') end subroutine set_coord_linear @@ -498,6 +550,7 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters ! Local variables real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. + real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3]. character(len=40) :: mdl = "set_coord_to_none" ! This subroutine's name. integer :: k, nz nz = GV%ke @@ -507,11 +560,20 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file) call get_param(param_file, mdl, "GFS" , g_fs, & "The reduced gravity at the free surface.", units="m s-2", & default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) + call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & + "The reference potential density used for layer 1.", & + units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) g_prime(1) = g_fs do k=2,nz ; g_prime(k) = 0. ; enddo - Rlay(1) = GV%Rho0 - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + Rlay(1) = Rlay_Ref + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo + else + do k=2,nz + Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) + enddo + endif call callTree_leave(trim(mdl)//'()') @@ -522,8 +584,8 @@ end subroutine set_coord_to_none subroutine write_vertgrid_file(GV, US, param_file, directory) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - character(len=*), intent(in) :: directory !< The directory into which to place the file. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + character(len=*), intent(in) :: directory !< The directory into which to place the file. ! Local variables character(len=240) :: filepath type(vardesc) :: vars(2) From 2f6e86e342608f07f5cc3fb478a22cc6089075b3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 21 Jul 2023 12:00:14 -0400 Subject: [PATCH 08/14] *Non-Boussinesq revision of MOM_thickness_diffuse Refactored thickness_diffuse when in non-Boussinesq mode to avoid any dependencies on the Boussinesq reference density, and to translate the volume streamfunction into the mass streamfunction using an appropriately defined in-situ density averaged to the interfaces at velocity points. This form follows the suggestions of Appendix A.3.2 of Griffies and Greatbatch (Ocean Modelling, 2012) when in non-Boussinesq mode. Thickness_diffuse_full was also revised to work properly in non-Boussinesq mode (and not depend on the Boussinesq reference density) when no equation of state is used. As a part of these changes, the code now uses thickness-based streamfunctions and other thickness-based internal calculations in MOM_thickness_diffuse. For example, the overturning streamfunctions with this change are now in m3/s in Boussinesq mode, but kg/s in non-Boussinesq mode. These changes use a call to thickness_to_dz to set up a separate variable with the vertical distance across layers, and in non-Boussinesq mode they use tv%SpV_avg to estimate in situ densities. Additional debugging checksums were added to thickness_diffuse. The code changes are extensive with 15 new or renamed internal variables, and changes to the units of 9 other internal variables and 3 arguments to the private routine streamfn_solver. After this change, GV%Rho, GV%Z_to_H and GV%H_to_Z are no longer used in any non-Boussinesq calculations (12 such instances having been elimated). Because some calculations have to be redone with the separate thickness and dz variables, this will be more expensive than the original version. No public interfaces are changed, and all answers are bitwise identical in Boussinesq or semiBoussinesq mode, but they will change in non-Boussinesq mode when the isopycnal height diffusion parameterization is used. --- .../lateral/MOM_thickness_diffuse.F90 | 286 ++++++++++++------ 1 file changed, 195 insertions(+), 91 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index a7ff2f1c0a..8617795e16 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -14,7 +14,7 @@ module MOM_thickness_diffuse use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_io, only : MOM_read_data, slasher -use MOM_interface_heights, only : find_eta +use MOM_interface_heights, only : find_eta, thickness_to_dz use MOM_isopycnal_slopes, only : vert_fill_TS use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type @@ -439,7 +439,6 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif endif - !$OMP do do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo !$OMP do @@ -458,6 +457,12 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%debug) then call uvchksum("Kh_[uv]", Kh_u, Kh_v, G%HI, haloshift=0, & scale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.) + call uvchksum("Kh_[uv]_CFL", Kh_u_CFL, Kh_v_CFL, G%HI, haloshift=0, & + scale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.) + if (Resoln_scaled) then + call uvchksum("Res_fn_[uv]", VarMix%Res_fn_u, VarMix%Res_fn_v, G%HI, haloshift=0, & + scale=1.0, scalar_pair=.true.) + endif call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, G%HI, haloshift=0) call hchksum(h, "thickness_diffuse_1 h", G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m) @@ -628,14 +633,17 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. h_frac ! The fraction of the mass in the column above the bottom ! interface of a layer that is within a layer [nondim]. 0 m] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: & Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [nondim] hN2_y_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency - ! at v-points [L2 Z-1 T-2 ~> m s-2], used for calculating PE release + ! at v-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], + ! used for calculating the potential energy release real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [nondim] hN2_x_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency - ! at u-points [L2 Z-1 T-2 ~> m s-2], used for calculating PE release + ! at u-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], + ! used for calculating the potential energy release real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & pres, & ! The pressure at an interface [R L2 T-2 ~> Pa]. h_avail_rsum ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -670,8 +678,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: Work_v(SZI_(G),SZJB_(G)) ! The work done by the isopycnal height diffusion ! integrated over v-point water columns [R Z L4 T-3 ~> W] real :: Work_h ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2]. - real :: PE_release_h ! The amount of potential energy released by GM averaged over an h-cell [L4 Z-1 T-3 ~> m3 s-3] - ! The calculation is equal to h * S^2 * N^2 * kappa_GM. + real :: PE_release_h ! The amount of potential energy released by GM averaged over an h-cell + ! [R Z L2 T-3 ~> W m-2]. The calculation equals rho0 * h * S^2 * N^2 * kappa_GM. real :: I4dt ! 1 / 4 dt [T-1 ~> s-1]. real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A) ! and below (B) the interface times the grid spacing [R ~> kg m-3]. @@ -686,57 +694,71 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! [Z R ~> kg m-2]. real :: hg2A, hg2B, hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4]. real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2]. - real :: dzaL, dzaR ! Temporary thicknesses [Z ~> m] + real :: dzg2A, dzg2B ! Squares of geometric mean vertical layer extents [Z2 ~> m2]. + real :: dzaA, dzaB ! Arithmetic mean vertical layer extents [Z ~> m]. + real :: dzaL, dzaR ! Temporary vertical layer extents [Z ~> m] real :: wtA, wtB ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6] real :: wtL, wtR ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5] real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4]. real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4]. - real :: h_harm ! Harmonic mean layer thickness [H ~> m or kg m-2]. - real :: c2_h_u(SZIB_(G),SZK_(GV)+1) ! Wave speed squared divided by h at u-points [L2 Z-1 T-2 ~> m s-2]. - real :: c2_h_v(SZI_(G),SZK_(GV)+1) ! Wave speed squared divided by h at v-points [L2 Z-1 T-2 ~> m s-2]. - real :: hN2_u(SZIB_(G),SZK_(GV)+1) ! Thickness in m times N2 at interfaces above u-points [L2 Z-1 T-2 ~> m s-2]. - real :: hN2_v(SZI_(G),SZK_(GV)+1) ! Thickness in m times N2 at interfaces above v-points [L2 Z-1 T-2 ~> m s-2]. + real :: dz_harm ! Harmonic mean layer vertical extent [Z ~> m]. + real :: c2_dz_u(SZIB_(G),SZK_(GV)+1) ! Wave speed squared divided by dz at u-points times rescaling + ! factors from depths to thicknesses [H2 L2 Z-3 T-2 ~> m s-2 or kg m-2 s-2] + real :: c2_dz_v(SZI_(G),SZK_(GV)+1) ! Wave speed squared divided by dz at v-points times rescaling + ! factors from depths to thicknesses [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2] + real :: dzN2_u(SZIB_(G),SZK_(GV)+1) ! Vertical extent times N2 at interfaces above u-points times + ! rescaling factors from vertical to horizontal distances [L2 Z-1 T-2 ~> m s-2] + real :: dzN2_v(SZI_(G),SZK_(GV)+1) ! Vertical extent times N2 at interfaces above v-points times + ! rescaling factors from vertical to horizontal distances [L2 Z-1 T-2 ~> m s-2] real :: Sfn_est ! A preliminary estimate (before limiting) of the overturning - ! streamfunction [Z L2 T-1 ~> m3 s-1]. - real :: Sfn_unlim_u(SZIB_(G),SZK_(GV)+1) ! Streamfunction for u-points [Z L2 T-1 ~> m3 s-1]. - real :: Sfn_unlim_v(SZI_(G),SZK_(GV)+1) ! Streamfunction for v-points [Z L2 T-1 ~> m3 s-1]. + ! streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1]. + real :: Sfn_unlim_u(SZIB_(G),SZK_(GV)+1) ! Volume streamfunction for u-points [Z L2 T-1 ~> m3 s-1] + real :: Sfn_unlim_v(SZI_(G),SZK_(GV)+1) ! Volume streamfunction for v-points [Z L2 T-1 ~> m3 s-1] real :: slope2_Ratio_u(SZIB_(G),SZK_(GV)+1) ! The ratio of the slope squared to slope_max squared [nondim] real :: slope2_Ratio_v(SZI_(G),SZK_(GV)+1) ! The ratio of the slope squared to slope_max squared [nondim] real :: Sfn_in_h ! The overturning streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] (note that ! the units are different from other Sfn vars). - real :: Sfn_safe ! The streamfunction that goes linearly back to 0 at the surface. This is a - ! good thing to use when the slope is so large as to be meaningless [Z L2 T-1 ~> m3 s-1]. + real :: Sfn_safe ! The streamfunction that goes linearly back to 0 at the surface + ! [H L2 T-1 ~> m3 s-1 or kg s-1]. This is a good value to use when the + ! slope is so large as to be meaningless, usually due to weak stratification. real :: Slope ! The slope of density surfaces, calculated in a way that is always ! between -1 and 1 after undoing dimensional scaling, [Z L-1 ~> nondim] real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8]. real :: I_slope_max2 ! The inverse of slope_max squared [L2 Z-2 ~> nondim]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. + real :: hn_2 ! Half of h_neglect [H ~> m or kg m-2]. real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]. real :: dz_neglect ! A thickness [Z ~> m], that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. + real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: G_scale ! The gravitational acceleration times a unit conversion ! factor [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. logical :: find_work ! If true, find the change in energy due to the fluxes. integer :: nk_linear ! The number of layers over which the streamfunction goes to 0. real :: G_rho0 ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2]. + real :: Rho_avg ! The in situ density averaged to an interface [R ~> kg m-3] real :: N2_floor ! A floor for N2 to avoid degeneracy in the elliptic solver ! times unit conversion factors [L2 Z-2 T-2 ~> s-2] + real :: N2_unlim ! An unlimited estimate of the buoyancy frequency + ! times unit conversion factors [L2 Z-2 T-2 ~> s-2] real :: Tl(5) ! copy of T in local stencil [C ~> degC] real :: mn_T ! mean of T in local stencil [C ~> degC] real :: mn_T2 ! mean of T**2 in local stencil [C2 ~> degC2] real :: hl(5) ! Copy of local stencil of H [H ~> m] real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1] + real :: Z_to_H ! A conversion factor from heights to thicknesses, perhaps based on + ! a spatially variable local density [H Z-1 ~> nondim or kg m-3] real :: Tsgs2(SZI_(G),SZJ_(G),SZK_(GV)) ! Sub-grid temperature variance [C2 ~> degC2] real :: diag_sfn_x(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! Diagnostic of the x-face streamfunction ! [H L2 T-1 ~> m3 s-1 or kg s-1] real :: diag_sfn_unlim_x(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! Diagnostic of the x-face streamfunction before - ! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1] + ! applying limiters [Z L2 T-1 ~> m3 s-1] real :: diag_sfn_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction ! [H L2 T-1 ~> m3 s-1 or kg s-1] real :: diag_sfn_unlim_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction before - ! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1] + ! applying limiters [Z L2 T-1 ~> m3 s-1] logical :: present_slope_x, present_slope_y, calc_derivatives integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of ! state calculations at u-points. @@ -753,10 +775,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV I_slope_max2 = 1.0 / (CS%slope_max**2) G_scale = GV%g_Earth * GV%H_to_Z - h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 - dz_neglect = GV%H_subroundoff*GV%H_to_Z + h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 ; hn_2 = 0.5*h_neglect + dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect**2 G_rho0 = GV%g_Earth / GV%Rho0 - N2_floor = CS%N2_floor*US%Z_to_L**2 + N2_floor = CS%N2_floor * US%Z_to_L**2 use_EOS = associated(tv%eqn_of_state) present_slope_x = PRESENT(slope_x) @@ -779,6 +801,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo, larger_h_denom=.true.) endif + ! Rescale the thicknesses, perhaps using the specific volume. + call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) + if (CS%use_FGNV_streamfn .and. .not. associated(cg1)) call MOM_error(FATAL, & "cg1 must be associated when using FGNV streamfunction.") @@ -824,20 +849,21 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV EOSdom_h1(:) = EOS_domain(G%HI, halo=1) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & - !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & - !$OMP h_neglect2,int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, & - !$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1,diag_sfn_x, & - !$OMP diag_sfn_unlim_x,N2_floor,EOSdom_u,EOSdom_h1,use_stanley,Tsgs2, & - !$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) & - !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & + !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz,dz_neglect,dz_neglect2, & + !$OMP h_neglect2,hn_2,I_slope_max2,int_slope_u,KH_u,uhtot, & + !$OMP h_frac,h_avail_rsum,uhD,h_avail,Work_u,CS,slope_x,cg1, & + !$OMP diag_sfn_x,diag_sfn_unlim_x,N2_floor,EOSdom_u,EOSdom_h1, & + !$OMP use_stanley,Tsgs2,present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) & + !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u,G_scale, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & - !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & + !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h,N2_unlim, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & - !$OMP drdx,mag_grad2,Slope,slope2_Ratio_u,hN2_u, & - !$OMP Sfn_unlim_u,drdi_u,drdkDe_u,h_harm,c2_h_u, & + !$OMP dzg2A,dzg2B,dzaA,dzaB,dz_harm,Z_to_H, & + !$OMP drdx,mag_grad2,Slope,slope2_Ratio_u,dzN2_u, & + !$OMP Sfn_unlim_u,Rho_avg,drdi_u,drdkDe_u,c2_dz_u, & !$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives) do j=js,je - do I=is-1,ie ; hN2_u(I,1) = 0. ; hN2_u(I,nz+1) = 0. ; enddo + do I=is-1,ie ; dzN2_u(I,1) = 0. ; dzN2_u(I,nz+1) = 0. ; enddo do K=nz,2,-1 if (find_work .and. .not.(use_EOS)) then drdiA = 0.0 ; drdiB = 0.0 @@ -907,9 +933,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect if (GV%Boussinesq) then dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z - else + elseif (GV%semi_Boussinesq) then dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect + else + dzaL = 0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect + dzaR = 0.5*(dz(i+1,j,k-1) + dz(i+1,j,k)) + dz_neglect endif ! Use the harmonic mean thicknesses to weight the horizontal gradients. ! These unnormalized weights have been rearranged to minimize divisions. @@ -924,10 +953,23 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV haA = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect haB = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect - ! hN2_u is used with the FGNV streamfunction formulation - hN2_u(I,K) = (0.5 * GV%H_to_Z * ( hg2A / haA + hg2B / haB )) * & - max(drdz*G_rho0, N2_floor) + if (GV%Boussinesq) then + N2_unlim = drdz*G_rho0 + else + N2_unlim = (GV%g_Earth*GV%RZ_to_H) * & + ((wtL * drdkL + wtR * drdkR) / (haL*wtL + haR*wtR)) + endif + + dzg2A = dz(i,j,k-1)*dz(i+1,j,k-1) + dz_neglect2 + dzg2B = dz(i,j,k)*dz(i+1,j,k) + dz_neglect2 + dzaA = 0.5*(dz(i,j,k-1) + dz(i+1,j,k-1)) + dz_neglect + dzaB = 0.5*(dz(i,j,k) + dz(i+1,j,k)) + dz_neglect + ! dzN2_u is used with the FGNV streamfunction formulation + dzN2_u(I,K) = (0.5 * ( dzg2A / dzaA + dzg2B / dzaB )) * max(N2_unlim, N2_floor) + if (find_work .and. CS%GM_src_alt) & + hN2_x_PE(I,j,k) = (0.5 * ( hg2A / haA + hg2B / haB )) * max(N2_unlim, N2_floor) endif + if (present_slope_x) then Slope = slope_x(I,j,k) slope2_Ratio_u(I,K) = Slope**2 * I_slope_max2 @@ -958,11 +1000,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max) - hN2_x_PE(I,j,k) = hN2_u(I,K) if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope - ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. - Sfn_unlim_u(I,K) = -((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) + ! Estimate the streamfunction at each interface [H L2 T-1 ~> m3 s-1 or kg s-1]. + Sfn_unlim_u(I,K) = -(KH_u(I,j,K)*G%dy_Cu(I,j))*Slope ! Avoid moving dense water upslope from below the level of ! the bottom on the receiving side. @@ -992,10 +1033,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) - hN2_u(I,K) = GV%g_prime(K) + dzN2_u(I,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) - hN2_u(I,K) = N2_floor * dz_neglect + dzN2_u(I,K) = N2_floor * dz_neglect Sfn_unlim_u(I,K) = 0. endif ! if (k > nk_linear) if (CS%id_sfn_unlim_x>0) diag_sfn_unlim_x(I,j,K) = Sfn_unlim_u(I,K) @@ -1004,10 +1045,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%use_FGNV_streamfn) then do k=1,nz ; do I=is-1,ie ; if (G%OBCmaskCu(I,j)>0.) then - h_harm = max( h_neglect, & - 2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) ) - c2_h_u(I,k) = CS%FGNV_scale * & - ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (GV%H_to_Z*h_harm) + dz_harm = max( dz_neglect, & + 2. * dz(i,j,k) * dz(i+1,j,k) / ( ( dz(i,j,k) + dz(i+1,j,k) ) + dz_neglect ) ) + c2_dz_u(I,k) = CS%FGNV_scale * ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / dz_harm endif ; enddo ; enddo ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. @@ -1016,7 +1056,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do K=2,nz Sfn_unlim_u(I,K) = (1. + CS%FGNV_scale) * Sfn_unlim_u(I,K) enddo - call streamfn_solver(nz, c2_h_u(I,:), hN2_u(I,:), Sfn_unlim_u(I,:)) + call streamfn_solver(nz, c2_dz_u(I,:), dzN2_u(I,:), Sfn_unlim_u(I,:)) else do K=2,nz Sfn_unlim_u(I,K) = 0. @@ -1027,25 +1067,36 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do K=nz,2,-1 do I=is-1,ie + + if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then + Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i+1,j,k) + h(i+1,j,k-1))) + 4.0*hn_2 ) / & + ( ((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k) + (h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1)) + & + ((h(i+1,j,k)+hn_2)*tv%SpV_avg(i+1,j,k) + (h(i+1,j,k-1)+hn_2)*tv%SpV_avg(i+1,j,k-1)) ) + ! Use an average density to convert the volume streamfunction estimate into a mass streamfunction. + Z_to_H = (GV%RZ_to_H*Rho_avg) + else + Z_to_H = GV%Z_to_H + endif + if (k > nk_linear) then if (use_EOS) then if (uhtot(I,j) <= 0.0) then ! The transport that must balance the transport below is positive. - Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i,j,k)) * GV%H_to_Z + Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i,j,k)) else ! (uhtot(I,j) > 0.0) - Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i+1,j,k)) * GV%H_to_Z + Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i+1,j,k)) endif - ! The actual streamfunction at each interface. - Sfn_est = (Sfn_unlim_u(I,K) + slope2_Ratio_u(I,K)*Sfn_safe) / (1.0 + slope2_Ratio_u(I,K)) - else ! With .not.use_EOS, the layers are constant density. - Sfn_est = Sfn_unlim_u(I,K) + ! Determine the actual streamfunction at each interface. + Sfn_est = (Z_to_H*Sfn_unlim_u(I,K) + slope2_Ratio_u(I,K)*Sfn_safe) / (1.0 + slope2_Ratio_u(I,K)) + else ! When use_EOS is false, the layers are constant density. + Sfn_est = Z_to_H*Sfn_unlim_u(I,K) endif ! Make sure that there is enough mass above to allow the streamfunction ! to satisfy the boundary condition of 0 at the surface. - Sfn_in_H = min(max(Sfn_est * GV%Z_to_H, -h_avail_rsum(i,j,K)), h_avail_rsum(i+1,j,K)) + Sfn_in_H = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i+1,j,K)) ! The actual transport is limited by the mass available in the two ! neighboring grid cells. @@ -1076,6 +1127,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! else ! sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i+1,j,k)) ! endif + endif uhtot(I,j) = uhtot(I,j) + uhD(I,j,k) @@ -1087,6 +1139,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! A second order centered estimate is used for the density transferred ! between water columns. + if (allocated(tv%SpV_avg)) then + G_scale = GV%H_to_RZ * GV%g_Earth / Rho_avg + else + G_scale = GV%g_Earth * GV%H_to_Z + endif + Work_u(I,j) = Work_u(I,j) + G_scale * & ( uhtot(I,j) * drdkDe_u(I,K) - & (uhD(I,j,k) * drdi_u(I,k)) * 0.25 * & @@ -1099,18 +1157,19 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Calculate the meridional fluxes and gradients. - !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & - !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & + !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S,dz, & + !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,dz_neglect2, & !$OMP h_neglect2,int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, & - !$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1,diag_sfn_y, & - !$OMP diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,Tsgs2, & - !$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & + !$OMP I_slope_max2,vhD,h_avail,Work_v,CS,slope_y,cg1,hn_2,& + !$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,& + !$OMP Tsgs2, present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v,S_h,S_hr, & - !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & - !$OMP drho_dT_dT_h,drho_dT_dT_hr, scrap,pres_h,T_h,T_hr, & + !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA,G_scale, & + !$OMP drho_dT_dT_h,drho_dT_dT_hr,scrap,pres_h,T_h,T_hr, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,pres_hr, & - !$OMP drdy,mag_grad2,Slope,slope2_Ratio_v,hN2_v, & - !$OMP Sfn_unlim_v,drdj_v,drdkDe_v,h_harm,c2_h_v, & + !$OMP dzg2A,dzg2B,dzaA,dzaB,dz_harm,Z_to_H, & + !$OMP drdy,mag_grad2,Slope,slope2_Ratio_v,dzN2_v,N2_unlim, & + !$OMP Sfn_unlim_v,Rho_avg,drdj_v,drdkDe_v,c2_dz_v, & !$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives) do J=js-1,je do K=nz,2,-1 @@ -1186,11 +1245,15 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2 haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect + if (GV%Boussinesq) then dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z - else + elseif (GV%semi_Boussinesq) then dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect + else + dzaL = 0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect + dzaR = 0.5*(dz(i,j+1,k-1) + dz(i,j+1,k)) + dz_neglect endif ! Use the harmonic mean thicknesses to weight the horizontal gradients. ! These unnormalized weights have been rearranged to minimize divisions. @@ -1205,9 +1268,22 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV haA = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect haB = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect - ! hN2_v is used with the FGNV streamfunction formulation - hN2_v(i,K) = (0.5 * GV%H_to_Z * ( hg2A / haA + hg2B / haB )) * & - max(drdz*G_rho0, N2_floor) + if (GV%Boussinesq) then + N2_unlim = drdz*G_rho0 + else + N2_unlim = (GV%g_Earth*GV%RZ_to_H) * & + ((wtL * drdkL + wtR * drdkR) / (haL*wtL + haR*wtR)) + endif + + dzg2A = dz(i,j,k-1)*dz(i,j+1,k-1) + dz_neglect2 + dzg2B = dz(i,j,k)*dz(i,j+1,k) + dz_neglect2 + dzaA = 0.5*(dz(i,j,k-1) + dz(i,j+1,k-1)) + dz_neglect + dzaB = 0.5*(dz(i,j,k) + dz(i,j+1,k)) + dz_neglect + + ! dzN2_v is used with the FGNV streamfunction formulation + dzN2_v(i,K) = (0.5*( dzg2A / dzaA + dzg2B / dzaB )) * max(N2_unlim, N2_floor) + if (find_work .and. CS%GM_src_alt) & + hN2_y_PE(i,J,k) = (0.5*( hg2A / haA + hg2B / haB )) * max(N2_unlim, N2_floor) endif if (present_slope_y) then Slope = slope_y(i,J,k) @@ -1239,10 +1315,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max) - hN2_y_PE(i,J,k) = hN2_v(i,K) if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope - ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) ! Avoid moving dense water upslope from below the level of @@ -1273,10 +1347,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) - hN2_v(i,K) = GV%g_prime(K) + dzN2_v(i,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) - hN2_v(i,K) = N2_floor * dz_neglect + dzN2_v(i,K) = N2_floor * dz_neglect Sfn_unlim_v(i,K) = 0. endif ! if (k > nk_linear) if (CS%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,J,K) = Sfn_unlim_v(i,K) @@ -1285,10 +1359,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%use_FGNV_streamfn) then do k=1,nz ; do i=is,ie ; if (G%OBCmaskCv(i,J)>0.) then - h_harm = max( h_neglect, & - 2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) ) - c2_h_v(i,k) = CS%FGNV_scale * & - ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (GV%H_to_Z*h_harm) + dz_harm = max( dz_neglect, & + 2. * dz(i,j,k) * dz(i,j+1,k) / ( ( dz(i,j,k) + dz(i,j+1,k) ) + dz_neglect ) ) + c2_dz_v(i,k) = CS%FGNV_scale * ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / dz_harm endif ; enddo ; enddo ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. @@ -1297,7 +1370,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do K=2,nz Sfn_unlim_v(i,K) = (1. + CS%FGNV_scale) * Sfn_unlim_v(i,K) enddo - call streamfn_solver(nz, c2_h_v(i,:), hN2_v(i,:), Sfn_unlim_v(i,:)) + call streamfn_solver(nz, c2_dz_v(i,:), dzN2_v(i,:), Sfn_unlim_v(i,:)) else do K=2,nz Sfn_unlim_v(i,K) = 0. @@ -1308,25 +1381,35 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do K=nz,2,-1 do i=is,ie + if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then + Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i,j+1,k) + h(i,j+1,k-1))) + 4.0*hn_2 ) / & + ( ((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k) + (h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1)) + & + ((h(i,j+1,k)+hn_2)*tv%SpV_avg(i,j+1,k) + (h(i,j+1,k-1)+hn_2)*tv%SpV_avg(i,j+1,k-1)) ) + ! Use an average density to convert the volume streamfunction estimate into a mass streamfunction. + Z_to_H = (GV%RZ_to_H*Rho_avg) + else + Z_to_H = GV%Z_to_H + endif + if (k > nk_linear) then if (use_EOS) then if (vhtot(i,J) <= 0.0) then ! The transport that must balance the transport below is positive. - Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j,k)) * GV%H_to_Z + Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j,k)) else ! (vhtot(I,j) > 0.0) - Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j+1,k)) * GV%H_to_Z + Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j+1,k)) endif - ! The actual streamfunction at each interface. - Sfn_est = (Sfn_unlim_v(i,K) + slope2_Ratio_v(i,K)*Sfn_safe) / (1.0 + slope2_Ratio_v(i,K)) - else ! With .not.use_EOS, the layers are constant density. - Sfn_est = Sfn_unlim_v(i,K) + ! Find the actual streamfunction at each interface. + Sfn_est = (Z_to_H*Sfn_unlim_v(i,K) + slope2_Ratio_v(i,K)*Sfn_safe) / (1.0 + slope2_Ratio_v(i,K)) + else ! When use_EOS is false, the layers are constant density. + Sfn_est = Z_to_H*Sfn_unlim_v(i,K) endif ! Make sure that there is enough mass above to allow the streamfunction ! to satisfy the boundary condition of 0 at the surface. - Sfn_in_H = min(max(Sfn_est * GV%Z_to_H, -h_avail_rsum(i,j,K)), h_avail_rsum(i,j+1,K)) + Sfn_in_H = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i,j+1,K)) ! The actual transport is limited by the mass available in the two ! neighboring grid cells. @@ -1367,6 +1450,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! A second order centered estimate is used for the density transferred ! between water columns. + if (allocated(tv%SpV_avg)) then + G_scale = GV%H_to_RZ * GV%g_Earth / Rho_avg + else + G_scale = GV%g_Earth * GV%H_to_Z + endif + Work_v(i,J) = Work_v(i,J) + G_scale * & ( vhtot(i,J) * drdkDe_v(i,K) - & (vhD(i,J,k) * drdj_v(i,k)) * 0.25 * & @@ -1383,7 +1472,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do J=js-1,je ; do i=is,ie ; vhD(i,J,1) = -vhtot(i,J) ; enddo ; enddo else EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) - !$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB) + !$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB,G_scale) do j=js,je if (use_EOS) then do I=is-1,ie @@ -1397,9 +1486,15 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do I=is-1,ie uhD(I,j,1) = -uhtot(I,j) + G_scale = GV%g_Earth * GV%H_to_Z if (use_EOS) then drdiB = drho_dT_u(I) * (T(i+1,j,1)-T(i,j,1)) + & drho_dS_u(I) * (S(i+1,j,1)-S(i,j,1)) + if (allocated(tv%SpV_avg)) then + G_scale = GV%H_to_RZ * GV%g_Earth * & + ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1) + (h(i+1,j,1)+hn_2) * tv%SpV_avg(i+1,j,1)) / & + ( (h(i,j,1) + h(i+1,j,1)) + 2.0*hn_2 ) ) + endif endif if (CS%use_GM_work_bug) then Work_u(I,j) = Work_u(I,j) + G_scale * & @@ -1414,7 +1509,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo EOSdom_v(:) = EOS_domain(G%HI) - !$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB) + !$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB,G_scale) do J=js-1,je if (use_EOS) then do i=is,ie @@ -1428,9 +1523,15 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV do i=is,ie vhD(i,J,1) = -vhtot(i,J) + G_scale = GV%g_Earth * GV%H_to_Z if (use_EOS) then drdjB = drho_dT_v(i) * (T(i,j+1,1)-T(i,j,1)) + & drho_dS_v(i) * (S(i,j+1,1)-S(i,j,1)) + if (allocated(tv%SpV_avg)) then + G_scale = GV%H_to_RZ * GV%g_Earth * & + ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1) + (h(i,j+1,1)+hn_2) * tv%SpV_avg(i,j+1,1)) / & + ( (h(i,j,1) + h(i,j+1,1)) + 2.0*hn_2 ) ) + endif endif Work_v(i,J) = Work_v(i,J) - G_scale * & ( (vhD(i,J,1) * drdjB) * 0.25 * & @@ -1451,11 +1552,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then do j=js,je ; do i=is,ie ; do k=nz,1,-1 - PE_release_h = -0.25*(KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & + PE_release_h = -0.25 * (GV%H_to_RZ*US%L_to_Z**2) * & + (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + & Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + & Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) - MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + US%L_to_Z**2 * GV%Rho0 * PE_release_h + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h enddo ; enddo ; enddo endif ; endif @@ -1471,16 +1573,18 @@ end subroutine thickness_diffuse_full !> Tridiagonal solver for streamfunction at interfaces subroutine streamfn_solver(nk, c2_h, hN2, sfn) integer, intent(in) :: nk !< Number of layers - real, dimension(nk), intent(in) :: c2_h !< Wave speed squared over thickness in layers [L2 Z-1 T-2 ~> m s-2] - real, dimension(nk+1), intent(in) :: hN2 !< Thickness times N2 at interfaces [L2 Z-1 T-2 ~> m s-2] - real, dimension(nk+1), intent(inout) :: sfn !< Streamfunction [Z L2 T-1 ~> m3 s-1] or arbitrary units + real, dimension(nk), intent(in) :: c2_h !< Wave speed squared over thickness in layers, rescaled to + !! [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2] + real, dimension(nk+1), intent(in) :: hN2 !< Thickness times N2 at interfaces times rescaling factors + !! [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2] + real, dimension(nk+1), intent(inout) :: sfn !< Streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] or arbitrary units !! On entry, equals diffusivity times slope. !! On exit, equals the streamfunction. ! Local variables real :: c1(nk) ! The dependence of the final streamfunction on the values below [nondim] real :: d1 ! The complement of c1(k) (i.e., 1 - c1(k)) [nondim] - real :: b_denom ! A term in the denominator of beta [L2 Z-1 T-2 ~> m s-2] - real :: beta ! The normalization for the pivot [Z T2 L-2 ~> s2 m-1] + real :: b_denom ! A term in the denominator of beta [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2] + real :: beta ! The normalization for the pivot [Z2 T2 H-1 L-2 ~> s2 m-1 or m2 s2 kg-1] integer :: k sfn(1) = 0. From 147ddf1d5707e79f3268c430d396f33e7d9956c3 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 25 Jul 2023 08:17:15 -0800 Subject: [PATCH 09/14] Brine plume (#401) * Salt data structures * First steps at brine plume: pass info from SIS2 * The brine plume parameterization, - including now passing the dimensional scaling tests. * Fix problem when running Tidal_bay case with gnu. * Avoiding visc_rem issues inside land mask. Tweaking the brine plume code. * Using the proper MLD in the brine plumes - it now works better on restart * Always including MLD in call to applyBoundary... - I could move it up and make it not optional. * Adding some OpenMP directives to brine plumes --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 13 ++- docs/zotero.bib | 9 ++ src/core/MOM_continuity_PPM.F90 | 8 +- src/core/MOM_forcing_type.F90 | 16 ++-- .../vertical/MOM_diabatic_aux.F90 | 82 +++++++++++++++++-- .../vertical/MOM_diabatic_driver.F90 | 8 +- .../vertical/MOM_set_viscosity.F90 | 34 ++++---- 7 files changed, 132 insertions(+), 38 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 251f37290d..a8398c3cc8 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -62,7 +62,7 @@ module MOM_surface_forcing_gfdl !! from MOM_domains) to indicate the staggering of !! the winds that are being provided in calls to !! update_ocean_model. - logical :: use_temperature !< If true, temp and saln used as state variables + logical :: use_temperature !< If true, temp and saln used as state variables. real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim]. real :: Rho0 !< Boussinesq reference density [R ~> kg m-3] @@ -175,6 +175,7 @@ module MOM_surface_forcing_gfdl real, pointer, dimension(:,:) :: t_flux =>NULL() !< sensible heat flux [W m-2] real, pointer, dimension(:,:) :: q_flux =>NULL() !< specific humidity flux [kg m-2 s-1] real, pointer, dimension(:,:) :: salt_flux =>NULL() !< salt flux [kg m-2 s-1] + real, pointer, dimension(:,:) :: excess_salt =>NULL() !< salt left behind by brine rejection [kg m-2 s-1] real, pointer, dimension(:,:) :: lw_flux =>NULL() !< long wave radiation [W m-2] real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation [W m-2] real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation [W m-2] @@ -304,6 +305,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + if (associated(IOB%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed) + do j=js-2,je+2 ; do i=is-2,ie+2 fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) @@ -576,6 +579,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call check_mask_val_consistency(IOB%salt_flux(i-i0,j-j0), G%mask2dT(i,j), i, j, 'salt_flux', G) enddo ; enddo endif + if (associated(IOB%excess_salt)) then + do j=js,je ; do i=is,ie + fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*(kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0)) + enddo ; enddo + endif !#CTRL# if (associated(CS%ctrl_forcing_CSp)) then !#CTRL# do j=js,je ; do i=is,ie @@ -1729,6 +1737,9 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) if (associated(iobt%mass_berg)) then chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks endif + if (associated(iobt%excess_salt)) then + chks = field_chksum( iobt%excess_salt ) ; if (root) write(outunit,100) 'iobt%excess_salt ', chks + endif 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/docs/zotero.bib b/docs/zotero.bib index c0c7ee3bd9..5acaee968a 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2738,3 +2738,12 @@ @article{kraus1967 journal = {Tellus} } +@article{Nguyen2009, + doi = {10.1029/2008JC005121}, + year = {2009}, + journal = {JGR Oceans}, + volume = {114}, + author = {A. T. Nguyen and D. Menemenlis and R. Kwok}, + title = {Improved modeling of the Arctic halocline with a subgrid-scale brine rejection parameterization}, + pages = {C11014} +} diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 090d1ee0fb..73c6503242 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -378,9 +378,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j)) else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif - if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)) & + if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) & du_max_CFL(I) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(I,k) - if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)) & + if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) & du_min_CFL(I) = -(dx_E*CFL_dt + u(I,j,k)) / visc_rem(I,k) enddo ; enddo endif @@ -1201,9 +1201,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j)) dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1)) else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif - if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)) & + if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) & dv_max_CFL(i) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(i,k) - if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)) & + if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) & dv_min_CFL(i) = -(dy_N*CFL_dt + v(i,J,k)) / visc_rem(i,k) enddo ; enddo endif diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 4897771100..9623256b88 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -135,8 +135,10 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & salt_flux => NULL(), & !< net salt flux into the ocean [R Z T-1 ~> kgSalt m-2 s-1] salt_flux_in => NULL(), & !< salt flux provided to the ocean from coupler [R Z T-1 ~> kgSalt m-2 s-1] - salt_flux_added => NULL() !< additional salt flux from restoring or flux adjustment before adjustment + salt_flux_added => NULL(), & !< additional salt flux from restoring or flux adjustment before adjustment !! to net zero [R Z T-1 ~> kgSalt m-2 s-1] + salt_left_behind => NULL() !< salt left in ocean at the surface from brine rejection + !! [R Z T-1 ~> kgSalt m-2 s-1] ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice) real, pointer, dimension(:,:) :: p_surf_full => NULL() @@ -746,15 +748,15 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & endif ! Salt fluxes - Net_salt(i) = 0.0 - if (do_NSR) Net_salt_rate(i) = 0.0 + net_salt(i) = 0.0 + if (do_NSR) net_salt_rate(i) = 0.0 ! Convert salt_flux from kg (salt)/(m^2 * s) to ! Boussinesq: (ppt * m) ! non-Bouss: (g/m^2) if (associated(fluxes%salt_flux)) then - Net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H + net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H !Repeat above code for 'rate' term - if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H + if (do_NSR) net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H endif ! Diagnostics follow... @@ -1947,6 +1949,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, diag%axesT1,Time,'Salt flux into ocean at surface due to restoring or flux adjustment', & units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) + handles%id_saltFluxAdded = register_diag_field('ocean_model', 'salt_left_behind', & + diag%axesT1,Time,'Salt left in ocean at surface due to ice formation', & + units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) + handles%id_saltFluxGlobalAdj = register_scalar_field('ocean_model', & 'salt_flux_global_restoring_adjustment', Time, diag, & 'Adjustment needed to balance net global salt flux into ocean at surface', & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index ba265af5e2..3742e93229 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -67,7 +67,10 @@ module MOM_diabatic_aux !! e-folding depth of incoming shortwave radiation. type(external_field) :: sbc_chl !< A handle used in time interpolation of !! chlorophyll read from a file. - logical :: chl_from_file !< If true, chl_a is read from a file. + logical :: chl_from_file !< If true, chl_a is read from a file. + logical :: do_brine_plume !< If true, insert salt flux below the surface according to + !! a parameterization by \cite Nguyen2009. + integer :: brine_plume_n !< The exponent in the brine plume parameterization. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output @@ -1034,7 +1037,7 @@ end subroutine diagnoseMLDbyEnergy subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, & aggregate_FW_forcing, evap_CFL_limit, & minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, & - SkinBuoyFlux ) + SkinBuoyFlux, MLD) type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1064,6 +1067,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. + real, pointer, dimension(:,:), optional :: MLD!< Mixed layer depth for brine plumes [Z ~> m] ! Local variables integer, parameter :: maxGroundings = 5 @@ -1102,7 +1106,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate) ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - netMassInOut_rate! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1] + netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1] + mixing_depth ! Mixed layer depth [Z -> m] real, dimension(SZI_(G), SZK_(GV)) :: & h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2] T2d, & ! A 2-d copy of the layer temperatures [C ~> degC] @@ -1132,6 +1137,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! and rejected brine are initially applied in vanishingly thin layers at the ! top of the layer before being mixed throughout the layer. logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux. + real, dimension(SZI_(G)) :: dK ! Depth [Z ~> m]. + real, dimension(SZI_(G)) :: A_brine ! Constant [Z-(n+1) ~> m-(n+1)]. + real :: fraction_left_brine ! Sum for keeping track of the fraction of brine so far (in depth) + real :: plume_fraction ! Sum for keeping track of the fraction of brine so far (in depth) + real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, k, nz, nb character(len=45) :: mesg @@ -1139,6 +1149,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Idt = 1.0 / dt + plume_flux = 0.0 calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) calculate_buoyancy = present(SkinBuoyFlux) @@ -1158,6 +1169,17 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 endif + if (CS%do_brine_plume .and. .not. associated(MLD)) then + call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& + "Brine plume parameterization requires a mixed-layer depth,\n"//& + "currently coming from the energetic PBL scheme.") + endif + if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then + call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& + "Brine plume parameterization requires DO_BRINE_PLUME\n"//& + "to be turned on in SIS2 as well as MOM6.") + endif + ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total ! depth of the ocean is vanishing. It does not (yet) handle a value of zero. ! To accommodate vanishing upper layers, we need to allow for an instantaneous @@ -1173,9 +1195,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,& !$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, & !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, & - !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, & + !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,& !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, & - !$OMP EnthalpyConst) & + !$OMP EnthalpyConst,MLD) & !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, & !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, & !$OMP IforcingDepthScale, & @@ -1183,7 +1205,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, & !$OMP netmassinout_rate,netheat_rate,netsalt_rate, & !$OMP drhodt,drhods,pen_sw_bnd_rate, & - !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst) & + !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, & + !$OMP mixing_depth,A_brine,fraction_left_brine, & + !$OMP plume_flux,plume_fraction,dK) & !$OMP firstprivate(SurfPressure) do j=js,je ! Work in vertical slices for efficiency @@ -1300,6 +1324,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! ocean (and corresponding outward heat content), and ignoring penetrative SW. ! B/ update mass, salt, temp from mass leaving ocean. ! C/ update temp due to penetrative SW + if (CS%do_brine_plume) then + do i=is,ie + mixing_depth(i) = max(MLD(i,j) - minimum_forcing_depth * GV%H_to_Z, minimum_forcing_depth * GV%H_to_Z) + mixing_depth(i) = min(mixing_depth(i), max(sum(h(i,j,:)), GV%angstrom_h) * GV%H_to_Z) + A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1)) + enddo + endif + do i=is,ie if (G%mask2dT(i,j) > 0.) then @@ -1372,8 +1404,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t enddo ! k=1,1 ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt. + fraction_left_brine = 1.0 do k=1,nz - ! Place forcing into this layer if this layer has nontrivial thickness. ! For layers thin relative to 1/IforcingDepthScale, then distribute ! forcing into deeper layers. @@ -1388,6 +1420,33 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i) endif + if (CS%do_brine_plume .and. associated(fluxes%salt_left_behind)) then + if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then + ! Place forcing into this layer by depth for brine plume parameterization. + if (k == 1) then + dK(i) = 0.5 * h(i,j,k) * GV%H_to_Z ! Depth of center of layer K + plume_flux = - (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + plume_fraction = 1.0 + else + dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * GV%H_to_Z ! Depth of center of layer K + plume_flux = 0.0 + endif + if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then + plume_fraction = min(fraction_left_brine, A_brine(i) * dK(i) ** CS%brine_plume_n & + * h(i,j,k) * GV%H_to_Z) + else + IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) ) + ! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale. + plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale) + endif + fraction_left_brine = fraction_left_brine - plume_fraction + plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) & + * GV%RZ_to_H + else + plume_flux = 0.0 + endif + endif + ! Change in state due to forcing dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) ) @@ -1432,7 +1491,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t endif Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness T2d(i,k) = (hOld*T2d(i,k) + dTemp)*Ithickness - tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt)*Ithickness + tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt + plume_flux)*Ithickness elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (h<0)') write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j) @@ -1702,6 +1761,13 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori CS%use_calving_heat_content = .false. endif + call get_param(param_file, mdl, "DO_BRINE_PLUME", CS%do_brine_plume, & + "If true, use a brine plume parameterization from "//& + "Nguyen et al., 2009.", default=.false.) + call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, & + "If using the brine plume parameterization, set the integer exponent.", & + default=5, do_not_log=.not.CS%do_brine_plume) + if (useALEalgorithm) then CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, & Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0c28c063ea..466ebbabca 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -815,7 +815,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) @@ -875,7 +875,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) endif ! endif for CS%use_energetic_PBL @@ -1360,7 +1360,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) @@ -1414,7 +1414,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) endif ! endif for CS%use_energetic_PBL diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 47d4dffef6..481aa5e9fc 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1027,22 +1027,24 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) endif endif - if (CS%body_force_drag .and. (h_bbl_drag(i) > 0.0)) then - ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force. - h_sum = 0.0 - I_hwtot = 1.0 / h_bbl_drag(i) - do k=nz,1,-1 - h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * I_hwtot - if (m==1) then - visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (CS%cdrag*US%L_to_Z*umag_avg(I)) * h_bbl_fr - else - visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (CS%cdrag*US%L_to_Z*umag_avg(i)) * h_bbl_fr - endif - h_sum = h_sum + h_at_vel(i,k) - if (h_sum >= h_bbl_drag(i)) exit ! The top of this layer is above the drag zone. - enddo - ! Do not enhance the near-bottom viscosity in this case. - Kv_bbl = CS%Kv_BBL_min + if (CS%body_force_drag) then + if (h_bbl_drag(i) > 0.0) then + ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force. + h_sum = 0.0 + I_hwtot = 1.0 / h_bbl_drag(i) + do k=nz,1,-1 + h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * I_hwtot + if (m==1) then + visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (CS%cdrag*US%L_to_Z*umag_avg(I)) * h_bbl_fr + else + visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (CS%cdrag*US%L_to_Z*umag_avg(i)) * h_bbl_fr + endif + h_sum = h_sum + h_at_vel(i,k) + if (h_sum >= h_bbl_drag(i)) exit ! The top of this layer is above the drag zone. + enddo + ! Do not enhance the near-bottom viscosity in this case. + Kv_bbl = CS%Kv_BBL_min + endif endif kv_bbl = max(CS%Kv_BBL_min, kv_bbl) From d5ba107af21757390fb82d123ae1bf9c236ec0e4 Mon Sep 17 00:00:00 2001 From: Spencer Jones <41342785+cspencerjones@users.noreply.github.com> Date: Tue, 25 Jul 2023 13:14:17 -0500 Subject: [PATCH 10/14] +add h to drifters interface (#408) This commit brings the drifters interface up-to-date with the current version of the drifters package, which requires h (layer thickness) to calculate the vertical movement of particles. The interfaces in the code and in config_src/external are updated to pass this information to the drifters package. --- .../external/drifters/MOM_particles.F90 | 44 +++++++++++++------ src/core/MOM.F90 | 17 +++++-- 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/config_src/external/drifters/MOM_particles.F90 b/config_src/external/drifters/MOM_particles.F90 index aad918e5a4..fa3840c6c2 100644 --- a/config_src/external/drifters/MOM_particles.F90 +++ b/config_src/external/drifters/MOM_particles.F90 @@ -11,19 +11,20 @@ module MOM_particles_mod implicit none ; private public particles, particles_run, particles_init, particles_save_restart, particles_end +public particles_to_k_space, particles_to_z_space contains !> Initializes particles container "parts" -subroutine particles_init(parts, Grid, Time, dt, u, v) +subroutine particles_init(parts, Grid, Time, dt, u, v, h) ! Arguments type(particles), pointer, intent(out) :: parts !< Container for all types and memory type(ocean_grid_type), target, intent(in) :: Grid !< Grid type from parent model type(time_type), intent(in) :: Time !< Time type from parent model - real, intent(in) :: dt !< particle timestep [s] - real, dimension(:,:,:), intent(in) :: u !< Zonal velocity field [m s-1] - real, dimension(:,:,:), intent(in) :: v !< Meridional velocity field [m s-1] - + real, intent(in) :: dt !< particle timestep in seconds [T ~> s] + real, dimension(:,:,:),intent(in) :: u !< Zonal velocity field [L T-1 ~> m s-1] + real, dimension(:,:,:),intent(in) :: v !< Meridional velocity field [L T-1 ~> m s-1] + real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2] end subroutine particles_init !> The main driver the steps updates particles @@ -31,8 +32,8 @@ subroutine particles_run(parts, time, uo, vo, ho, tv, stagger) ! Arguments type(particles), pointer :: parts !< Container for all types and memory type(time_type), intent(in) :: time !< Model time - real, dimension(:,:,:), intent(in) :: uo !< Ocean zonal velocity [m s-1] - real, dimension(:,:,:), intent(in) :: vo !< Ocean meridional velocity [m s-1] + real, dimension(:,:,:), intent(in) :: uo !< Ocean zonal velocity [L T-1 ~>m s-1] + real, dimension(:,:,:), intent(in) :: vo !< Ocean meridional velocity [L T-1~> m s-1] real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< structure containing pointers to available thermodynamic fields integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered @@ -41,21 +42,38 @@ end subroutine particles_run !>Save particle locations (and sometimes other vars) to restart file -subroutine particles_save_restart(parts, temp, salt) +subroutine particles_save_restart(parts, h, temp, salt) ! Arguments type(particles), pointer :: parts !< Container for all types and memory - real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature - real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity + real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2] + real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature [C ~> degC] + real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity [S ~> ppt] end subroutine particles_save_restart !> Deallocate all memory and disassociated pointer -subroutine particles_end(parts, temp, salt) +subroutine particles_end(parts, h, temp, salt) ! Arguments type(particles), pointer :: parts !< Container for all types and memory - real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature - real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity + real, dimension(:,:,:),intent(in) :: h !< Thickness of each layer [H ~> m or kg m-2] + real, dimension(:,:,:), optional, intent(in) :: temp !< Optional container for temperature [C ~> degC] + real, dimension(:,:,:), optional, intent(in) :: salt !< Optional container for salinity [S ~> ppt] end subroutine particles_end +subroutine particles_to_k_space(parts, h) + ! Arguments + type(particles), pointer :: parts !< Container for all types and memory + real, dimension(:,:,:),intent(in) :: h !< Thickness of layers [H ~> m or kg m-2] + +end subroutine particles_to_k_space + + +subroutine particles_to_z_space(parts, h) + ! Arguments + type(particles), pointer :: parts !< Container for all types and memory + real, dimension(:,:,:),intent(in) :: h !< Thickness of layers [H ~> m or kg m-2] + +end subroutine particles_to_z_space + end module MOM_particles_mod diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index df3a308e85..d2df26524d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -162,7 +162,7 @@ module MOM use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end - +use MOM_particles_mod, only : particles_to_k_space, particles_to_z_space implicit none ; private #include @@ -1541,6 +1541,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call preAle_tracer_diagnostics(CS%tracer_Reg, G, GV) + if (CS%use_particles) then + call particles_to_z_space(CS%particles, h) + endif + if (CS%debug) then call MOM_state_chksum("Pre-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US, omit_corners=.true.) call hchksum(tv%T,"Pre-ALE T", G%HI, haloshift=1, omit_corners=.true., scale=US%C_to_degC) @@ -1588,6 +1592,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_end(id_clock_ALE) endif ! endif for the block "if ( CS%use_ALE_algorithm )" + + if (CS%use_particles) then + call particles_to_k_space(CS%particles, h) + endif + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) call create_group_pass(pass_uv_T_S_h, u, v, G%Domain, halo=dynamics_stencil) if (associated(tv%T)) & @@ -3232,7 +3241,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) G => CS%G ; GV => CS%GV ; US => CS%US if (CS%use_particles) then - call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v) + call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v, CS%h) endif ! Write initial conditions @@ -3935,7 +3944,7 @@ subroutine save_MOM6_internal_state(CS, dirs, time, stamp_time) ! Could call save_restart(CS%restart_CSp) here - if (CS%use_particles) call particles_save_restart(CS%particles) + if (CS%use_particles) call particles_save_restart(CS%particles, CS%h) end subroutine save_MOM6_internal_state @@ -3978,7 +3987,7 @@ subroutine MOM_end(CS) endif if (CS%use_particles) then - call particles_end(CS%particles) + call particles_end(CS%particles, CS%h) deallocate(CS%particles) endif From 249a0788d8a041cd7a4cef497aa1c5b2cf75931d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 11 Jul 2023 17:40:32 -0400 Subject: [PATCH 11/14] +*Revise units of arguments to vert_fill_TS Pass dt_kappa_smooth to calc_isoneutral_slopes and vert_fill_TS in units of [H Z ~> m2 or kg m-1] instead of [Z2 ~> m2] for consistency with the units of other diffusivities in the code and to reduce the depenency on the Boussinesq reference density in non-Boussinesq configurations. In addition to the changes to the units of these two arguments, there is a new unit_scale_type argument to vert_fill_TS and MOM_calc_varT and a new verticalGrid_type argument to MOM_stoch_eos_init. The units of 4 vertical diffusivities in the control structures in 4 different modules are also changed accordingly. All answers are bitwise identical in Boussinesq mode, but they can change for some non-Boussinesq configurations. There are new mandatory arguments to three publicly visible routines. --- src/core/MOM.F90 | 4 +-- src/core/MOM_isopycnal_slopes.F90 | 34 +++++++++++-------- src/core/MOM_stoch_eos.F90 | 12 ++++--- src/parameterizations/lateral/MOM_MEKE.F90 | 3 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 4 +-- .../lateral/MOM_thickness_diffuse.F90 | 10 +++--- .../vertical/MOM_internal_tide_input.F90 | 6 ++-- .../vertical/MOM_set_diffusivity.F90 | 12 +++---- 8 files changed, 46 insertions(+), 39 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index d2df26524d..d7e2d74735 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1090,7 +1090,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call cpu_clock_end(id_clock_stoch) call cpu_clock_begin(id_clock_varT) if (CS%use_stochastic_EOS) then - call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt) + call MOM_calc_varT(G, GV, US, h, CS%tv, CS%stoch_eos_CS, dt) if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1) endif call cpu_clock_end(id_clock_varT) @@ -3022,7 +3022,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & new_sim = is_new_run(restart_CSp) if (use_temperature) then - CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp) + CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS%stoch_eos_CS, restart_CSp) else CS%use_stochastic_EOS = .false. endif diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 07dd19b0a6..73ae8a9816 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -36,8 +36,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables - real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity - !! times a smoothing timescale [Z2 ~> m2]. + real, intent(in) :: dt_kappa_smooth !< A smoothing vertical + !! diffusivity times a smoothing + !! timescale [H Z ~> m2 or kg m-1] logical, intent(in) :: use_stanley !< turn on stanley param in slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] @@ -142,7 +143,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 - dz_neglect = GV%H_subroundoff * GV%H_to_Z + dz_neglect = GV%dZ_subroundoff local_open_u_BC = .false. local_open_v_BC = .false. @@ -195,9 +196,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan if (use_EOS) then if (present(halo)) then - call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, halo+1) + call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, halo+1) else - call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, 1) + call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, 1) endif endif @@ -341,9 +342,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan slope = slope * max(g%mask2dT(i,j),g%mask2dT(i+1,j)) endif slope_x(I,j,K) = slope - if (present(dzSxN)) dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & - + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N - * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2 + if (present(dzSxN)) & + dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & + + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N + * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2 enddo ! I enddo ; enddo ! end of j-loop @@ -477,9 +479,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan slope = slope * max(g%mask2dT(i,j),g%mask2dT(i,j+1)) endif slope_y(i,J,K) = slope - if (present(dzSyN)) dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & - + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N - * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2 + if (present(dzSyN)) & + dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & + + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N + * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2 enddo ! i enddo ; enddo ! end of j-loop @@ -488,14 +491,15 @@ end subroutine calc_isoneutral_slopes !> Returns tracer arrays (nominally T and S) with massless layers filled with !! sensible values, by diffusing vertically with a small but constant diffusivity. -subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom) +subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, larger_h_denom) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T_in !< Input temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S_in !< Input salinity [S ~> ppt] real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing - !! times a smoothing timescale [Z2 ~> m2]. + !! times a smoothing timescale [H Z ~> m2 or kg m-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T_f !< Filled temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S_f !< Filled salinity [S ~> ppt] integer, optional, intent(in) :: halo_here !< Number of halo points to work on, @@ -525,10 +529,10 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke h_neglect = GV%H_subroundoff - kap_dt_x2 = (2.0*kappa_dt)*GV%Z_to_H**2 + kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H. h0 = h_neglect if (present(larger_h_denom)) then - if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*GV%Z_to_H + if (larger_h_denom) h0 = 1.0e-16*sqrt(0.5*kap_dt_x2) endif if (kap_dt_x2 <= 0.0) then diff --git a/src/core/MOM_stoch_eos.F90 b/src/core/MOM_stoch_eos.F90 index deb878e99c..2bd742be6d 100644 --- a/src/core/MOM_stoch_eos.F90 +++ b/src/core/MOM_stoch_eos.F90 @@ -40,7 +40,7 @@ module MOM_stoch_eos real :: stanley_coeff !< Coefficient correlating the temperature gradient !! and SGS T variance [nondim]; if <0, turn off scheme in all codes real :: stanley_a !< a in exp(aX) in stochastic coefficient [nondim] - real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] + real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] !>@{ Diagnostic IDs integer :: id_stoch_eos = -1, id_stoch_phi = -1, id_tvar_sgs = -1 @@ -51,9 +51,10 @@ module MOM_stoch_eos contains !> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used. -logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_CS) +logical function MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS, restart_CS) type(time_type), intent(in) :: Time !< Time for stochastic process type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse type(diag_ctrl), target, intent(inout) :: diag !< Structure used to control diagnostics @@ -80,7 +81,7 @@ logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_C call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s, & + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T, & do_not_log=(CS%stanley_coeff<0.0)) ! Don't run anything if STANLEY_COEFF < 0 @@ -193,9 +194,10 @@ subroutine post_stoch_EOS_diags(CS, tv, diag) end subroutine post_stoch_EOS_diags !> Computes a parameterization of the SGS temperature variance -subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) +subroutine MOM_calc_varT(G, GV, US, h, tv, CS, dt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure @@ -219,7 +221,7 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) ! extreme gradients along layers which are vanished against topography. It is ! still a poor approximation in the interior when coordinates are strongly tilted. if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0) - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo_here=1, larger_h_denom=.true.) do k=1,G%ke do j=G%jsc,G%jec diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 02338fab96..0ef261a956 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -1578,7 +1578,8 @@ subroutine ML_MEKE_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, f h_v(i,J,k) = 0.5*(h(i,j,k)*G%mask2dT(i,j) + h(i,j+1,k)*G%mask2dT(i,j+1)) + GV%Angstrom_H enddo; enddo; enddo; call find_eta(h, tv, G, GV, US, e, halo_size=2) - call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7, .false., slope_x, slope_y) + ! Note the hard-coded dimenisional constant in the following line. + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7*GV%m2_s_to_HZ_T, .false., slope_x, slope_y) call pass_vector(slope_x, slope_y, G%Domain) do j=js-1,je+1; do i=is-1,ie+1 slope_x_vert_avg(I,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), GV%H_subroundoff) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 8f0aa02b12..e74e48055a 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -135,7 +135,7 @@ module MOM_lateral_mixing_coeffs !! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power) real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m] real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim] - real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] + real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] integer :: Res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any !! positive integer power may be used, but even powers !! and especially 2 are coded to be more efficient. @@ -1265,7 +1265,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) endif if (CS%calculate_Eady_growth_rate) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8617795e16..f24f790d06 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -44,9 +44,9 @@ module MOM_thickness_diffuse real :: Kh_eta_bg !< Background isopycnal height diffusivity [L2 T-1 ~> m2 s-1] real :: Kh_eta_vel !< Velocity scale that is multiplied by the grid spacing to give !! the isopycnal height diffusivity [L T-1 ~> m s-1] - real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim]. - real :: kappa_smooth !< Vertical diffusivity used to interpolate more - !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim] + real :: kappa_smooth !< Vertical diffusivity used to interpolate more sensible values + !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] logical :: thickness_diffuse !< If true, interfaces heights are diffused. logical :: use_FGNV_streamfn !< If true, use the streamfunction formulation of !! Ferrari et al., 2010, which effectively emphasizes @@ -798,7 +798,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (use_EOS) then halo = 1 ! Default halo to fill is 1 - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo, larger_h_denom=.true.) endif ! Rescale the thicknesses, perhaps using the specific volume. @@ -2191,7 +2191,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", CS%use_FGNV_streamfn, & "If true, use the streamfunction formulation of "//& "Ferrari et al., 2010, which effectively emphasizes "//& diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 7ec612f141..95e33929df 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -41,7 +41,7 @@ module MOM_int_tide_input real :: TKE_itide_max !< Maximum Internal tide conversion !! available to mix above the BBL [R Z3 T-3 ~> W m-2] real :: kappa_fill !< Vertical diffusivity used to interpolate sensible values - !! of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real, allocatable, dimension(:,:) :: TKE_itidal_coef !< The time-invariant field that enters the TKE_itidal input calculation [R Z3 T-2 ~> J m-2]. @@ -118,7 +118,7 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) ! Smooth the properties through massless layers. if (use_EOS) then - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_fill*dt, T_f, S_f, G, GV, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_fill*dt, T_f, S_f, G, GV, US, larger_h_denom=.true.) endif call find_N2_bottom(h, tv, T_f, S_f, itide%h2, fluxes, G, GV, US, N2_bot) @@ -352,7 +352,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_fill, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9dc7b81c46..dfd264c92a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -87,7 +87,7 @@ module MOM_set_diffusivity real :: Kd_add !< uniform diffusivity added everywhere without !! filtering or scaling [Z2 T-1 ~> m2 s-1]. real :: Kd_smooth !< Vertical diffusivity used to interpolate more - !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + !! sensible values of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostic output timing logical :: limit_dissipation !< If enabled, dissipation is limited to be larger @@ -274,7 +274,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i integer :: i, j, k, is, ie, js, je, nz, isd, ied, jsd, jed - real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [Z2 ~> m2] + real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [H Z ~> m2 or kg m-1] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -289,7 +289,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (CS%answer_date < 20190101) then ! These hard-coded dimensional parameters are being replaced. - kappa_dt_fill = US%m_to_Z**2 * 1.e-3 * 7200. + kappa_dt_fill = 1.e-3*GV%m2_s_to_HZ_T * 7200.*US%s_to_T else kappa_dt_fill = CS%Kd_smooth * dt endif @@ -340,7 +340,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then call full_convection(G, GV, US, h, tv, T_f, S_f, fluxes%p_surf, & - (GV%Z_to_H**2)*kappa_dt_fill, halo=1) + GV%Z_to_H*kappa_dt_fill, halo=1) call calc_kappa_shear_vertex(u, v, h, T_f, S_f, tv, fluxes%p_surf, visc%Kd_shear, & visc%TKE_turb, visc%Kv_shear_Bu, dt, G, GV, US, CS%kappaShear_CSp) @@ -380,7 +380,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call hchksum(tv%S, "before vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt) call hchksum(h, "before vert_fill_TS h",G%HI, scale=GV%H_to_m) endif - call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, US, larger_h_denom=.true.) if (CS%debug) then call hchksum(tv%T, "after vert_fill_TS tv%T", G%HI, scale=US%C_to_degC) call hchksum(tv%S, "after vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt) @@ -2212,7 +2212,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ call get_param(param_file, mdl, "KD_SMOOTH", CS%Kd_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & From e465b1f6f913c97905b9d8aa9d9aec1110154a01 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 26 Jul 2023 04:55:19 -0400 Subject: [PATCH 12/14] Add comment justifying rescaling in vert_fill_TS Added a comment justifying the use of a fixed rescaling factor for the diffusivity used in vert_fill_TS. All answers and output are identical. --- src/core/MOM_isopycnal_slopes.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 73ae8a9816..29c547148d 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -529,6 +529,11 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke h_neglect = GV%H_subroundoff + ! The use of the fixed rescaling factor in the next line avoids an extra call to thickness_to_dz() + ! and the use of an extra 3-d array of vertical distnaces across layers (dz). This would be more + ! physically consistent, but it would also be more expensive, and given that this routine applies + ! a small (but arbitrary) amount of mixing to clean up the properties of nearly massless layers, + ! the added expense is hard to justify. kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H. h0 = h_neglect if (present(larger_h_denom)) then From 636d6109b9c14f7c5479dc5514913159d7cd97e3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 18 Jul 2023 05:27:01 -0400 Subject: [PATCH 13/14] +*Add and use find_ustar Added the new public interface find_ustar to extract the friction velocity from either a forcing type argument, or a mech_forcing_type argument, either directly or from tau_mag, and in non-Boussinesq mode by using the time-evolving surface specific volume. Find_ustar is an overloaded interface to find_ustar_fluxes or find_ustar_mech_forcing, which are the same but for the type of one of their arguments. For now, the subroutines bulkmixedlayer, mixedlayer_restrajt_OM4, mixedlayer_restrat_Bodner and mixedlayer_restrat_BML are calling find_ustar to avoid code duplication during the transition to work in fully non-Boussinesq mode, but it will eventually be used in about another half dozen other places. All Boussinesq answers are bitwise identical, but non-Boussinesq answers will change and become less dependent on the Boussinesq reference density, and there is a new publicly visible interface wrapping two subroutines. --- src/core/MOM_forcing_type.F90 | 141 +++++++++++++++++- .../lateral/MOM_mixed_layer_restrat.F90 | 36 +++-- .../vertical/MOM_bulk_mixed_layer.F90 | 19 ++- 3 files changed, 182 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9623256b88..a6d35903ee 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -29,7 +29,7 @@ module MOM_forcing_type public extractFluxes1d, extractFluxes2d, optics_type public MOM_forcing_chksum, MOM_mech_forcing_chksum -public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d +public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar public forcing_accumulate, fluxes_accumulate public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type @@ -53,6 +53,12 @@ module MOM_forcing_type module procedure allocate_mech_forcing_from_ref end interface allocate_mech_forcing +!> Determine the friction velocity from a forcing type or a mechanical forcing type. +interface find_ustar + module procedure find_ustar_fluxes + module procedure find_ustar_mech_forcing +end interface find_ustar + ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units @@ -1077,6 +1083,139 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv, end subroutine calculateBuoyancyFlux2d +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(forcing), intent(in) :: fluxes !< Surface fluxes container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] or + !! [H T-1 ~> m s-1 or kg m-2 s-1], depending on H_T_units. + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) & + call MOM_error(FATAL, "find_ustar_fluxes requires that either ustar or tau_mag be associated.") + + if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = fluxes%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * fluxes%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_fluxes + + +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(mech_forcing), intent(in) :: forces !< Surface forces container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(forces%ustar) .or. associated(forces%tau_mag))) & + call MOM_error(FATAL, "find_ustar_mech requires that either ustar or tau_mag be associated.") + + if (associated(forces%ustar) .and. (GV%Boussinesq .or. .not.associated(forces%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = forces%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * forces%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_mech_forcing + + !> Write out chksums for thermodynamic fluxes. subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) character(len=*), intent(in) :: mesg !< message diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 206773ecb0..5b7ec60dee 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -11,7 +11,7 @@ module MOM_mixed_layer_restrat use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_forcing_type, only : mech_forcing +use MOM_forcing_type, only : mech_forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_lateral_mixing_coeffs, only : VarMix_CS @@ -184,6 +184,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! The wind friction velocity, calculated using + ! the Boussinesq reference density or the time-evolving surface density + ! in non-Boussinesq mode [Z T-1 ~> m s-1] MLD_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2] htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av_fast, & ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] @@ -254,6 +257,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// & "The resolution argument, Rd/dx, was not associated.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD. !! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA pRef_MLD(:) = 0. @@ -408,7 +414,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, if (CS%debug) then call hchksum(h,'mixed_layer_restrat: h', G%HI, haloshift=1, scale=GV%H_to_m) - call hchksum(forces%ustar,'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(MLD_fast,'mixed_layer_restrat: MLD', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(Rml_av_fast,'mixed_layer_restrat: rml', G%HI, haloshift=1, & scale=US%m_to_Z*US%L_T_to_m_s**2) @@ -421,7 +427,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! U - Component !$OMP do do j=js,je ; do I=is-1,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -508,7 +514,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! V- component !$OMP do do J=js-1,je ; do i=is,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -711,6 +717,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d wpup ! Turbulent vertical momentum [ ????? ~> m2 s-2] real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] + real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity, calculated using the Boussinesq + ! reference density or the time-evolving surface density in non-Boussinesq + ! mode [Z T-1 ~> m s-1] real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [degC ppt] real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [ppt2] real :: dmu(SZK_(GV)) ! Change in mu(z) across layer k [nondim] @@ -762,12 +771,15 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d call pass_var(bflux, G%domain, halo=1) + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%debug) then call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(BLD, 'mle_Bodner: BLD in', G%HI, haloshift=1, scale=US%Z_to_m) if (associated(bflux)) & call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, scale=US%Z_to_m**2*US%s_to_T**3) - call hchksum(forces%ustar,'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', & G%HI, haloshift=1, scale=US%Z_to_m) call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', & @@ -793,7 +805,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d do j = js-1, je+1 ; do i = is-1, ie+1 w_star3 = max(0., -bflux(i,j)) * BLD(i,j) & ! (this line in Z3 T-3 ~> m3 s-3) * ( ( US%Z_to_m * US%s_to_T )**3 ) ! m3 s-3 - u_star3 = ( US%Z_to_m * US%s_to_T * forces%ustar(i,j) )**3 ! m3 s-3 + u_star3 = ( US%Z_to_m * US%s_to_T * U_star_2d(i,j) )**3 ! m3 s-3 wpup(i,j) = max( CS%min_wstar2, & ! The max() avoids division by zero later ( CS%mstar * u_star3 + CS%nstar * w_star3 )**two_thirds ) & ! (this line m2 s-2) * ( ( US%m_to_Z * US%T_to_s )**2 ) ! Z2 T-2 ~> m2 s-2 @@ -965,7 +977,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d ! Offer diagnostic fields for averaging. if (query_averaging_enabled(CS%diag)) then - if (CS%id_ustar > 0) call post_data(CS%id_ustar, forces%ustar, CS%diag) + if (CS%id_ustar > 0) call post_data(CS%id_ustar, U_star_2d, CS%diag) if (CS%id_bflux > 0) call post_data(CS%id_bflux, bflux, CS%diag) if (CS%id_wpup > 0) call post_data(CS%id_wpup, wpup, CS%diag) if (CS%id_Rml > 0) call post_data(CS%id_Rml, buoy_av, CS%diag) @@ -1053,6 +1065,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! The wind friction velocity, calculated using + ! the Boussinesq reference density or the time-evolving surface density + ! in non-Boussinesq mode [Z T-1 ~> m s-1] htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] @@ -1110,6 +1125,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) if (CS%use_Stanley_ML) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// & "The Stanley parameterization is not available with the BML.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + ! Fix this later for nkml >= 3. p0(:) = 0.0 @@ -1145,7 +1163,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do j=js,je ; do I=is-1,ie h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) @@ -1196,7 +1214,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do J=js-1,je ; do i=is,ie h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 66e2dfa6b2..ceba8dad1a 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -9,7 +9,7 @@ module MOM_bulk_mixed_layer use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_forcing_type, only : extractFluxes1d, forcing +use MOM_forcing_type, only : extractFluxes1d, forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_opacity, only : absorbRemainingSW, optics_type, extract_optics_slice use MOM_unit_scaling, only : unit_scale_type @@ -235,6 +235,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ksort ! The sorted k-index that each original layer goes to. real, dimension(SZI_(G),SZJ_(G)) :: & h_miss ! The summed absolute mismatch [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d ! The wind friction velocity, calculated using the Boussinesq reference density or + ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1] real, dimension(SZI_(G)) :: & TKE, & ! The turbulent kinetic energy available for mixing over a ! time step [Z L2 T-2 ~> m3 s-2]. @@ -412,6 +415,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C max_BL_det(:) = -1 EOSdom(:) = EOS_domain(G%HI) + ! Extract the friction velocity from the forcing type. + call find_ustar(fluxes, tv, U_star_2d, G, GV, US) + !$OMP parallel default(shared) firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) & !$OMP private(h,u,v,h_orig,eps,T,S,opacity_band,d_ea,d_eb,R0,Rcv,ksort, & !$OMP dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,htot,Ttot,Stot,TKE,Conv_en, & @@ -513,7 +519,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! First the TKE at the depth of free convection that is available ! to drive mixing is calculated. - call find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & + call find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) @@ -1252,7 +1258,7 @@ end subroutine mixedlayer_convection !> This subroutine determines the TKE available at the depth of free !! convection to drive mechanical entrainment. -subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & +subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -1265,6 +1271,10 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, type(forcing), intent(in) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields !! have NULL pointers. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: U_star_2d !< The wind friction velocity, calculated + !! using the Boussinesq reference density or + !! the time-evolving surface density in + !! non-Boussinesq mode [Z T-1 ~> m s-1] real, dimension(SZI_(G)), intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy source !! due to free convection [Z L2 T-2 ~> m3 s-2]. real, dimension(SZI_(G)), intent(in) :: dKE_FC !< The vertically integrated change in @@ -1325,7 +1335,8 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega do i=is,ie - U_star = fluxes%ustar(i,j) + U_star = U_star_2d(i,j) + if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & U_star = (1.0 - fluxes%frac_shelf_h(i,j)) * U_star + & From 878fd1ef671456a804c598df606e4cf7c803c203 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 24 Jul 2023 14:14:29 -0400 Subject: [PATCH 14/14] +wave_speed arg mono_N2_depth in thickness units Changed the units of the optional mono_N2_depth argument to wave_speed, wave_speed_init and wave_speed_set_param in thickness units instead of height units. Accordingly, the units of one element each in the diagnostics_CS and wave_speed_CS and a local variable in VarMix_init are also changed to thickness units. The unit descriptions of some comments describing diagnostics were also amended to also describe the non-Boussinesq versions. Because this is essentially just changing when the unit conversion occurs, all answers are bitwise identical, but there are changes to the units of an optional argument in 3 publicly visible routines. --- src/diagnostics/MOM_diagnostics.F90 | 24 +++---- src/diagnostics/MOM_wave_speed.F90 | 64 ++++++++++--------- .../lateral/MOM_lateral_mixing_coeffs.F90 | 4 +- 3 files changed, 49 insertions(+), 43 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index cf8b042c14..157c7268bf 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -56,7 +56,7 @@ module MOM_diagnostics !! monotonic for the purposes of calculating the equivalent !! barotropic wave speed [nondim]. real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of - !! calculating the equivalent barotropic wave speed [Z ~> m]. + !! calculating the equivalent barotropic wave speed [H ~> m or kg m-2]. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. @@ -984,7 +984,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_dKEdt > 0) then - ! Calculate the time derivative of the layer KE [H L2 T-3 ~> m3 s-3]. + ! Calculate the time derivative of the layer KE [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * CS%du_dt(I,j,k) @@ -1006,7 +1006,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_PE_to_KE > 0) then - ! Calculate the potential energy to KE term [H L2 T-3 ~> m3 s-3]. + ! Calculate the potential energy to KE term [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu(I,j,k) @@ -1025,7 +1025,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_BT > 0) then - ! Calculate the barotropic contribution to KE term [H L2 T-3 ~> m3 s-3]. + ! Calculate the barotropic contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt(I,j,k) @@ -1044,7 +1044,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_Coradv > 0) then - ! Calculate the KE source from the combined Coriolis and advection terms [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from the combined Coriolis and advection terms [H L2 T-3 ~> m3 s-3 or W m-2]. ! The Coriolis source should be zero, but is not due to truncation errors. There should be ! near-cancellation of the global integral of this spurious Coriolis source. do k=1,nz @@ -1069,7 +1069,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_adv > 0) then - ! Calculate the KE source from along-layer advection [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from along-layer advection [H L2 T-3 ~> m3 s-3 or W m-2]. ! NOTE: All terms in KE_adv are multiplied by -1, which can easily produce ! negative zeros and may signal a reproducibility issue over land. ! We resolve this by re-initializing and only evaluating over water points. @@ -1098,7 +1098,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_visc > 0) then - ! Calculate the KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc(I,j,k) @@ -1117,7 +1117,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_visc_gl90 > 0) then - ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_gl90(I,j,k) @@ -1136,7 +1136,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_stress > 0) then - ! Calculate the KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_str(I,j,k) @@ -1155,7 +1155,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_horvisc > 0) then - ! Calculate the KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu(I,j,k) @@ -1174,7 +1174,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS endif if (CS%id_KE_dia > 0) then - ! Calculate the KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]. + ! Calculate the KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3 or W m-2]. do k=1,nz do j=js,je ; do I=Isq,Ieq KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_dia(I,j,k) @@ -1594,7 +1594,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", CS%mono_N2_depth, & "The depth below which N2 is limited as monotonic for the "// & "purposes of calculating the equivalent barotropic wave speed.", & - units='m', scale=US%m_to_Z, default=-1.) + units='m', scale=GV%m_to_H, default=-1.) call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, & "The fractional tolerance for finding the wave speeds.", & units="nondim", default=0.001) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 5757e25cd1..c2b671f1c6 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -38,7 +38,7 @@ module MOM_wave_speed !! wave speed [nondim]. This parameter controls the default behavior of !! wave_speed() which can be overridden by optional arguments. real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of - !! calculating the equivalent barotropic wave speed [Z ~> m]. + !! calculating the equivalent barotropic wave speed [H ~> m or kg m-2]. !! If this parameter is negative, this limiting does not occur. !! This parameter controls the default behavior of wave_speed() which !! can be overridden by optional arguments. @@ -81,7 +81,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ !! for the purposes of calculating vertical modal structure [nondim]. real, optional, intent(in) :: mono_N2_depth !< A depth below which N2 is limited as !! monotonic for the purposes of calculating vertical - !! modal structure [Z ~> m]. + !! modal structure [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: modal_structure !< Normalized model structure [nondim] @@ -157,9 +157,11 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ real :: sum_hc ! The sum of the layer thicknesses [Z ~> m] real :: gp ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2] real :: N2min ! A minimum buoyancy frequency, including a slope rescaling factor [L2 Z-2 T-2 ~> s-2] + logical :: below_mono_N2_frac ! True if an interface is below the fractional depth where N2 should not increase. + logical :: below_mono_N2_depth ! True if an interface is below the absolute depth where N2 should not increase. logical :: l_use_ebt_mode, calc_modal_structure real :: l_mono_N2_column_fraction ! A local value of mono_N2_column_fraction [nondim] - real :: l_mono_N2_depth ! A local value of mono_N2_column_depth [Z ~> m] + real :: l_mono_N2_depth ! A local value of mono_N2_column_depth [H ~> m or kg m-2] real :: mode_struct(SZK_(GV)) ! The mode structure [nondim], but it is also temporarily ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6. real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2] @@ -214,18 +216,11 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ c2_scale = US%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results. min_h_frac = tol_Hfrac / real(nz) -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,tv,& -!$OMP calc_modal_structure,l_use_ebt_mode,modal_structure, & -!$OMP l_mono_N2_column_fraction,l_mono_N2_depth,CS, & -!$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale, & -!$OMP better_est,cg1_min2,tol_merge,tol_solve,c2_scale) & -!$OMP private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, & -!$OMP Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT,drho_dS, & -!$OMP drxh_sum,kc,Hc,Hc_H,Tc,Sc,I_Hnew,gprime,& -!$OMP Rc,speed2_tot,Igl,Igu,lam0,lam,lam_it,dlam, & -!$OMP mode_struct,sum_hc,N2min,gp,hw, & -!$OMP ms_min,ms_max,ms_sq,H_top,H_bot,I_Htot,merge, & -!$OMP det,ddet,det_it,ddet_it) + !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,tv,use_EOS, & + !$OMP CS,min_h_frac,calc_modal_structure,l_use_ebt_mode, & + !$OMP modal_structure,l_mono_N2_column_fraction,l_mono_N2_depth, & + !$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale,cg1_min2, & + !$OMP better_est,tol_solve,tol_merge,c2_scale) do j=js,je ! First merge very thin layers with the one above (or below if they are ! at the top). This also transposes the row order so that columns can @@ -335,7 +330,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * max(0.0,Rf(k,i)-Rf(k-1,i)) enddo endif - endif + endif ! use_EOS ! Find gprime across each internal interface, taking care of convective instabilities by ! merging layers. If the estimated wave speed is too small, simply return zero. @@ -452,24 +447,34 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ Igu(1) = 0. ! Neumann condition for pressure modes sum_hc = Hc(1) N2min = gprime(2)/Hc(1) + + below_mono_N2_frac = .false. + below_mono_N2_depth = .false. do k=2,kc hw = 0.5*(Hc(k-1)+Hc(k)) gp = gprime(K) if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then - !### Change to: if ( ((htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) .or. & ) ) - if ( (((G%bathyT(i,j)+G%Z_ref) - sum_hc < l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) .or. & - ((l_mono_N2_depth >= 0.) .and. (sum_hc > l_mono_N2_depth))) .and. & - (gp > N2min*hw) ) then - ! Filters out regions where N2 increases with depth but only in a lower fraction + ! Determine whether N2 estimates should not be allowed to increase with depth. + if (l_mono_N2_column_fraction>0.) then + !### Change to: (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) + below_mono_N2_frac = ((G%bathyT(i,j)+G%Z_ref) - GV%H_to_Z*sum_hc < & + l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) + endif + if (l_mono_N2_depth >= 0.) below_mono_N2_depth = (sum_hc > GV%H_to_Z*l_mono_N2_depth) + + if ( (gp > N2min*hw) .and. (below_mono_N2_frac .or. below_mono_N2_depth) ) then + ! Filters out regions where N2 increases with depth, but only in a lower fraction ! of the water column or below a certain depth. gp = N2min * hw else N2min = gp / hw endif endif + Igu(k) = 1.0/(gp*Hc(k)) Igl(k-1) = 1.0/(gp*Hc(k-1)) sum_hc = sum_hc + Hc(k) + if (better_est) then ! Estimate that the ebt_mode is sqrt(2) times the speed of the flat bottom modes. speed2_tot = speed2_tot + 2.0 * gprime(K)*((H_top(K) * H_bot(K)) * I_Htot) @@ -690,7 +695,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s H_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m] H_bot, & ! The distance of each filtered interface from the bottom [Z ~> m] gprime, & ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. - N2 ! The Brunt Vaissalla freqency squared [T-2 ~> s-2] + N2 ! The buoyancy freqency squared [T-2 ~> s-2] real, dimension(SZK_(GV),SZI_(G)) :: & Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] Tf, & ! Layer temperatures after very thin layers are combined [C ~> degC] @@ -704,7 +709,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s Sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt] Rc, & ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2] - real :: I_Htot ! The inverse of the total filtered thicknesses [Z ~> m] + real :: I_Htot ! The inverse of the total filtered thicknesses [Z-1 ~> m-1] real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and its ! derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim]. ! The exact value should not matter for the final result if it is an even power of 2. @@ -797,6 +802,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! Simplifying the following could change answers at roundoff. Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth) use_EOS = associated(tv%eqn_of_state) + if (CS%c1_thresh < 0.0) & call MOM_error(FATAL, "INTERNAL_WAVE_CG1_THRESH must be set to a non-negative "//& "value via wave_speed_init for wave_speeds to be used.") @@ -978,7 +984,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s I_Hnew = 1.0 / (Hc(kc) + Hc(kc-1)) Tc(kc-1) = (Hc(kc)*Tc(kc) + Hc(kc-1)*Tc(kc-1)) * I_Hnew Sc(kc-1) = (Hc(kc)*Sc(kc) + Hc(kc-1)*Sc(kc-1)) * I_Hnew - Hc(kc-1) = (Hc(kc) + Hc(kc-1)) + Hc(kc-1) = Hc(kc) + Hc(kc-1) kc = kc - 1 else ; exit ; endif enddo @@ -1006,7 +1012,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s if (merge) then ! Merge this layer with the one above and backtrack. Rc(kc) = (Hc(kc)*Rc(kc) + Hf(k,i)*Rf(k,i)) / (Hc(kc) + Hf(k,i)) - Hc(kc) = (Hc(kc) + Hf(k,i)) + Hc(kc) = Hc(kc) + Hf(k,i) ! Backtrack to remove any convective instabilities above... Note ! that the tolerance is a factor of two larger, to avoid limit how ! far back we go. @@ -1019,7 +1025,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s if (merge) then ! Merge the two bottommost layers. At this point kc = k2. Rc(kc-1) = (Hc(kc)*Rc(kc) + Hc(kc-1)*Rc(kc-1)) / (Hc(kc) + Hc(kc-1)) - Hc(kc-1) = (Hc(kc) + Hc(kc-1)) + Hc(kc-1) = Hc(kc) + Hc(kc-1) kc = kc - 1 else ; exit ; endif enddo @@ -1109,7 +1115,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s do k=1,kc w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) ![Z L4 T-4] enddo - renorm = sqrt(htot(i)*a_int/w2avg) ![L-2 T-2] + renorm = sqrt(htot(i)*a_int/w2avg) ! [T2 L-2] do K=1,kc+1 ; mode_struct(K) = renorm * mode_struct(K) ; enddo ! after renorm, mode_struct is again [nondim] @@ -1437,7 +1443,7 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de !! calculating the vertical modal structure [nondim]. real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited !! as monotonic for the purposes of calculating the - !! vertical modal structure [Z ~> m]. + !! vertical modal structure [H ~> m or kg m-2]. logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions !! that recover the remapping answers from 2018. Otherwise !! use more robust but mathematically equivalent expressions. @@ -1489,7 +1495,7 @@ subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_ !! calculating the vertical modal structure [nondim]. real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited !! as monotonic for the purposes of calculating the - !! vertical modal structure [Z ~> m]. + !! vertical modal structure [H ~> m or kg m-2]. logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions !! that recover the remapping answers from 2018. Otherwise !! use more robust but mathematically equivalent expressions. diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e74e48055a..df26f3f6a4 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1103,7 +1103,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) real :: oneOrTwo ! A variable that may be 1 or 2, depending on which form ! of the equatorial deformation radius us used [nondim] real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when - ! calculating the first-mode wave speed [Z ~> m] + ! calculating the first-mode wave speed [H ~> m or kg m-2] real :: KhTr_passivity_coeff ! Coefficient setting the ratio between along-isopycnal tracer ! mixing and interface height mixing [nondim] real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The @@ -1241,7 +1241,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "The depth below which N2 is monotonized to avoid stratification "//& "artifacts from altering the equivalent barotropic mode structure. "//& "This monotonzization is disabled if this parameter is negative.", & - units="m", default=-1.0, scale=US%m_to_Z) + units="m", default=-1.0, scale=GV%m_to_H) allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif