From 2577c8aa246f0933984d1f2a037092c8f6e1272f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 12 Sep 2016 17:48:23 -0400 Subject: [PATCH 1/4] +Revamped MOM_diapyc_energy_req.F90 Updated the test-code for calculating diapycnal mixing energy requirements, and added calls (controlled by the undocumented run-time parameter DEBUG_ENERGY_REQ) within MOM_diabatic_driver to call this code. The energy requiments can now be applied for successive addition of diffusivities. This checkin involves changes to some internal interfaces, but all solutions are bitwise identical. --- .../vertical/MOM_diabatic_driver.F90 | 17 +- .../vertical/MOM_diapyc_energy_req.F90 | 908 +++++++++++++++--- 2 files changed, 804 insertions(+), 121 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 93e2bfb4b6..b6024426dc 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -17,6 +17,8 @@ module MOM_diabatic_driver use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_target_grids use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags +use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end +use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS use MOM_diffConvection, only : diffConvection_CS, diffConvection_init use MOM_diffConvection, only : diffConvection_calculate, diffConvection_end use MOM_domains, only : pass_var, To_West, To_South @@ -78,7 +80,7 @@ module MOM_diabatic_driver !! in the surface boundary layer. logical :: use_kappa_shear !< If true, use the kappa_shear module to find the !! shear-driven diapycnal diffusivity. - logical :: use_cvmix_shear !< If true, use the CVMix module to find the + logical :: use_cvmix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal diffusivity. logical :: use_sponge !< If true, sponges may be applied anywhere in the !! domain. The exact location and properties of @@ -139,6 +141,7 @@ module MOM_diabatic_driver !! is statically unstable. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debugConservation !< If true, monitor conservation and extrema. + logical :: debug_energy_req ! If true, test the mixing energy requirement code. type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output real :: MLDdensityDifference !< Density difference used to determine MLD_user integer :: nsw !< SW_NBANDS @@ -196,6 +199,7 @@ module MOM_diabatic_driver type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(KPP_CS), pointer :: KPP_CSp => NULL() type(diffConvection_CS), pointer :: Conv_CSp => NULL() + type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass @@ -358,6 +362,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) endif if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, G) + if (CS%debug_energy_req) & + call diapyc_energy_req_test(h, dt, tv, G, GV, CS%diapyc_en_rec_CSp) + + call cpu_clock_begin(id_clock_set_diffusivity) call set_BBL_TKE(u, v, h, fluxes, visc, G, GV, CS%set_diff_CSp) call cpu_clock_end(id_clock_set_diffusivity) @@ -1842,6 +1850,9 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, write out verbose debugging data.", default=.false.) call get_param(param_file, mod, "DEBUG_CONSERVATION", CS%debugConservation, & "If true, monitor conservation and extrema.", default=.false.) + + call get_param(param_file, mod, "DEBUG_ENERGY_REQ", CS%debug_energy_req, & + "If true, debug the energy requirements.", default=.false., do_not_log=.true.) call get_param(param_file, mod, "MIX_BOUNDARY_TRACERS", CS%mix_boundary_tracers, & "If true, mix the passive tracers in massless layers at \n"//& "the bottom into the interior as though a diffusivity of \n"//& @@ -2182,6 +2193,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, call regularize_layers_init(Time, G, param_file, diag, CS%regularize_layers_CSp) + if (CS%debug_energy_req) & + call diapyc_energy_req_init(G, param_file, CS%diapyc_en_rec_CSp) ! obtain information about the number of bands for penetrative shortwave if (use_temperature) then @@ -2219,6 +2232,8 @@ subroutine diabatic_driver_end(CS) if (CS%useConvection) call diffConvection_end(CS%Conv_CSp) if (CS%use_energetic_PBL) & call energetic_PBL_end(CS%energetic_PBL_CSp) + if (CS%debug_energy_req) & + call diapyc_energy_req_end(CS%diapyc_en_rec_CSp) if (associated(CS%optics)) then call opacity_end(CS%opacity_CSp, CS%optics) diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index 2713a99a96..bbdeb197c2 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -49,18 +49,23 @@ module MOM_diapyc_energy_req contains -subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV) +subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, CS, Kd_int) type(ocean_grid_type), intent(in) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_3d type(thermo_var_ptrs), intent(inout) :: tv real, intent(in) :: dt + type(diapyc_energy_req_CS), pointer :: CS + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: Kd_int + ! Arguments: h_3d - Layer thickness before entrainment, in m or kg m-2. -! (in/out) tv - A structure containing pointers to any available +! (in) tv - A structure containing pointers to any available ! thermodynamic fields. Absent fields have NULL ptrs. ! (in) dt - The amount of time covered by this call, in s. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. +! (in) CS - This module's control structure +! (in,opt) Kd_int - Interface diffusivities. real, dimension(SZK_(G)) :: & T0, S0, & ! T0 & S0 are columns of initial temperatures and salinities, in degC and g/kg. @@ -75,8 +80,14 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV) logical :: surface_BL, bottom_BL is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + if (.not. associated(CS)) call MOM_error(FATAL, "diapyc_energy_req_test: "// & + "Module must be initialized before it is used.") + !$OMP do do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + if (present(Kd_int)) then + do k=1,nz+1 ; Kd(K) = Kd_int(i,j,K) ; enddo + else htot = 0.0 ; h_top(1) = 0.0 do k=1,nz T0(k) = tv%T(i,j,k) ; S0(k) = tv%S(i,j,k) @@ -95,8 +106,9 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV) Kd(1) = 0.0 ; Kd(nz+1) = 0.0 do K=2,nz tmp1 = h_top(K) * h_bot(K) * GV%H_to_m - Kd(k) = ustar * 0.41 * (tmp1*ustar) / (absf*tmp1 + htot*ustar) + Kd(K) = ustar * 0.41 * (tmp1*ustar) / (absf*tmp1 + htot*ustar) enddo + endif call diapyc_energy_req_calc(h_col, T0, S0, Kd, energy_Kd, dt, tv, G, GV) endif ; enddo ; enddo @@ -133,13 +145,24 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV p_lay, & ! Average pressure of a layer, in Pa. dSV_dT, & ! Partial derivative of specific volume with temperature, in m3 kg-1 K-1. dSV_dS, & ! Partial derivative of specific volume with salinity, in m3 kg-1 / (g kg-1). - T0, S0, & - Te, Se, & + T0, S0, & ! Initial temperatures and salinities. + Te, Se, & ! Running incomplete estimates of the new temperatures and salinities. + Te_a, Se_a, & ! Running incomplete estimates of the new temperatures and salinities. + Te_b, Se_b, & ! Running incomplete estimates of the new temperatures and salinities. + Tf, Sf, & ! New final values of the temperatures and salinities. dTe, dSe, & ! Running (1-way) estimates of temperature and salinity change. + dTe_a, dSe_a, & ! Running (1-way) estimates of temperature and salinity change. + dTe_b, dSe_b, & ! Running (1-way) estimates of temperature and salinity change. + Te_b1_b, Te_b1_a, & + Se_b1_b, Se_b1_a, & dT_to_dPE, & ! Partial derivative of column potential energy with the temperature dS_to_dPE, & ! and salinity changes within a layer, in J m-2 K-1 and J m-2 / (g kg-1). - b_denom_1, & ! The first term in the denominator of b1, in m or kg m-2. - c1, & ! c1 is used by the tridiagonal solver, ND. + b_den_1a, & ! The first term in the denominator of b1 in a downward-oriented + ! tridiagonal solver, in m or kg m-2. + b_den_1b, & ! The first term in the denominator of b1 in an upward-oriented + ! tridiagonal solver, in m or kg m-2. + c1_a, & ! c1_a is used by a downward-oriented tridiagonal solver, ND. + c1_b, & ! c1_b is used by an upward-oriented tridiagonal solver, ND. h_tr ! h_tr is h at tracer points with a h_neglect added to ! ensure positive definiteness, in m or kg m-2. real, dimension(SZK_(G)+1) :: & @@ -148,17 +171,40 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV dS_to_dPE_a, & dT_to_dPE_b, & dS_to_dPE_b, & - PE_chg_k, & ! The integrated potential energy change within a timestep due - ! to the diffusivity at interface K, in J m-2. - PEchg, & - Kddt_h ! The diapycnal diffusivity times a timestep divided by the +! PEchg, & + Kddt_h_a, & ! + Kddt_h_b, & ! + Kddt_h, & ! The diapycnal diffusivity times a timestep divided by the ! average thicknesses around a layer, in m or kg m-2. + Kd_so_far ! The value of Kddt_h that has been applied already in + ! calculating the energy changes, in m or kg m-2. + real, dimension(SZK_(G)+1,4) :: & + PE_chg_k ! The integrated potential energy change within a timestep due + ! to the diffusivity at interface K for 3 different orders of + ! accumulating the diffusivities, in J m-2. real :: & b1 ! b1 is used by the tridiagonal solver, in m-1 or m2 kg-1. + real :: & + I_b1 ! The inverse of b1, in m or kg m-2. +! real :: & +! b1_a, b1_b ! b1_a & b1_b are used by the tridiagonal solver, in m-1 or m2 kg-1. + real :: b_dens, bdt1 + real :: dT_c , dS_c + real :: PEc_core + real :: Kd0, dKd +! real :: d0, dn +! real :: I_d0, I_dn real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected, in m. real :: dTe_term, dSe_term real :: Kddt_h_guess + real :: dMass ! The mass per unit area within a layer, in kg m-2. + real :: dPres ! The hydrostatic pressure change across a layer, in Pa. + real :: dMKE_max ! The maximum amount of mean kinetic energy that could be + ! converted to turbulent kinetic energy if the velocity in + ! the layer below an interface were homogenized with all of + ! the water above the interface, in J m-2 = kg s-2. + real :: PE_change real :: htot real :: dT_k, dT_km1, dS_k, dS_km1 ! Temporary arrays real :: b1Kd ! Temporary arrays @@ -167,25 +213,30 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV ! The following are a bunch of diagnostic arrays for debugging purposes. real, dimension(SZK_(G)) :: & - Ta, Sa + Ta, Sa, Tb, Sb real, dimension(SZK_(G)+1) :: & dPEa_dKd, dPEa_dKd_est, dPEa_dKd_err, dPEa_dKd_trunc, dPEa_dKd_err_norm, & dPEb_dKd, dPEb_dKd_est, dPEb_dKd_err, dPEb_dKd_trunc, dPEb_dKd_err_norm real :: PE_chg_tot1A, PE_chg_tot2A, T_chg_totA real :: PE_chg_tot1B, PE_chg_tot2B, T_chg_totB + real :: PE_chg_tot1C, PE_chg_tot2C, T_chg_totC + real :: PE_chg_tot1D, PE_chg_tot2D, T_chg_totD real, dimension(SZK_(G)+1) :: dPEchg_dKd real :: PE_chg(6) real, dimension(6) :: dT_k_itt, dS_k_itt, dT_km1_itt, dS_km1_itt real :: I_G_Earth real :: dKd_rat_dKd, ddT_k_dKd, ddS_k_dKd, ddT_km1_dKd, ddS_km1_dKd - integer :: k, nz, itt, max_itt - logical :: surface_BL, bottom_BL, debug + integer :: k, nz, itt, max_itt, k_cent + logical :: surface_BL, bottom_BL, central, halves, debug + logical :: old_PE_calc nz = G%ke h_neglect = GV%H_subroundoff I_G_Earth = 1.0 / GV%g_Earth - surface_BL = .true. ; bottom_BL = .true. ; debug = .true. + surface_BL = .true. ; bottom_BL = .true. ; halves = .true. + central = .true. ; K_cent = nz/2 + debug = .true. dPEa_dKd(:) = 0.0 ; dPEa_dKd_est(:) = 0.0 ; dPEa_dKd_err(:) = 0.0 ; dPEa_dKd_err_norm(:) = 0.0 ; dPEa_dKd_trunc(:) = 0.0 dPEb_dKd(:) = 0.0 ; dPEb_dKd_est(:) = 0.0 ; dPEb_dKd_err(:) = 0.0 ; dPEb_dKd_err_norm(:) = 0.0 ; dPEb_dKd_trunc(:) = 0.0 @@ -214,55 +265,94 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV call calculate_specific_vol_derivs(T0, S0, p_lay, dSV_dT, dSV_dS, 1, nz, tv%eqn_of_state) do k=1,nz - dT_to_dPE(k) = 0.5 * I_G_Earth *(pres(K+1)**2 - pres(K)**2) * dSV_dT(k) - dS_to_dPE(k) = 0.5 * I_G_Earth *(pres(K+1)**2 - pres(K)**2) * dSV_dS(k) + dMass = GV%H_to_kg_m2 * h_tr(k) + dPres = GV%g_Earth * dMass + dT_to_dPE(k) = (dMass * (pres(K) + 0.5*dPres)) * dSV_dT(k) + dS_to_dPE(k) = (dMass * (pres(K) + 0.5*dPres)) * dSV_dS(k) +! dT_to_dColHt(k) = dMass * dSV_dT(k) +! dS_to_dColHt(k) = dMass * dSV_dS(k) enddo - PE_chg_k(1) = 0.0 ; PE_chg_k(nz+1) = 0.0 - PEchg(:) = 0.0 ; dPEchg_dKd(:) = 0.0 +! PE_chg_k(1) = 0.0 ; PE_chg_k(nz+1) = 0.0 + ! PEchg(:) = 0.0 + PE_chg_k(:,:) = 0.0 + dPEchg_dKd(:) = 0.0 if (surface_BL) then ! This version is appropriate for a surface boundary layer. + old_PE_calc = .false. + + ! Set up values appropriate for no diffusivity. + do k=1,nz + b_den_1a(k) = h_tr(k) ; b_den_1b(k) = h_tr(k) + dT_to_dPE_a(k) = dT_to_dPE(k) ; dS_to_dPE_a(k) = dS_to_dPE(k) + dT_to_dPE_b(k) = dT_to_dPE(k) ; dS_to_dPE_b(k) = dS_to_dPE(k) + enddo - b_denom_1(1) = h_tr(1) - dT_to_dPE_a(1) = dT_to_dPE(1) - dS_to_dPE_a(1) = dS_to_dPE(1) do K=2,nz ! At this point, Kddt_h(K) will be unknown because its value may depend ! on how much energy is available. ! Precalculate some temporary expressions that are independent of Kddt_h_guess. - if (K==2) then - dT_km1_t2 = (T0(k)-T0(k-1)) - dS_km1_t2 = (S0(k)-S0(k-1)) - dTe_t2 = 0.0 ; dSe_t2 = 0.0 + if (old_PE_calc) then + if (K==2) then + dT_km1_t2 = (T0(k)-T0(k-1)) + dS_km1_t2 = (S0(k)-S0(k-1)) + dTe_t2 = 0.0 ; dSe_t2 = 0.0 + else + dTe_t2 = Kddt_h(K-1) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + dSe_t2 = Kddt_h(K-1) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + dT_km1_t2 = (T0(k)-T0(k-1)) - & + (Kddt_h(K-1) / b_den_1a(k-1)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + dS_km1_t2 = (S0(k)-S0(k-1)) - & + (Kddt_h(K-1) / b_den_1a(k-1)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + endif + dTe_term = dTe_t2 + b_den_1a(k-1) * (T0(k-1)-T0(k)) + dSe_term = dSe_t2 + b_den_1a(k-1) * (S0(k-1)-S0(k)) else - dTe_t2 = Kddt_h(K-1) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) - dSe_t2 = Kddt_h(K-1) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) - dT_km1_t2 = (T0(k)-T0(k-1)) - & - (Kddt_h(K-1) / b_denom_1(k-1)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) - dS_km1_t2 = (S0(k)-S0(k-1)) - & - (Kddt_h(K-1) / b_denom_1(k-1)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + if (K<=2) then + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + else + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2) + Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2) + endif + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) endif - dTe_term = dTe_t2 + b_denom_1(k-1) * (T0(k-1)-T0(k)) - dSe_term = dSe_t2 + b_denom_1(k-1) * (S0(k-1)-S0(k)) ! Find the energy change due to a guess at the strength of diffusion at interface K. Kddt_h_guess = Kddt_h(K) - call find_PE_chg(Kddt_h_guess, h_tr(k), b_denom_1(k-1), & - dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, & - dT_to_dPE(k), dS_to_dPE(k), dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & - PE_chg_k(k), dPEa_dKd(k)) + if (old_PE_calc) then + call find_PE_chg_orig(Kddt_h_guess, h_tr(k), b_den_1a(k-1), & + dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, & + dT_to_dPE(k), dS_to_dPE(k), dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg_k(k,1), dPEa_dKd(k)) + else + call find_PE_chg(0.0, Kddt_h_guess, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg_k(K,1), dPEc_dKd=dPEa_dKd(K)) + endif if (debug) then do itt=1,5 Kddt_h_guess = (1.0+0.01*(itt-3))*Kddt_h(K) - call find_PE_chg(Kddt_h_guess, h_tr(k), b_denom_1(k-1), & - dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, & - dT_to_dPE(k), dS_to_dPE(k), dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & - PE_chg=PE_chg(itt)) + if (old_PE_calc) then + call find_PE_chg_orig(Kddt_h_guess, h_tr(k), b_den_1a(k-1), & + dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, & + dT_to_dPE(k), dS_to_dPE(k), dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg(itt)) + else + call find_PE_chg(0.0, Kddt_h_guess, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg(itt)) + endif enddo ! Compare with a 4th-order finite difference estimate. dPEa_dKd_est(k) = (4.0*(PE_chg(4)-Pe_chg(2))/(0.02*Kddt_h(K)) - & @@ -277,8 +367,8 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV ! At this point, the final value of Kddt_h(K) is known, so the estimated ! properties for layer k-1 can be calculated. - b1 = 1.0 / (b_denom_1(k-1) + Kddt_h(K)) - c1(K) = Kddt_h(K) * b1 + b1 = 1.0 / (b_den_1a(k-1) + Kddt_h(K)) + c1_a(K) = Kddt_h(K) * b1 if (k==2) then Te(1) = b1*(h_tr(1)*T0(1)) Se(1) = b1*(h_tr(1)*S0(1)) @@ -286,84 +376,120 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV Te(k-1) = b1 * (h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2)) Se(k-1) = b1 * (h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2)) endif - dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) - dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) + if (old_PE_calc) then + dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) + dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) + endif - b_denom_1(k) = h_tr(k) + (b_denom_1(k-1) * b1) * Kddt_h(K) - dT_to_dPE_a(k) = dT_to_dPE(k) + c1(K)*dT_to_dPE_a(k-1) - dS_to_dPE_a(k) = dS_to_dPE(k) + c1(K)*dS_to_dPE_a(k-1) + b_den_1a(k) = h_tr(k) + (b_den_1a(k-1) * b1) * Kddt_h(K) + dT_to_dPE_a(k) = dT_to_dPE(k) + c1_a(K)*dT_to_dPE_a(k-1) + dS_to_dPE_a(k) = dS_to_dPE(k) + c1_a(K)*dS_to_dPE_a(k-1) enddo - b1 = 1.0 / (b_denom_1(nz)) - Te(nz) = b1 * (h_tr(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) - Se(nz) = b1 * (h_tr(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) - dTe(nz) = b1 * Kddt_h(nz) * ((T0(nz-1)-T0(nz)) + dTe(nz-1)) - dSe(nz) = b1 * Kddt_h(nz) * ((S0(nz-1)-S0(nz)) + dSe(nz-1)) + b1 = 1.0 / (b_den_1a(nz)) + Tf(nz) = b1 * (h_tr(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) + Sf(nz) = b1 * (h_tr(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) + if (old_PE_calc) then + dTe(nz) = b1 * Kddt_h(nz) * ((T0(nz-1)-T0(nz)) + dTe(nz-1)) + dSe(nz) = b1 * Kddt_h(nz) * ((S0(nz-1)-S0(nz)) + dSe(nz-1)) + endif do k=nz-1,1,-1 - Te(k) = Te(k) + c1(K+1)*Te(k+1) - Se(k) = Se(k) + c1(K+1)*Se(k+1) + Tf(k) = Te(k) + c1_a(K+1)*Tf(k+1) + Sf(k) = Se(k) + c1_a(K+1)*Sf(k+1) enddo if (debug) then - do k=1,nz ; Ta(k) = Te(k) ; Sa(k) = Se(k) ; enddo + do k=1,nz ; Ta(k) = Tf(k) ; Sa(k) = Sf(k) ; enddo PE_chg_tot1A = 0.0 ; PE_chg_tot2A = 0.0 ; T_chg_totA = 0.0 do k=1,nz - PE_chg_tot1A = PE_chg_tot1A + (dT_to_dPE(k) * (Te(k) - T0(k)) + & - dS_to_dPE(k) * (Se(k) - S0(k))) - T_chg_totA = T_chg_totA + h_tr(k) * (Te(k) - T0(k)) + PE_chg_tot1A = PE_chg_tot1A + (dT_to_dPE(k) * (Tf(k) - T0(k)) + & + dS_to_dPE(k) * (Sf(k) - S0(k))) + T_chg_totA = T_chg_totA + h_tr(k) * (Tf(k) - T0(k)) enddo do K=2,nz - PE_chg_tot2A = PE_chg_tot2A + PE_chg_k(K) + PE_chg_tot2A = PE_chg_tot2A + PE_chg_k(K,1) enddo endif endif if (bottom_BL) then ! This version is appropriate for a bottom boundary layer. + old_PE_calc = .false. + + ! Set up values appropriate for no diffusivity. + do k=1,nz + b_den_1a(k) = h_tr(k) ; b_den_1b(k) = h_tr(k) + dT_to_dPE_a(k) = dT_to_dPE(k) ; dS_to_dPE_a(k) = dS_to_dPE(k) + dT_to_dPE_b(k) = dT_to_dPE(k) ; dS_to_dPE_b(k) = dS_to_dPE(k) + enddo - b_denom_1(nz) = h_tr(nz) - dT_to_dPE_b(nz) = dT_to_dPE(nz) - dS_to_dPE_b(nz) = dS_to_dPE(nz) do K=nz,2,-1 ! Loop over interior interfaces. ! At this point, Kddt_h(K) will be unknown because its value may depend ! on how much energy is available. ! Precalculate some temporary expressions that are independent of Kddt_h_guess. - if (K==nz) then - dT_k_t2 = (T0(k-1)-T0(k)) - dS_k_t2 = (S0(k-1)-S0(k)) - dTe_t2 = 0.0 ; dSe_t2 = 0.0 + if (old_PE_calc) then + if (K==nz) then + dT_k_t2 = (T0(k-1)-T0(k)) + dS_k_t2 = (S0(k-1)-S0(k)) + dTe_t2 = 0.0 ; dSe_t2 = 0.0 + else + dTe_t2 = Kddt_h(K+1) * ((T0(k+1) - T0(k)) + dTe(k+1)) + dSe_t2 = Kddt_h(K+1) * ((S0(k+1) - S0(k)) + dSe(k+1)) + dT_k_t2 = (T0(k-1)-T0(k)) - & + (Kddt_h(k+1)/ b_den_1b(k)) * ((T0(k+1) - T0(k)) + dTe(k+1)) + dS_k_t2 = (S0(k-1)-S0(k)) - & + (Kddt_h(k+1)/ b_den_1b(k)) * ((S0(k+1) - S0(k)) + dSe(k+1)) + endif + dTe_term = dTe_t2 + b_den_1b(k) * (T0(k)-T0(k-1)) + dSe_term = dSe_t2 + b_den_1b(k) * (S0(k)-S0(k-1)) else - dTe_t2 = Kddt_h(K+1) * ((T0(k+1) - T0(k)) + dTe(k+1)) - dSe_t2 = Kddt_h(K+1) * ((S0(k+1) - S0(k)) + dSe(k+1)) - dT_k_t2 = (T0(k-1)-T0(k)) - & - (Kddt_h(k+1)/ b_denom_1(k)) * ((T0(k+1) - T0(k)) + dTe(k+1)) - dS_k_t2 = (S0(k-1)-S0(k)) - & - (Kddt_h(k+1)/ b_denom_1(k)) * ((S0(k+1) - S0(k)) + dSe(k+1)) + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + if (K>=nz) then + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + else + Te_b1_b(k) = h_tr(k) * T0(k) + Kddt_h(K+1) * Te(k+1) + Se_b1_b(k) = h_tr(k) * S0(k) + Kddt_h(k+1) * Se(k+1) + endif endif - dTe_term = dTe_t2 + b_denom_1(k) * (T0(k)-T0(k-1)) - dSe_term = dSe_t2 + b_denom_1(k) * (S0(k)-S0(k-1)) ! Find the energy change due to a guess at the strength of diffusion at interface K. - Kddt_h_guess = Kddt_h(K) - call find_PE_chg(Kddt_h_guess, h_tr(k-1), b_denom_1(k), & - dTe_term, dSe_term, dT_k_t2, dS_k_t2, & - dT_to_dPE(k-1), dS_to_dPE(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & - PE_chg_k(K), dPEb_dKd(K)) + if (old_PE_calc) then + call find_PE_chg_orig(Kddt_h_guess, h_tr(k-1), b_den_1b(k), & + dTe_term, dSe_term, dT_k_t2, dS_k_t2, & + dT_to_dPE(k-1), dS_to_dPE(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg_k(K,2), dPEc_dKd=dPEb_dKd(K)) + else + call find_PE_chg(0.0, Kddt_h_guess, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg_k(K,2), dPEc_dKd=dPEb_dKd(K)) + endif if (debug) then ! Compare with a 4th-order finite difference estimate. do itt=1,5 Kddt_h_guess = (1.0+0.01*(itt-3))*Kddt_h(K) - call find_PE_chg(Kddt_h_guess, h_tr(k-1), b_denom_1(k), & + if (old_PE_calc) then + call find_PE_chg_orig(Kddt_h_guess, h_tr(k-1), b_den_1b(k), & dTe_term, dSe_term, dT_k_t2, dS_k_t2, & dT_to_dPE(k-1), dS_to_dPE(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & PE_chg=PE_chg(itt)) + else + call find_PE_chg(0.0, Kddt_h_guess, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_chg(itt)) + endif enddo dPEb_dKd_est(k) = (4.0*(PE_chg(4)-Pe_chg(2))/(0.02*Kddt_h(K)) - & @@ -378,8 +504,8 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV ! At this point, the final value of Kddt_h(K) is known, so the estimated ! properties for layer k can be calculated. - b1 = 1.0 / (b_denom_1(k) + Kddt_h(K)) - c1(K) = Kddt_h(K) * b1 + b1 = 1.0 / (b_den_1b(k) + Kddt_h(K)) + c1_b(K) = Kddt_h(K) * b1 if (k==nz) then Te(nz) = b1* (h_tr(nz)*T0(nz)) Se(nz) = b1* (h_tr(nz)*S0(nz)) @@ -387,59 +513,551 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV Te(k) = b1 * (h_tr(k) * T0(k) + Kddt_h(K+1) * Te(k+1)) Se(k) = b1 * (h_tr(k) * S0(k) + Kddt_h(k+1) * Se(k+1)) endif - dTe(k) = b1 * ( Kddt_h(K)*(T0(k-1)-T0(k)) + dTe_t2 ) - dSe(k) = b1 * ( Kddt_h(K)*(S0(k-1)-S0(k)) + dSe_t2 ) + if (old_PE_calc) then + dTe(k) = b1 * ( Kddt_h(K)*(T0(k-1)-T0(k)) + dTe_t2 ) + dSe(k) = b1 * ( Kddt_h(K)*(S0(k-1)-S0(k)) + dSe_t2 ) + endif - b_denom_1(k-1) = h_tr(k-1) + (b_denom_1(k) * b1) * Kddt_h(K) - dT_to_dPE_b(k-1) = dT_to_dPE(k-1) + c1(K)*dT_to_dPE_b(k) - dS_to_dPE_b(k-1) = dS_to_dPE(k-1) + c1(K)*dS_to_dPE_b(k) + b_den_1b(k-1) = h_tr(k-1) + (b_den_1b(k) * b1) * Kddt_h(K) + dT_to_dPE_b(k-1) = dT_to_dPE(k-1) + c1_b(K)*dT_to_dPE_b(k) + dS_to_dPE_b(k-1) = dS_to_dPE(k-1) + c1_b(K)*dS_to_dPE_b(k) enddo - b1 = 1.0 / (b_denom_1(1)) - Te(1) = b1 * (h_tr(1) * T0(1) + Kddt_h(2) * Te(2)) - Se(1) = b1 * (h_tr(1) * S0(1) + Kddt_h(2) * Se(2)) - dTe(1) = b1 * Kddt_h(2) * ((T0(2)-T0(1)) + dTe(2)) - dSe(1) = b1 * Kddt_h(2) * ((S0(2)-S0(1)) + dSe(2)) + b1 = 1.0 / (b_den_1b(1)) + Tf(1) = b1 * (h_tr(1) * T0(1) + Kddt_h(2) * Te(2)) + Sf(1) = b1 * (h_tr(1) * S0(1) + Kddt_h(2) * Se(2)) + if (old_PE_calc) then + dTe(1) = b1 * Kddt_h(2) * ((T0(2)-T0(1)) + dTe(2)) + dSe(1) = b1 * Kddt_h(2) * ((S0(2)-S0(1)) + dSe(2)) + endif do k=2,nz - Te(k) = Te(k) + c1(K)*Te(k-1) - Se(k) = Se(k) + c1(K)*Se(k-1) + Tf(k) = Te(k) + c1_b(K)*Tf(k-1) + Sf(k) = Se(k) + c1_b(K)*Sf(k-1) enddo if (debug) then + do k=1,nz ; Tb(k) = Tf(k) ; Sb(k) = Sf(k) ; enddo PE_chg_tot1B = 0.0 ; PE_chg_tot2B = 0.0 ; T_chg_totB = 0.0 do k=1,nz - PE_chg_tot1B = PE_chg_tot1B + (dT_to_dPE(k) * (Te(k) - T0(k)) + & - dS_to_dPE(k) * (Se(k) - S0(k))) - T_chg_totB = T_chg_totB + h_tr(k) * (Te(k) - T0(k)) + PE_chg_tot1B = PE_chg_tot1B + (dT_to_dPE(k) * (Tf(k) - T0(k)) + & + dS_to_dPE(k) * (Sf(k) - S0(k))) + T_chg_totB = T_chg_totB + h_tr(k) * (Tf(k) - T0(k)) + enddo + do K=2,nz + PE_chg_tot2B = PE_chg_tot2B + PE_chg_k(K,2) + enddo + endif + + endif + + if (central) then + + ! Set up values appropriate for no diffusivity. + do k=1,nz + b_den_1a(k) = h_tr(k) ; b_den_1b(k) = h_tr(k) + dT_to_dPE_a(k) = dT_to_dPE(k) ; dS_to_dPE_a(k) = dS_to_dPE(k) + dT_to_dPE_b(k) = dT_to_dPE(k) ; dS_to_dPE_b(k) = dS_to_dPE(k) + enddo + + ! Calculate the dependencies on layers above. + Kddt_h_a(1) = 0.0 + do K=2,nz ! Loop over interior interfaces. + ! First calculate some terms that are independent of the change in Kddt_h(K). + Kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K). + if (K<=2) then + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + else + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te_a(k-2) + Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se_a(k-2) + endif + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + + Kddt_h_a(K) = 0.0 + if (K < K_cent) Kddt_h_a(K) = Kddt_h(K) + + dKd = Kddt_h_a(K) + + call find_PE_chg(Kd0, dKd, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_change) + PE_chg_k(K,3) = PE_change + + b1 = 1.0 / (b_den_1a(k-1) + Kddt_h_a(K)) + c1_a(K) = Kddt_h_a(K) * b1 + if (k==2) then + Te_a(1) = b1*(h_tr(1)*T0(1)) + Se_a(1) = b1*(h_tr(1)*S0(1)) + else + Te_a(k-1) = b1 * (h_tr(k-1) * T0(k-1) + Kddt_h_a(K-1) * Te_a(k-2)) + Se_a(k-1) = b1 * (h_tr(k-1) * S0(k-1) + Kddt_h_a(K-1) * Se_a(k-2)) + endif + + b_den_1a(k) = h_tr(k) + (b_den_1a(k-1) * b1) * Kddt_h_a(K) + dT_to_dPE_a(k) = dT_to_dPE(k) + c1_a(K)*dT_to_dPE_a(k-1) + dS_to_dPE_a(k) = dS_to_dPE(k) + c1_a(K)*dS_to_dPE_a(k-1) + enddo + + ! Calculate the dependencies on layers below. + Kddt_h_b(nz+1) = 0.0 + do K=nz,2,-1 ! Loop over interior interfaces. + ! First calculate some terms that are independent of the change in Kddt_h(K). + Kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K). +! if (K<=2) then + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) +! else +! Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te_a(k-2) +! Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se_a(k-2) +! endif + if (K>=nz) then + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + else + Te_b1_b(k) = h_tr(k) * T0(k) + Kddt_h(K+1) * Te_b(k+1) + Se_b1_b(k) = h_tr(k) * S0(k) + Kddt_h(k+1) * Se_b(k+1) + endif + + Kddt_h_b(K) = 0.0 ; if (K > K_cent) Kddt_h_b(K) = Kddt_h(K) + dKd = Kddt_h_b(K) + + call find_PE_chg(Kd0, dKd, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_change) + PE_chg_k(K,3) = PE_chg_k(K,3) + PE_change + + b1 = 1.0 / (b_den_1b(k) + Kddt_h_b(K)) + c1_b(K) = Kddt_h_b(K) * b1 + if (k==nz) then + Te_b(k) = b1 * (h_tr(k)*T0(k)) + Se_b(k) = b1 * (h_tr(k)*S0(k)) + else + Te_b(k) = b1 * (h_tr(k) * T0(k) + Kddt_h_b(K+1) * Te_b(k+1)) + Se_b(k) = b1 * (h_tr(k) * S0(k) + Kddt_h_b(k+1) * Se_b(k+1)) + endif + + b_den_1b(k-1) = h_tr(k-1) + (b_den_1b(k) * b1) * Kddt_h_b(K) + dT_to_dPE_b(k-1) = dT_to_dPE(k-1) + c1_b(K)*dT_to_dPE_b(k) + dS_to_dPE_b(k-1) = dS_to_dPE(k-1) + c1_b(K)*dS_to_dPE_b(k) + + enddo + + ! Calculate the final solution for the layers surrounding interface K_cent + K = K_cent + + ! First calculate some terms that are independent of the change in Kddt_h(K). + Kd0 = 0.0 ! This might need to be changed - it is the already applied value of Kddt_h(K). + if (K<=2) then + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + else + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) + Kddt_h(K-1) * Te_a(k-2) + Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + Kddt_h(K-1) * Se_a(k-2) + endif + if (K>=nz) then + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + else + Te_b1_b(k) = h_tr(k) * T0(k) + Kddt_h(K+1) * Te_b(k+1) + Se_b1_b(k) = h_tr(k) * S0(k) + Kddt_h(k+1) * Se_b(k+1) + endif + + dKd = Kddt_h(K) + + call find_PE_chg(Kd0, dKd, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_change) + PE_chg_k(K,3) = PE_chg_k(K,3) + PE_change + + + ! To derive the following, first to a partial update for the estimated + ! temperatures and salinities in the layers around this interface: + ! b1_a = 1.0 / (b_den_1a(k-1) + Kddt_h(K)) + ! b1_b = 1.0 / (b_den_1b(k) + Kddt_h(K)) + ! Te_up(k) = Te_b1_b(k) * b1_b ; Se_up(k) = Se_b1_b(k) * b1_b + ! Te_up(k-1) = Te_b1_a(k-1) * b1_a ; Se_up(k-1) = Se_b1_a(k-1) * b1_a + ! Find the final values of T & S for both layers around K_cent, using that + ! c1_a(K) = Kddt_h(K) * b1_a ; c1_b(K) = Kddt_h(K) * b1_b + ! Tf(K_cent-1) = Te_up(K_cent-1) + c1_a(K_cent)*Tf(K_cent) + ! Tf(K_cent) = Te_up(K_cent) + c1_b(K_cent)*Tf(K_cent-1) + ! but further exploiting the expressions for c1_a and c1_b to avoid + ! subtraction in the denominator, and use only a single division. + b1 = 1.0 / (b_den_1a(k-1)*b_den_1b(k) + Kddt_h(K)*(b_den_1a(k-1) + b_den_1b(k))) + Tf(k-1) = ((b_den_1b(k) + Kddt_h(K)) * Te_b1_a(k-1) + & + Kddt_h(K) * Te_b1_b(k) ) * b1 + Sf(k-1) = ((b_den_1b(k) + Kddt_h(K)) * Se_b1_a(k-1) + & + Kddt_h(K) * Se_b1_b(k) ) * b1 + Tf(k) = (Kddt_h(K) * Te_b1_a(k-1) + & + (b_den_1a(k-1) + Kddt_h(K)) * Te_b1_b(k) ) * b1 + Sf(k) = (Kddt_h(K) * Se_b1_a(k-1) + & + (b_den_1a(k-1) + Kddt_h(K)) * Se_b1_b(k) ) * b1 + + c1_a(K) = Kddt_h(K) / (b_den_1a(k-1) + Kddt_h(K)) + c1_b(K) = Kddt_h(K) / (b_den_1b(k) + Kddt_h(K)) + + ! Now update the other layer working outward from k_cent to determine the final + ! temperatures and salinities. + do k=K_cent-2,1,-1 + Tf(k) = Te_a(k) + c1_a(K+1)*Tf(k+1) + Sf(k) = Se_a(k) + c1_a(K+1)*Sf(k+1) + enddo + do k=K_cent+1,nz + Tf(k) = Te_b(k) + c1_b(K)*Tf(k-1) + Sf(k) = Se_b(k) + c1_b(K)*Sf(k-1) + enddo + + if (debug) then + PE_chg_tot1C = 0.0 ; PE_chg_tot2C = 0.0 ; T_chg_totC = 0.0 + do k=1,nz + PE_chg_tot1C = PE_chg_tot1C + (dT_to_dPE(k) * (Tf(k) - T0(k)) + & + dS_to_dPE(k) * (Sf(k) - S0(k))) + T_chg_totC = T_chg_totC + h_tr(k) * (Tf(k) - T0(k)) enddo do K=2,nz - PE_chg_tot2B = PE_chg_tot2B + PE_chg_k(K) + PE_chg_tot2C = PE_chg_tot2C + PE_chg_k(K,3) enddo endif endif - energy_Kd = 0.0 ; do K=2,nz ; energy_Kd = energy_Kd + PE_chg_k(K) ; enddo + if (halves) then + + ! Set up values appropriate for no diffusivity. + do k=1,nz + b_den_1a(k) = h_tr(k) ; b_den_1b(k) = h_tr(k) + dT_to_dPE_a(k) = dT_to_dPE(k) ; dS_to_dPE_a(k) = dS_to_dPE(k) + dT_to_dPE_b(k) = dT_to_dPE(k) ; dS_to_dPE_b(k) = dS_to_dPE(k) + enddo + do K=1,nz+1 + Kd_so_far(K) = 0.0 + enddo + + ! Calculate the dependencies on layers above. + Kddt_h_a(1) = 0.0 + do K=2,nz ! Loop over interior interfaces. + ! First calculate some terms that are independent of the change in Kddt_h(K). + Kd0 = Kd_so_far(K) + if (K<=2) then + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) ; Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + else + Te_b1_a(k-1) = h_tr(k-1) * T0(k-1) + Kd_so_far(K-1) * Te(k-2) + Se_b1_a(k-1) = h_tr(k-1) * S0(k-1) + Kd_so_far(K-1) * Se(k-2) + endif + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + + dKd = 0.5 * Kddt_h(K) - Kd_so_far(K) + + call find_PE_chg(Kd0, dKd, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_change) + + PE_chg_k(K,4) = PE_change + + Kd_so_far(K) = Kd_so_far(K) + dKd + + b1 = 1.0 / (b_den_1a(k-1) + Kd_so_far(K)) + c1_a(K) = Kd_so_far(K) * b1 + if (k==2) then + Te(1) = b1*(h_tr(1)*T0(1)) + Se(1) = b1*(h_tr(1)*S0(1)) + else + Te(k-1) = b1 * (h_tr(k-1) * T0(k-1) + Kd_so_far(K-1) * Te(k-2)) + Se(k-1) = b1 * (h_tr(k-1) * S0(k-1) + Kd_so_far(K-1) * Se(k-2)) + endif + + b_den_1a(k) = h_tr(k) + (b_den_1a(k-1) * b1) * Kd_so_far(K) + dT_to_dPE_a(k) = dT_to_dPE(k) + c1_a(K)*dT_to_dPE_a(k-1) + dS_to_dPE_a(k) = dS_to_dPE(k) + c1_a(K)*dS_to_dPE_a(k-1) + enddo + + ! Calculate the dependencies on layers below. + do K=nz,2,-1 ! Loop over interior interfaces. + ! First calculate some terms that are independent of the change in Kddt_h(K). + Kd0 = Kd_so_far(K) + if (K>=nz) then + Te_b1_b(k) = h_tr(k) * T0(k) ; Se_b1_b(k) = h_tr(k) * S0(k) + else + Te_b1_b(k) = h_tr(k) * T0(k) + Kd_so_far(K+1) * Te(k+1) + Se_b1_b(k) = h_tr(k) * S0(k) + Kd_so_far(k+1) * Se(k+1) + endif + + dKd = Kddt_h(K) - Kd_so_far(K) + + call find_PE_chg(Kd0, dKd, b_den_1a(k-1), b_den_1b(k), & + Te_b1_a(k-1), Se_b1_a(k-1), Te_b1_b(k), Se_b1_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + 0.0, 0.0, 0.0, 0.0, 0.0, & + PE_chg=PE_change) + + PE_chg_k(K,4) = PE_chg_k(K,4) + PE_change + + Kd_so_far(K) = Kd_so_far(K) + dKd + + b1 = 1.0 / (b_den_1b(k) + Kd_so_far(K)) + c1_b(K) = Kd_so_far(K) * b1 + if (k==nz) then + Te(k) = b1 * (h_tr(k)*T0(k)) + Se(k) = b1 * (h_tr(k)*S0(k)) + else + Te(k) = b1 * (h_tr(k) * T0(k) + Kd_so_far(K+1) * Te(k+1)) + Se(k) = b1 * (h_tr(k) * S0(k) + Kd_so_far(k+1) * Se(k+1)) + endif + + b_den_1b(k-1) = h_tr(k-1) + (b_den_1b(k) * b1) * Kd_so_far(K) + dT_to_dPE_b(k-1) = dT_to_dPE(k-1) + c1_b(K)*dT_to_dPE_b(k) + dS_to_dPE_b(k-1) = dS_to_dPE(k-1) + c1_b(K)*dS_to_dPE_b(k) + + enddo + + ! Now update the other layer working down from the top to determine the + ! final temperatures and salinities. + b1 = 1.0 / (b_den_1b(1)) + Tf(1) = b1 * (h_tr(1) * T0(1) + Kddt_h(2) * Te(2)) + Sf(1) = b1 * (h_tr(1) * S0(1) + Kddt_h(2) * Se(2)) + do k=2,nz + Tf(k) = Te(k) + c1_b(K)*Tf(k-1) + Sf(k) = Se(k) + c1_b(K)*Sf(k-1) + enddo + + if (debug) then + PE_chg_tot1D = 0.0 ; PE_chg_tot2D = 0.0 ; T_chg_totD = 0.0 + do k=1,nz + PE_chg_tot1D = PE_chg_tot1D + (dT_to_dPE(k) * (Tf(k) - T0(k)) + & + dS_to_dPE(k) * (Sf(k) - S0(k))) + T_chg_totD = T_chg_totD + h_tr(k) * (Tf(k) - T0(k)) + enddo + do K=2,nz + PE_chg_tot2D = PE_chg_tot2D + PE_chg_k(K,4) + enddo + endif + + endif + + energy_Kd = 0.0 ; do K=2,nz ; energy_Kd = energy_Kd + PE_chg_k(K,1) ; enddo energy_Kd = energy_Kd / dt K=nz end subroutine diapyc_energy_req_calc -subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interfaces's diapycnal diffusivity times a timestep. +subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & + dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, & + pres, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & + PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0) + real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface, in units of H (m or kg-2). + real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface, in units of H (m or kg-2). + real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers, in K H. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers, in K H. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers, in K H. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers, in K H. + real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers above, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers above, in J m-2 ppt-1. + real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers below, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers below, in J m-2 ppt-1. + real, intent(in) :: pres !< The hydrostatic interface pressure, which is used to relate + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing, in Pa. + real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above, in m K-1. + real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above, in m ppt-1. + real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below, in m K-1. + real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below, in m ppt-1. + + real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying + !! Kddt_h at the present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h, + !! in units of J m-2 H-1. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realizedd by applying a huge value of Kddt_h at the + !! present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the + !! limit where Kddt_h = 0, in J m-2 H-1. + + real :: hps ! The sum of the two effective pivot thicknesses, in H. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term, in H2. + real :: dT_c ! The core term in the expressions for the temperature changes, in K H2. + real :: dS_c ! The core term in the expressions for the salinity changes, in psu H2. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes, J m-3. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes, J m-3. + real :: ColHt_chg ! The change in the column height, in m. + real :: y1 ! A local temporary term, in units of H-3 or H-4 in various contexts. + + ! The expression for the change in potential energy used here is derived + ! from the expression for the final estimates of the changes in temperature + ! and salinities, and then extensively manipulated to get it into its most + ! succint form. The derivation is not necessarily obvious, but it demonstably + ! works by comparison with separate calculations of the energy changes after + ! the tridiagonal solver for the final changes in temperature and salinity are + ! applied. + + hps = hp_a + hp_b + bdt1 = hp_a * hp_b + Kddt_h0 * hps + dT_c = hp_a * Th_b - hp_b * Th_a + dS_c = hp_a * Sh_b - hp_b * Sh_a + PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & + hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) + ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & + hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) + + if (present(PE_chg)) then + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKddt_h. + y1 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + PE_chg = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) PE_chg = PE_chg - pres * ColHt_chg + endif + + if (present(dPEc_dKd)) then + ! Find the derivative of the potential energy change with dKddt_h. + y1 = 1.0 / (bdt1 + dKddt_h * hps)**2 + dPEc_dKd = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres * ColHt_chg + endif + + if (present(dPE_max)) then + ! This expression is the limit of PE_chg for infinite dKddt_h. + y1 = 1.0 / (bdt1 * hps) + dPE_max = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPE_max = dPE_max - pres * ColHt_chg + endif + + if (present(dPEc_dKd_0)) then + ! This expression is the limit of dPEc_dKd for dKddt_h = 0. + y1 = 1.0 / bdt1**2 + dPEc_dKd_0 = PEc_core * y1 + ColHt_chg = ColHt_core * y1 + if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres * ColHt_chg + endif + +end subroutine find_PE_chg + + +subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, & - dT_to_dPEa, dS_to_dPEa, PE_chg, dPE_chg_dKd) - real, intent(in) :: Kddt_h, h_k, b_den_1, dTe_term, dSe_term - real, intent(in) :: dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k - real, intent(in) :: dT_to_dPEa, dS_to_dPEa - real, optional, intent(out) :: PE_chg - real, optional, intent(out) :: dPE_chg_dKd + dT_to_dPEa, dS_to_dPEa, pres, dT_to_dColHt_k, & + dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, & + PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0) + real, intent(in) :: Kddt_h !< The diffusivity at an interface times the time step and + !! divided by the average of the thicknesses around the + !! interface, in units of H (m or kg-2). + real, intent(in) :: h_k !< The thickness of the layer below the interface, in H. + real, intent(in) :: b_den_1 !< The first term in the denominator of the pivot + !! for the tridiagonal solver, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above, in H. + real, intent(in) :: dTe_term !< A diffusivity-independent term related to the + !! temperature change in the layer below the interface, in K H. + real, intent(in) :: dSe_term !< A diffusivity-independent term related to the + !! salinity change in the layer below the interface, in ppt H. + real, intent(in) :: dT_km1_t2 !< A diffusivity-independent term related to the + !! temperature change in the layer above the interface, in K. + real, intent(in) :: dS_km1_t2 !< A diffusivity-independent term related to the + !! salinity change in the layer above the interface, in ppt. + real, intent(in) :: pres !< The hydrostatic interface pressure, which is used to relate + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing, in Pa. + real, intent(in) :: dT_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers below, in J m-2 K-1. + real, intent(in) :: dS_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers below, in J m-2 ppt-1. + real, intent(in) :: dT_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the temperatures of all the layers above, in J m-2 K-1. + real, intent(in) :: dS_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column + !! potential energy, including all implicit diffusive changes + !! in the salinities of all the layers above, in J m-2 ppt-1. + real, intent(in) :: dT_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below, in m K-1. + real, intent(in) :: dS_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below, in m ppt-1. + real, intent(in) :: dT_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above, in m K-1. + real, intent(in) :: dS_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above, in m ppt-1. + + real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying + !! Kddt_h at the present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h, + !! in units of J m-2 H-1. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realizedd by applying a huge value of Kddt_h at the + !! present interface, in J m-2. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the + !! limit where Kddt_h = 0, in J m-2 H-1. ! This subroutine determines the total potential energy change due to mixing ! at an interface, including all of the implicit effects of the prescribed ! mixing at interfaces above. Everything here is derived by careful manipulation -! of the robust tridiagonal solvers used for tracers by MOM6. +! of the robust tridiagonal solvers used for tracers by MOM6. The results are +! positive for mixing in a stably stratified environment. ! The comments describing these arguments are for a downward mixing pass, but ! this routine can also be used for an upward pass with the sense of direction ! reversed. @@ -460,28 +1078,52 @@ subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & ! temperature change in the layer above the interface, in K. ! (in) dS_km1_t2 - A diffusivity-independent term related to the ! salinity change in the layer above the interface, in ppt. -! (in) dT_to_dPE_k - A factor (1/2*pres*mass_lay*dSpec_vol/dT) relating +! (in) dT_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating ! a layer's temperature change to the change in column ! potential energy, in J m-2 K-1. -! (in) dS_to_dPE_k - A factor (1/2*pres*mass_lay*dSpec_vol/dS) relating +! (in) dS_to_dPE_k - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating ! a layer's salinity change to the change in column ! potential energy, in J m-2 ppt-1. -! (in) dT_to_dPEa - A factor (1/2*pres*mass_lay*dSpec_vol/dT) relating +! (in) dT_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dT) relating ! a layer's temperature change to the change in column ! potential energy, including all implicit diffusive changes ! in the temperatures of all the layers above, in J m-2 K-1. -! (in) dS_to_dPEa - A factor (1/2*pres*mass_lay*dSpec_vol/dS) relating +! (in) dS_to_dPEa - A factor (pres_lay*mass_lay*dSpec_vol/dS) relating ! a layer's salinity change to the change in column ! potential energy, including all implicit diffusive changes ! in the salinities of all the layers above, in J m-2 ppt-1. +! (in) pres - The hydrostatic interface pressure, which is used to relate +! the changes in column thickness to the energy that is radiated +! as gravity waves and unavailable to drive mixing, in Pa. +! (in) dT_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dT) relating +! a layer's temperature change to the change in column +! height, in m K-1. +! (in) dS_to_dColHt_k - A factor (mass_lay*dSColHtc_vol/dS) relating +! a layer's salinity change to the change in column +! height, in m ppt-1. +! (in) dT_to_dColHta - A factor (mass_lay*dSColHtc_vol/dT) relating +! a layer's temperature change to the change in column +! height, including all implicit diffusive changes +! in the temperatures of all the layers above, in m K-1. +! (in) dS_to_dColHta - A factor (mass_lay*dSColHtc_vol/dS) relating +! a layer's salinity change to the change in column +! height, including all implicit diffusive changes +! in the salinities of all the layers above, in m ppt-1. ! (out,opt) PE_chg - The change in column potential energy from applying -! Kddt_h at the present interface, in J m-2. +! Kddt_h at the present interface, in J m-2. ! (out,opt) dPE_chg_dKd - The partial derivative of PE_chg with Kddt_h, -! in units of J m-2 H-1. - +! in units of J m-2 H-1. +! (out,opt) PE_max - The maximum change in column potential energy that could +! be realizedd by applying a huge value of Kddt_h at the +! present interface, in J m-2. +! (out,opt) dPEc_dKc0 - The partial derivative of PE_chg with Kddt_h in the +! limit where Kddt_h = 0, in J m-2 H-1. ! b_den_1 - The first term in the denominator of b1, in m or kg m-2. real :: b1 ! b1 is used by the tridiagonal solver, in H-1. real :: b1Kd ! Temporary array (nondim.) + real :: ColHt_chg ! The change in column thickness in m. + real :: dColHt_max ! The change in column thickess for infinite diffusivity, in m. + real :: dColHt_dKd ! The partial derivative of column thickess with diffusivity, in s m-1. real :: dT_k, dT_km1 ! Temporary arrays in K. real :: dS_k, dS_km1 ! Temporary arrays in ppt. real :: I_Kr_denom, dKr_dKd ! Temporary arrays in H-2 and nondim. @@ -509,9 +1151,12 @@ subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & PE_chg = (dT_to_dPE_k * dT_k + dT_to_dPEa * dT_km1) + & (dS_to_dPE_k * dS_k + dS_to_dPEa * dS_km1) + ColHt_chg = (dT_to_dColHt_k * dT_k + dT_to_dColHta * dT_km1) + & + (dS_to_dColHt_k * dS_k + dS_to_dColHta * dS_km1) + if (ColHt_chg < 0.0) PE_chg = PE_chg - pres * ColHt_chg endif - if (present(dPE_chg_dKd)) then + if (present(dPEc_dKd)) then ! Find the derivatives of the temperature and salinity changes with Kddt_h. dKr_dKd = (h_k*b_den_1) * I_Kr_denom**2 @@ -521,11 +1166,34 @@ subroutine find_PE_chg(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & ddS_km1_dKd = (b1**2 * b_den_1) * ( dS_k + dS_km1_t2 ) + b1Kd * ddS_k_dKd ! Calculate the partial derivative of Pe_chg with Kddt_h. - dPE_chg_dKd = (dT_to_dPE_k * ddT_k_dKd + dT_to_dPEa * ddT_km1_dKd) + & - (dS_to_dPE_k * ddS_k_dKd + dS_to_dPEa * ddS_km1_dKd) + dPEc_dKd = (dT_to_dPE_k * ddT_k_dKd + dT_to_dPEa * ddT_km1_dKd) + & + (dS_to_dPE_k * ddS_k_dKd + dS_to_dPEa * ddS_km1_dKd) + dColHt_dKd = (dT_to_dColHt_k * ddT_k_dKd + dT_to_dColHta * ddT_km1_dKd) + & + (dS_to_dColHt_k * ddS_k_dKd + dS_to_dColHta * ddS_km1_dKd) + if (dColHt_dKd < 0.0) dPEc_dKd = dPEc_dKd - pres * dColHt_dKd endif -end subroutine find_PE_chg + if (present(dPE_max)) then + ! This expression is the limit of PE_chg for infinite Kddt_h. + dPE_max = (dT_to_dPEa * dT_km1_t2 + dS_to_dPEa * dS_km1_t2) + & + ((dT_to_dPE_k + dT_to_dPEa) * dTe_term + & + (dS_to_dPE_k + dS_to_dPEa) * dSe_term) / (b_den_1 + h_k) + dColHt_max = (dT_to_dColHta * dT_km1_t2 + dS_to_dColHta * dS_km1_t2) + & + ((dT_to_dColHt_k + dT_to_dColHta) * dTe_term + & + (dS_to_dColHt_k + dS_to_dColHta) * dSe_term) / (b_den_1 + h_k) + if (dColHt_max < 0.0) dPE_max = dPE_max - pres*dColHt_max + endif + + if (present(dPEc_dKd_0)) then + ! This expression is the limit of dPEc_dKd for Kddt_h = 0. + dPEc_dKd_0 = (dT_to_dPEa * dT_km1_t2 + dS_to_dPEa * dS_km1_t2) / (b_den_1) + & + (dT_to_dPE_k * dTe_term + dS_to_dPE_k * dSe_term) / (h_k*b_den_1) + dColHt_dKd = (dT_to_dColHta * dT_km1_t2 + dS_to_dColHta * dS_km1_t2) / (b_den_1) + & + (dT_to_dColHt_k * dTe_term + dS_to_dColHt_k * dSe_term) / (h_k*b_den_1) + if (dColHt_dKd < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres*dColHt_dKd + endif + +end subroutine find_PE_chg_orig subroutine diapyc_energy_req_init(G, param_file, CS) type(ocean_grid_type), intent(in) :: G From 36275aed478d21de3c9051cf7155a2a0a50ea2b7 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 13 Sep 2016 14:47:36 -0400 Subject: [PATCH 2/4] *Altered initial conditions for circle_obcs - The circle_obcs experiment ostensibly tests the Flather open boundary conditions except it was modified to use 10 layers and also test the Orlanski boundary conditions too. The initialization code appears to have been written for a 2-layer system and essentially only forces a very high-mode. This change produces a mixed mode initial condition (thickness anomaly linear with z) but with most energy in the first baroclinic mode. - Changes circle_obcs experiment only. --- src/user/circle_obcs_initialization.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index 80e326674f..e16cb73ae5 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -93,8 +93,16 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file) ! if (rad <= 6.*diskrad) h(i,j,k) = h(i,j,k)+10.0*exp( -0.5*( rad**2 ) ) rad = min( rad, 1. ) ! Flatten outside radius of diskrad rad = rad*(2.*asin(1.)) ! Map 0-1 to 0-pi - h(i,j,k) = h(i,j,k) + 10.0*0.5*(1.+cos(rad)) ! cosine bell - h(i,j,k-1) = h(i,j,k-1) - 10.0*0.5*(1.+cos(rad)) ! cosine bell + if (Nz==1) then + ! The model is barotropic + h(i,j,k) = h(i,j,k) + 1.0*0.5*(1.+cos(rad)) ! cosine bell + else + ! The model is baroclinic + do k = 1, Nz + h(i,j,k) = h(i,j,k) - 0.5*(1.+cos(rad)) & ! cosine bell + * 5.0 * real( 2*k-Nz ) + enddo + endif enddo ; enddo end subroutine circle_obcs_initialize_thickness From 0aef287ab102d1d047372ea9d72c79dc6f3417d7 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 13 Sep 2016 10:56:27 -0400 Subject: [PATCH 3/4] Fixed extract_word() for n>number of words - extract_word("a,b,c",",",5) was returning "c" when it should have returned "". It turns out we had no instances where we used this or tested for it. - Fixed bug. - Added use case in unit_tests. - No answer changes. --- src/framework/MOM_string_functions.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index e5e5df681b..64c14ad213 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -255,7 +255,7 @@ end function extractWord endif endif enddo - if (b<=ns) extract_word = trim(string(b:ns)) + if (b<=ns .and. nw==n-1) extract_word = trim(string(b:ns)) end function extract_word !> Returns string with all spaces removed. @@ -298,6 +298,8 @@ logical function string_functions_unit_tests() call localTest(extractWord("One Two,Three",3),"Three") call localTest(extractWord("One Two, Three",3),"Three") call localTest(extractWord(" One Two,Three",1),"One") + call localTest(extract_word("One,Two,Three",",",3),"Three") + call localTest(extract_word("One,Two,Three",",",4),"") write(*,*) '==========================================================' contains subroutine localTest(str1,str2) From fcdb84edcdc7d25e627113c69e1117b7da2a8cb6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 14 Sep 2016 14:34:16 -0400 Subject: [PATCH 4/4] Fixes needed to "not" extend OBCs beyond domain. - Added %Flather logical to segment type to differentiate Flather (for barotropic) from Orlanski (baroclinic) boundary conditions. - Bugfix: Barotropic-direction south had wrong sign for setting the tangential ubt. Does not appear to change answers? - Changed syntax of OBC_SEGMENT_%%% string to allow multiple boundary conditions type per segment, e.g. FLATHER,ORLANSKI - zonal_mass_flux() and meridional_mass_flux() in continuity_PPM now use %Flather instead of %radiation just before reconciling fluxes with the barotropic fluxes. - Fixed need to extend boundaries beyond domain in run-time parameters. Extended domains used to be needed to avoid masking the corner vorticity points. Added code to detect corner joins of segments and handle explicitly. - Removed code to interpret I=-1:N+1 (extend boundaries no longer needed). - Use "hack" in setup_*_point_obc() to extend segments by one point (in code) to recover "corner" behavior of old OBCs. - Added back setting of tangential velocity components in continuity_PPM. - set_Flather_data() has commented out code to set "h" across corners when segments are not extended. - Cleaned up setting of "direction" in setup_u_point_obc() and setup_v_point_obc(). - Simplied using of directions/flags in continuity_PPM. - No answer changes with consistent switch to SEGMENT parameters. TODO: - Allow coexistence of specified and radiative BCs within zonal_mass_flux() and meridional_mass_flux(). - Understand why the above bugfix did not change answers? - Figure out why extending segments to recover corner behavior is necessary: LAPLACIAN=BIHARMONIC=False still failed, suggesting PV term needs examination. --- src/core/MOM_barotropic.F90 | 10 +- src/core/MOM_continuity_PPM.F90 | 72 +++-- src/core/MOM_legacy_barotropic.F90 | 10 +- src/core/MOM_open_boundary.F90 | 408 ++++++++++++++++------------- 4 files changed, 283 insertions(+), 217 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index dd209fdb9f..abddf27738 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2378,7 +2378,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, uhbt(I,j) = BT_OBC%uhbt(I,j) ubt(I,j) = BT_OBC%ubt_outer(I,j) vel_trans = ubt(I,j) - elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 @@ -2411,7 +2411,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif vel_trans = ubt(I,j) elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S) then - if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then + if ((vbt(i,J)+vbt(i+1,J)) < 0.0) then ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) else ubt(I,j) = BT_OBC%ubt_outer(I,j) @@ -2438,7 +2438,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, vhbt(i,J) = BT_OBC%vhbt(i,J) vbt(i,J) = BT_OBC%vbt_outer(i,J) vel_trans = vbt(i,J) - elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(I,j) ! CFL v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl<1 @@ -2537,7 +2537,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) .and. & associated(BT_OBC%OBC_mask_u)) then do j=js,je ; do I=is-1,ie ; if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt(I-1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 @@ -2564,7 +2564,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south) .and. & associated(BT_OBC%OBC_mask_v)) then do J=js-1,je ; do i=is,ie ; if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL v_inlet = cfl*vbt(i,J-1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 91a544cd03..662fa6608a 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -191,22 +191,26 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, call cpu_clock_end(id_clock_update) if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west) then - do k=1,nz ; do j=LB%jsh,LB%jeh - do I=LB%ish,LB%ieh+1 + do k=1,nz + do j=LB%jsh,LB%jeh ; do I=LB%ish,LB%ieh+1 if (OBC%OBC_segment_u(I-1,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I-1,j))%radiation & - .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & + if (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E) & h(i,j,k) = h_input(i-1,j,k) endif enddo do i=LB%ish-1,LB%ieh if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation & - .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) & h(i,j,k) = h_input(i+1,j,k) endif - enddo - enddo ; enddo + enddo ; enddo + do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E)) & + v(i,J,k) = v(i-1,J,k) + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W)) & + v(i,J,k) = v(i+1,J,k) + enddo ; enddo + enddo endif LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec @@ -228,18 +232,22 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, do k=1,nz do J=LB%jsh,LB%jeh+1 ; do i=LB%ish-1,LB%ieh+1 if (OBC%OBC_segment_v(i,J-1) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J-1))%radiation & - .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & + if (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N) & h(i,j,k) = h_input(i,j-1,k) endif enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation & - .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) & h(i,j,k) = h_input(i,j+1,k) endif enddo ; enddo + do j=LB%jsh-1,LB%jeh+1 ; do I=LB%ish-1,LB%ieh + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N)) & + u(I,j,k) = u(I,j-1,k) + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S)) & + u(I,j,k) = u(I,j+1,k) + enddo ; enddo enddo endif else ! .not. x_first @@ -261,18 +269,22 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, do k=1,nz do J=LB%jsh,LB%jeh+1 ; do i=LB%ish-1,LB%ieh+1 if (OBC%OBC_segment_v(i,J-1) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J-1))%radiation & - .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & + if (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N) & h(i,j,k) = h_input(i,j-1,k) endif enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation & - .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) & h(i,j,k) = h_input(i,j+1,k) endif enddo ; enddo + do j=LB%jsh-1,LB%jeh+1 ; do I=LB%ish-1,LB%ieh + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N)) & + u(I,j,k) = u(I,j-1,k) + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S)) & + u(I,j,k) = u(I,j+1,k) + enddo ; enddo enddo endif @@ -292,22 +304,26 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, call cpu_clock_end(id_clock_update) if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west) then - do k=1,nz ; do j=LB%jsh,LB%jeh - do I=LB%ish,LB%ieh+1 + do k=1,nz + do j=LB%jsh,LB%jeh ; do I=LB%ish,LB%ieh+1 if (OBC%OBC_segment_u(I-1,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I-1,j))%radiation & - .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & + if (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E) & h(i,j,k) = h_input(i-1,j,k) endif enddo do i=LB%ish-1,LB%ieh if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation & - .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) & h(i,j,k) = h_input(i+1,j,k) endif - enddo - enddo ; enddo + enddo ; enddo + do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E)) & + v(i,J,k) = v(i-1,J,k) + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W)) & + v(i,J,k) = v(i+1,J,k) + enddo ; enddo + enddo endif endif @@ -524,8 +540,10 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) if (.not.do_i(I)) any_simple_OBC = .true. enddo ; else if (apply_OBC_flather) then ; do I=ish-1,ieh + ! This is a tangential condition and is needed for unknown reasons and + ! probably implies that we made a calculation elsewhere that we should not have. do_i(I) = .not.(OBC%OBC_mask_u(I,j) .and. & - (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather) .and. & ((OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N) .or. & (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S))) enddo ; else ; do I=ish-1,ieh @@ -1290,8 +1308,10 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) if (.not.do_i(i)) any_simple_OBC = .true. enddo ; else if (apply_OBC_flather) then ; do i=ish,ieh + ! This is a tangential condition and is needed for unknown reasons and + ! probably implies that we made a calculation elsewhere that we should not have. do_i(i) = .not.(OBC%OBC_mask_v(i,J) .and. & - (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather) .and. & ((OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E) .or. & (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W))) enddo ; else ; do i=ish,ieh diff --git a/src/core/MOM_legacy_barotropic.F90 b/src/core/MOM_legacy_barotropic.F90 index 4cd5b21ccc..5b70381459 100644 --- a/src/core/MOM_legacy_barotropic.F90 +++ b/src/core/MOM_legacy_barotropic.F90 @@ -2236,7 +2236,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, uhbt(I,j) = BT_OBC%uhbt(I,j) ubt(I,j) = BT_OBC%ubt_outer(I,j) vel_trans = ubt(I,j) - elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 @@ -2269,7 +2269,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif vel_trans = ubt(I,j) elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S) then - if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then + if ((vbt(i,J)+vbt(i+1,J)) < 0.0) then ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) else ubt(I,j) = BT_OBC%ubt_outer(I,j) @@ -2296,7 +2296,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, vhbt(i,J) = BT_OBC%vhbt(i,J) vbt(i,J) = BT_OBC%vbt_outer(i,J) vel_trans = vbt(i,J) - elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(I,j) ! CFL v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl<1 @@ -2395,7 +2395,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) .and. & associated(BT_OBC%OBC_mask_u)) then do j=js,je ; do I=is-1,ie ; if (BT_OBC%OBC_mask_u(I,j)) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%Flather) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt(I-1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 @@ -2422,7 +2422,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south) .and. & associated(BT_OBC%OBC_mask_v)) then do J=js-1,je ; do i=is,ie ; if (BT_OBC%OBC_mask_v(i,J)) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%Flather) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL v_inlet = cfl*vbt(i,J-1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 31cddec1ad..c38b442092 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1,4 +1,3 @@ -! This file is part of MOM6. See LICENSE.md for the license. !> Controls where open boundary conditions are applied module MOM_open_boundary @@ -43,7 +42,9 @@ module MOM_open_boundary !> Open boundary segment type - we'll have one for each open segment !! to describe that segment. type, public :: OBC_segment_type - logical :: radiation !< Radiation boundary. + logical :: Flather !< If true, applies Flather + Chapman radiation of barotropic gravity waves. + logical :: radiation !< If true, 1D Orlanksi radiation boundary conditions are applied. + !! If False, a gradient condition is applied. logical :: radiation2D !< Oblique waves supported at radiation boundary. logical :: nudged !< Optional supplement to radiation boundary. logical :: specified !< Boundary fixed to external value. @@ -193,6 +194,7 @@ subroutine open_boundary_config(G, param_file, OBC) ! Allocate everything allocate(OBC%OBC_segment_list(0:OBC%number_of_segments)) do l=0,OBC%number_of_segments + OBC%OBC_segment_list(l)%Flather = .false. OBC%OBC_segment_list(l)%radiation = .false. OBC%OBC_segment_list(l)%radiation2D = .false. OBC%OBC_segment_list(l)%nudged = .false. @@ -250,93 +252,102 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) integer, intent(in) :: l_seg !< which segment is this? ! Local variables integer :: I_obc, Js_obc, Je_obc ! Position of segment in global index space - integer :: j, this_kind - character(len=32) :: action_str + integer :: j, this_kind, a_loop + character(len=32) :: action_str(5) ! This returns the global indices for the segment call parse_segment_str(G%ieg, G%jeg, segment_str, I_obc, Js_obc, Je_obc, action_str ) I_obc = I_obc - G%idg_offset ! Convert to local tile indices on this tile Js_obc = Js_obc - G%jdg_offset ! Convert to local tile indices on this tile Je_obc = Je_obc - G%jdg_offset ! Convert to local tile indices on this tile + this_kind = OBC_NONE -! if (Je_obc>Js_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E -! if (Je_obcJs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E - if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA - if (Je_obcJs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E - if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA - if (Je_obcG%HI%IedB) return ! Boundary is not on tile - if (max(Js_obc,Je_obc)G%HI%JedB) return ! Segment is not on tile - - ! These four lines extend the open boundary into the halo region of tiles on the edge of the physical - ! domain. They are used to reproduce the checksums of the circle_obcs test case and will be removed - ! in the fullness of time. -AJA -! These were causing grief in the supercritical problem. - KSH -! if (Js_obc == G%HI%JscB) Js_obc = G%HI%jsd-1 -! if (Js_obc == G%HI%JecB) Js_obc = G%HI%jed -! if (Je_obc == G%HI%JscB) Je_obc = G%HI%jsd-1 -! if (Je_obc == G%HI%JecB) Je_obc = G%HI%jed - - do j=G%HI%jsd, G%HI%jed - if (j>min(Js_obc,Je_obc) .and. j<=max(Js_obc,Je_obc)) then - OBC%OBC_mask_u(I_obc,j) = .true. - OBC%OBC_segment_u(I_obc,j) = l_seg - if (Je_obc>Js_obc) then ! East is outward - if (this_kind == OBC_FLATHER) then - OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_E ! We only use direction for Flather (maybe) - ! Set v points outside segment - OBC%OBC_mask_v(i_obc+1,J) = .true. - if (OBC%OBC_segment_v(i_obc+1,J) == OBC_NONE) then - OBC%OBC_direction_v(i_obc+1,J) = OBC_DIRECTION_E - OBC%OBC_segment_v(i_obc+1,J) = l_seg - endif - OBC%OBC_mask_v(i_obc+1,J-1) = .true. - if (OBC%OBC_segment_v(i_obc+1,J-1) == OBC_NONE) then - OBC%OBC_direction_v(i_obc+1,J-1) = OBC_DIRECTION_E - OBC%OBC_segment_v(i_obc+1,J-1) = l_seg - endif - endif - else ! West is outward - if (this_kind == OBC_FLATHER) then - OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_W ! We only use direction for Flather - ! Set v points outside segment - OBC%OBC_mask_v(i_obc,J) = .true. - if (OBC%OBC_segment_v(i_obc,J) == OBC_NONE) then - OBC%OBC_direction_v(i_obc,J) = OBC_DIRECTION_W - OBC%OBC_segment_v(i_obc,J) = l_seg + if (Je_obc>Js_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E + if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA + if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA + if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA + if (Je_obcG%HI%IedB) return ! Boundary is not on tile + if (max(Js_obc,Je_obc)G%HI%JedB) return ! Segment is not on tile + + do j=G%HI%jsd, G%HI%jed + if (j>min(Js_obc,Je_obc) .and. j<=max(Js_obc,Je_obc)) then + OBC%OBC_mask_u(I_obc,j) = .true. + OBC%OBC_segment_u(I_obc,j) = l_seg + if (Je_obc>Js_obc) then ! East is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_E ! We only use direction for Flather (maybe) + ! Set v points outside segment + OBC%OBC_mask_v(i_obc+1,J) = .true. + if (OBC%OBC_segment_v(i_obc+1,J) == OBC_NONE) then + OBC%OBC_direction_v(i_obc+1,J) = OBC_DIRECTION_E + OBC%OBC_segment_v(i_obc+1,J) = l_seg + endif + OBC%OBC_mask_v(i_obc+1,J-1) = .true. + if (OBC%OBC_segment_v(i_obc+1,J-1) == OBC_NONE) then + OBC%OBC_direction_v(i_obc+1,J-1) = OBC_DIRECTION_E + OBC%OBC_segment_v(i_obc+1,J-1) = l_seg + endif endif - OBC%OBC_mask_v(i_obc,J-1) = .true. - if (OBC%OBC_segment_v(i_obc,J-1) == OBC_NONE) then - OBC%OBC_direction_v(i_obc,J-1) = OBC_DIRECTION_W - OBC%OBC_segment_v(i_obc,J-1) = l_seg + else ! West is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_W ! We only use direction for Flather + ! Set v points outside segment + OBC%OBC_mask_v(i_obc,J) = .true. + if (OBC%OBC_segment_v(i_obc,J) == OBC_NONE) then + OBC%OBC_direction_v(i_obc,J) = OBC_DIRECTION_W + OBC%OBC_segment_v(i_obc,J) = l_seg + endif + OBC%OBC_mask_v(i_obc,J-1) = .true. + if (OBC%OBC_segment_v(i_obc,J-1) == OBC_NONE) then + OBC%OBC_direction_v(i_obc,J-1) = OBC_DIRECTION_W + OBC%OBC_segment_v(i_obc,J-1) = l_seg + endif endif endif endif - endif - enddo + enddo + enddo ! a_loop end subroutine setup_u_point_obc @@ -348,90 +359,102 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) integer, intent(in) :: l_seg !< which segment is this? ! Local variables integer :: J_obc, Is_obc, Ie_obc ! Position of segment in global index space - integer :: i, this_kind - character(len=32) :: action_str + integer :: i, this_kind, a_loop + character(len=32) :: action_str(5) ! This returns the global indices for the segment call parse_segment_str(G%ieg, G%jeg, segment_str, J_obc, Is_obc, Ie_obc, action_str ) J_obc = J_obc - G%jdg_offset ! Convert to local tile indices on this tile Is_obc = Is_obc - G%idg_offset ! Convert to local tile indices on this tile Ie_obc = Ie_obc - G%idg_offset ! Convert to local tile indices on this tile + this_kind = OBC_NONE - if (trim(action_str) == 'FLATHER') then - this_kind = OBC_FLATHER - if (Ie_obc>Is_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_S - if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA - if (Ie_obcIs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_S - if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA - if (Ie_obcG%HI%JedB) return ! Boundary is not on tile - if (max(Is_obc,Ie_obc)G%HI%IedB) return ! Segment is not on tile - - ! These four lines extend the open boundary into the halo region of tiles on the edge of the physical - ! domain. They are used to reproduce the checksums of the circle_obcs test case and will be removed - ! in the fullness of time. -AJA -! These cause trouble with DOME -! if (Is_obc == G%HI%IscB) Is_obc = G%HI%isd-1 -! if (Is_obc == G%HI%IecB) Is_obc = G%HI%ied -! if (Ie_obc == G%HI%IscB) Ie_obc = G%HI%isd-1 -! if (Ie_obc == G%HI%IecB) Ie_obc = G%HI%ied - - do i=G%HI%isd, G%HI%ied - if (i>min(Is_obc,Ie_obc) .and. i<=max(Is_obc,Ie_obc)) then - OBC%OBC_mask_v(i,J_obc) = .true. - OBC%OBC_segment_v(i,J_obc) = l_seg - if (Is_obc>Ie_obc) then ! North is outward - if (this_kind == OBC_FLATHER) then - OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_N ! We only use direction for Flather - ! Set u points outside segment - OBC%OBC_mask_u(I,j_obc+1) = .true. - if (OBC%OBC_segment_u(I,j_obc+1) == OBC_NONE) then - OBC%OBC_direction_u(I,j_obc+1) = OBC_DIRECTION_N - OBC%OBC_segment_u(I,j_obc+1) = l_seg - endif - OBC%OBC_mask_u(I-1,j_obc+1) = .true. - if (OBC%OBC_segment_u(I-1,j_obc+1) == OBC_NONE) then - OBC%OBC_direction_u(I-1,j_obc+1) = OBC_DIRECTION_N - OBC%OBC_segment_u(I-1,j_obc+1) = l_seg - endif - endif - else ! South is outward - if (this_kind == OBC_FLATHER) then - OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_S ! We only use direction for Flather - ! Set u points outside segment - OBC%OBC_mask_u(I,j_obc) = .true. - if (OBC%OBC_segment_u(I,j_obc) == OBC_NONE) then - OBC%OBC_direction_u(I,j_obc) = OBC_DIRECTION_S - OBC%OBC_segment_u(I,j_obc) = l_seg + if (Ie_obc>Is_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_S + if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA + if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA + if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA + if (Ie_obcG%HI%JedB) return ! Boundary is not on tile + if (max(Is_obc,Ie_obc)G%HI%IedB) return ! Segment is not on tile + + do i=G%HI%isd, G%HI%ied + if (i>min(Is_obc,Ie_obc) .and. i<=max(Is_obc,Ie_obc)) then + OBC%OBC_mask_v(i,J_obc) = .true. + OBC%OBC_segment_v(i,J_obc) = l_seg + if (Is_obc>Ie_obc) then ! North is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_N ! We only use direction for Flather + ! Set u points outside segment + OBC%OBC_mask_u(I,j_obc+1) = .true. + if (OBC%OBC_segment_u(I,j_obc+1) == OBC_NONE) then + OBC%OBC_direction_u(I,j_obc+1) = OBC_DIRECTION_N + OBC%OBC_segment_u(I,j_obc+1) = l_seg + endif + OBC%OBC_mask_u(I-1,j_obc+1) = .true. + if (OBC%OBC_segment_u(I-1,j_obc+1) == OBC_NONE) then + OBC%OBC_direction_u(I-1,j_obc+1) = OBC_DIRECTION_N + OBC%OBC_segment_u(I-1,j_obc+1) = l_seg + endif endif - OBC%OBC_mask_u(I-1,j_obc) = .true. - if (OBC%OBC_segment_u(I-1,j_obc) == OBC_NONE) then - OBC%OBC_direction_u(I-1,j_obc) = OBC_DIRECTION_S - OBC%OBC_segment_u(I-1,j_obc) = l_seg + else ! South is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_S ! We only use direction for Flather + ! Set u points outside segment + OBC%OBC_mask_u(I,j_obc) = .true. + if (OBC%OBC_segment_u(I,j_obc) == OBC_NONE) then + OBC%OBC_direction_u(I,j_obc) = OBC_DIRECTION_S + OBC%OBC_segment_u(I,j_obc) = l_seg + endif + OBC%OBC_mask_u(I-1,j_obc) = .true. + if (OBC%OBC_segment_u(I-1,j_obc) == OBC_NONE) then + OBC%OBC_direction_u(I-1,j_obc) = OBC_DIRECTION_S + OBC%OBC_segment_u(I-1,j_obc) = l_seg + endif endif endif endif - endif - enddo + enddo + enddo ! a_loop end subroutine setup_v_point_obc @@ -443,11 +466,12 @@ subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_ integer, intent(out) :: l !< The value of I=l, if segment_str begins with I=l, or the value of J=l integer, intent(out) :: m !< The value of J=m, if segment_str begins with I=, or the value of I=m integer, intent(out) :: n !< The value of J=n, if segment_str begins with I=, or the value of I=n - character(len=*), intent(out) :: action_str !< The "string" part of segment_str + character(len=*), intent(out) :: action_str(:) !< The "string" part of segment_str ! Local variables character(len=24) :: word1, word2, m_word, n_word !< Words delineated by commas in a string in form of "I=%,J=%:%,string" integer :: l_max !< Either ni_global or nj_global, depending on whether segment_str begins with "I=" or "J=" integer :: mn_max !< Either nj_global or ni_global, depending on whether segment_str begins with "I=" or "J=" + integer :: j ! Process first word which will started with either 'I=' or 'J=' word1 = extract_word(segment_str,',',1) @@ -496,7 +520,9 @@ subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_ endif ! Type of open boundary condition - action_str = extract_word(segment_str,',',3) + do j = 1, size(action_str) + action_str(j) = extract_word(segment_str,',',2+j) + enddo contains @@ -513,14 +539,8 @@ integer function interpret_int_expr(string, imax) if (len_trim(string)==1 .and. string(1:1)=='N') then interpret_int_expr = imax elseif (string(1:1)=='N') then - read(string(3:slen),*,err=911) interpret_int_expr - if (string(2:2)=='-') then - interpret_int_expr = imax - interpret_int_expr - elseif (string(2:2)=='+') then - interpret_int_expr = imax + interpret_int_expr - else - goto 911 - endif + read(string(2:slen),*,err=911) interpret_int_expr + interpret_int_expr = imax - interpret_int_expr else read(string(1:slen),*,err=911) interpret_int_expr endif @@ -617,34 +637,30 @@ subroutine open_boundary_impose_normal_slope(OBC, G, depth) real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points ! Local variables integer :: i, j + logical :: bc_north, bc_south, bc_east, bc_west if (.not.associated(OBC)) return - if (associated(OBC%OBC_segment_u)) then - do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 - if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == & - OBC_DIRECTION_E) depth(i+1,j) = depth(i,j) - if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == & - OBC_DIRECTION_W) depth(i,j) = depth(i+1,j) - endif -! if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) depth(i+1,j) = depth(i,j) -! if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) depth(i,j) = depth(i+1,j) - enddo ; enddo - endif - - if (associated(OBC%OBC_segment_v)) then - do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied - if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == & - OBC_DIRECTION_N) depth(i,j+1) = depth(i,j) - if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == & - OBC_DIRECTION_S) depth(i,j) = depth(i,j+1) - endif -! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) depth(i,j+1) = depth(i,j) -! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) depth(i,j) = depth(i,j+1) - enddo ; enddo - endif + do J=G%jsd+1,G%jed-1 ; do i=G%isd+1,G%ied-1 + bc_north = .false. ; bc_south = .false. ; bc_east = .false. ; bc_west = .false. + if (associated(OBC%OBC_segment_u)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_E) bc_east = .true. + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I-1,j))%direction == OBC_DIRECTION_W) bc_west = .true. + endif + if (associated(OBC%OBC_segment_v)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_N) bc_north = .true. + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J-1))%direction == OBC_DIRECTION_S) bc_south = .true. + endif + if (bc_north) depth(i,j+1) = depth(i,j) + if (bc_south) depth(i,j-1) = depth(i,j) + if (bc_east) depth(i+1,j) = depth(i,j) + if (bc_west) depth(i-1,j) = depth(i,j) + ! Convex corner cases + if (bc_north.and.bc_east) depth(i+1,j+1) = depth(i,j) + if (bc_north.and.bc_west) depth(i-1,j+1) = depth(i,j) + if (bc_south.and.bc_east) depth(i+1,j-1) = depth(i,j) + if (bc_south.and.bc_west) depth(i-1,j-1) = depth(i,j) + enddo ; enddo end subroutine open_boundary_impose_normal_slope @@ -703,7 +719,7 @@ subroutine open_boundary_impose_land_mask(OBC, G) end subroutine open_boundary_impose_land_mask -!> Diagnose radiation conditions at open boundaries +!> Apply radiation conditions to 3D u,v (,h) at open boundaries subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & h_new, h_old, G) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure @@ -1205,6 +1221,36 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) h(i,j+1,k) = h(i,j,k) if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) h(i,j,k) = h(i,j+1,k) enddo ; enddo ; enddo +! When we do not extend segments, this commented block was needed to +! get the same'ish h's. +! do k=1,nz ; do j=jsd,jed-1 ; do i=isd,ied +! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) h(i,j+1,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd+1,jed ; do i=isd,ied +! if (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_S) h(i,j-1,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd,jed ; do i=isd,ied-1 +! if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) h(i+1,j,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd,jed ; do i=isd+1,ied +! if (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_W) h(i-1,j,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd,jed-1 ; do i=isd,ied-1 +! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N .and. & +! OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) h(i+1,j+1,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd,jed-1 ; do i=isd+1,ied +! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N .and. & +! OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_W) h(i-1,j+1,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd+1,jed ; do i=isd,ied-1 +! if (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_S .and. & +! OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) h(i+1,j-1,k) = h(i,j,k) +! enddo ; enddo ; enddo +! do k=1,nz ; do j=jsd+1,jed ; do i=isd+1,ied +! if (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_S .and. & +! OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_W) h(i-1,j-1,k) = h(i,j,k) +! enddo ; enddo ; enddo end subroutine set_Flather_data