diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9a3bf2d062..52f31014e2 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1103,24 +1103,14 @@ subroutine build_grid_HyCOM1( G, GV, h, tv, dzInterface, remapCS, CS ) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position type(remapping_CS), intent(in) :: remapCS !< Remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure + ! Local variables integer :: i, j, k, nz real, dimension(SZK_(GV)) :: T_col, S_col, p_col, rho_col, h_col_new ! Layer quantities real, dimension(SZK_(GV)+1) :: z_col, z_col_new ! Interface positions relative to the surface in H units (m or kg m-2) real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col in H units (m or kg m-2) - real, dimension(CS%nk,2) :: ppoly_i_E ! Edge value of polynomial - real, dimension(CS%nk,2) :: ppoly_i_S ! Edge slope of polynomial - real, dimension(CS%nk,CS%degree_i+1) :: ppoly_i_coefficients ! Coefficients of polynomial - real :: nominal_z ! Nominal depth of interface is using z* (m or Pa) - real :: stretching ! z* stretching, converts z* to z. - real :: hNew - logical :: maximum_depths_set ! If true, the maximum depths of interface have been set. - logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set. - integer :: ppoly_degree nz = GV%ke - maximum_depths_set = allocated(CS%max_interface_depths) - maximum_h_set = allocated(CS%max_layer_thickness) if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_HyCOM1 : "//& "Target densities must be set before build_grid_HyCOM1 is called.") @@ -1134,49 +1124,8 @@ subroutine build_grid_HyCOM1( G, GV, h, tv, dzInterface, remapCS, CS ) ! once the final grid has been determined). T_col(:) = tv%T(i,j,:) S_col(:) = tv%S(i,j,:) - z_col(1) = 0. ! Work downward rather than bottom up - do K = 1, nz - z_col(K+1) = z_col(K) + h(i,j,k) ! Work in units of h (m or Pa) - p_col(k) = CS%ref_pressure + CS%compressibility_fraction * & - ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure ) - enddo - ! Work bottom recording potential density - call calculate_density(T_col, S_col, p_col, & - rho_col, 1, nz, tv%eqn_of_state ) - ! This ensures the potential density profile is monotonic - ! although not necessarily single valued. - do k = nz-1, 1, -1 - rho_col(k) = min( rho_col(k), rho_col(k+1) ) - enddo - - ! Interpolates for the target interface position with the rho_col profile - call regridding_set_ppolys(rho_col, CS, nz, h(i,j,:), ppoly_i_E, ppoly_i_S, & - ppoly_i_coefficients, ppoly_degree) - ! Based on global density profile, interpolate to generate a new grid - call interpolate_grid(nz, h(i,j,:), z_col, ppoly_i_E, ppoly_i_coefficients, & - CS%target_density, ppoly_degree, nz, h_col_new, z_col_new) - - ! Sweep down the interfaces and make sure that the interface is at least - ! as deep as a nominal target z* grid - nominal_z = 0. - stretching = z_col(nz+1) / G%bathyT(i,j) * GV%m_to_H ! Stretches z* to z - do k = 2, nz+1 - nominal_z = nominal_z + CS%coordinateResolution(k-1) * stretching - z_col_new(k) = max( z_col_new(k), nominal_z ) - z_col_new(k) = min( z_col_new(k), z_col(nz+1) ) - enddo - - if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz - ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. - ! Recall that z_col_new is positive downward. - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), & - z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; elseif (maximum_depths_set) then ; do K=2,nz - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K)) - enddo ; elseif (maximum_h_set) then ; do k=2,nz - z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; endif + build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T_col, S_col, z_col_new) ! Calculate the final change in grid position after blending new and old grids call filtered_grid_motion( CS, nz, z_col, z_col_new, dz_col ) @@ -1192,6 +1141,79 @@ subroutine build_grid_HyCOM1( G, GV, h, tv, dzInterface, remapCS, CS ) end subroutine build_grid_HyCOM1 +subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, z_col_new) + type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options + type(EOS_type), pointer :: eqn_of_state !< Equation of state structure + integer, intent(in) :: nz !< Number of levels + real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, dimension(nz), intent(in) :: T, S !< T and S for column + real, dimension(nz), intent(in) :: h !< Layer thicknesses, in m + real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces + + ! Local variables + integer :: i, j, k, nz + real, dimension(SZK_(GV)) :: p_col, rho_col, h_col_new ! Layer quantities + real, dimension(SZK_(GV)+1) :: z_col ! Interface positions relative to the surface in H units (m or kg m-2) + real, dimension(CS%nk,2) :: ppoly_i_E ! Edge value of polynomial + real, dimension(CS%nk,2) :: ppoly_i_S ! Edge slope of polynomial + real, dimension(CS%nk,CS%degree_i+1) :: ppoly_i_coefficients ! Coefficients of polynomial + real :: stretching ! z* stretching, converts z* to z. + real :: nominal_z ! Nominal depth of interface is using z* (m or Pa) + real :: hNew + logical :: maximum_depths_set ! If true, the maximum depths of interface have been set. + logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set. + integer :: ppoly_degree + + maximum_depths_set = allocated(CS%max_interface_depths) + maximum_h_set = allocated(CS%max_layer_thickness) + + z_col(1) = 0. ! Work downward rather than bottom up + do K = 1, nz + z_col(K+1) = z_col(K) + h(i,j,k) ! Work in units of h (m or Pa) + p_col(k) = CS%ref_pressure + CS%compressibility_fraction * & + ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure ) + enddo + + ! Work bottom recording potential density + call calculate_density(T, S, p_col, rho_col, 1, nz, eqn_of_state) + ! This ensures the potential density profile is monotonic + ! although not necessarily single valued. + do k = nz-1, 1, -1 + rho_col(k) = min( rho_col(k), rho_col(k+1) ) + enddo + + ! Interpolates for the target interface position with the rho_col profile + call regridding_set_ppolys(rho_col, CS, nz, h(i,j,:), ppoly_i_E, ppoly_i_S, & + ppoly_i_coefficients, ppoly_degree) + ! Based on global density profile, interpolate to generate a new grid + call interpolate_grid(nz, h(i,j,:), z_col, ppoly_i_E, ppoly_i_coefficients, & + CS%target_density, ppoly_degree, nz, h_col_new, z_col_new) + + ! Sweep down the interfaces and make sure that the interface is at least + ! as deep as a nominal target z* grid + nominal_z = 0. + stretching = z_col(nz+1) / G%bathyT(i,j) * GV%m_to_H ! Stretches z* to z + do k = 2, nz+1 + nominal_z = nominal_z + CS%coordinateResolution(k-1) * stretching + z_col_new(k) = max( z_col_new(k), nominal_z ) + z_col_new(k) = min( z_col_new(k), z_col(nz+1) ) + enddo + + if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz + ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. + ! Recall that z_col_new is positive downward. + z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), & + z_col_new(K-1) + CS%max_layer_thickness(k-1)) + enddo ; elseif (maximum_depths_set) then ; do K=2,nz + z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K)) + enddo ; elseif (maximum_h_set) then ; do k=2,nz + z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1)) + enddo ; endif + +end subroutine build_hycom1_column + + !> Builds a grid that tracks density interfaces for water that is denser than !! the surface density plus an increment of some number of layers, and uses all !! lighter layers uniformly above this location. Note that this amounts to @@ -1209,6 +1231,7 @@ subroutine build_grid_SLight( G, GV, h, tv, dzInterface, remapCS, CS ) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position type(remapping_CS), intent(in) :: remapCS !< Remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure + ! Local variables real, dimension(SZK_(GV)) :: T_col, S_col, p_col, rho_col ! Layer quantities real, dimension(SZK_(GV)) :: h_col ! A column of layer thicknesses on the original grid in H units. @@ -1268,268 +1291,320 @@ subroutine build_grid_SLight( G, GV, h, tv, dzInterface, remapCS, CS ) h_col(k) = h(i,j,k) enddo - z_col(1) = 0. ! Work downward rather than bottom up - do K=1,nz - z_col(K+1) = z_col(K) + h_col(k) ! Work in units of h (m or Pa) - p_col(k) = CS%ref_pressure + CS%compressibility_fraction * & - ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure ) - enddo + call build_slight_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, z_col_new) - if (z_col(nz+1) - z_col(1) < nz*CS%min_thickness) then - ! This is a nearly massless total depth, so distribute the water evenly. - dz = (z_col(nz+1) - z_col(1)) / real(nz) - do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo - else + ! Calculate the final change in grid position after blending new and old grids + call filtered_grid_motion( CS, nz, z_col, z_col_new, dz_col ) + do K=1,nz+1 ; dzInterface(i,j,K) = -dz_col(K) ; enddo +#ifdef __DO_SAFETY_CHECKS__ + if (dzInterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!' + if (dzInterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!' +#endif - call calculate_density(T_col, S_col, p_col, rho_col, 1, nz, & - tv%eqn_of_state) + ! This adjusts things robust to round-off errors + call adjust_interface_motion( nz, CS%min_thickness, h(i,j,:), dzInterface(i,j,:) ) - ! Find the locations of the target potential densities, flagging - ! locations in apparently unstable regions as not reliable. - call rho_interfaces_col(rho_col, h_col, z_col, CS%target_density, nz, & - z_col_new, CS, reliable, debug=.true.) + else ! on land + dzInterface(i,j,:) = 0. + endif ! mask2dT + enddo; enddo ! i,j - ! Ensure that the interfaces are at least CS%min_thickness apart. - if (CS%min_thickness > 0.0) then - ! Move down interfaces below overly thin layers. - do K=2,nz ; if (z_col_new(K) < z_col_new(K-1) + CS%min_thickness) then - z_col_new(K) = z_col_new(K-1) + CS%min_thickness - endif ; enddo - ! Now move up any interfaces that are too close to the bottom. - do K=nz,2,-1 ; if (z_col_new(K) > z_col_new(K+1) - CS%min_thickness) then - z_col_new(K) = z_col_new(K+1) - CS%min_thickness - else - exit ! No more interfaces can be too close to the bottom. - endif ; enddo - endif +end subroutine build_grid_SLight - ! Fix up the unreliable regions. - kur_ss = 2 ! reliable(1) and reliable(nz+1) must always be true. - do - ! Search for the uppermost unreliable interface postion. - kur1 = nz+2 - do K=kur_ss,nz ; if (.not.reliable(K)) then - kur1 = K ; exit - endif ; enddo - if (kur1 > nz) exit ! Everything is now reliable. +subroutine build_slight_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, z_col_new) + type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options + type(EOS_type), pointer :: eqn_of_state !< Equation of state structure + integer, intent(in) :: nz !< Number of levels + real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, dimension(nz), intent(in) :: T, S !< T and S for column + real, dimension(nz), intent(in) :: h !< Layer thicknesses, in m + real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces - kur2 = kur1-1 ! For error checking. - do K=kur1+1,nz+1 ; if (reliable(K)) then - kur2 = K-1 ; kur_ss = K ; exit - endif ; enddo - if (kur2 < kur1) call MOM_error(FATAL, "Bad unreliable range.") - - dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1) - ! drho = CS%target_density(kur2+1) - CS%target_density(kur1-1) - ! Perhaps reset the wgt and cowgt depending on how bad the old interface - ! locations were. - wgt = 1.0 ; cowgt = 0.0 ! = 1.0-wgt - do K=kur1,kur2 - z_col_new(K) = cowgt*z_col_new(K) + & - wgt * (z_col_new(kur1-1) + dz_ur*(K - (kur1-1)) / ((kur2 - kur1) + 2)) - enddo - enddo + ! Local variables + real, dimension(SZK_(GV)) :: T_col, S_col, p_col, rho_col ! Layer quantities + real, dimension(SZK_(GV)) :: h_col ! A column of layer thicknesses on the original grid in H units. + real, dimension(SZK_(GV)) :: T_f, S_f ! Filtered ayer quantities + real, dimension(SZK_(GV)+1) :: z_col, z_col_new ! Interface positions relative to the surface in H units (m or kg m-2) + real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col in H units (m or kg m-2) + logical, dimension(SZK_(GV)+1) :: reliable ! If true, this interface is in a reliable position. + real, dimension(SZK_(GV)+1) :: T_int, S_int ! Temperature and salinity interpolated to interfaces. + real, dimension(SZK_(GV)+1) :: rho_tmp, drho_dp, p_IS, p_R + real, dimension(SZK_(GV)+1) :: drhoIS_dT, drhoIS_dS + real, dimension(SZK_(GV)+1) :: drhoR_dT, drhoR_dS + real, dimension(SZK_(GV)+1) :: strat_rat + real :: H_to_cPa + real :: drIS, drR, Fn_now, I_HStol, Fn_zero_val + real :: z_int_unst + real :: dz ! A uniform layer thickness in very shallow water, in H. + real :: dz_ur ! The total thickness of an unstable region, in H. + real :: wgt, cowgt ! A weight and its complement, nondim. + real :: rho_ml_av ! The average potential density in a near-surface region, in kg m-3. + real :: H_ml_av ! A thickness to try to use in taking the near-surface average, in H. + real :: rho_x_z ! A cumulative integral of a density, in kg m-3 H. + real :: z_wt ! The thickness actually used in taking the near-surface average, in H. + real :: k_interior ! The (real) value of k where the interior grid starts. + real :: k_int2 ! The (real) value of k where the interior grid starts. + real :: z_interior ! The depth where the interior grid starts, in H. + real :: z_ml_fix ! The depth at which the fixed-thickness near-surface layers end, in H. + real :: dz_dk ! The thickness of layers between the fixed-thickness + ! near-surface layars and the interior, in H. + real :: Lfilt ! A filtering lengthscale, in H. + logical :: maximum_depths_set ! If true, the maximum depths of interface have been set. + logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set. + real :: k2_used, k2here, dz_sum, z_max + integer :: k2 + real :: h_tr, b_denom_1, b1, d1 ! Temporary variables used by the tridiagonal solver. + real, dimension(SZK_(GV)) :: c1 ! Temporary variables used by the tridiagonal solver. + integer :: kur1, kur2 ! The indicies at the top and bottom of an unreliable region. + integer :: kur_ss ! The index to start with in the search for the next unstable region. + integer :: i, j, k, nz, nkml - ! Determine which interfaces are in the s-space region and the depth extent - ! of this region. - z_wt = 0.0 ; rho_x_z = 0.0 - H_ml_av = GV%m_to_H*CS%Rho_ml_avg_depth - do k=1,nz - if (z_wt + h_col(k) >= H_ml_av) then - rho_x_z = rho_x_z + rho_col(k) * (H_ml_av - z_wt) - z_wt = H_ml_av - exit - else - rho_x_z = rho_x_z + rho_col(k) * h_col(k) - z_wt = z_wt + h_col(k) - endif - enddo - if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt - - nkml = CS%nz_fixed_surface - ! Find the interface that matches rho_ml_av. - if (rho_ml_av <= CS%target_density(nkml)) then - k_interior = CS%nlay_ml_offset + real(nkml) - elseif (rho_ml_av > CS%target_density(nz+1)) then - k_interior = real(nz+1) - else ; do K=nkml,nz - if ((rho_ml_av >= CS%target_density(K)) .and. & - (rho_ml_av < CS%target_density(K+1))) then - k_interior = (CS%nlay_ml_offset + K) + & - (rho_ml_av - CS%target_density(K)) / & - (CS%target_density(K+1) - CS%target_density(K)) - exit - endif - enddo ; endif - if (k_interior > real(nz+1)) k_interior = real(nz+1) + z_col(1) = 0. ! Work downward rather than bottom up + do K=1,nz + z_col(K+1) = z_col(K) + h_col(k) ! Work in units of h (m or Pa) + p_col(k) = CS%ref_pressure + CS%compressibility_fraction * & + ( 0.5 * ( z_col(K) + z_col(K+1) ) * GV%H_to_Pa - CS%ref_pressure ) + enddo + + if (z_col(nz+1) - z_col(1) < nz*CS%min_thickness) then + ! This is a nearly massless total depth, so distribute the water evenly. + dz = (z_col(nz+1) - z_col(1)) / real(nz) + do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo + else + + call calculate_density(T_col, S_col, p_col, rho_col, 1, nz, & + tv%eqn_of_state) + + ! Find the locations of the target potential densities, flagging + ! locations in apparently unstable regions as not reliable. + call rho_interfaces_col(rho_col, h_col, z_col, CS%target_density, nz, & + z_col_new, CS, reliable, debug=.true.) + + ! Ensure that the interfaces are at least CS%min_thickness apart. + if (CS%min_thickness > 0.0) then + ! Move down interfaces below overly thin layers. + do K=2,nz ; if (z_col_new(K) < z_col_new(K-1) + CS%min_thickness) then + z_col_new(K) = z_col_new(K-1) + CS%min_thickness + endif ; enddo + ! Now move up any interfaces that are too close to the bottom. + do K=nz,2,-1 ; if (z_col_new(K) > z_col_new(K+1) - CS%min_thickness) then + z_col_new(K) = z_col_new(K+1) - CS%min_thickness + else + exit ! No more interfaces can be too close to the bottom. + endif ; enddo + endif + + ! Fix up the unreliable regions. + kur_ss = 2 ! reliable(1) and reliable(nz+1) must always be true. + do + ! Search for the uppermost unreliable interface postion. + kur1 = nz+2 + do K=kur_ss,nz ; if (.not.reliable(K)) then + kur1 = K ; exit + endif ; enddo + if (kur1 > nz) exit ! Everything is now reliable. + + kur2 = kur1-1 ! For error checking. + do K=kur1+1,nz+1 ; if (reliable(K)) then + kur2 = K-1 ; kur_ss = K ; exit + endif ; enddo + if (kur2 < kur1) call MOM_error(FATAL, "Bad unreliable range.") + + dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1) +! drho = CS%target_density(kur2+1) - CS%target_density(kur1-1) + ! Perhaps reset the wgt and cowgt depending on how bad the old interface + ! locations were. + wgt = 1.0 ; cowgt = 0.0 ! = 1.0-wgt + do K=kur1,kur2 + z_col_new(K) = cowgt*z_col_new(K) + & + wgt * (z_col_new(kur1-1) + dz_ur*(K - (kur1-1)) / ((kur2 - kur1) + 2)) + enddo + enddo + + ! Determine which interfaces are in the s-space region and the depth extent + ! of this region. + z_wt = 0.0 ; rho_x_z = 0.0 + H_ml_av = GV%m_to_H*CS%Rho_ml_avg_depth + do k=1,nz + if (z_wt + h_col(k) >= H_ml_av) then + rho_x_z = rho_x_z + rho_col(k) * (H_ml_av - z_wt) + z_wt = H_ml_av + exit + else + rho_x_z = rho_x_z + rho_col(k) * h_col(k) + z_wt = z_wt + h_col(k) + endif + enddo + if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt + + nkml = CS%nz_fixed_surface + ! Find the interface that matches rho_ml_av. + if (rho_ml_av <= CS%target_density(nkml)) then + k_interior = CS%nlay_ml_offset + real(nkml) + elseif (rho_ml_av > CS%target_density(nz+1)) then + k_interior = real(nz+1) + else ; do K=nkml,nz + if ((rho_ml_av >= CS%target_density(K)) .and. & + (rho_ml_av < CS%target_density(K+1))) then + k_interior = (CS%nlay_ml_offset + K) + & + (rho_ml_av - CS%target_density(K)) / & + (CS%target_density(K+1) - CS%target_density(K)) + exit + endif + enddo ; endif + if (k_interior > real(nz+1)) k_interior = real(nz+1) - ! Linearly interpolate to find z_interior. This could be made more sophisticated. - K = int(ceiling(k_interior)) - z_interior = (K-k_interior)*z_col_new(K-1) + (1.0+(k_interior-K))*z_col_new(K) + ! Linearly interpolate to find z_interior. This could be made more sophisticated. + K = int(ceiling(k_interior)) + z_interior = (K-k_interior)*z_col_new(K-1) + (1.0+(k_interior-K))*z_col_new(K) - if (CS%fix_haloclines) then + if (CS%fix_haloclines) then ! ! Identify regions above the reference pressure where the chosen ! ! potential density significantly underestimates the actual ! ! stratification, and use these to find a second estimate of ! ! z_int_unst and k_interior. - if (CS%halocline_filter_length > 0.0) then - Lfilt = CS%halocline_filter_length*GV%m_to_H - - ! Filter the temperature and salnity with a fixed lengthscale. - h_tr = h_col(1) + GV%H_subroundoff - b1 = 1.0 / (h_tr + Lfilt) ; d1 = h_tr * b1 - T_f(1) = (b1*h_tr)*T_col(1) ; S_f(1) = (b1*h_tr)*S_col(1) - do k=2,nz - c1(k) = Lfilt * b1 - h_tr = h_col(k) + GV%H_subroundoff ; b_denom_1 = h_tr + d1*Lfilt - b1 = 1.0 / (b_denom_1 + Lfilt) ; d1 = b_denom_1 * b1 - T_f(k) = b1 * (h_tr*T_col(k) + Lfilt*T_f(k-1)) - S_f(k) = b1 * (h_tr*S_col(k) + Lfilt*S_f(k-1)) - enddo - do k=nz-1,1,-1 - T_f(k) = T_f(k) + c1(k+1)*T_f(k+1) ; S_f(k) = S_f(k) + c1(k+1)*S_f(k+1) - enddo - else - do k=1,nz ; T_f(k) = T_col(k) ; S_f(k) = S_col(k) ; enddo - endif + if (CS%halocline_filter_length > 0.0) then + Lfilt = CS%halocline_filter_length*GV%m_to_H + + ! Filter the temperature and salnity with a fixed lengthscale. + h_tr = h_col(1) + GV%H_subroundoff + b1 = 1.0 / (h_tr + Lfilt) ; d1 = h_tr * b1 + T_f(1) = (b1*h_tr)*T_col(1) ; S_f(1) = (b1*h_tr)*S_col(1) + do k=2,nz + c1(k) = Lfilt * b1 + h_tr = h_col(k) + GV%H_subroundoff ; b_denom_1 = h_tr + d1*Lfilt + b1 = 1.0 / (b_denom_1 + Lfilt) ; d1 = b_denom_1 * b1 + T_f(k) = b1 * (h_tr*T_col(k) + Lfilt*T_f(k-1)) + S_f(k) = b1 * (h_tr*S_col(k) + Lfilt*S_f(k-1)) + enddo + do k=nz-1,1,-1 + T_f(k) = T_f(k) + c1(k+1)*T_f(k+1) ; S_f(k) = S_f(k) + c1(k+1)*S_f(k+1) + enddo + else + do k=1,nz ; T_f(k) = T_col(k) ; S_f(k) = S_col(k) ; enddo + endif - T_int(1) = T_f(1) ; S_int(1) = S_f(1) - do K=2,nz - T_int(K) = 0.5*(T_f(k-1) + T_f(k)) ; S_int(K) = 0.5*(S_f(k-1) + S_f(k)) - p_IS(K) = z_col(K) * GV%H_to_Pa - p_R(K) = CS%ref_pressure + CS%compressibility_fraction * ( p_IS(K) - CS%ref_pressure ) - enddo - T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz) - p_IS(nz+1) = z_col(nz+1) * GV%H_to_Pa - call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, 2, nz-1, & - tv%eqn_of_state) - call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, & - tv%eqn_of_state) - if (CS%compressibility_fraction > 0.0) then - call calculate_compress(T_int, S_int, p_R, rho_tmp, drho_dp, 2, nz-1, & - tv%eqn_of_state) - else - do K=2,nz ; drho_dp(K) = 0.0 ; enddo - endif - - H_to_cPa = CS%compressibility_fraction*GV%H_to_Pa - strat_rat(1) = 1.0 - do K=2,nz - drIS = drhoIS_dT(K) * (T_f(k) - T_f(k-1)) + & - drhoIS_dS(K) * (S_f(k) - S_f(k-1)) - drR = (drhoR_dT(K) * (T_f(k) - T_f(k-1)) + & - drhoR_dS(K) * (S_f(k) - S_f(k-1))) + & - drho_dp(K) * (H_to_cPa*0.5*(h_col(k) + h_col(k-1))) - - if (drIS <= 0.0) then - strat_rat(K) = 2.0 ! Maybe do this? => ; if (drR < 0.0) strat_rat(K) = -2.0 + T_int(1) = T_f(1) ; S_int(1) = S_f(1) + do K=2,nz + T_int(K) = 0.5*(T_f(k-1) + T_f(k)) ; S_int(K) = 0.5*(S_f(k-1) + S_f(k)) + p_IS(K) = z_col(K) * GV%H_to_Pa + p_R(K) = CS%ref_pressure + CS%compressibility_fraction * ( p_IS(K) - CS%ref_pressure ) + enddo + T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz) + p_IS(nz+1) = z_col(nz+1) * GV%H_to_Pa + call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, 2, nz-1, & + tv%eqn_of_state) + call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, & + tv%eqn_of_state) + if (CS%compressibility_fraction > 0.0) then + call calculate_compress(T_int, S_int, p_R, rho_tmp, drho_dp, 2, nz-1, & + tv%eqn_of_state) + else + do K=2,nz ; drho_dp(K) = 0.0 ; enddo + endif + + H_to_cPa = CS%compressibility_fraction*GV%H_to_Pa + strat_rat(1) = 1.0 + do K=2,nz + drIS = drhoIS_dT(K) * (T_f(k) - T_f(k-1)) + & + drhoIS_dS(K) * (S_f(k) - S_f(k-1)) + drR = (drhoR_dT(K) * (T_f(k) - T_f(k-1)) + & + drhoR_dS(K) * (S_f(k) - S_f(k-1))) + & + drho_dp(K) * (H_to_cPa*0.5*(h_col(k) + h_col(k-1))) + + if (drIS <= 0.0) then + strat_rat(K) = 2.0 ! Maybe do this? => ; if (drR < 0.0) strat_rat(K) = -2.0 + else + strat_rat(K) = 2.0*max(drR,0.0) / (drIS + abs(drR)) + endif + enddo + strat_rat(nz+1) = 1.0 + + z_int_unst = 0.0 ; Fn_now = 0.0 + Fn_zero_val = min(2.0*CS%halocline_strat_tol, & + 0.5*(1.0 + CS%halocline_strat_tol)) + if (CS%halocline_strat_tol > 0.0) then + ! Use Adcroft's reciprocal rule. + I_HStol = 0.0 ; if (Fn_zero_val - CS%halocline_strat_tol > 0.0) & + I_HStol = 1.0 / (Fn_zero_val - CS%halocline_strat_tol) + do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then + z_int_unst = z_int_unst + Fn_now * h_col(k) + if (strat_rat(K) <= Fn_zero_val) then + if (strat_rat(K) <= CS%halocline_strat_tol) then ; Fn_now = 1.0 else - strat_rat(K) = 2.0*max(drR,0.0) / (drIS + abs(drR)) + Fn_now = max(Fn_now, (Fn_zero_val - strat_rat(K)) * I_HStol) endif - enddo - strat_rat(nz+1) = 1.0 - - z_int_unst = 0.0 ; Fn_now = 0.0 - Fn_zero_val = min(2.0*CS%halocline_strat_tol, & - 0.5*(1.0 + CS%halocline_strat_tol)) - if (CS%halocline_strat_tol > 0.0) then - ! Use Adcroft's reciprocal rule. - I_HStol = 0.0 ; if (Fn_zero_val - CS%halocline_strat_tol > 0.0) & - I_HStol = 1.0 / (Fn_zero_val - CS%halocline_strat_tol) - do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then - z_int_unst = z_int_unst + Fn_now * h_col(k) - if (strat_rat(K) <= Fn_zero_val) then - if (strat_rat(K) <= CS%halocline_strat_tol) then ; Fn_now = 1.0 - else - Fn_now = max(Fn_now, (Fn_zero_val - strat_rat(K)) * I_HStol) - endif - endif - endif ; enddo - else - do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then - z_int_unst = z_int_unst + Fn_now * h_col(k) - if (strat_rat(K) <= CS%halocline_strat_tol) Fn_now = 1.0 - endif ; enddo endif + endif ; enddo + else + do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then + z_int_unst = z_int_unst + Fn_now * h_col(k) + if (strat_rat(K) <= CS%halocline_strat_tol) Fn_now = 1.0 + endif ; enddo + endif - if (z_interior < z_int_unst) then - ! Find a second estimate of the extent of the s-coordinate region. - kur1 = max(int(ceiling(k_interior)),2) - if (z_col_new(kur1-1) < z_interior) then - k_int2 = kur1 - do K = kur1,nz+1 ; if (z_col_new(K) >= z_int_unst) then - ! This is linear interpolation again. - if (z_col_new(K-1) >= z_int_unst) & - call MOM_error(FATAL,"build_grid_SLight, bad halocline structure.") - k_int2 = real(K-1) + (z_int_unst - z_col_new(K-1)) / & - (z_col_new(K) - z_col_new(K-1)) - exit - endif ; enddo - if (z_col_new(nz+1) < z_int_unst) then - ! This should be unnecessary. - z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1) - endif - - ! Now take the larger values. - if (k_int2 > k_interior) then - k_interior = k_int2 ; z_interior = z_int_unst - endif - endif + if (z_interior < z_int_unst) then + ! Find a second estimate of the extent of the s-coordinate region. + kur1 = max(int(ceiling(k_interior)),2) + if (z_col_new(kur1-1) < z_interior) then + k_int2 = kur1 + do K = kur1,nz+1 ; if (z_col_new(K) >= z_int_unst) then + ! This is linear interpolation again. + if (z_col_new(K-1) >= z_int_unst) & + call MOM_error(FATAL,"build_grid_SLight, bad halocline structure.") + k_int2 = real(K-1) + (z_int_unst - z_col_new(K-1)) / & + (z_col_new(K) - z_col_new(K-1)) + exit + endif ; enddo + if (z_col_new(nz+1) < z_int_unst) then + ! This should be unnecessary. + z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1) + endif + + ! Now take the larger values. + if (k_int2 > k_interior) then + k_interior = k_int2 ; z_interior = z_int_unst endif - - endif ! fix_haloclines - - z_col_new(1) = 0.0 - do K=2,nkml+1 - z_col_new(K) = min((K-1)*CS%dz_ml_min, & - z_col_new(nz+1) - CS%min_thickness*(nz+1-K)) - enddo - z_ml_fix = z_col_new(nkml+1) - if (z_interior > z_ml_fix) then - dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1)) - do K=nkml+2,int(floor(k_interior)) - z_col_new(K) = z_ml_fix + dz_dk * (K - (nkml+1)) - enddo - else ! The fixed-thickness z-region penetrates into the interior. - do K=nkml+2,nz - if (z_col_new(K) <= z_col_new(CS%nz_fixed_surface+1)) then - z_col_new(K) = z_col_new(CS%nz_fixed_surface+1) - else ; exit ; endif - enddo endif + endif - if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz - ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. - ! Recall that z_col_new is positive downward. - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), & - z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; elseif (maximum_depths_set) then ; do K=2,nz - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K)) - enddo ; elseif (maximum_h_set) then ; do k=2,nz - z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; endif - - endif ! Total thickness exceeds nz*CS%min_thickness. - + endif ! fix_haloclines - ! Calculate the final change in grid position after blending new and old grids - call filtered_grid_motion( CS, nz, z_col, z_col_new, dz_col ) - do K=1,nz+1 ; dzInterface(i,j,K) = -dz_col(K) ; enddo -#ifdef __DO_SAFETY_CHECKS__ - if (dzInterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!' - if (dzInterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!' -#endif + z_col_new(1) = 0.0 + do K=2,nkml+1 + z_col_new(K) = min((K-1)*CS%dz_ml_min, & + z_col_new(nz+1) - CS%min_thickness*(nz+1-K)) + enddo + z_ml_fix = z_col_new(nkml+1) + if (z_interior > z_ml_fix) then + dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1)) + do K=nkml+2,int(floor(k_interior)) + z_col_new(K) = z_ml_fix + dz_dk * (K - (nkml+1)) + enddo + else ! The fixed-thickness z-region penetrates into the interior. + do K=nkml+2,nz + if (z_col_new(K) <= z_col_new(CS%nz_fixed_surface+1)) then + z_col_new(K) = z_col_new(CS%nz_fixed_surface+1) + else ; exit ; endif + enddo + endif - ! This adjusts things robust to round-off errors - call adjust_interface_motion( nz, CS%min_thickness, h(i,j,:), dzInterface(i,j,:) ) + if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz + ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. + ! Recall that z_col_new is positive downward. + z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), & + z_col_new(K-1) + CS%max_layer_thickness(k-1)) + enddo ; elseif (maximum_depths_set) then ; do K=2,nz + z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K)) + enddo ; elseif (maximum_h_set) then ; do k=2,nz + z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1)) + enddo ; endif - else ! on land - dzInterface(i,j,:) = 0. - endif ! mask2dT - enddo; enddo ! i,j + endif ! Total thickness exceeds nz*CS%min_thickness. -end subroutine build_grid_SLight +end subroutine build_slight_column !> Finds the new interface locations in a column of water that match the !! prescribed target densities.