diff --git a/src/tracer/MOM_generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 similarity index 100% rename from src/tracer/MOM_generic_tracer.F90 rename to config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7ee90d746a..e5d1f46086 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -986,7 +986,11 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS do j=js,je ; do i=is,ie CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j) enddo ; enddo - if (CS%IDs%id_ssh_inst > 0) call post_data(CS%IDs%id_ssh_inst, ssh, CS%diag) + if (CS%IDs%id_ssh_inst > 0) then + call enable_averages(dt, Time_local, CS%diag) + call post_data(CS%IDs%id_ssh_inst, ssh, CS%diag) + call disable_averaging(CS%diag) + endif call cpu_clock_end(id_clock_dynamics) endif @@ -4014,7 +4018,7 @@ subroutine extract_surface_state(CS, sfc_state_in) numberOfErrors=0 ! count number of errors do j=js,je ; do i=is,ie if (G%mask2dT(i,j)>0.) then - localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) - G%Z_ref & + localError = sfc_state%sea_lev(i,j) < -G%bathyT(i,j) - G%Z_ref & .or. sfc_state%sea_lev(i,j) >= CS%bad_val_ssh_max & .or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max & .or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick @@ -4041,7 +4045,7 @@ subroutine extract_surface_state(CS, sfc_state_in) write(msg(1:240),'(2(a,i4,1x),4(a,f8.3,1x),6(a,es11.4))') & 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & - 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), & + 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), & 'D=',US%Z_to_m*(G%bathyT(i,j)+G%Z_ref), 'SSH=',US%Z_to_m*sfc_state%sea_lev(i,j), & 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 21d484b9b6..70e0a738c4 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4959,15 +4959,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0) then call MOM_read_data(wave_drag_file, wave_drag_u, CS%lin_drag_u, G%Domain, & - position=EAST_FACE, scale=GV%m_to_H*US%T_to_s) - call pass_var(CS%lin_drag_u, G%Domain) - CS%lin_drag_u(:,:) = wave_drag_scale * CS%lin_drag_u(:,:) - + position=EAST_FACE, scale=wave_drag_scale*GV%m_to_H*US%T_to_s) call MOM_read_data(wave_drag_file, wave_drag_v, CS%lin_drag_v, G%Domain, & - position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s) - call pass_var(CS%lin_drag_v, G%Domain) - CS%lin_drag_v(:,:) = wave_drag_scale * CS%lin_drag_v(:,:) - + position=NORTH_FACE, scale=wave_drag_scale*GV%m_to_H*US%T_to_s) + call pass_vector(CS%lin_drag_u, CS%lin_drag_v, G%domain, direction=To_All+SCALAR_PAIR) else allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index abce27909b..de4e9190df 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -332,9 +332,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_masso > 0) then mass_cell(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - mass_cell(i,j) = mass_cell(i,j) + (GV%H_to_kg_m2*h(i,j,k)) * US%L_to_m**2*G%areaT(i,j) + mass_cell(i,j) = mass_cell(i,j) + (GV%H_to_RZ*h(i,j,k)) * G%areaT(i,j) enddo ; enddo ; enddo - masso = reproducing_sum(mass_cell) + masso = reproducing_sum(mass_cell, unscale=US%RZL2_to_kg) call post_data(CS%id_masso, masso, CS%diag) endif @@ -1644,8 +1644,9 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag Time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2', conversion=GV%H_to_kg_m2, & standard_name='sea_water_mass_per_unit_area', v_extensive=.true.) - CS%id_masso = register_scalar_field('ocean_model', 'masso', Time, & - diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass') + CS%id_masso = register_scalar_field('ocean_model', 'masso', Time, diag, & + 'Mass of liquid ocean', units='kg', conversion=US%RZL2_to_kg, & + standard_name='sea_water_mass') CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, & long_name='Cell Thickness', standard_name='cell_thickness', & diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index f5ff19630b..85029c5152 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -340,8 +340,8 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci real :: areaTm(SZI_(G),SZJ_(G)) ! A masked version of areaT [L2 ~> m2]. real :: KE(SZK_(GV)) ! The total kinetic energy of a layer [J]. real :: PE(SZK_(GV)+1)! The available potential energy of an interface [J]. - real :: KE_tot ! The total kinetic energy [J]. - real :: PE_tot ! The total available potential energy [J]. + real :: KE_tot ! The total kinetic energy [R Z L4 T-2 ~> J]. + real :: PE_tot ! The total available potential energy [R Z L4 T-2 ~> J]. real :: Z_0APE(SZK_(GV)+1) ! The uniform depth which overlies the same ! volume as is below an interface [Z ~> m]. real :: H_0APE(SZK_(GV)+1) ! A version of Z_0APE, converted to m, usually positive [m]. @@ -351,9 +351,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci ! the total mass of the ocean [m2 s-2]. real :: vol_lay(SZK_(GV)) ! The volume of fluid in a layer [Z L2 ~> m3]. real :: volbelow ! The volume of all layers beneath an interface [Z L2 ~> m3]. - real :: mass_lay(SZK_(GV)) ! The mass of fluid in a layer [kg]. + real :: mass_lay(SZK_(GV)) ! The mass of fluid in a layer [R Z L2 ~> kg]. + real :: mass_tot_RZL2 ! The total mass of the ocean [R Z L2 ~> kg]. real :: mass_tot ! The total mass of the ocean [kg]. - real :: vol_tot ! The total ocean volume [m3]. + real :: vol_tot ! The total ocean volume [Z L2 ~> m3]. real :: mass_chg ! The change in total ocean mass of fresh water since ! the last call to this subroutine [kg]. real :: mass_anom ! The change in fresh water that cannot be accounted for @@ -399,14 +400,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & tmp1 ! A temporary array used in reproducing sums [various] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & - PE_pt ! The potential energy at each point [J]. + PE_pt ! The potential energy at each point [R Z L4 T-2 ~> J]. real, dimension(SZI_(G),SZJ_(G)) :: & - Temp_int, Salt_int ! Layer and cell integrated heat and salt [J] and [g Salt]. - real :: HL2_to_kg ! A conversion factor from a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1] - real :: KE_scale_factor ! The combination of unit rescaling factors in the kinetic energy - ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or 1] - real :: PE_scale_factor ! The combination of unit rescaling factors in the potential energy - ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> 1] + Temp_int, Salt_int ! Layer and cell integrated heat and salt [Q R Z L2 ~> J] and [g Salt]. + real :: RZL4_to_J ! The combination of unit rescaling factors to convert the spatially integrated + ! kinetic or potential energies into mks units [T2 kg m2 R-1 Z-1 L-4 s-2 ~> 1] integer :: num_nc_fields ! The number of fields that will actually go into ! the NetCDF file. integer :: i, j, k, is, ie, js, je, nz, m, Isq, Ieq, Jsq, Jeq, isr, ier, jsr, jer @@ -532,7 +530,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) - HL2_to_kg = GV%H_to_kg_m2*US%L_to_m**2 + RZL4_to_J = US%RZL2_to_kg*US%L_T_to_m_s**2 ! Used to unscale energies if (.not.associated(CS)) call MOM_error(FATAL, & "write_energy: Module must be initialized before it is used.") @@ -544,28 +542,21 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci areaTm(i,j) = G%mask2dT(i,j)*G%areaT(i,j) enddo ; enddo - if (GV%Boussinesq) then - tmp1(:,:,:) = 0.0 - do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = h(i,j,k) * (HL2_to_kg*areaTm(i,j)) - enddo ; enddo ; enddo + tmp1(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do i=is,ie + tmp1(i,j,k) = h(i,j,k) * (GV%H_to_RZ*areaTm(i,j)) + enddo ; enddo ; enddo + mass_tot_RZL2 = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP, unscale=US%RZL2_to_kg) - mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) - do k=1,nz ; vol_lay(k) = (US%m_to_L**2*GV%H_to_Z/GV%H_to_kg_m2)*mass_lay(k) ; enddo + if (GV%Boussinesq) then + do k=1,nz ; vol_lay(k) = (1.0 / GV%Rho0) * mass_lay(k) ; enddo else - tmp1(:,:,:) = 0.0 - do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j) - enddo ; enddo ; enddo - mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) - if (CS%do_APE_calc) then call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = US%Z_to_m*US%L_to_m**2*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) + tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo - vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay) - do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2 * vol_lay(k) ; enddo + vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay, unscale=US%Z_to_m*US%L_to_m**2) endif endif ! Boussinesq @@ -692,7 +683,6 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci ! Calculate the Available Potential Energy integrated over each interface. With a nonlinear ! equation of state or with a bulk mixed layer this calculation is only approximate. ! With an ALE model this does not make sense and should be revisited. - PE_scale_factor = US%RZ_to_kg_m2*US%L_to_m**2*US%L_T_to_m_s**2 PE_pt(:,:,:) = 0.0 if (GV%Boussinesq) then do j=js,je ; do i=is,ie @@ -702,7 +692,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref)) hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref) hbot = (hbot + ABS(hbot)) * 0.5 - PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j)) * (GV%Rho0*GV%g_prime(K)) * & + PE_pt(i,j,K) = (0.5 * areaTm(i,j)) * (GV%Rho0*GV%g_prime(K)) * & (hint * hint - hbot * hbot) enddo enddo ; enddo @@ -711,7 +701,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci do K=nz,1,-1 hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) - PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * & + PE_pt(i,j,K) = (0.5 * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * & (hint * hint - hbot * hbot) enddo enddo ; enddo @@ -720,47 +710,44 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci do K=nz,2,-1 hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) - PE_pt(i,j,K) = (0.25 * PE_scale_factor * areaTm(i,j) * & + PE_pt(i,j,K) = (0.25 * areaTm(i,j) * & ((GV%Rlay(k)+GV%Rlay(k-1))*GV%g_prime(K))) * & (hint * hint - hbot * hbot) enddo hint = Z_0APE(1) + eta(i,j,1) ! eta and H_0 have opposite signs. hbot = max(Z_0APE(1) - (G%bathyT(i,j) + G%Z_ref), 0.0) - PE_pt(i,j,1) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * & + PE_pt(i,j,1) = (0.5 * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * & (hint * hint - hbot * hbot) enddo ; enddo endif - PE_tot = reproducing_sum(PE_pt, isr, ier, jsr, jer, sums=PE) + PE_tot = reproducing_sum(PE_pt, isr, ier, jsr, jer, sums=PE, unscale=RZL4_to_J) do k=1,nz+1 ; H_0APE(K) = US%Z_to_m*Z_0APE(K) ; enddo else PE_tot = 0.0 do k=1,nz+1 ; PE(K) = 0.0 ; H_0APE(K) = 0.0 ; enddo endif -! Calculate the Kinetic Energy integrated over each layer. - KE_scale_factor = HL2_to_kg*US%L_T_to_m_s**2 + ! Calculate the Kinetic Energy integrated over each layer. tmp1(:,:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = (0.25 * KE_scale_factor * (areaTm(i,j) * h(i,j,k))) * & + tmp1(i,j,k) = (0.25 * GV%H_to_RZ*(areaTm(i,j) * h(i,j,k))) * & (((u(I-1,j,k)**2) + (u(I,j,k)**2)) + ((v(i,J-1,k)**2) + (v(i,J,k)**2))) enddo ; enddo ; enddo - KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=KE) + KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=KE, unscale=RZL4_to_J) - toten = KE_tot + PE_tot - - Salt = 0.0 ; Heat = 0.0 + ! Use reproducing sums to do global integrals relate to the heat, salinity and water budgets. if (CS%use_temperature) then Temp_int(:,:) = 0.0 ; Salt_int(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - Salt_int(i,j) = Salt_int(i,j) + US%S_to_ppt*tv%S(i,j,k) * & - (h(i,j,k)*(HL2_to_kg * areaTm(i,j))) - Temp_int(i,j) = Temp_int(i,j) + (US%Q_to_J_kg*tv%C_p * tv%T(i,j,k)) * & - (h(i,j,k)*(HL2_to_kg * areaTm(i,j))) + Salt_int(i,j) = Salt_int(i,j) + tv%S(i,j,k) * (h(i,j,k)*(GV%H_to_RZ * areaTm(i,j))) + Temp_int(i,j) = Temp_int(i,j) + (tv%C_p * tv%T(i,j,k)) * (h(i,j,k)*(GV%H_to_RZ * areaTm(i,j))) enddo ; enddo ; enddo - salt_EFP = reproducing_sum_EFP(Salt_int, isr, ier, jsr, jer, only_on_PE=.true.) - heat_EFP = reproducing_sum_EFP(Temp_int, isr, ier, jsr, jer, only_on_PE=.true.) + salt_EFP = reproducing_sum_EFP(Salt_int, isr, ier, jsr, jer, only_on_PE=.true., & + unscale=US%RZL2_to_kg*US%S_to_ppt) + heat_EFP = reproducing_sum_EFP(Temp_int, isr, ier, jsr, jer, only_on_PE=.true., & + unscale=US%RZL2_to_kg*US%Q_to_J_kg) ! Combining the sums avoids multiple blocking all-PE updates. EFP_list(1) = salt_EFP ; EFP_list(2) = heat_EFP ; EFP_list(3) = CS%fresh_water_in_EFP @@ -770,13 +757,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci salt_EFP = EFP_list(1) ; heat_EFP = EFP_list(2) ; CS%fresh_water_in_EFP = EFP_list(3) CS%net_salt_in_EFP = EFP_list(4) ; CS%net_heat_in_EFP = EFP_list(5) - Salt = EFP_to_real(salt_EFP) - Heat = EFP_to_real(heat_EFP) else call EFP_sum_across_PEs(CS%fresh_water_in_EFP) endif -! Calculate the maximum CFL numbers. + ! Calculate the maximum CFL numbers. max_CFL(1:2) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq CFL_Iarea = G%IareaT(i,j) @@ -803,7 +788,12 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci call max_across_PEs(max_CFL, 2) + ! From this point onward, many of the calculations set or use variables in unscaled (mks) units. + + Salt = 0.0 ; Heat = 0.0 if (CS%use_temperature) then + Salt = EFP_to_real(salt_EFP) + Heat = EFP_to_real(heat_EFP) if (CS%previous_calls == 0) then CS%salt_prev_EFP = salt_EFP ; CS%heat_prev_EFP = heat_EFP endif @@ -826,6 +816,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci endif mass_chg = EFP_to_real(mass_chg_EFP) + mass_tot = US%RZL2_to_kg * mass_tot_RZL2 if (CS%use_temperature) then salin = Salt / mass_tot salin_anom = Salt_anom / mass_tot @@ -833,6 +824,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci temp = heat / (mass_tot*US%Q_to_J_kg*US%degC_to_C*tv%C_p) temp_anom = Heat_anom / (mass_tot*US%Q_to_J_kg*US%degC_to_C*tv%C_p) endif + toten = RZL4_to_J * (KE_tot + PE_tot) En_mass = toten / mass_tot call get_time(day, start_of_day, num_days) @@ -952,9 +944,12 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci call CS%fileenergy_nc%write_field(CS%fields(1), real(CS%ntrunc), reday) call CS%fileenergy_nc%write_field(CS%fields(2), toten, reday) + do k=1,nz+1 ; PE(K) = RZL4_to_J*PE(K) ; enddo call CS%fileenergy_nc%write_field(CS%fields(3), PE, reday) + do k=1,nz ; KE(k) = RZL4_to_J*KE(k) ; enddo call CS%fileenergy_nc%write_field(CS%fields(4), KE, reday) call CS%fileenergy_nc%write_field(CS%fields(5), H_0APE, reday) + do k=1,nz ; mass_lay(k) = US%RZL2_to_kg*mass_lay(k) ; enddo call CS%fileenergy_nc%write_field(CS%fields(6), mass_lay, reday) call CS%fileenergy_nc%write_field(CS%fields(7), mass_tot, reday) @@ -1018,19 +1013,17 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) !! to MOM_sum_output_init. ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & - FW_in, & ! The net fresh water input, integrated over a timestep [kg]. + FW_in, & ! The net fresh water input, integrated over a timestep [R Z L2 ~> kg]. salt_in, & ! The total salt added by surface fluxes, integrated ! over a time step [ppt kg]. heat_in ! The total heat added by surface fluxes, integrated - ! over a time step [J]. - real :: RZL2_to_kg ! A combination of scaling factors for mass [kg R-1 Z-1 L-2 ~> 1] - real :: QRZL2_to_J ! A combination of scaling factors for heat [J Q-1 R-1 Z-1 L-2 ~> 1] + ! over a time step [Q R Z L2 ~> J]. type(EFP_type) :: & FW_in_EFP, & ! The net fresh water input, integrated over a timestep ! and summed over space [kg]. salt_in_EFP, & ! The total salt added by surface fluxes, integrated - ! over a time step and summed over space [ppt kg]. + ! over a time step and summed over space [R Z L2 ~> ppt kg]. heat_in_EFP ! The total heat added by boundary fluxes, integrated ! over a time step and summed over space [J]. @@ -1038,14 +1031,11 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - RZL2_to_kg = US%L_to_m**2*US%RZ_to_kg_m2 - QRZL2_to_J = RZL2_to_kg*US%Q_to_J_kg - FW_in(:,:) = 0.0 if (associated(fluxes%evap)) then if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then do j=js,je ; do i=is,ie - FW_in(i,j) = RZL2_to_kg * dt*G%areaT(i,j)*(fluxes%evap(i,j) + & + FW_in(i,j) = dt*G%areaT(i,j)*(fluxes%evap(i,j) + & (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + & (fluxes%fprec(i,j) + fluxes%frunoff(i,j)))) enddo ; enddo @@ -1056,27 +1046,26 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) endif if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie - FW_in(i,j) = FW_in(i,j) + RZL2_to_kg*dt * & - G%areaT(i,j) * fluxes%seaice_melt(i,j) + FW_in(i,j) = FW_in(i,j) + dt * G%areaT(i,j) * fluxes%seaice_melt(i,j) enddo ; enddo ; endif salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0 if (CS%use_temperature) then if (associated(fluxes%sw)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt * QRZL2_to_J * G%areaT(i,j) * (fluxes%sw(i,j) + & + heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * (fluxes%sw(i,j) + & (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j)))) enddo ; enddo ; endif if (associated(fluxes%seaice_melt_heat)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt * QRZL2_to_J * G%areaT(i,j) * & + heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & fluxes%seaice_melt_heat(i,j) enddo ; enddo ; endif ! smg: new code ! include heat content from water transport across ocean surface ! if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie -! heat_in(i,j) = heat_in(i,j) + dt * QRZL2_to_J * G%areaT(i,j) * & +! heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & ! (fluxes%heat_content_lprec(i,j) + (fluxes%heat_content_fprec(i,j) & ! + (fluxes%heat_content_lrunoff(i,j) + (fluxes%heat_content_frunoff(i,j) & ! + (fluxes%heat_content_cond(i,j) + (fluxes%heat_content_vprec(i,j) & @@ -1086,41 +1075,41 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) ! smg: old code if (associated(fluxes%heat_content_evap)) then do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt * QRZL2_to_J * G%areaT(i,j) * & + heat_in(i,j) = heat_in(i,j) + dt * G%areaT(i,j) * & (fluxes%heat_content_evap(i,j) + fluxes%heat_content_lprec(i,j) + & fluxes%heat_content_cond(i,j) + fluxes%heat_content_fprec(i,j) + & fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j)) enddo ; enddo elseif (associated(tv%TempxPmE)) then do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + (tv%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%TempxPmE(i,j) + heat_in(i,j) = heat_in(i,j) + (tv%C_p * G%areaT(i,j)) * tv%TempxPmE(i,j) enddo ; enddo elseif (associated(fluxes%evap)) then do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + (US%Q_to_J_kg*tv%C_p * sfc_state%SST(i,j)) * FW_in(i,j) + heat_in(i,j) = heat_in(i,j) + (tv%C_p * sfc_state%SST(i,j)) * FW_in(i,j) enddo ; enddo endif ! The following heat sources may or may not be used. if (associated(tv%internal_heat)) then do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + (tv%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%internal_heat(i,j) + heat_in(i,j) = heat_in(i,j) + (tv%C_p * G%areaT(i,j)) * tv%internal_heat(i,j) enddo ; enddo endif if (associated(tv%frazil)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + QRZL2_to_J * G%areaT(i,j) * tv%frazil(i,j) + heat_in(i,j) = heat_in(i,j) + G%areaT(i,j) * tv%frazil(i,j) enddo ; enddo ; endif if (associated(fluxes%heat_added)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + QRZL2_to_J * dt*G%areaT(i,j) * fluxes%heat_added(i,j) + heat_in(i,j) = heat_in(i,j) + dt*G%areaT(i,j) * fluxes%heat_added(i,j) enddo ; enddo ; endif ! if (associated(sfc_state%sw_lost)) then ; do j=js,je ; do i=is,ie -! heat_in(i,j) = heat_in(i,j) - US%L_to_m**2*G%areaT(i,j) * sfc_state%sw_lost(i,j) +! sfc_state%sw_lost must be in units of [Q R Z ~> J m-2] +! heat_in(i,j) = heat_in(i,j) - G%areaT(i,j) * sfc_state%sw_lost(i,j) ! enddo ; enddo ; endif if (associated(fluxes%salt_flux)) then ; do j=js,je ; do i=is,ie ! integrate salt_flux in [R Z T-1 ~> kgSalt m-2 s-1] to give [ppt kg] - salt_in(i,j) = RZL2_to_kg * dt * & - G%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j)) + salt_in(i,j) = dt * G%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j)) enddo ; enddo ; endif endif @@ -1129,9 +1118,12 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) ! The on-PE sums are stored here, but the sums across PEs are deferred to ! the next call to write_energy to avoid extra barriers. isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) - FW_in_EFP = reproducing_sum_EFP(FW_in, isr, ier, jsr, jer, only_on_PE=.true.) - heat_in_EFP = reproducing_sum_EFP(heat_in, isr, ier, jsr, jer, only_on_PE=.true.) - salt_in_EFP = reproducing_sum_EFP(salt_in, isr, ier, jsr, jer, only_on_PE=.true.) + FW_in_EFP = reproducing_sum_EFP(FW_in, isr, ier, jsr, jer, only_on_PE=.true., & + unscale=US%RZL2_to_kg) + heat_in_EFP = reproducing_sum_EFP(heat_in, isr, ier, jsr, jer, only_on_PE=.true., & + unscale=US%RZL2_to_kg*US%Q_to_J_kg) + salt_in_EFP = reproducing_sum_EFP(salt_in, isr, ier, jsr, jer, only_on_PE=.true., & + unscale=US%RZL2_to_kg) CS%fresh_water_in_EFP = CS%fresh_water_in_EFP + FW_in_EFP CS%net_salt_in_EFP = CS%net_salt_in_EFP + salt_in_EFP diff --git a/src/framework/MOM_coms.F90 b/src/framework/MOM_coms.F90 index e4f5235da8..8471dc0b29 100644 --- a/src/framework/MOM_coms.F90 +++ b/src/framework/MOM_coms.F90 @@ -51,7 +51,8 @@ module MOM_coms logical :: NaN_error = .false. !< This becomes true if a NaN is encountered. logical :: debug = .false. !< Making this true enables debugging output. -!> Find an accurate and order-invariant sum of a distributed 2d or 3d field +!> Find an accurate and order-invariant sum of a distributed 2d or 3d field, in some cases after +!! undoing the scaling of the input array and restoring that scaling in the returned value interface reproducing_sum module procedure reproducing_sum_2d, reproducing_sum_3d end interface reproducing_sum @@ -91,8 +92,9 @@ module MOM_coms !! the result returned as an extended fixed point type that can be converted back to a real number !! using EFP_to_real. This technique is described in Hallberg & Adcroft, 2014, Parallel Computing, !! doi:10.1016/j.parco.2014.04.007. -function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, only_on_PE) result(EFP_sum) - real, dimension(:,:), intent(in) :: array !< The array to be summed in arbitrary units [a] +function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, only_on_PE, unscale) result(EFP_sum) + real, dimension(:,:), intent(in) :: array !< The array to be summed in arbitrary units [a], or in + !! arbitrary scaled units [A ~> a] if unscale is present integer, optional, intent(in) :: isr !< The starting i-index of the sum, noting !! that the array indices starts at 1 integer, optional, intent(in) :: ier !< The ending i-index of the sum, noting @@ -102,15 +104,17 @@ function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, integer, optional, intent(in) :: jer !< The ending j-index of the sum, noting !! that the array indices starts at 1 logical, optional, intent(in) :: overflow_check !< If present and false, disable - !! checking for overflows in incremental results. - !! This can speed up calculations if the number - !! of values being summed is small enough - integer, optional, intent(out) :: err !< If present, return an error code instead of - !! triggering any fatal errors directly from - !! this routine. + !! checking for overflows in incremental results. + !! This can speed up calculations if the number + !! of values being summed is small enough + integer, optional, intent(out) :: err !< If present, return an error code instead of + !! triggering any fatal errors directly from + !! this routine. logical, optional, intent(in) :: only_on_PE !< If present and true, do not do the sum - !! across processors, only reporting the local sum - type(EFP_type) :: EFP_sum !< The result in extended fixed point format + !! across processors, only reporting the local sum + real, optional, intent(in) :: unscale !< A factor that is used to undo scaling of array before it is + !! summed, often to compensate for the scaling in [a A-1 ~> 1] + type(EFP_type) :: EFP_sum !< The result in extended fixed point format ! This subroutine uses a conversion to an integer representation ! of real numbers to give order-invariant sums that will reproduce @@ -120,7 +124,8 @@ function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, integer(kind=int64) :: ival, prec_error real :: rs ! The remaining value to add, in arbitrary units [a] real :: max_mag_term ! A running maximum magnitude of the values in arbitrary units [a] - logical :: over_check, do_sum_across_PEs + real :: descale ! A local copy of unscale if it is present [a A-1 ~> 1] or 1 + logical :: over_check, do_sum_across_PEs, do_unscale character(len=256) :: mesg integer :: i, j, n, is, ie, js, je, sgn @@ -130,7 +135,7 @@ function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, prec_error = ((2_int64)**62 + ((2_int64)**62 - 1)) / num_PEs() - is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2 ) + is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2) if (present(isr)) then if (isr < is) call MOM_error(FATAL, "Value of isr too small in reproducing_EFP_sum_2d.") is = isr @@ -150,32 +155,40 @@ function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, over_check = .true. ; if (present(overflow_check)) over_check = overflow_check do_sum_across_PEs = .true. ; if (present(only_on_PE)) do_sum_across_PEs = .not.only_on_PE + do_unscale = .false. ; if (present(unscale)) do_unscale = (unscale /= 1.0) + descale = 1.0 ; if (do_unscale) descale = unscale overflow_error = .false. ; NaN_error = .false. ; max_mag_term = 0.0 ints_sum(:) = 0 if (over_check) then if ((je+1-js)*(ie+1-is) < max_count_prec) then - do j=js,je ; do i=is,ie - call increment_ints_faster(ints_sum, array(i,j), max_mag_term) - enddo ; enddo + ! This is the most common case, so handle the do_unscale case separately for efficiency. + if (do_unscale) then + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sum, unscale*array(i,j), max_mag_term) + enddo ; enddo + else + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sum, array(i,j), max_mag_term) + enddo ; enddo + endif call carry_overflow(ints_sum, prec_error) elseif ((ie+1-is) < max_count_prec) then do j=js,je do i=is,ie - call increment_ints_faster(ints_sum, array(i,j), max_mag_term) + call increment_ints_faster(ints_sum, descale*array(i,j), max_mag_term) enddo call carry_overflow(ints_sum, prec_error) enddo else do j=js,je ; do i=is,ie - call increment_ints(ints_sum, real_to_ints(array(i,j), prec_error), & - prec_error) + call increment_ints(ints_sum, real_to_ints(descale*array(i,j), prec_error), prec_error) enddo ; enddo endif else do j=js,je ; do i=is,ie sgn = 1 ; if (array(i,j)<0.0) sgn = -1 - rs = abs(array(i,j)) + rs = abs(descale*array(i,j)) do n=1,ni ival = int(rs*I_pr(n), kind=int64) rs = rs - ival*pr(n) @@ -213,13 +226,15 @@ function reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, end function reproducing_EFP_sum_2d + !> This subroutine uses a conversion to an integer representation of real numbers to give an !! order-invariant sum of distributed 2-D arrays that reproduces across domain decomposition. !! This technique is described in Hallberg & Adcroft, 2014, Parallel Computing, !! doi:10.1016/j.parco.2014.04.007. function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & - overflow_check, err, only_on_PE) result(sum) - real, dimension(:,:), intent(in) :: array !< The array to be summed in arbitrary units [a] + overflow_check, err, only_on_PE, unscale) result(sum) + real, dimension(:,:), intent(in) :: array !< The array to be summed in arbitrary units [a], or in + !! arbitrary scaled units [A ~> a] if unscale is present integer, optional, intent(in) :: isr !< The starting i-index of the sum, noting !! that the array indices starts at 1 integer, optional, intent(in) :: ier !< The ending i-index of the sum, noting @@ -228,7 +243,7 @@ function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & !! that the array indices starts at 1 integer, optional, intent(in) :: jer !< The ending j-index of the sum, noting !! that the array indices starts at 1 - type(EFP_type), optional, intent(out) :: EFP_sum !< The result in extended fixed point format + type(EFP_type), optional, intent(out) :: EFP_sum !< The result in extended fixed point format logical, optional, intent(in) :: reproducing !< If present and false, do the sum !! using the naive non-reproducing approach logical, optional, intent(in) :: overflow_check !< If present and false, disable @@ -240,16 +255,18 @@ function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & !! this routine. logical, optional, intent(in) :: only_on_PE !< If present and true, do not do the sum !! across processors, only reporting the local sum - real :: sum !< The sum of the values in array in arbitrary units [a] - - ! This subroutine uses a conversion to an integer representation - ! of real numbers to give order-invariant sums that will reproduce - ! across PE count. This idea comes from R. Hallberg and A. Adcroft. + real, optional, intent(in) :: unscale !< A factor that is used to undo scaling of array before it is + !! summed, often to compensate for the scaling in [a A-1 ~> 1] + real :: sum !< The sum of the values in array in the same + !! arbitrary units as array [a] or [A ~> a] + ! Local variables integer(kind=int64), dimension(ni) :: ints_sum integer(kind=int64) :: prec_error - real :: rsum(1) ! The running sum, in arbitrary units [a] - logical :: repro, do_sum_across_PEs + real :: rsum(1) ! The running sum, in arbitrary units [a] + real :: descale ! A local copy of unscale if it is present [a A-1 ~> 1] or 1 + real :: I_unscale ! The reciprocal of unscale [A a-1 ~> 1] + logical :: repro, do_sum_across_PEs, do_unscale character(len=256) :: mesg type(EFP_type) :: EFP_val ! An extended fixed point version of the sum integer :: i, j, is, ie, js, je @@ -260,7 +277,7 @@ function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & prec_error = ((2_int64)**62 + ((2_int64)**62 - 1)) / num_PEs() - is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2 ) + is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2) if (present(isr)) then if (isr < is) call MOM_error(FATAL, "Value of isr too small in reproducing_sum_2d.") is = isr @@ -280,19 +297,25 @@ function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & repro = .true. ; if (present(reproducing)) repro = reproducing do_sum_across_PEs = .true. ; if (present(only_on_PE)) do_sum_across_PEs = .not.only_on_PE + do_unscale = .false. ; if (present(unscale)) do_unscale = (unscale /= 1.0) + descale = 1.0 ; I_unscale = 1.0 + if (do_unscale) then + descale = unscale + if (abs(unscale) > 0.0) I_unscale = 1.0 / unscale + endif if (repro) then - EFP_val = reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, only_on_PE) - sum = ints_to_real(EFP_val%v) + EFP_val = reproducing_EFP_sum_2d(array, isr, ier, jsr, jer, overflow_check, err, only_on_PE, unscale) + sum = ints_to_real(EFP_val%v) * I_unscale if (present(EFP_sum)) EFP_sum = EFP_val if (debug) ints_sum(:) = EFP_sum%v(:) else rsum(1) = 0.0 do j=js,je ; do i=is,ie - rsum(1) = rsum(1) + array(i,j) + rsum(1) = rsum(1) + descale*array(i,j) enddo ; enddo if (do_sum_across_PEs) call sum_across_PEs(rsum,1) - sum = rsum(1) + sum = rsum(1) * I_unscale if (present(err)) then ; err = 0 ; endif @@ -312,7 +335,7 @@ function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, & endif if (debug) then - write(mesg,'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni) + write(mesg,'("2d RS: ", ES24.16, 6 Z17.16)') sum*descale, ints_sum(1:ni) call MOM_mesg(mesg, 3) endif @@ -322,9 +345,10 @@ end function reproducing_sum_2d !! order-invariant sum of distributed 3-D arrays that reproduces across domain decomposition. !! This technique is described in Hallberg & Adcroft, 2014, Parallel Computing, !! doi:10.1016/j.parco.2014.04.007. -function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_sums, err, only_on_PE) & +function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_sums, err, only_on_PE, unscale) & result(sum) - real, dimension(:,:,:), intent(in) :: array !< The array to be summed in arbitrary units [a] + real, dimension(:,:,:), intent(in) :: array !< The array to be summed in arbitrary units [a], or in + !! arbitrary scaled units [A ~> a] if unscale is present integer, optional, intent(in) :: isr !< The starting i-index of the sum, noting !! that the array indices starts at 1 integer, optional, intent(in) :: ier !< The ending i-index of the sum, noting @@ -333,28 +357,31 @@ function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_su !! that the array indices starts at 1 integer, optional, intent(in) :: jer !< The ending j-index of the sum, noting !! that the array indices starts at 1 - real, dimension(:), optional, intent(out) :: sums !< The sums by vertical layer in abitrary units [a] + real, dimension(:), optional, intent(out) :: sums !< The sums by vertical layer in the same + !! abitrary units as array [a] or [A ~> a] type(EFP_type), optional, intent(out) :: EFP_sum !< The result in extended fixed point format type(EFP_type), dimension(:), & optional, intent(out) :: EFP_lay_sums !< The sums by vertical layer in EFP format - integer, optional, intent(out) :: err !< If present, return an error code instead of - !! triggering any fatal errors directly from - !! this routine. + integer, optional, intent(out) :: err !< If present, return an error code instead of + !! triggering any fatal errors directly from + !! this routine. logical, optional, intent(in) :: only_on_PE !< If present and true, do not do the sum - !! across processors, only reporting the local sum - real :: sum !< The sum of the values in array in arbitrary units [a] - - ! This subroutine uses a conversion to an integer representation - ! of real numbers to give order-invariant sums that will reproduce - ! across PE count. This idea comes from R. Hallberg and A. Adcroft. + !! across processors, only reporting the local sum + real, optional, intent(in) :: unscale !< A factor that is used to undo scaling of array before it is + !! summed, often to compensate for the scaling in [a A-1 ~> 1] + real :: sum !< The sum of the values in array in the same + !! arbitrary units as array [a] or [A ~> a] + ! Local variables real :: val ! The real number that is extracted in arbitrary units [a] real :: max_mag_term ! A running maximum magnitude of the val's in arbitrary units [a] + real :: descale ! A local copy of unscale if it is present [a A-1 ~> 1] or 1 + real :: I_unscale ! The Adcroft reciprocal of unscale [A a-1 ~> 1] integer(kind=int64), dimension(ni) :: ints_sum integer(kind=int64), dimension(ni,size(array,3)) :: ints_sums integer(kind=int64) :: prec_error character(len=256) :: mesg - logical :: do_sum_across_PEs + logical :: do_sum_across_PEs, do_unscale integer :: i, j, k, is, ie, js, je, ke, isz, jsz, n if (num_PEs() > max_count_prec) call MOM_error(FATAL, & @@ -384,6 +411,7 @@ function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_su jsz = je+1-js; isz = ie+1-is do_sum_across_PEs = .true. ; if (present(only_on_PE)) do_sum_across_PEs = .not.only_on_PE + do_unscale = .false. ; if (present(unscale)) do_unscale = (unscale /= 1.0) if (present(sums) .or. present(EFP_lay_sums)) then if (present(sums)) then ; if (size(sums) < ke) then @@ -396,22 +424,28 @@ function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_su overflow_error = .false. ; NaN_error = .false. ; max_mag_term = 0.0 if (jsz*isz < max_count_prec) then do k=1,ke - do j=js,je ; do i=is,ie - call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term) - enddo ; enddo + if (do_unscale) then + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sums(:,k), unscale*array(i,j,k), max_mag_term) + enddo ; enddo + else + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term) + enddo ; enddo + endif call carry_overflow(ints_sums(:,k), prec_error) enddo elseif (isz < max_count_prec) then do k=1,ke ; do j=js,je do i=is,ie - call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term) + call increment_ints_faster(ints_sums(:,k), descale*array(i,j,k), max_mag_term) enddo call carry_overflow(ints_sums(:,k), prec_error) enddo ; enddo else do k=1,ke ; do j=js,je ; do i=is,ie call increment_ints(ints_sums(:,k), & - real_to_ints(array(i,j,k), prec_error), prec_error) + real_to_ints(descale*array(i,j,k), prec_error), prec_error) enddo ; enddo ; enddo endif if (present(err)) then @@ -458,21 +492,27 @@ function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_su overflow_error = .false. ; NaN_error = .false. ; max_mag_term = 0.0 if (jsz*isz < max_count_prec) then do k=1,ke - do j=js,je ; do i=is,ie - call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term) - enddo ; enddo + if (do_unscale) then + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sum, unscale*array(i,j,k), max_mag_term) + enddo ; enddo + else + do j=js,je ; do i=is,ie + call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term) + enddo ; enddo + endif call carry_overflow(ints_sum, prec_error) enddo elseif (isz < max_count_prec) then do k=1,ke ; do j=js,je do i=is,ie - call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term) + call increment_ints_faster(ints_sum, descale*array(i,j,k), max_mag_term) enddo call carry_overflow(ints_sum, prec_error) enddo ; enddo else do k=1,ke ; do j=js,je ; do i=is,ie - call increment_ints(ints_sum, real_to_ints(array(i,j,k), prec_error), & + call increment_ints(ints_sum, real_to_ints(descale*array(i,j,k), prec_error), & prec_error) enddo ; enddo ; enddo endif @@ -504,6 +544,15 @@ function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, EFP_lay_su endif endif + if (do_unscale) then + ! Revise the sum to restore the scaling of input array before it is returned + I_unscale = 0.0 ; if (abs(unscale) > 0.0) I_unscale = 1.0 / unscale + sum = sum * I_unscale + if (present(sums)) then + do k=1,ke ; sums(k) = sums(k) * I_unscale ; enddo + endif + endif + end function reproducing_sum_3d !> Convert a real number into the array of integers constitute its extended-fixed-point representation @@ -516,9 +565,11 @@ function real_to_ints(r, prec_error, overflow) result(ints) logical, optional, intent(inout) :: overflow !< Returns true if the conversion is being !! done on a value that is too large to be represented integer(kind=int64), dimension(ni) :: ints + ! This subroutine converts a real number to an equivalent representation ! using several long integers. + ! Local variables real :: rs ! The remaining value to add, in arbitrary units [a] character(len=80) :: mesg integer(kind=int64) :: ival, prec_err diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90 index 868352102e..8b4f9266a8 100644 --- a/src/framework/MOM_unit_scaling.F90 +++ b/src/framework/MOM_unit_scaling.F90 @@ -47,6 +47,7 @@ module MOM_unit_scaling real :: QRZ_T_to_W_m2 !< Convert heat fluxes from Q R Z T-1 to W m-2 [W T Q-1 R-1 Z-1 m-2 ~> 1] ! Not used enough: real :: kg_m2_to_RZ !< Convert mass loads from kg m-2 to R Z [R Z m2 kg-1 ~> 1] real :: RZ_to_kg_m2 !< Convert mass loads from R Z to kg m-2 [kg R-1 Z-1 m-2 ~> 1] + real :: RZL2_to_kg !< Convert masses from R Z L2 to kg [kg R-1 Z-1 L-2 ~> 1] real :: kg_m2s_to_RZ_T !< Convert mass fluxes from kg m-2 s-1 to R Z T-1 [R Z m2 s T-1 kg-1 ~> 1] real :: RZ_T_to_kg_m2s !< Convert mass fluxes from R Z T-1 to kg m-2 s-1 [T kg R-1 Z-1 m-2 s-1 ~> 1] real :: RZ3_T3_to_W_m2 !< Convert turbulent kinetic energy fluxes from R Z3 T-3 to W m-2 [W T3 R-1 Z-3 m-2 ~> 1] @@ -224,6 +225,8 @@ subroutine set_unit_scaling_combos(US) ! Wind stresses: US%RLZ_T2_to_Pa = US%R_to_kg_m3 * US%L_T_to_m_s**2 * US%Z_to_L US%Pa_to_RLZ_T2 = US%kg_m3_to_R * US%m_s_to_L_T**2 * US%L_to_Z + ! Masses: + US%RZL2_to_kg = US%R_to_kg_m3 * US%Z_to_m * US%L_to_m**2 end subroutine set_unit_scaling_combos diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index eabd376512..0a257b6ca1 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1314,17 +1314,14 @@ subroutine compute_global_grid_integrals(G, US) ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming ! Masked and unscaled cell areas [m2] - real :: area_scale ! A scaling factor for area into MKS units [m2 L-2 ~> 1] - integer :: i,j - - area_scale = US%L_to_m**2 + integer :: i, j tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) + tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j) enddo ; enddo - G%areaT_global = reproducing_sum(tmpForSumming) + G%areaT_global = US%L_to_m**2 * reproducing_sum(tmpForSumming, unscale=US%L_to_m**2) if (G%areaT_global == 0.0) & call MOM_error(FATAL, "compute_global_grid_integrals: zero ocean area (check topography?)") diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 5c6bd75da5..b5a36ba9d2 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1296,7 +1296,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl bfsfc=surfBuoyFlux, & ! surface buoyancy flux [m2 s-3] uStar=surfFricVel, & ! surface friction velocity [m s-1] CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - CS%Vt2(i,j,:) = US%m_to_Z*US%T_to_s * Vt2_1d(:) + CS%Vt2(i,j,:) = US%m_to_Z**2*US%T_to_s**2 * Vt2_1d(:) endif ! recompute wscale for diagnostics, now that we in fact know boundary layer depth diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e984d5831c..819b06ee50 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -858,7 +858,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, & CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US) @@ -1410,7 +1410,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, & CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f10e2f445d..e32bc8de2d 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -5,6 +5,7 @@ module MOM_energetic_PBL use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_coms, only : EFP_type, real_to_EFP, EFP_to_real, operator(+), assignment(=), EFP_sum_across_PEs +use MOM_debugging, only : hchksum use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -16,7 +17,7 @@ module MOM_energetic_PBL use MOM_intrinsic_functions, only : cuberoot use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number use MOM_stochastics, only : stochastic_CS @@ -59,6 +60,8 @@ module MOM_energetic_PBL !! self-consistent mixed layer depth. Otherwise use the false position !! after a maximum and minimum bound have been evaluated and the !! returned value from the previous guess or bisection before this. + logical :: MLD_iter_bug !< If true use buggy logic that gives the wrong bounds for the next + !! iteration when successive guesses increase by exactly EPBL_MLD_TOLERANCE. integer :: max_MLD_its !< The maximum number of iterations that can be used to find a !! self-consistent mixed layer depth with Use_MLD_iteration. real :: MixLenExponent !< Exponent in the mixing length shape-function [nondim]. @@ -67,6 +70,10 @@ module MOM_energetic_PBL real :: MKE_to_TKE_effic !< The efficiency with which mean kinetic energy released by !! mechanically forced entrainment of the mixed layer is converted to !! TKE [nondim]. + logical :: direct_calc !< If true and there is no conversion from mean kinetic energy to ePBL + !! turbulent kinetic energy, use a direct calculation of the + !! diffusivity that is supported by a given energy input instead of the + !! more general but slower iterative solver. real :: ustar_min !< A minimum value of ustar to avoid numerical problems [Z T-1 ~> m s-1]. !! If the value is small enough, this should not affect the solution. real :: Ekman_scale_coef !< A nondimensional scaling factor controlling the inhibition of the @@ -92,9 +99,8 @@ module MOM_energetic_PBL !! Making this larger increases the diffusivity. real :: vstar_surf_fac !< If (wT_scheme == wT_from_RH18) this is the proportionality coefficient between !! ustar and the surface mechanical contribution to vstar [nondim] - real :: vstar_scale_fac !< An overall nondimensional scaling factor for vstar times a unit - !! conversion factor [Z s T-1 m-1 ~> nondim]. Making this larger increases - !! the diffusivity. + real :: vstar_scale_fac !< An overall nondimensional scaling factor for vstar [nondim]. Making + !! this larger increases the diffusivity. !mstar related options integer :: mstar_scheme !< An encoded integer to determine which formula is used to set mstar @@ -155,6 +161,43 @@ module MOM_energetic_PBL !! the Ekman depth over the Obukhov depth with destabilizing forcing [nondim]. real :: Max_Enhance_M = 5. !< The maximum allowed LT enhancement to the mixing [nondim]. + !/ Bottom boundary layer mixing related options + real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL [nondim] + logical :: Use_BBLD_iteration !< If true, use the proximity to the top of the actively turbulent + !! bottom boundary layer to constrain the mixing lengths. + real :: TKE_decay_BBL !< The ratio of the natural Ekman depth to the TKE decay scale for + !! bottom boundary layer mixing [nondim] + real :: min_BBL_mix_len !< The minimum mixing length scale that will be used by ePBL in the bottom + !! boundary layer mixing [Z ~> m]. The default (0) does not set a minimum. + real :: MixLenExponent_BBL !< Exponent in the bottom boundary layer mixing length shape-function [nondim]. + !! 1 is law-of-the-wall at top and bottom, + !! 2 is more KPP like. + real :: BBLD_tol !< The tolerance for the iteratively determined bottom boundary layer depth [Z ~> m]. + !! This is only used with USE_MLD_ITERATION. + integer :: max_BBLD_its !< The maximum number of iterations that can be used to find a self-consistent + !! bottom boundary layer depth. + integer :: wT_scheme_BBL !< An enumerated value indicating the method for finding the bottom boundary + !! layer turbulent velocity scale. There are currently two options: + !! wT_mwT_from_cRoot_TKE is the original (TKE_remaining)^1/3 + !! wT_from_RH18 is the version described by Reichl and Hallberg, 2018 + real :: vstar_scale_fac_BBL !< An overall nondimensional scaling factor for wT in the bottom boundary layer [nondim]. + !! Making this larger increases the bottom boundary layer diffusivity.", & + real :: vstar_surf_fac_BBL !< If (wT_scheme_BBL == wT_from_RH18) this is the proportionality coefficient between + !! ustar and the bottom boundayer layer mechanical contribution to vstar [nondim] + real :: Ekman_scale_coef_BBL !< A nondimensional scaling factor controlling the inhibition of the + !! diffusive length scale by rotation in the bottom boundary layer [nondim]. + !! Making this larger decreases the bottom boundary layer diffusivity. + logical :: decay_adjusted_BBL_TKE !< If true, include an adjustment factor in the bottom boundary layer + !! energetics that accounts for an exponential decay of TKE from a + !! near-bottom source and an assumed piecewise linear linear profile + !! of the buoyancy flux response to a change in a diffusivity. + + !/ Options for documenting differences from parameter choices + integer :: options_diff !< If positive, this is a coded integer indicating a pair of + !! settings whose differences are diagnosed in a passive diagnostic mode + !! via extra calls to ePBL_column. If this is 0 or negative no extra + !! calls occur. + !/ Others type(time_type), pointer :: Time=>NULL() !< A pointer to the ocean model's clock. @@ -170,38 +213,30 @@ module MOM_energetic_PBL !! potential energy change code. Otherwise, it uses a newer version !! that can work with successive increments to the diffusivity in !! upward or downward passes. + logical :: debug !< If true, write verbose checksums for debugging purposes. type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the !! timing of diagnostic output. real, allocatable, dimension(:,:) :: & - ML_depth !< The mixed layer depth determined by active mixing in ePBL [H ~> m or kg m-2] - ! These are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3]. + ML_depth !< The mixed layer depth determined by active mixing in ePBL, which may + !! be used for the first guess in the next time step [H ~> m or kg m-2] real, allocatable, dimension(:,:) :: & - diag_TKE_wind, & !< The wind source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_MKE, & !< The resolved KE source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_conv, & !< The convective source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_forcing, & !< The TKE sink required to mix surface penetrating shortwave heating - !! [R Z3 T-3 ~> W m-2]. - diag_TKE_mech_decay, & !< The decay of mechanical TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_conv_decay, & !< The decay of convective TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_mixing, & !< The work done by TKE to deepen the mixed layer [R Z3 T-3 ~> W m-2]. - ! These additional diagnostics are also 2d. - MSTAR_MIX, & !< Mstar used in EPBL [nondim] - MSTAR_LT, & !< Mstar due to Langmuir turbulence [nondim] - LA, & !< Langmuir number [nondim] - LA_MOD !< Modified Langmuir number [nondim] + BBL_depth !< The bottom boundary layer depth determined by active mixing in ePBL [H ~> m or kg m-2] type(EFP_type), dimension(2) :: sum_its !< The total number of iterations and columns worked on + type(EFP_type), dimension(2) :: sum_its_BBL !< The total number of iterations and columns worked on - real, allocatable, dimension(:,:,:) :: & - Velocity_Scale, & !< The velocity scale used in getting Kd [Z T-1 ~> m s-1] - Mixing_Length !< The length scale used in getting Kd [Z ~> m] !>@{ Diagnostic IDs integer :: id_ML_depth = -1, id_hML_depth = -1, id_TKE_wind = -1, id_TKE_mixing = -1 integer :: id_TKE_MKE = -1, id_TKE_conv = -1, id_TKE_forcing = -1 integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1 integer :: id_Mixing_Length = -1, id_Velocity_Scale = -1 + integer :: id_Kd_BBL = -1, id_BBL_Mix_Length = -1, id_BBL_Vel_Scale = -1 + integer :: id_TKE_BBL = -1, id_TKE_BBL_mixing = -1, id_TKE_BBL_decay = -1 + integer :: id_ustar_BBL = -1, id_BBL_decay_scale = -1, id_BBL_depth = -1 integer :: id_MSTAR_mix = -1, id_LA_mod = -1, id_LA = -1, id_MSTAR_LT = -1 + ! The next options are used when passively diagnosing sensitivities from parameter choices + integer :: id_opt_diff_Kd_ePBL = -1, id_opt_maxdiff_Kd_ePBL = -1, id_opt_diff_hML_depth = -1 !>@} end type energetic_PBL_CS @@ -236,11 +271,14 @@ module MOM_energetic_PBL !>@{ Local column copies of energy change diagnostics, all in [R Z3 T-3 ~> W m-2]. real :: dTKE_conv, dTKE_forcing, dTKE_wind, dTKE_mixing ! Local column diagnostics [R Z3 T-3 ~> W m-2] real :: dTKE_MKE, dTKE_mech_decay, dTKE_conv_decay ! Local column diagnostics [R Z3 T-3 ~> W m-2] + real :: dTKE_BBL, dTKE_BBL_decay, dTKE_BBL_mixing ! Local column diagnostics [R Z3 T-3 ~> W m-2] !>@} real :: LA !< The value of the Langmuir number [nondim] real :: LAmod !< The modified Langmuir number by convection [nondim] real :: mstar !< The value of mstar used in ePBL [nondim] real :: mstar_LT !< The portion of mstar due to Langmuir turbulence [nondim] + integer :: OBL_its !< The number of iterations used to find a self-consistent surface boundary layer depth + integer :: BBL_its !< The number of iterations used to find a self-consistent bottom boundary layer depth end type ePBL_column_diags contains @@ -249,7 +287,7 @@ module MOM_energetic_PBL !! mixed layer model. It assumes that heating, cooling and freshwater fluxes !! have already been applied. All calculations are done implicitly, and there !! is no stability limit on the time step. -subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, US, CS, & stoch_CS, dSV_dT, dSV_dS, TKE_forced, buoy_flux, Waves ) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -279,6 +317,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields have !! NULL ptrs. + type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields real, intent(in) :: dt !< Time increment [T ~> s]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(out) :: Kd_int !< The diagnosed diffusivities at interfaces @@ -287,7 +327,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]. type(wave_parameters_CS), pointer :: Waves !< Waves control structure for Langmuir turbulence - type(stochastic_CS), pointer :: stoch_CS !< The control structure returned by a previous + type(stochastic_CS), pointer :: stoch_CS !< The control structure returned by a previous ! This subroutine determines the diffusivities from the integrated energetics ! mixed layer model. It assumes that heating, cooling and freshwater fluxes @@ -335,12 +375,17 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS u, & ! The zonal velocity [L T-1 ~> m s-1]. v ! The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZK_(GV)+1) :: & - Kd, & ! The diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + Kd, & ! The diapycnal diffusivity due to ePBL [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. mixvel, & ! A turbulent mixing velocity [Z T-1 ~> m s-1]. mixlen, & ! A turbulent mixing length [Z ~> m]. - SpV_dt ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) + mixvel_BBL, & ! A bottom boundary layer turbulent mixing velocity [Z T-1 ~> m s-1]. + mixlen_BBL, & ! A bottom boundary layer turbulent mixing length [Z ~> m]. + Kd_BBL, & ! The bottom boundary layer diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + SpV_dt, & ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0), + ! in [R-1 T-1 ~> m3 kg-1 s-1], used to convert local TKE into a turbulence velocity cubed. + SpV_dt_cf ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) ! times conversion factors for answer dates before 20240101 in - ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the convsersion factors for + ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the conversion factors for ! answer dates of 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], used to ! convert local TKE into a turbulence velocity cubed. real :: h_neglect ! A thickness that is so small it is usually lost @@ -351,6 +396,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS real :: U_Star_Mean ! The surface friction without gustiness [Z T-1 ~> m s-1]. real :: mech_TKE ! The mechanically generated turbulent kinetic energy available for mixing over a ! timestep before the application of the efficiency in mstar [R Z3 T-2 ~> J m-2] + real :: u_star_BBL ! The bottom boundary layer friction velocity [H T-1 ~> m s-1 or kg m-2 s-1]. + real :: BBL_TKE ! The mechanically generated turbulent kinetic energy available for bottom + ! boundary layer mixing within a timestep [R Z3 T-2 ~> J m-2] real :: I_rho ! The inverse of the Boussinesq reference density times a ratio of scaling ! factors [Z L-1 R-1 ~> m3 kg-1] real :: I_dt ! The Adcroft reciprocal of the timestep [T-1 ~> s-1] @@ -358,9 +406,63 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! step [R-1 T-1 ~> m3 kg-1 s-1] real :: B_Flux ! The surface buoyancy flux [Z2 T-3 ~> m2 s-3] real :: MLD_io ! The mixed layer depth found by ePBL_column [Z ~> m] + real :: BBLD_io ! The bottom boundary layer thickness found by ePBL_BBL_column [Z ~> m] + real :: MLD_in ! The first guess at the mixed layer depth [Z ~> m] + real :: BBLD_in ! The first guess at the bottom boundary layer thickness [Z ~> m] type(ePBL_column_diags) :: eCD ! A container for passing around diagnostics. + ! The following variables are used for diagnostics + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + diag_Velocity_Scale, & ! The velocity scale used in getting Kd [Z T-1 ~> m s-1] + diag_Mixing_Length, & ! The length scale used in getting Kd [Z ~> m] + Kd_BBL_3d, & ! The bottom boundary layer diffusivities [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + BBL_Vel_Scale, & ! The velocity scale used in getting the BBL part of Kd [Z T-1 ~> m s-1] + BBL_Mix_Length ! The length scale used in getting the BBL part of Kd [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)) :: & + ! The next 7 diagnostics are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3]. + diag_TKE_wind, & ! The wind source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_MKE, & ! The resolved KE source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_conv, & ! The convective source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_forcing, & ! The TKE sink required to mix surface penetrating shortwave heating [R Z3 T-3 ~> W m-2] + diag_TKE_mech_decay, & ! The decay of mechanical TKE [R Z3 T-3 ~> W m-2] + diag_TKE_conv_decay, & ! The decay of convective TKE [R Z3 T-3 ~> W m-2] + diag_TKE_mixing, & ! The work done by TKE to deepen the mixed layer [R Z3 T-3 ~> W m-2] + diag_TKE_BBL, & ! The source of TKE to the bottom boundary layer [R Z3 T-3 ~> W m-2]. + diag_TKE_BBL_mixing, & ! The work done by TKE to thicken the bottom boundary layer [R Z3 T-3 ~> W m-2]. + diag_TKE_BBL_decay, & ! The work lost to decy of mechanical TKE in the bottom boundary + ! layer [R Z3 T-3 ~> W m-2]. + diag_ustar_BBL, & ! The bottom boundary layer friction velocity [H T-1 ~> m s-1 or kg m-2 s-1] + diag_BBL_decay_scale, & ! The bottom boundary layer TKE decay length scale [H ~> m] + + diag_mStar_MIX, & ! Mstar used in EPBL [nondim] + diag_mStar_LT, & ! Mstar due to Langmuir turbulence [nondim] + diag_LA, & ! Langmuir number [nondim] + diag_LA_MOD ! Modified Langmuir number [nondim] + + ! The following variables are only used for diagnosing sensitivities to ePBL settings + real, dimension(SZK_(GV)+1) :: & + Kd_1, Kd_2 ! Diapycnal diffusivities found with different ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: diff_Kd(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The change in diapycnal diffusivities found with different + ! ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: max_abs_diff_Kd(SZI_(G),SZJ_(G)) ! The column maximum magnitude of the change in diapycnal + ! diffusivities found with different ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: diff_hML_depth(SZI_(G),SZJ_(G)) ! The change in diagnosed active mixing layer depth with + ! different ePBL options [Z ~> m] + real :: BLD_1, BLD_2 ! Surface or bottom boundary layer depths found with different ePBL_column options [Z ~> m] + real :: SpV_scale1 ! A factor that accounts for the varying scaling of SpV_dt with answer date + ! [nondim] or [T3 m3 Z-3 s-3 ~> 1] + real :: SpV_scale2 ! A factor that accounts for the varying scaling of SpV_dt with answer date + ! [nondim] or [Z3 s3 T-3 m-3 ~> 1] + real :: SpV_dt_tmp(SZK_(GV)+1) ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) + ! times conversion factors for answer dates before 20240101 in + ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the conversion factors for + ! answer dates of 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], used to + ! convert local TKE into a turbulence velocity cubed. + type(ePBL_column_diags) :: eCD_tmp ! A container for not passing around diagnostics. + type(energetic_PBL_CS) :: CS_tmp1, CS_tmp2 ! Copies of the energetic PBL control structure that + ! can be modified to test for sensitivities + integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -386,16 +488,60 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Zero out diagnostics before accumulation. if (CS%TKE_diagnostics) then - !!OMP parallel do default(none) shared(is,ie,js,je,CS) + !!OMP parallel do default(shared) do j=js,je ; do i=is,ie - CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_MKE(i,j) = 0.0 - CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_forcing(i,j) = 0.0 - CS%diag_TKE_mixing(i,j) = 0.0 ; CS%diag_TKE_mech_decay(i,j) = 0.0 - CS%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced(i,j) = 0.0 + diag_TKE_wind(i,j) = 0.0 ; diag_TKE_MKE(i,j) = 0.0 + diag_TKE_conv(i,j) = 0.0 ; diag_TKE_forcing(i,j) = 0.0 + diag_TKE_mixing(i,j) = 0.0 ; diag_TKE_mech_decay(i,j) = 0.0 + diag_TKE_conv_decay(i,j) = 0.0 !; diag_TKE_unbalanced(i,j) = 0.0 enddo ; enddo + if (CS%ePBL_BBL_effic > 0.0) then + !!OMP parallel do default(shared) + do j=js,je ; do i=is,ie + diag_TKE_BBL(i,j) = 0.0 ; diag_TKE_BBL_mixing(i,j) = 0.0 + diag_TKE_BBL_decay(i,j) = 0.0 + enddo ; enddo + endif + endif + if (CS%debug .or. (CS%id_Mixing_Length>0)) diag_Mixing_Length(:,:,:) = 0.0 + if (CS%debug .or. (CS%id_Velocity_Scale>0)) diag_Velocity_Scale(:,:,:) = 0.0 + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) BBL_Mix_Length(:,:,:) = 0.0 + if (CS%debug .or. (CS%id_BBL_Vel_Scale>0)) BBL_Vel_Scale(:,:,:) = 0.0 + if (CS%id_Kd_BBL > 0) Kd_BBL_3d(:,:,:) = 0.0 + if (CS%id_ustar_BBL > 0) diag_ustar_BBL(:,:) = 0.0 + if (CS%id_BBL_decay_scale > 0) diag_BBL_decay_scale(:,:) = 0.0 + endif + + ! CS_tmp is used to test sensitivity to parameter setting changes. + if (CS%options_diff > 0) then + CS_tmp1 = CS ; CS_tmp2 = CS + SpV_scale1 = 1.0 ; SpV_scale2 = 1.0 + + if (CS%options_diff == 1) then + CS_tmp1%orig_PE_calc = .true. ; CS_tmp2%orig_PE_calc = .false. + elseif (CS%options_diff == 2) then + CS_tmp1%answer_date = 20181231 ; CS_tmp2%answer_date = 20240101 + elseif (CS%options_diff == 3) then + CS_tmp1%direct_calc = .true. ; CS_tmp2%direct_calc = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%orig_PE_calc = .false. ; CS_tmp2%orig_PE_calc = .false. + elseif (CS%options_diff == 4) then + CS_tmp1%direct_calc = .true. ; CS_tmp2%direct_calc = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%ePBL_BBL_effic = 0.2 ; CS_tmp2%ePBL_BBL_effic = 0.2 + elseif (CS%options_diff == 5) then + CS_tmp1%decay_adjusted_BBL_TKE = .true. ; CS_tmp2%decay_adjusted_BBL_TKE = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%ePBL_BBL_effic = 0.2 ; CS_tmp2%ePBL_BBL_effic = 0.2 + endif + ! This logic is needed because the scaling of SpV_dt changes with answer date. + if (CS_tmp1%answer_date < 20240101) SpV_scale1 = US%m_to_Z**3 * US%T_to_s**3 + if (CS_tmp2%answer_date < 20240101) SpV_scale2 = US%m_to_Z**3 * US%T_to_s**3 + if (CS%id_opt_diff_Kd_ePBL > 0) diff_Kd(:,:,:) = 0.0 + if (CS%id_opt_maxdiff_Kd_ePBL > 0) max_abs_diff_Kd(:,:) = 0.0 + if (CS%id_opt_diff_hML_depth > 0) diff_hML_depth(:,:) = 0.0 endif - ! if (CS%id_Mixing_Length>0) CS%Mixing_Length(:,:,:) = 0.0 - ! if (CS%id_Velocity_Scale>0) CS%Velocity_Scale(:,:,:) = 0.0 !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,I_dt, & !!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int) @@ -414,7 +560,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS if ((dt > 0.0) .and. GV%Boussinesq .or. .not.allocated(tv%SpV_avg)) then if (CS%answer_date < 20240101) then do K=1,nz+1 - SpV_dt(K) = (US%Z_to_m**3*US%s_to_T**3) / (dt*GV%Rho0) + SpV_dt(K) = 1.0 / (dt*GV%Rho0) enddo else do K=1,nz+1 @@ -457,19 +603,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS endif if (allocated(tv%SpV_avg) .and. .not.GV%Boussinesq) then - if (CS%answer_date < 20240101) then - SpV_dt(1) = (US%Z_to_m**3*US%s_to_T**3) * tv%SpV_avg(i,j,1) * I_dt - do K=2,nz - SpV_dt(K) = (US%Z_to_m**3*US%s_to_T**3) * 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt - enddo - SpV_dt(nz+1) = (US%Z_to_m**3*US%s_to_T**3) * tv%SpV_avg(i,j,nz) * I_dt - else - SpV_dt(1) = tv%SpV_avg(i,j,1) * I_dt - do K=2,nz - SpV_dt(K) = 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt - enddo - SpV_dt(nz+1) = tv%SpV_avg(i,j,nz) * I_dt - endif + SpV_dt(1) = tv%SpV_avg(i,j,1) * I_dt + do K=2,nz + SpV_dt(K) = 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt + enddo + SpV_dt(nz+1) = tv%SpV_avg(i,j,nz) * I_dt endif B_flux = buoy_flux(i,j) @@ -491,77 +629,186 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Perhaps provide a first guess for MLD based on a stored previous value. MLD_io = -1.0 if (CS%MLD_iteration_guess .and. (CS%ML_depth(i,j) > 0.0)) MLD_io = CS%ML_depth(i,j) + BBLD_io = 0.0 + + ! Store the initial guesses at the boundary layer depths for testing sensitivities. + MLD_in = MLD_io + if (CS%answer_date < 20240101) then + do K=1,nz+1 ; SpV_dt_cf(K) = (US%Z_to_m**3*US%s_to_T**3) * SpV_dt(K) ; enddo + else + do K=1,nz+1 ; SpV_dt_cf(K) = SpV_dt(K) ; enddo + endif if (stoch_CS%pert_epbl) then ! stochastics are active - call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_cf, TKE_forcing, B_flux, absf, & u_star, u_star_mean, mech_TKE, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j, & TKE_gen_stoch=stoch_CS%epbl1_wts(i,j), TKE_diss_stoch=stoch_CS%epbl2_wts(i,j)) else - call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_cf, TKE_forcing, B_flux, absf, & u_star, u_star_mean, mech_TKE, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j) endif + ! Add the diffusivity due to bottom boundary layer mixing, if there is energy to drive this mixing. + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j) + BBLD_in = BBLD_io + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) + + do K=1,nz+1 ; Kd(K) = Kd(K) + Kd_BBL(K) ; enddo + if (CS%id_Kd_BBL > 0) then ; do K=1,nz+1 + Kd_BBL_3d(i,j,K) = Kd_BBL(K) + enddo ; endif + if (CS%id_ustar_BBL > 0) diag_ustar_BBL(i,j) = u_star_BBL + if ((CS%id_BBL_decay_scale > 0) .and. (CS%TKE_decay * absf > 0)) & + diag_BBL_decay_scale(i,j) = u_star_BBL / (CS%TKE_decay * absf) + endif + ! Copy the diffusivities to a 2-d array. do K=1,nz+1 Kd_2d(i,K) = Kd(K) enddo CS%ML_depth(i,j) = MLD_io + CS%BBL_depth(i,j) = BBLD_io if (CS%TKE_diagnostics) then - CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + eCD%dTKE_MKE - CS%diag_TKE_conv(i,j) = CS%diag_TKE_conv(i,j) + eCD%dTKE_conv - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + eCD%dTKE_forcing - CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + eCD%dTKE_wind - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) + eCD%dTKE_mixing - CS%diag_TKE_mech_decay(i,j) = CS%diag_TKE_mech_decay(i,j) + eCD%dTKE_mech_decay - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + eCD%dTKE_conv_decay - ! CS%diag_TKE_unbalanced(i,j) = CS%diag_TKE_unbalanced(i,j) + eCD%dTKE_unbalanced + diag_TKE_MKE(i,j) = diag_TKE_MKE(i,j) + eCD%dTKE_MKE + diag_TKE_conv(i,j) = diag_TKE_conv(i,j) + eCD%dTKE_conv + diag_TKE_forcing(i,j) = diag_TKE_forcing(i,j) + eCD%dTKE_forcing + diag_TKE_wind(i,j) = diag_TKE_wind(i,j) + eCD%dTKE_wind + diag_TKE_mixing(i,j) = diag_TKE_mixing(i,j) + eCD%dTKE_mixing + diag_TKE_mech_decay(i,j) = diag_TKE_mech_decay(i,j) + eCD%dTKE_mech_decay + diag_TKE_conv_decay(i,j) = diag_TKE_conv_decay(i,j) + eCD%dTKE_conv_decay + ! diag_TKE_unbalanced(i,j) = diag_TKE_unbalanced(i,j) + eCD%dTKE_unbalanced endif - ! Write to 3-D for outputting Mixing length and velocity scale. - if (CS%id_Mixing_Length>0) then ; do k=1,nz - CS%Mixing_Length(i,j,k) = mixlen(k) + ! Write mixing length and velocity scale to 3-D arrays for diagnostic output + if (CS%debug .or. (CS%id_Mixing_Length > 0)) then ; do K=1,nz+1 + diag_Mixing_Length(i,j,K) = mixlen(K) enddo ; endif - if (CS%id_Velocity_Scale>0) then ; do k=1,nz - CS%Velocity_Scale(i,j,k) = mixvel(k) + if (CS%debug .or. (CS%id_Velocity_Scale > 0)) then ; do K=1,nz+1 + diag_Velocity_Scale(i,j,K) = mixvel(K) enddo ; endif - if (allocated(CS%mstar_mix)) CS%mstar_mix(i,j) = eCD%mstar - if (allocated(CS%mstar_lt)) CS%mstar_lt(i,j) = eCD%mstar_LT - if (allocated(CS%La)) CS%La(i,j) = eCD%LA - if (allocated(CS%La_mod)) CS%La_mod(i,j) = eCD%LAmod + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) then ; do k=1,nz + BBL_Mix_Length(i,j,k) = mixlen_BBL(k) + enddo ; endif + if (CS%debug .or. (CS%id_BBL_Vel_Scale>0)) then ; do k=1,nz + BBL_Vel_Scale(i,j,k) = mixvel_BBL(k) + enddo ; endif + if (CS%id_TKE_BBL>0) & + diag_TKE_BBL(i,j) = diag_TKE_BBL(i,j) + BBL_TKE + endif + if (CS%id_MSTAR_MIX > 0) diag_mStar_mix(i,j) = eCD%mstar + if (CS%id_MSTAR_LT > 0) diag_mStar_lt(i,j) = eCD%mstar_LT + if (CS%id_LA > 0) diag_LA(i,j) = eCD%LA + if (CS%id_LA_MOD > 0) diag_LA_mod(i,j) = eCD%LAmod + if (report_avg_its) then + CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(eCD%OBL_its)) + CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) + if (CS%ePBL_BBL_effic > 0.0) then + CS%sum_its_BBL(1) = CS%sum_its_BBL(1) + real_to_EFP(real(eCD%BBL_its)) + CS%sum_its_BBL(2) = CS%sum_its_BBL(2) + real_to_EFP(1.0) + endif + endif + + if (CS%options_diff > 0) then + ! Call ePBL_column of ePBL_BBL_column with different parameter settings to diagnose sensitivities. + ! These do not change the model state, and are only used for diagnostic purposes. + if (CS%options_diff < 4) then + BLD_1 = MLD_in ; BLD_2 = MLD_in + do K=1,nz+1 ; SpV_dt_tmp(K) = SpV_scale1 * SpV_dt(K) ; enddo + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_tmp, TKE_forcing, & + B_flux, absf, u_star, u_star_mean, mech_TKE, dt, BLD_1, Kd_1, & + mixvel, mixlen, GV, US, CS_tmp1, eCD_tmp, Waves, G, i, j) + do K=1,nz+1 ; SpV_dt_tmp(K) = SpV_scale2 * SpV_dt(K) ; enddo + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_tmp, TKE_forcing, & + B_flux, absf, u_star, u_star_mean, mech_TKE, dt, BLD_2, Kd_2, & + mixvel, mixlen, GV, US, CS_tmp2, eCD_tmp, Waves, G, i, j) + else + BLD_1 = BBLD_in ; BLD_2 = BBLD_in + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_1, BLD_1, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp1, eCD_tmp) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_2, BLD_2, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp2, eCD_tmp) + endif + + if (CS%id_opt_diff_Kd_ePBL > 0) then + do K=1,nz+1 ; diff_Kd(i,j,K) = Kd_1(K) - Kd_2(K) ; enddo + endif + if (CS%id_opt_maxdiff_Kd_ePBL > 0) then + max_abs_diff_Kd(i,j) = 0.0 + do K=1,nz+1 ; max_abs_diff_Kd(i,j) = max(max_abs_diff_Kd(i,j), abs(Kd_1(K) - Kd_2(K))) ; enddo + endif + if (CS%id_opt_diff_hML_depth > 0) diff_hML_depth(i,j) = BLD_1 - BLD_2 + endif + else ! End of the ocean-point part of the i-loop ! For masked points, Kd_int must still be set (to 0) because it has intent out. do K=1,nz+1 ; Kd_2d(i,K) = 0. ; enddo CS%ML_depth(i,j) = 0.0 - endif ; enddo ! Close of i-loop - Note unusual loop order! + CS%BBL_depth(i,j) = 0.0 + endif ; enddo ! Close of i-loop - Note the unusual loop order, with k-loops inside i-loops. do K=1,nz+1 ; do i=is,ie ; Kd_int(i,j,K) = Kd_2d(i,K) ; enddo ; enddo enddo ! j-loop + if (CS%debug .and. (CS%ePBL_BBL_effic > 0.0)) then + call hchksum(visc%TKE_BBL, "ePBL visc%TKE_BBL", G%HI, unscale=GV%H_to_MKS*US%Z_to_m**2*US%s_to_T**3) + call hchksum(visc%ustar_BBL, "ePBL visc%ustar_BBL", G%HI, unscale=GV%H_to_MKS*US%s_to_T) + call hchksum(Kd_int, "End of ePBL Kd_int", G%HI, unscale=GV%H_to_MKS*US%Z_to_m*US%s_to_T) + call hchksum(diag_Velocity_Scale, "ePBL Velocity_Scale", G%HI, unscale=US%Z_to_m*US%s_to_T) + call hchksum(diag_Mixing_Length, "ePBL Mixing_Length", G%HI, unscale=US%Z_to_m) + call hchksum(BBL_Vel_Scale, "ePBL BBL_Vel_Scale", G%HI, unscale=US%Z_to_m*US%s_to_T) + call hchksum(BBL_Mix_Length, "ePBL BBL_Mix_Length", G%HI, unscale=US%Z_to_m) + endif + if (CS%id_ML_depth > 0) call post_data(CS%id_ML_depth, CS%ML_depth, CS%diag) if (CS%id_hML_depth > 0) call post_data(CS%id_hML_depth, CS%ML_depth, CS%diag) - if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, CS%diag_TKE_wind, CS%diag) - if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, CS%diag_TKE_MKE, CS%diag) - if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, CS%diag_TKE_conv, CS%diag) - if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, CS%diag_TKE_forcing, CS%diag) - if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, CS%diag_TKE_mixing, CS%diag) + if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, diag_TKE_wind, CS%diag) + if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, diag_TKE_MKE, CS%diag) + if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, diag_TKE_conv, CS%diag) + if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, diag_TKE_forcing, CS%diag) + if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, diag_TKE_mixing, CS%diag) if (CS%id_TKE_mech_decay > 0) & - call post_data(CS%id_TKE_mech_decay, CS%diag_TKE_mech_decay, CS%diag) + call post_data(CS%id_TKE_mech_decay, diag_TKE_mech_decay, CS%diag) if (CS%id_TKE_conv_decay > 0) & - call post_data(CS%id_TKE_conv_decay, CS%diag_TKE_conv_decay, CS%diag) - if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, CS%Mixing_Length, CS%diag) - if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, CS%Velocity_Scale, CS%diag) - if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, CS%MSTAR_MIX, CS%diag) - if (CS%id_LA > 0) call post_data(CS%id_LA, CS%LA, CS%diag) - if (CS%id_LA_MOD > 0) call post_data(CS%id_LA_MOD, CS%LA_MOD, CS%diag) - if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, CS%MSTAR_LT, CS%diag) + call post_data(CS%id_TKE_conv_decay, diag_TKE_conv_decay, CS%diag) + if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, diag_Mixing_Length, CS%diag) + if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, diag_Velocity_Scale, CS%diag) + if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, diag_mStar_MIX, CS%diag) + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, Kd_BBL_3d, CS%diag) + if (CS%id_BBL_Mix_Length > 0) call post_data(CS%id_BBL_Mix_Length, BBL_Mix_Length, CS%diag) + if (CS%id_BBL_Vel_Scale > 0) call post_data(CS%id_BBL_Vel_Scale, BBL_Vel_Scale, CS%diag) + if (CS%id_ustar_BBL > 0) call post_data(CS%id_ustar_BBL, diag_ustar_BBL, CS%diag) + if (CS%id_BBL_decay_scale > 0) call post_data(CS%id_BBL_decay_scale, diag_BBL_decay_scale, CS%diag) + if (CS%id_TKE_BBL > 0) call post_data(CS%id_TKE_BBL, diag_TKE_BBL, CS%diag) + if (CS%id_TKE_BBL_mixing > 0) call post_data(CS%id_TKE_BBL_mixing, diag_TKE_BBL_mixing, CS%diag) + if (CS%id_TKE_BBL_decay > 0) call post_data(CS%id_TKE_BBL_decay, diag_TKE_BBL_decay, CS%diag) + if (CS%id_BBL_depth > 0) call post_data(CS%id_BBL_depth, CS%BBL_depth, CS%diag) + endif + if (CS%id_LA > 0) call post_data(CS%id_LA, diag_LA, CS%diag) + if (CS%id_LA_MOD > 0) call post_data(CS%id_LA_MOD, diag_LA_MOD, CS%diag) + if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, diag_mStar_LT, CS%diag) if (stoch_CS%pert_epbl) then if (stoch_CS%id_epbl1_wts > 0) call post_data(stoch_CS%id_epbl1_wts, stoch_CS%epbl1_wts, CS%diag) if (stoch_CS%id_epbl2_wts > 0) call post_data(stoch_CS%id_epbl2_wts, stoch_CS%epbl2_wts, CS%diag) endif + if (CS%options_diff > 0) then + ! These diagnostics are only for determining sensitivities to different ePBL settings. + if (CS%id_opt_diff_Kd_ePBL > 0) call post_data(CS%id_opt_diff_Kd_ePBL, diff_Kd, CS%diag) + if (CS%id_opt_maxdiff_Kd_ePBL > 0) call post_data(CS%id_opt_maxdiff_Kd_ePBL, max_abs_diff_Kd, CS%diag) + if (CS%id_opt_diff_hML_depth > 0) call post_data(CS%id_opt_diff_hML_depth, diff_hML_depth, CS%diag) + endif + end subroutine energetic_PBL @@ -591,7 +838,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !! divided by dt or 1.0 / (dt * Rho0), times conversion !! factors for answer dates before 20240101 in !! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without - !! the convsersion factors for answer dates of + !! the conversion factors for answer dates of !! 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], !! used to convert local TKE into a turbulence !! velocity cubed. @@ -618,14 +865,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !! [Z T-1 ~> m s-1]. real, dimension(SZK_(GV)+1), & intent(out) :: mixlen !< The mixing length scale used in Kd [Z ~> m]. - type(energetic_PBL_CS), intent(inout) :: CS !< Energetic PBL control structure + type(energetic_PBL_CS), intent(in) :: CS !< Energetic PBL control structure type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. type(wave_parameters_CS), pointer :: Waves !< Waves control structure for Langmuir turbulence - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + integer, intent(in) :: i !< The i-index to work on (used for Waves) + integer, intent(in) :: j !< The j-index to work on (used for Waves) real, optional, intent(in) :: TKE_gen_stoch !< random factor used to perturb TKE generation [nondim] real, optional, intent(in) :: TKE_diss_stoch !< random factor used to perturb TKE dissipation [nondim] - integer, intent(in) :: i !< The i-index to work on (used for Waves) - integer, intent(in) :: j !< The i-index to work on (used for Waves) ! This subroutine determines the diffusivities in a single column from the integrated energetics ! planetary boundary layer (ePBL) model. It assumes that heating, cooling and freshwater fluxes @@ -684,6 +931,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, Se, & ! Estimated final values of S in the column [S ~> ppt]. dTe, & ! Running (1-way) estimates of temperature change [C ~> degC]. dSe, & ! Running (1-way) estimates of salinity change [S ~> ppt]. + hp_a, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers above [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in a downward-oriented tridiagonal solver. Th_a, & ! An effective temperature times a thickness in the layer above, including implicit ! mixing effects with other yet higher layers [C H ~> degC m or degC kg m-2]. Sh_a, & ! An effective salinity times a thickness in the layer above, including implicit @@ -698,12 +948,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! the boundary layer [nondim]. h_dz_int, & ! The ratio of the layer thicknesses over the vertical distances ! across the layers surrounding an interface [H Z-1 ~> nondim or kg m-3] - Kddt_h ! The diapycnal diffusivity times a timestep divided by the - ! average thicknesses around a layer [H ~> m or kg m-2]. + Kddt_h ! The total diapycnal diffusivity at an interface times a timestep divided by the + ! average thicknesses around an interface [H ~> m or kg m-2]. real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. - real :: hp_a ! An effective pivot thickness of the layer including the effects - ! of coupling with layers above [H ~> m or kg m-2]. This is the first term - ! in the denominator of b1 in a downward-oriented tridiagonal solver. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -757,8 +1004,8 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, real :: dMKE_src_dK ! The partial derivative of MKE_src with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2] - real :: Kddt_h_g0 ! The first guess diapycnal diffusivity times a timestep divided - ! by the average thicknesses around a layer [H ~> m or kg m-2]. + real :: Kddt_h_g0 ! The first guess of the change in diapycnal diffusivity times a timestep + ! divided by the average thicknesses around an interface [H ~> m or kg m-2]. real :: PE_chg_max ! The maximum PE change for very large values of Kddt_h(K) [R Z3 T-2 ~> J m-2]. real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K) ! for very small values of Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. @@ -806,12 +1053,18 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! can improve this. real :: dMLD_min ! The change in diagnosed mixed layer depth when the guess is min_MLD [Z ~> m] real :: dMLD_max ! The change in diagnosed mixed layer depth when the guess is max_MLD [Z ~> m] - logical :: OBL_converged ! Flag for convergence of MLD integer :: OBL_it ! Iteration counter + real :: TKE_used ! The TKE used to support mixing at an interface [R Z3 T-2 ~> J m-2]. + ! real :: Kd_add ! The additional diffusivity at an interface [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: frac_in_BL ! The fraction of the energy required to support dKd_max that is suppiled by + ! max_PE_chg, used here to determine a fractional layer contribution to the + ! boundary layer thickness [nondim] real :: Surface_Scale ! Surface decay scale for vstar [nondim] logical :: calc_Te ! If true calculate the expected final temperature and salinity values. + logical :: no_MKE_conversion ! If true, there is no conversion from MKE to TKE, so a simpler solver can be used. logical :: debug ! This is used as a hard-coded value for debugging. + logical :: convectively_unstable ! If true, there is convective instability at an interface. ! The following arrays are used only for debugging purposes. real :: dPE_debug ! An estimate of the potential energy change [R Z3 T-2 ~> J m-2] @@ -837,6 +1090,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, debug = .false. ! Change this hard-coded value for debugging. calc_Te = (debug .or. (.not.CS%orig_PE_calc)) + no_MKE_conversion = ((CS%direct_calc) .and. (CS%MKE_to_TKE_effic == 0.0)) h_neglect = GV%H_subroundoff dz_neglect = GV%dZ_subroundoff @@ -905,13 +1159,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif ! Iterate to determine a converged EPBL depth. - OBL_converged = .false. do OBL_it=1,CS%Max_MLD_Its - if (.not. OBL_converged) then - ! If not using MLD_Iteration flag loop to only execute once. - if (.not.CS%Use_MLD_iteration) OBL_converged = .true. - + !### The old indenting is being retained for now to simplify the review of some pending changes. if (debug) then ; mech_TKE_k(:) = 0.0 ; conv_PErel_k(:) = 0.0 ; endif ! Reset ML_depth @@ -923,7 +1173,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess), u_star_mean, i, j, dz, Waves, & U_H=u, V_H=v) call find_mstar(CS, US, B_flux, u_star, MLD_guess, absf, & - MStar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,& + mstar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,& mstar_LT=mstar_LT) else call find_mstar(CS, US, B_flux, u_star, MLD_guess, absf, mstar_total) @@ -931,10 +1181,10 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !/ Apply MStar to get mech_TKE if ((CS%answer_date < 20190101) .and. (CS%mstar_scheme==Use_Fixed_MStar)) then - mech_TKE = (dt*MSTAR_total*GV%Rho0) * u_star**3 + mech_TKE = (dt*mstar_total*GV%Rho0) * u_star**3 else - mech_TKE = MSTAR_total * mech_TKE_in - ! mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) + mech_TKE = mstar_total * mech_TKE_in + ! mech_TKE = mstar_total * (dt*GV%Rho0* u_star**3) endif ! stochastically perturb mech_TKE in the UFS if (present(TKE_gen_stoch)) mech_TKE = mech_TKE*TKE_gen_stoch @@ -996,7 +1246,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif Kd(1) = 0.0 ; Kddt_h(1) = 0.0 - hp_a = h(1) + hp_a(1) = h(1) dT_to_dPE_a(1) = dT_to_dPE(1) ; dT_to_dColHt_a(1) = dT_to_dColHt(1) dS_to_dPE_a(1) = dS_to_dPE(1) ; dS_to_dColHt_a(1) = dS_to_dColHt(1) @@ -1121,14 +1371,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! tridiagonal solver for the whole column to be completed for debugging ! purposes, and also allows for something akin to convective adjustment ! in unstable interior regions? - b1 = 1.0 / hp_a + b1 = 1.0 / hp_a(k-1) c1(K) = 0.0 if (CS%orig_PE_calc) then dTe(k-1) = b1 * ( dTe_t2 ) dSe(k-1) = b1 * ( dSe_t2 ) endif - hp_a = h(k) + hp_a(k) = h(k) dT_to_dPE_a(k) = dT_to_dPE(k) dS_to_dPE_a(k) = dS_to_dPE(k) dT_to_dColHt_a(k) = dT_to_dColHt(k) @@ -1145,12 +1395,12 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, dS_km1_t2 = (S0(k)-S0(k-1)) else dT_km1_t2 = (T0(k)-T0(k-1)) - & - (Kddt_h(K-1) / hp_a) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + (Kddt_h(K-1) / hp_a(k-1)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) dS_km1_t2 = (S0(k)-S0(k-1)) - & - (Kddt_h(K-1) / hp_a) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + (Kddt_h(K-1) / hp_a(k-1)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) endif - dTe_term = dTe_t2 + hp_a * (T0(k-1)-T0(k)) - dSe_term = dSe_t2 + hp_a * (S0(k-1)-S0(k)) + dTe_term = dTe_t2 + hp_a(k-1) * (T0(k-1)-T0(k)) + dSe_term = dSe_t2 + hp_a(k-1) * (S0(k-1)-S0(k)) else if (K<=2) then Th_a(k-1) = h(k-1) * T0(k-1) ; Sh_a(k-1) = h(k-1) * S0(k-1) @@ -1205,7 +1455,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif endif hbs_here = min(hb_hs(K), MixLen_shape(K)) - mixlen(K) = MAX(CS%min_mix_len, ((dz_tt*hbs_here)*vstar) / & + mixlen(K) = max(CS%min_mix_len, ((dz_tt*hbs_here)*vstar) / & ((CS%Ekman_scale_coef * absf) * (dz_tt*hbs_here) + vstar)) !Note setting Kd_guess0 to vstar * CS%vonKar * mixlen(K) here will ! change the answers. Therefore, skipping that. @@ -1221,26 +1471,40 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, mixvel(K) = vstar ! Track vstar Kddt_h_g0 = Kd_guess0 * dt_h - if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kddt_h_g0, h(k), hp_a, dTe_term, dSe_term, & + if (no_MKE_conversion) then + ! Without conversion from MKE to TKE, the updated diffusivity can be determined directly. + ! Replace h(k) with hp_b(k) = h(k), and dT_to_dPE with dT_to_dPE_b, etc., for a 2-direction solver. + call find_Kd_from_PE_chg(0.0, Kd_guess0, dt_h, tot_TKE, hp_a(k-1), h(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt(k), dS_to_dColHt(k), Kd_add=Kd(K), PE_chg=TKE_used, & + dPE_max=PE_chg_max, frac_dKd_max_PE=frac_in_BL) + convectively_unstable = (PE_chg_max < 0.0) + PE_chg_g0 = TKE_used ! This is only used in the convectively unstable limit. + MKE_src = 0.0 + elseif (CS%orig_PE_calc) then + call find_PE_chg_orig(Kddt_h_g0, h(k), hp_a(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), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) + convectively_unstable = (PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0)) + MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) else - call find_PE_chg(0.0, Kddt_h_g0, hp_a, h(k), & + call find_PE_chg(0.0, Kddt_h_g0, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & dT_to_dColHt(k), dS_to_dColHt(k), & PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) + convectively_unstable = (PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0)) + MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) endif - MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) - ! This block checks out different cases to determine Kd at the present interface. - if ((PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0))) then + if (convectively_unstable) then ! This column is convectively unstable. if (PE_chg_max <= 0.0) then ! Does MKE_src need to be included in the calculation of vstar here? @@ -1280,14 +1544,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, mixvel(K) = vstar if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kd(K)*dt_h, h(k), hp_a, dTe_term, dSe_term, & + call find_PE_chg_orig(Kd(K)*dt_h, h(k), hp_a(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), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=dPE_conv) else - call find_PE_chg(0.0, Kd(K)*dt_h, hp_a, h(k), & + call find_PE_chg(0.0, Kd(K)*dt_h, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & @@ -1316,6 +1580,31 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif Kddt_h(K) = Kd(K) * dt_h + + elseif (no_MKE_conversion) then ! (PE_chg_max >= 0.0) and use the diffusivity from find_Kd_from_PE_chg. + ! Kd(K) and TKE_used were already set by find_Kd_from_PE_chg. + + ! frac_in_BL = min((TKE_used / PE_chg_g0), 1.0) + if (sfc_connected) MLD_output = MLD_output + frac_in_BL*dz(k) + if (frac_in_BL < 1.0) sfc_disconnect = .true. + + ! Reduce the mechanical and convective TKE proportionately. + TKE_reduc = 0.0 ! tot_TKE could be 0 if Convectively_stable is false. + if ((tot_TKE > 0.0) .and. (tot_TKE > TKE_used)) TKE_reduc = (tot_TKE - TKE_used) / tot_TKE + + ! All TKE should have been consumed. + if (CS%TKE_diagnostics) then + eCD%dTKE_mixing = eCD%dTKE_mixing - TKE_used * I_dtdiag + eCD%dTKE_conv_decay = eCD%dTKE_conv_decay + & + (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel * I_dtdiag + endif + + tot_TKE = tot_TKE - TKE_used + mech_TKE = TKE_reduc*mech_TKE + conv_PErel = TKE_reduc*conv_PErel + + Kddt_h(K) = Kd(K) * dt_h + elseif (tot_TKE + (MKE_src - PE_chg_g0) >= 0.0) then ! This column is convectively stable and there is energy to support the suggested ! mixing. Keep that estimate. @@ -1366,19 +1655,19 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif do itt=1,max_itt if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kddt_h_guess, h(k), hp_a, dTe_term, dSe_term, & + call find_PE_chg_orig(Kddt_h_guess, h(k), hp_a(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), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=PE_chg, dPEc_dKd=dPEc_dKd ) else - call find_PE_chg(0.0, Kddt_h_guess, hp_a, h(k), & + call find_PE_chg(0.0, Kddt_h_guess, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & dT_to_dColHt(k), dS_to_dColHt(k), & - PE_chg=dPE_conv, dPEc_dKd=dPEc_dKd) + PE_chg=PE_chg, dPEc_dKd=dPEc_dKd) endif MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) @@ -1445,14 +1734,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, Kddt_h(K) = Kd(K) * dt_h ! At this point, the final value of Kddt_h(K) is known, so the ! estimated properties for layer k-1 can be calculated. - b1 = 1.0 / (hp_a + Kddt_h(K)) + b1 = 1.0 / (hp_a(k-1) + Kddt_h(K)) c1(K) = Kddt_h(K) * b1 if (CS%orig_PE_calc) then dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) endif - hp_a = h(k) + (hp_a * b1) * Kddt_h(K) + hp_a(k) = h(k) + (hp_a(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) dT_to_dColHt_a(k) = dT_to_dColHt(k) + c1(K)*dT_to_dColHt_a(k-1) @@ -1476,7 +1765,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif if (calc_Te) then - if (k==2) then + if (K==2) then Te(1) = b1*(h(1)*T0(1)) Se(1) = b1*(h(1)*S0(1)) else @@ -1489,7 +1778,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, if (debug) then ! Complete the tridiagonal solve for Te. - b1 = 1.0 / hp_a + b1 = 1.0 / hp_a(nz) Te(nz) = b1 * (h(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) Se(nz) = b1 * (h(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) dT_expect(nz) = Te(nz) - T0(nz) ; dS_expect(nz) = Se(nz) - S0(nz) @@ -1508,25 +1797,39 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, enddo mixing_debug = dPE_debug * I_dtdiag endif - k = nz ! This is here to allow a breakpoint to be set. - !/BGR - ! The following lines are used for the iteration - ! note the iteration has been altered to use the value predicted by - ! the TKE threshold (ML_DEPTH). This is because the MSTAR - ! is now dependent on the ML, and therefore the ML needs to be estimated - ! more precisely than the grid spacing. - - !New method uses ML_DEPTH as computed in ePBL routine + + if (OBL_it >= CS%Max_MLD_Its) exit + + ! The following lines are used for the iteration. Note the iteration has been altered + ! to use the value predicted by the TKE threshold (ML_DEPTH). This is because the MSTAR + ! is now dependent on the ML, and therefore the ML needs to be estimated more precisely + ! than the grid spacing. + + ! New method uses ML_DEPTH as computed in ePBL routine MLD_found = MLD_output - if (MLD_found - MLD_guess > CS%MLD_tol) then - min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess - elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then - OBL_converged = .true. ! Break convergence loop - else ! We know this guess was too deep - max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + + ! Find out whether to do another iteration and the new bounds on it. + if (CS%MLD_iter_bug) then + ! There is a bug in the logic here if (MLD_found - MLD_guess == CS%MLD_tol). + if (MLD_found - MLD_guess > CS%MLD_tol) then + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then + exit ! Break the MLD convergence loop + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + endif + else + if (abs(MLD_found - MLD_guess) < CS%MLD_tol) then + exit ! Break the MLD convergence loop + elseif (MLD_found > MLD_guess) then ! This guess was too shallow + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + endif endif - if (.not.OBL_converged) then ; if (CS%MLD_bisection) then + if (OBL_it < CS%Max_MLD_Its) then + if (CS%MLD_bisection) then ! For the next pass, guess the average of the minimum and maximum values. MLD_guess = 0.5*(min_MLD + max_MLD) else ! Try using the false position method or the returned value instead of simple bisection. @@ -1535,22 +1838,18 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! Both bounds have valid change estimates and are probably in the range of possible outputs. MLD_guess = (dMLD_min*max_MLD - dMLD_max*min_MLD) / (dMLD_min - dMLD_max) elseif ((MLD_found > min_MLD) .and. (MLD_found < max_MLD)) then - ! The output MLD_found is an interesting guess, as it likely to bracket the true solution + ! The output MLD_found is an interesting guess, as it is likely to bracket the true solution ! along with the previous value of MLD_guess and to be close to the solution. MLD_guess = MLD_found else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. MLD_guess = 0.5*(min_MLD + max_MLD) endif - endif ; endif - endif - if ((OBL_converged) .or. (OBL_it==CS%Max_MLD_Its)) then - if (report_avg_its) then - CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(OBL_it)) - CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) endif - exit endif + enddo ! Iteration loop for converged boundary layer thickness. + + eCD%OBL_its = min(OBL_it, CS%Max_MLD_Its) if (CS%Use_LT) then eCD%LA = LA ; eCD%LAmod = LAmod ; eCD%mstar = mstar_total ; eCD%mstar_LT = mstar_LT else @@ -1561,150 +1860,1159 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, end subroutine ePBL_column -!> This subroutine calculates the change in potential energy and or derivatives -!! for several changes in an interface'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_Z, 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, PE_ColHt_cor) - real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times - !! the time step and divided by the average of the - !! thicknesses around the interface [H ~> m or kg m-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 [H ~> m or kg m-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 [H ~> m or kg m-2]. - 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 [H ~> m or kg m-2]. - real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer - !! above, including implicit mixing effects with other - !! yet higher layers [C H ~> degC m or degC kg m-2]. - real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer - !! above, including implicit mixing effects with other - !! yet higher layers [S H ~> ppt m or ppt kg m-2]. - real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer - !! below, including implicit mixing effects with other - !! yet lower layers [C H ~> degC m or degC kg m-2]. - real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer - !! below, including implicit mixing effects with other - !! yet lower layers [S H ~> ppt m or ppt kg m-2]. - 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> J m-2 ppt-1]. - real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates - !! the changes in column thickness to the energy that is radiated - !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. - 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 [Z C-1 ~> m degC-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 [Z S-1 ~> 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 [Z C-1 ~> m degC-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 [Z S-1 ~> m ppt-1]. - - real, intent(out) :: PE_chg !< The change in column potential energy from applying - !! Kddt_h at the present interface [R Z3 T-2 ~> J m-2]. - real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h - !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. - real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could - !! be realized by applying a huge value of Kddt_h at the - !! present interface [R Z3 T-2 ~> 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 [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. - real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net - !! change in the column height [R Z3 T-2 ~> J m-2]. - - ! Local variables - real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. - real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. - real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. - real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. - real :: PEc_core ! The diffusivity-independent core term in the expressions - ! for the potential energy changes [R Z2 T-2 ~> J m-3]. - real :: ColHt_core ! The diffusivity-independent core term in the expressions - ! for the column height changes [H Z ~> m2 or kg m-1]. - real :: ColHt_chg ! The change in the column height [H ~> m or kg m-2]. - real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3]. - real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4]. - - ! 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 - ! succinct form. The derivation is not necessarily obvious, but it demonstrably - ! works by comparison with separate calculations of the energy changes after - ! the tridiagonal solver for the final changes in temperature and salinity are - ! applied. - - hps = hp_a + hp_b - bdt1 = hp_a * hp_b + Kddt_h0 * hps - dT_c = hp_a * Th_b - hp_b * Th_a - dS_c = hp_a * Sh_b - hp_b * Sh_a - PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & - hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) - ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & - hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) - - ! Find the change in column potential energy due to the change in the - ! diffusivity at this interface by dKddt_h. - y1_3 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) - PE_chg = PEc_core * y1_3 - ColHt_chg = ColHt_core * y1_3 - if (ColHt_chg < 0.0) PE_chg = PE_chg - pres_Z * ColHt_chg - if (present(PE_ColHt_cor)) PE_ColHt_cor = -pres_Z * min(ColHt_chg, 0.0) +!> This subroutine determines the diffusivities from a bottom boundary layer version of +!! the integrated energetics mixed layer model for a single column of water. +subroutine ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, absf, & + dt, Kd, BBL_TKE_in, u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]. + real, dimension(SZK_(GV)), intent(in) :: u !< Zonal velocities interpolated to h points + !! [L T-1 ~> m s-1]. + real, dimension(SZK_(GV)), intent(in) :: v !< Zonal velocities interpolated to h points + !! [L T-1 ~> m s-1]. + real, dimension(SZK_(GV)), intent(in) :: T0 !< The initial layer temperatures [C ~> degC]. + real, dimension(SZK_(GV)), intent(in) :: S0 !< The initial layer salinities [S ~> ppt]. - if (present(dPEc_dKd)) then - ! Find the derivative of the potential energy change with dKddt_h. - y1_4 = 1.0 / (bdt1 + dKddt_h * hps)**2 - dPEc_dKd = PEc_core * y1_4 - ColHt_chg = ColHt_core * y1_4 - if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres_Z * ColHt_chg - endif + real, dimension(SZK_(GV)), intent(in) :: dSV_dT !< The partial derivative of in-situ specific + !! volume with potential temperature + !! [R-1 C-1 ~> m3 kg-1 degC-1]. + real, dimension(SZK_(GV)), intent(in) :: dSV_dS !< The partial derivative of in-situ specific + !! volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. + real, dimension(SZK_(GV)+1), intent(in) :: SpV_dt !< Specific volume interpolated to interfaces + !! divided by dt (if non-Boussinesq) or + !! 1.0 / (dt * Rho0), in [R-1 T-1 ~> m3 kg-1 s-1], + !! used to convert local TKE into a turbulence + !! velocity cubed. + real, intent(in) :: absf !< The absolute value of the Coriolis parameter [T-1 ~> s-1]. + real, intent(in) :: dt !< Time increment [T ~> s]. + real, dimension(SZK_(GV)+1), & + intent(in) :: Kd !< The diffusivities at interfaces due to previously + !! applied mixing processes [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: BBL_TKE_in !< The mechanically generated turbulent + !! kinetic energy available for bottom boundary + !! layer mixing within a time step [R Z3 T-2 ~> J m-2]. + real, intent(in) :: u_star_BBL !< The bottom boundary layer friction velocity + !! in thickuness flux units [H T-1 ~> m s-1 or kg m-2 s-1] + real, dimension(SZK_(GV)+1), & + intent(out) :: Kd_BBL !< The bottom boundary layer contribution to diffusivities + !! at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(inout) :: BBLD_io !< A first guess at the bottom boundary layer depth on input, and + !! the calculated bottom boundary layer depth on output [Z ~> m] + real, dimension(SZK_(GV)+1), & + intent(out) :: mixvel_BBL !< The profile of boundary layer turbulent mixing + !! velocities [Z T-1 ~> m s-1] + real, dimension(SZK_(GV)+1), & + intent(out) :: mixlen_BBL !< The profile of bottom boundary layer turbulent + !! mixing lengths [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(energetic_PBL_CS), intent(in) :: CS !< Energetic PBL control structure + type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. - if (present(dPE_max)) then - ! This expression is the limit of PE_chg for infinite dKddt_h. - y1_3 = 1.0 / (bdt1 * hps) - dPE_max = PEc_core * y1_3 - ColHt_chg = ColHt_core * y1_3 - if (ColHt_chg < 0.0) dPE_max = dPE_max - pres_Z * ColHt_chg - endif +! This subroutine determines the contributions from diffusivities in a single column from a +! bottom-boundary layer adaptation of the integrated energetics planetary boundary layer (ePBL) +! model. It accounts for the possibility that the surface boundary diffusivities have already +! been determined. All calculations are done implicitly, and there is no stability limit on the +! time step. Only mechanical mixing in the bottom boundary layer is considered. (Geothermal heat +! fluxes are addressed elsewhere in the MOM6 code, and convection throughout the water column is +! handled by the surface version of ePBL.) There is no conversion of released mean kinetic energy +! into bottom boundary layer turbulent kinetic energy (at least for now), apart from the explicit +! energy that is supplied as an argument to this routine. - if (present(dPEc_dKd_0)) then - ! This expression is the limit of dPEc_dKd for dKddt_h = 0. - y1_4 = 1.0 / bdt1**2 - dPEc_dKd_0 = PEc_core * y1_4 - ColHt_chg = ColHt_core * y1_4 - if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres_Z * ColHt_chg - endif + ! Local variables + real, dimension(SZK_(GV)+1) :: & + pres_Z, & ! Interface pressures with a rescaling factor to convert interface height + ! movements into changes in column potential energy [R Z2 T-2 ~> kg m-1 s-2]. + dztop_dztot ! The distance from the surface divided by the thickness of the + ! water column [nondim]. + real :: mech_BBL_TKE ! The mechanically generated turbulent kinetic energy available for + ! bottom boundary layer mixing within a time step [R Z3 T-2 ~> J m-2]. + real :: TKE_eff_avail ! The turbulent kinetic energy that is effectively available to drive mixing + ! after any effects of exponentially decay have been taken into account + ! [R Z3 T-2 ~> J m-2] + real :: TKE_eff_used ! The amount of TKE_eff_avail that has been used to drive mixing [R Z3 T-2 ~> J m-2] + real :: htot ! The total thickness of the layers above an interface [H ~> m or kg m-2]. + real :: dztot ! The total depth of the layers above an interface [Z ~> m]. + real :: uhtot ! The depth integrated zonal velocities in the layers above [H L T-1 ~> m2 s-1 or kg m-1 s-1] + real :: vhtot ! The depth integrated meridional velocities in the layers above [H L T-1 ~> m2 s-1 or kg m-1 s-1] + real :: Idecay_len_TKE ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1]. + real :: dz_sum ! The total thickness of the water column [Z ~> m]. + + real, dimension(SZK_(GV)) :: & + dT_to_dColHt, & ! Partial derivative of the total column height with the temperature changes + ! within a layer [Z C-1 ~> m degC-1]. + dS_to_dColHt, & ! Partial derivative of the total column height with the salinity changes + ! within a layer [Z S-1 ~> m ppt-1]. + dT_to_dPE, & ! Partial derivatives of column potential energy with the temperature + ! changes within a layer, in [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE, & ! Partial derivatives of column potential energy with the salinity changes + ! within a layer, in [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + dT_to_dColHt_a, & ! Partial derivative of the total column height with the temperature changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [Z C-1 ~> m degC-1]. + dS_to_dColHt_a, & ! Partial derivative of the total column height with the salinity changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [Z S-1 ~> m ppt-1]. + dT_to_dPE_a, & ! Partial derivatives of column potential energy with the temperature changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE_a, & ! Partial derivative of column potential energy with the salinity changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + dT_to_dColHt_b, & ! Partial derivative of the total column height with the temperature changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [Z C-1 ~> m degC-1]. + dS_to_dColHt_b, & ! Partial derivative of the total column height with the salinity changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [Z S-1 ~> m ppt-1]. + dT_to_dPE_b, & ! Partial derivatives of column potential energy with the temperature changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE_b, & ! Partial derivative of column potential energy with the salinity changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + c1, & ! c1 is used by the tridiagonal solver [nondim]. + Te, & ! Estimated final values of T in the column [C ~> degC]. + Se, & ! Estimated final values of S in the column [S ~> ppt]. + dTe, & ! Running (1-way) estimates of temperature change [C ~> degC]. + dSe, & ! Running (1-way) estimates of salinity change [S ~> ppt]. + hp_a, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers above [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in a downward-oriented tridiagonal solver. + hp_b, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers below [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in an upward-oriented tridiagonal solver. + Th_a, & ! An effective temperature times a thickness in the layer above, including implicit + ! mixing effects with other yet higher layers [C H ~> degC m or degC kg m-2]. + Sh_a, & ! An effective salinity times a thickness in the layer above, including implicit + ! mixing effects with other yet higher layers [S H ~> ppt m or ppt kg m-2]. + Th_b, & ! An effective temperature times a thickness in the layer below, including implicit + ! mixing effects with other yet lower layers [C H ~> degC m or degC kg m-2]. + Sh_b ! An effective salinity times a thickness in the layer below, including implicit + ! mixing effects with other yet lower layers [S H ~> ppt m or ppt kg m-2]. + real, dimension(SZK_(GV)+1) :: & + MixLen_shape, & ! A nondimensional shape factor for the mixing length that + ! gives it an appropriate asymptotic value at the bottom of + ! the boundary layer [nondim]. + h_dz_int, & ! The ratio of the layer thicknesses over the vertical distances + ! across the layers surrounding an interface [H Z-1 ~> nondim or kg m-3] + Kddt_h ! The total diapycnal diffusivity at an interface times a timestep divided by the + ! average thicknesses around an interface [H ~> m or kg m-2]. + real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2]. + real :: dz_neglect ! A vertical distance that is so small it is usually lost + ! in roundoff and can be neglected [Z ~> m]. + real :: dMass ! The mass per unit area within a layer [Z R ~> kg m-2]. + real :: dPres ! The hydrostatic pressure change across a layer [R Z2 T-2 ~> Pa = J m-3]. + + real :: dt_h ! The timestep divided by the averages of the vertical distances around + ! a layer [T Z-1 ~> s m-1]. + real :: dz_top ! The distance from the surface [Z ~> m]. + real :: dz_rsum ! The distance from the seafloor [Z ~> m]. + real :: I_dzsum ! The inverse of dz_sum [Z-1 ~> m-1]. + real :: I_BBLD ! The inverse of the current value of BBLD [Z-1 ~> m-1]. + real :: dz_tt ! The distance from the surface or up to the next interface + ! that did not exhibit turbulent mixing from this scheme plus + ! a surface mixing roughness length given by dz_tt_min [Z ~> m]. + real :: dz_tt_min ! A surface roughness length [Z ~> m]. + real :: C1_3 ! = 1/3 [nondim] + real :: vstar ! An in-situ turbulent velocity [Z T-1 ~> m s-1]. + real :: BBLD_output ! The bottom boundary layer depth output from this routine [Z ~> m] + real :: hbs_here ! The local minimum of dztop_dztot and MixLen_shape [nondim] + real :: TKE_used ! The TKE used to support mixing at an interface [R Z3 T-2 ~> J m-2]. + real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2] + real :: Kddt_h_g0 ! The first guess of the change in diapycnal diffusivity times a timestep + ! divided by the average thicknesses around an interface [H ~> m or kg m-2]. + real :: Kddt_h_prev ! The diapycnal diffusivity before modification times a timestep divided + ! by the average thicknesses around an interface [H ~> m or kg m-2]. + real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K) + ! for very small values of Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real :: PE_chg ! The change in potential energy due to mixing at an + ! interface [R Z3 T-2 ~> J m-2], positive for the column increasing + ! in potential energy (i.e., consuming TKE). + real :: TKE_left ! The amount of turbulent kinetic energy left for the most + ! recent guess at Kddt_h(K) [R Z3 T-2 ~> J m-2]. + real :: dPEc_dKd ! The partial derivative of PE_chg with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real :: TKE_left_min, TKE_left_max ! Maximum and minimum values of TKE_left [R Z3 T-2 ~> J m-2]. + real :: Kddt_h_max, Kddt_h_min ! Maximum and minimum values of Kddt_h(K) [H ~> m or kg m-2]. + real :: Kddt_h_guess ! A guess at the value of Kddt_h(K) [H ~> m or kg m-2]. + real :: Kddt_h_next ! The next guess at the value of Kddt_h(K) [H ~> m or kg m-2]. + real :: dKddt_h ! The change between guesses at Kddt_h(K) [H ~> m or kg m-2]. + real :: dKddt_h_Newt ! The change between guesses at Kddt_h(K) with Newton's method [H ~> m or kg m-2]. + real :: Kddt_h_newt ! The Newton's method next guess for Kddt_h(K) [H ~> m or kg m-2]. + real :: exp_kh ! The nondimensional decay of TKE across a layer [nondim]. + real :: frac_in_BL ! The fraction of the energy required to support dKd_max that is suppiled by + ! max_PE_chg, used here to determine a fractional layer contribution to the + ! boundary layer thickness [nondim] + real :: TKE_rescale ! The effective fractional increase in energy available to + ! mixing at this interface once its exponential decay is accounted for [nondim] + logical :: use_Newt ! Use Newton's method for the next guess at Kddt_h(K). + logical :: convectively_stable ! If true the water column is convectively stable at this interface. + logical :: bot_connected ! If true the ocean is actively turbulent from the present + ! interface all the way down to the seafloor. + logical :: bot_disconnect ! If true, any turbulence has become disconnected + ! from the bottom. + + ! The following is only used for diagnostics. + real :: I_dtdiag ! = 1.0 / dt [T-1 ~> s-1]. + + real :: BBLD_guess, BBLD_found ! Bottom boundary layer depth guessed/found for iteration [Z ~> m] + real :: min_BBLD, max_BBLD ! Iteration bounds on BBLD [Z ~> m], which are adjusted at each step + real :: dBBLD_min ! The change in diagnosed mixed layer depth when the guess is min_BLD [Z ~> m] + real :: dBBLD_max ! The change in diagnosed mixed layer depth when the guess is max_BLD [Z ~> m] + logical :: BBL_converged ! Flag for convergence of BBLD + integer :: BBL_it ! Iteration counter + + real :: Surface_Scale ! Surface decay scale for vstar [nondim] + logical :: debug ! This is used as a hard-coded value for debugging. + logical :: no_MKE_conversion ! If true, there is conversion of MKE to TKE in this routine. + + ! The following arrays are used only for debugging purposes. + real :: dPE_debug ! An estimate of the potential energy change [R Z3 T-2 ~> J m-2] + real :: mixing_debug ! An estimate of the rate of change of potential energy due to mixing [R Z3 T-3 ~> W m-2] + real, dimension(20) :: TKE_left_itt ! The value of TKE_left after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(20) :: PE_chg_itt ! The value of PE_chg after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(20) :: Kddt_h_itt ! The value of Kddt_h_guess after each iteration [H ~> m or kg m-2] + real, dimension(20) :: dPEa_dKd_itt ! The value of dPEc_dKd after each iteration [R Z3 T-2 H-1 ~> J m-3 or J kg-1] +! real, dimension(20) :: MKE_src_itt ! The value of MKE_src after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(SZK_(GV)) :: dT_expect !< Expected temperature changes [C ~> degC] + real, dimension(SZK_(GV)) :: dS_expect !< Expected salinity changes [S ~> ppt] + real, dimension(SZK_(GV)) :: mech_BBL_TKE_k ! The mechanically generated turbulent kinetic energy + ! available for bottom boundary mixing over a time step for each layer [R Z3 T-2 ~> J m-2]. + integer, dimension(SZK_(GV)) :: num_itts + + integer :: k, nz, itt, max_itt + + nz = GV%ke + + debug = .false. ! Change this hard-coded value for debugging. + no_MKE_conversion = ((CS%direct_calc) ) ! .and. (CS%MKE_to_TKE_effic == 0.0)) + + ! Add bottom boundary layer mixing if there is energy to support it. + if ((CS%ePBL_BBL_effic <= 0.0) .or. (BBL_TKE_in <= 0.0)) then + ! There is no added bottom boundary layer mixing. + BBLD_io = 0.0 + Kd_BBL(:) = 0.0 + mixvel_BBL(:) = 0.0 ; mixlen_BBL(:) = 0.0 + eCD%BBL_its = 0 + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = 0.0 ; eCD%dTKE_BBL_decay = 0.0 ; eCD%dTKE_BBL = 0.0 + ! eCD%dTKE_BBL_MKE = 0.0 + endif + return + else + ! There will be added bottom boundary layer mixing. + + h_neglect = GV%H_subroundoff + dz_neglect = GV%dZ_subroundoff + + C1_3 = 1.0 / 3.0 + I_dtdiag = 1.0 / dt + max_itt = 20 + dz_tt_min = 0.0 + + ! The next two blocks of code could be shared with ePBL_column. + + ! Set up fields relating a layer's temperature and salinity changes to potential energy changes. + pres_Z(1) = 0.0 + do k=1,nz + dMass = GV%H_to_RZ * h(k) + dPres = US%L_to_Z**2 * GV%g_Earth * dMass + dT_to_dPE(k) = (dMass * (pres_Z(K) + 0.5*dPres)) * dSV_dT(k) + dS_to_dPE(k) = (dMass * (pres_Z(K) + 0.5*dPres)) * dSV_dS(k) + dT_to_dColHt(k) = dMass * dSV_dT(k) + dS_to_dColHt(k) = dMass * dSV_dS(k) + + pres_Z(K+1) = pres_Z(K) + dPres + enddo + + if (GV%Boussinesq) then + do K=1,nz+1 ; h_dz_int(K) = GV%Z_to_H ; enddo + else + h_dz_int(1) = (h(1) + h_neglect) / (dz(1) + dz_neglect) + do K=2,nz + h_dz_int(K) = (h(k-1) + h(k) + h_neglect) / (dz(k-1) + dz(k) + dz_neglect) + enddo + h_dz_int(nz+1) = (h(nz) + h_neglect) / (dz(nz) + dz_neglect) + endif + ! The two previous blocks of code could be shared with ePBL_column. + + ! Determine the total thickness (dz_sum) and the fractional distance from the top (dztop_dztot). + dz_sum = 0.0 ; do k=1,nz ; dz_sum = dz_sum + dz(k) ; enddo + I_dzsum = 0.0 ; if (dz_sum > 0.0) I_dzsum = 1.0 / dz_sum + dz_top = 0.0 + dztop_dztot(nz+1) = 0.0 + do k=1,nz + dz_top = dz_top + dz(k) + dztop_dztot(K) = dz_top * I_dzsum + enddo + + ! Set terms from a tridiagonal solver based on the previously determined diffusivities. + Kddt_h(1) = 0.0 + hp_a(1) = h(1) + dT_to_dPE_a(1) = dT_to_dPE(1) ; dT_to_dColHt_a(1) = dT_to_dColHt(1) + dS_to_dPE_a(1) = dS_to_dPE(1) ; dS_to_dColHt_a(1) = dS_to_dColHt(1) + do K=2,nz + dt_h = dt / max(0.5*(dz(k-1)+dz(k)), 1e-15*dz_sum) + Kddt_h(K) = Kd(K) * dt_h + b1 = 1.0 / (hp_a(k-1) + Kddt_h(K)) + c1(K) = Kddt_h(K) * b1 + hp_a(k) = h(k) + (hp_a(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) + dT_to_dColHt_a(k) = dT_to_dColHt(k) + c1(K)*dT_to_dColHt_a(k-1) + dS_to_dColHt_a(k) = dS_to_dColHt(k) + c1(K)*dS_to_dColHt_a(k-1) + if (K<=2) then + Te(k-1) = b1*(h(k-1)*T0(k-1)) ; Se(k-1) = b1*(h(k-1)*S0(k-1)) + Th_a(k-1) = h(k-1) * T0(k-1) ; Sh_a(k-1) = h(k-1) * S0(k-1) + else + Te(k-1) = b1 * (h(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2)) + Se(k-1) = b1 * (h(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2)) + Th_a(k-1) = h(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2) + Sh_a(k-1) = h(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2) + endif + enddo + Kddt_h(nz+1) = 0.0 + if (debug) then + ! Complete the tridiagonal solve for Te and Se, which may be useful for debugging. + b1 = 1.0 / hp_a(nz) + Te(nz) = b1 * (h(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) + Se(nz) = b1 * (h(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) + do k=nz-1,1,-1 + Te(k) = Te(k) + c1(K+1)*Te(k+1) + Se(k) = Se(k) + c1(K+1)*Se(k+1) + enddo + endif + + BBLD_guess = BBLD_io + + !/The following lines are for the iteration over BBLD + ! max_BBLD will initialized as ocean bottom depth + max_BBLD = 0.0 ; do k=1,nz ; max_BBLD = max_BBLD + dz(k) ; enddo + ! min_BBLD will be initialized to 0. + min_BBLD = 0.0 + ! Set values of the wrong signs to indicate that these changes are not based on valid estimates + dBBLD_min = -1.0*US%m_to_Z ; dBBLD_max = 1.0*US%m_to_Z + + ! If no first guess is provided for BBLD, try the middle of the water column + if (BBLD_guess <= min_BBLD) BBLD_guess = 0.5 * (min_BBLD + max_BBLD) + + ! Iterate to determine a converged EPBL bottom boundary layer depth. + do BBL_it=1,CS%max_BBLD_its + + if (debug) then ; mech_BBL_TKE_k(:) = 0.0 ; endif + + ! Reset BBL_depth + BBLD_output = dz(nz) + bot_connected = .true. + + mech_BBL_TKE = BBL_TKE_in + + if (CS%TKE_diagnostics) then + ! eCD%dTKE_BBL_MKE = 0.0 + eCD%dTKE_BBL_mixing = 0.0 + eCD%dTKE_BBL_decay = 0.0 + eCD%dTKE_BBL = mech_BBL_TKE * I_dtdiag + endif + + ! Store in 1D arrays for output. + do K=1,nz+1 ; mixvel_BBL(K) = 0.0 ; mixlen_BBL(K) = 0.0 ; enddo + + ! Determine the mixing shape function MixLen_shape. + if ((.not.CS%Use_BBLD_iteration) .or. & + (CS%transLay_scale >= 1.0) .or. (CS%transLay_scale < 0.0) ) then + do K=1,nz+1 + MixLen_shape(K) = 1.0 + enddo + elseif (BBLD_guess <= 0.0) then + if (CS%transLay_scale > 0.0) then ; do K=1,nz+1 + MixLen_shape(K) = CS%transLay_scale + enddo ; else ; do K=1,nz+1 + MixLen_shape(K) = 1.0 + enddo ; endif + else + ! Reduce the mixing length based on BBLD, with a quadratic + ! expression that follows KPP. + I_BBLD = 1.0 / BBLD_guess + dz_rsum = 0.0 + MixLen_shape(nz+1) = 1.0 + if (CS%MixLenExponent_BBL==2.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) )**2 + enddo + elseif (CS%MixLenExponent_BBL==1.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) ) + enddo + else ! (CS%MixLenExponent_BBL /= 2.0 or 1.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) )**CS%MixLenExponent_BBL + enddo + endif + endif + + Kd_BBL(nz+1) = 0.0 ; Kddt_h(nz+1) = 0.0 + hp_b(nz) = h(nz) + dT_to_dPE_b(nz) = dT_to_dPE(nz) ; dT_to_dColHt_b(nz) = dT_to_dColHt(nz) + dS_to_dPE_b(nz) = dS_to_dPE(nz) ; dS_to_dColHt_b(nz) = dS_to_dColHt(nz) + + htot = h(nz) ; dztot = dz(nz) ; uhtot = u(nz)*h(nz) ; vhtot = v(nz)*h(nz) + + if (debug) then + mech_BBL_TKE_k(nz) = mech_BBL_TKE + num_itts(:) = -1 + endif + + Idecay_len_TKE = (CS%TKE_decay_BBL * absf) / u_star_BBL + do K=nz,2,-1 + ! Apply dissipation to the TKE, here applied as an exponential decay + ! due to 3-d turbulent energy being lost to inefficient rotational modes. + ! The following form is often used for mechanical stirring from the surface. + ! There could be several different "flavors" of TKE that decay at different rates. + exp_kh = 1.0 + if (Idecay_len_TKE > 0.0) exp_kh = exp(-h(k)*Idecay_len_TKE) + if (CS%TKE_diagnostics) & + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay + (exp_kh-1.0) * mech_BBL_TKE * I_dtdiag + mech_BBL_TKE = mech_BBL_TKE * exp_kh + + if (debug) then + mech_BBL_TKE_k(K) = mech_BBL_TKE + endif + + ! Precalculate some temporary expressions that are independent of Kddt_h(K). + dt_h = dt / max(0.5*(dz(k-1)+dz(k)), 1e-15*dz_sum) + + ! This tests whether the layers above and below this interface are in + ! a convectively stable configuration, without considering any effects of + ! mixing at higher interfaces. It is an approximation to the more + ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of + ! mixing across interface K+1. The dT_to_dColHt here are effectively + ! mass-weighted estimates of dSV_dT. + Convectively_stable = ( 0.0 <= & + ( (dT_to_dColHt(k) + dT_to_dColHt(k-1) ) * (T0(k-1)-T0(k)) + & + (dS_to_dColHt(k) + dS_to_dColHt(k-1) ) * (S0(k-1)-S0(k)) ) ) + + if ((mech_BBL_TKE <= 0.0) .and. Convectively_stable) then + ! Energy is already exhausted, so set Kd_BBL = 0 and cycle or exit? + mech_BBL_TKE = 0.0 + Kd_BBL(K) = 0.0 ; Kddt_h(K) = Kd(K) * dt_h + bot_disconnect = .true. + ! if (.not.debug) exit + + else ! mech_BBL_TKE > 0.0 or this is a potentially convectively unstable profile. + bot_disconnect = .false. + + ! Precalculate some more temporary expressions that are independent of Kddt_h(K). + if (K>=nz) then + Th_b(k) = h(k) * T0(k) ; Sh_b(k) = h(k) * S0(k) + else + Th_b(k) = h(k) * T0(k) + Kddt_h(K+1) * Te(k+1) + Sh_b(k) = h(k) * S0(k) + Kddt_h(K+1) * Se(k+1) + endif + + ! Using Pr=1 and the diffusivity at the upper interface (once it is + ! known), determine how much resolved mean kinetic energy (MKE) will be + ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of + ! this to the mTKE budget available for mixing in the next layer. + ! This is not enabled yet for BBL mixing. + ! if ((CS%MKE_to_TKE_effic > 0.0) .and. (htot*h(k-1) > 0.0)) then + ! ! This is the energy that would be available from homogenizing the + ! ! velocities between layer k-1 and the layers below. + ! dMKE_max = (US%L_to_Z**2*GV%H_to_RZ * CS%MKE_to_TKE_effic) * 0.5 * & + ! (h(k-1) / ((htot + h(k-1))*htot)) * & + ! ((uhtot-u(k-1)*htot)**2 + (vhtot-v(k-1)*htot)**2) + ! ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be + ! ! extracted by mixing with a finite viscosity. + ! MKE2_Hharm = (htot + h(k-1) + 2.0*h_neglect) / & + ! ((htot+h_neglect) * (h(k-1)+h_neglect)) + ! else + ! dMKE_max = 0.0 + ! MKE2_Hharm = 0.0 + ! endif + + ! At this point, Kddt_h(K) will be unknown because its value may depend + ! on how much energy is available. + dz_tt = dztot + dz_tt_min + if (mech_BBL_TKE > 0.0) then + if (CS%wT_scheme_BBL==wT_from_cRoot_TKE) then + vstar = CS%vstar_scale_fac_BBL * cuberoot(SpV_dt(K)*mech_BBL_TKE) + elseif (CS%wT_scheme_BBL==wT_from_RH18) then + Surface_Scale = max(0.05, 1.0 - dztot / BBLD_guess) + vstar = (CS%vstar_scale_fac_BBL * Surface_Scale) * ( CS%vstar_surf_fac_BBL*u_star_BBL/h_dz_int(K) ) + endif + hbs_here = min(dztop_dztot(K), MixLen_shape(K)) + mixlen_BBL(K) = max(CS%min_BBL_mix_len, ((dz_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef_BBL * absf) * (dz_tt*hbs_here) + vstar)) + Kd_guess0 = (h_dz_int(K)*vstar) * CS%vonKar * mixlen_BBL(K) + else + vstar = 0.0 ; Kd_guess0 = 0.0 + endif + mixvel_BBL(K) = vstar ! Track vstar + + TKE_rescale = 1.0 + if (CS%decay_adjusted_BBL_TKE) then + ! Add a scaling factor that accounts for the exponential decay of TKE from a + ! near-bottom source and the assumption that an increase in the diffusivity at an + ! interface causes a linearly increasing buoyancy flux going from 0 at the bottom + ! to a peak at the interface, and then going back to 0 atop the layer above. + TKE_rescale = exp_decay_TKE_adjust(htot, h(k-1), Idecay_len_TKE) + endif + + TKE_eff_avail = TKE_rescale*mech_BBL_TKE + + if (no_MKE_conversion) then + ! Without conversion from MKE to TKE, the updated diffusivity can be determined directly. + call find_Kd_from_PE_chg(Kd(K), Kd_guess0, dt_h, TKE_eff_avail, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), Kd_add=Kd_BBL(K), PE_chg=TKE_eff_used, & + frac_dKd_max_PE=frac_in_BL) + + ! Do not add energy if the column is convectively unstable. This was handled previously + ! for mixing from the surface. + if (TKE_eff_used < 0.0) TKE_eff_used = 0.0 + + ! Convert back to the TKE that has actually been used. + if (CS%decay_adjusted_BBL_TKE) then + if (TKE_rescale == 0.0) then ! This probably never occurs, even at roundoff. + TKE_used = mech_BBL_TKE ! All the energy was dissipated before it could mix. + else + TKE_used = TKE_eff_used / TKE_rescale + endif + else + TKE_used = TKE_eff_used + endif + + if (bot_connected) BBLD_output = BBLD_output + frac_in_BL*dz(k-1) + if (frac_in_BL < 1.0) bot_disconnect = .true. + + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - TKE_eff_used * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (TKE_used-TKE_eff_used) * I_dtdiag + endif + + mech_BBL_TKE = mech_BBL_TKE - TKE_used + + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + + else + Kddt_h_prev = Kd(K) * dt_h + Kddt_h_g0 = Kd_guess0 * dt_h + ! Find the change in PE with the guess at the added bottom boundary layer mixing. + call find_PE_chg(Kddt_h_prev, Kddt_h_g0, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), & + PE_chg=PE_chg_g0, dPEc_dKd_0=dPEc_dKd_Kd0 ) + + ! MKE_src = 0.0 ! Enable later?: = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) + + ! Do not add energy if the column is convectively unstable. This was handled previously + ! for mixing from the surface. + if (PE_chg_g0 < 0.0) PE_chg_g0 = 0.0 + + ! This block checks out different cases to determine Kd at the present interface. + ! if (mech_BBL_TKE*TKE_rescale + (MKE_src - PE_chg_g0) >= 0.0) then + if (TKE_eff_avail - PE_chg_g0 >= 0.0) then + ! This column is convectively stable and there is energy to support the suggested + ! mixing, or it is convectively unstable. Keep this first estimate of Kd. + Kd_BBL(K) = Kd_guess0 + Kddt_h(K) = Kddt_h_prev + Kddt_h_g0 + + TKE_used = PE_chg_g0 / TKE_rescale + + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - PE_chg_g0 * I_dtdiag +! eCD%dTKE_MKE = eCD%dTKE_MKE + MKE_src * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (TKE_used - PE_chg_g0) * I_dtdiag + endif + + ! mech_BBL_TKE = mech_BBL_TKE + MKE_src - TKE_used + mech_BBL_TKE = mech_BBL_TKE - TKE_used + if (bot_connected) then + BBLD_output = BBLD_output + dz(k-1) + endif + + elseif (TKE_eff_avail == 0.0) then + ! This can arise if there is no energy input to drive mixing or if there + ! is such strong decay that the mech_BBL_TKE becomes 0 via an underflow. + Kd_BBL(K) = 0.0 ; Kddt_h(K) = Kddt_h_prev + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - mech_BBL_TKE * I_dtdiag + endif + mech_BBL_TKE = 0.0 + bot_disconnect = .true. + else + ! There is not enough energy to support the mixing, so reduce the + ! diffusivity to what can be supported. + Kddt_h_max = Kddt_h_g0 ; Kddt_h_min = 0.0 + ! TKE_left_max = TKE_eff_avail + (MKE_src - PE_chg_g0) + TKE_left_max = TKE_eff_avail - PE_chg_g0 + TKE_left_min = TKE_eff_avail + + ! As a starting guess, take the minimum of a false position estimate + ! and a Newton's method estimate starting from dKddt_h = 0.0 + ! Enable conversion from MKE to TKE in the bottom boundary layer later? + ! Kddt_h_guess = TKE_eff_avail * Kddt_h_max / max( PE_chg_g0 - MKE_src, & + ! Kddt_h_max * (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + Kddt_h_guess = TKE_eff_avail * Kddt_h_max / max( PE_chg_g0, Kddt_h_max * dPEc_dKd_Kd0 ) + ! The above expression is mathematically the same as the following + ! except it is not susceptible to division by zero when + ! dPEc_dKd_Kd0 = dMKE_max = 0 . + ! Kddt_h_guess = TKE_eff_avail * min( Kddt_h_max / (PE_chg_g0 - MKE_src), & + ! 1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + if (debug) then + TKE_left_itt(:) = 0.0 ; dPEa_dKd_itt(:) = 0.0 ; PE_chg_itt(:) = 0.0 + Kddt_h_itt(:) = 0.0 ! ; MKE_src_itt(:) = 0.0 + endif + do itt=1,max_itt + call find_PE_chg(Kddt_h_prev, Kddt_h_guess, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), & + PE_chg=PE_chg, dPEc_dKd=dPEc_dKd) + ! Enable conversion from MKE to TKE in the bottom boundary layer later? + ! MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) + ! dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) + + ! TKE_left = TKE_eff_avail + (MKE_src - PE_chg) + TKE_left = TKE_eff_avail - PE_chg + if (debug .and. itt<=20) then + Kddt_h_itt(itt) = Kddt_h_guess ! ; MKE_src_itt(itt) = MKE_src + PE_chg_itt(itt) = PE_chg ; dPEa_dKd_itt(itt) = dPEc_dKd + TKE_left_itt(itt) = TKE_left + endif + ! Store the new bounding values, bearing in mind that min and max + ! here refer to Kddt_h and dTKE_left/dKddt_h < 0: + if (TKE_left >= 0.0) then + Kddt_h_min = Kddt_h_guess ; TKE_left_min = TKE_left + else + Kddt_h_max = Kddt_h_guess ; TKE_left_max = TKE_left + endif + + ! Try to use Newton's method, but if it would go outside the bracketed + ! values use the false-position method instead. + use_Newt = .true. + ! if (dPEc_dKd - dMKE_src_dK <= 0.0) then + if (dPEc_dKd <= 0.0) then + use_Newt = .false. + else + ! dKddt_h_Newt = TKE_left / (dPEc_dKd - dMKE_src_dK) + dKddt_h_Newt = TKE_left / dPEc_dKd + Kddt_h_Newt = Kddt_h_guess + dKddt_h_Newt + if ((Kddt_h_Newt > Kddt_h_max) .or. (Kddt_h_Newt < Kddt_h_min)) & + use_Newt = .false. + endif + + if (use_Newt) then + Kddt_h_next = Kddt_h_guess + dKddt_h_Newt + dKddt_h = dKddt_h_Newt + else + Kddt_h_next = (TKE_left_max * Kddt_h_min - Kddt_h_max * TKE_left_min) / & + (TKE_left_max - TKE_left_min) + dKddt_h = Kddt_h_next - Kddt_h_guess + endif + + if ((abs(dKddt_h) < 1e-9*Kddt_h_guess) .or. (itt==max_itt)) then + ! Use the old value so that the energy calculation does not need to be repeated. + if (debug) num_itts(K) = itt + exit + else + Kddt_h_guess = Kddt_h_next + endif + enddo ! Inner iteration loop on itt. + Kd_BBL(K) = Kddt_h_guess / dt_h + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + + ! All TKE should have been consumed. + if (CS%TKE_diagnostics) then + ! eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - (TKE_eff_avail + MKE_src) * I_dtdiag + ! eCD%dTKE_BBL_MKE = eCD%dTKE_BBL_MKE + MKE_src * I_dtdiag + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - TKE_eff_avail * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (mech_BBL_TKE-TKE_eff_avail) * I_dtdiag + endif + + if (bot_connected) BBLD_output = BBLD_output + (PE_chg / PE_chg_g0) * dz(k-1) + + mech_BBL_TKE = 0.0 + bot_disconnect = .true. + endif ! End of convective or forced mixing cases to determine Kd. + endif + + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + endif ! tot_TKT > 0.0 branch. Kddt_h(K) has been set. + + ! 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 / (hp_b(k) + Kddt_h(K)) + c1(K) = Kddt_h(K) * b1 + + hp_b(k-1) = h(k-1) + (hp_b(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) + dT_to_dColHt_b(k-1) = dT_to_dColHt(k-1) + c1(K)*dT_to_dColHt_b(k) + dS_to_dColHt_b(k-1) = dS_to_dColHt(k-1) + c1(K)*dS_to_dColHt_b(k) + + ! Store integrated velocities and thicknesses for MKE conversion calculations. + if (bot_disconnect) then + ! There is no turbulence at this interface, so restart the running sums. + uhtot = u(k-1)*h(k-1) + vhtot = v(k-1)*h(k-1) + htot = h(k-1) + dztot = dz(k-1) + bot_connected = .false. + else + uhtot = uhtot + u(k-1)*h(k-1) + vhtot = vhtot + v(k-1)*h(k-1) + htot = htot + h(k-1) + dztot = dztot + dz(k-1) + endif + + if (K==nz) then + Te(k) = b1*(h(k)*T0(k)) + Se(k) = b1*(h(k)*S0(k)) + else + Te(k) = b1 * (h(k) * T0(k) + Kddt_h(K+1) * Te(k+1)) + Se(k) = b1 * (h(k) * S0(k) + Kddt_h(K+1) * Se(k+1)) + endif + enddo + Kd_BBL(1) = 0.0 + + if (debug) then + ! Complete the tridiagonal solve for Te with a downward pass. + b1 = 1.0 / hp_b(1) + Te(1) = b1 * (h(1) * T0(1) + Kddt_h(2) * Te(2)) + Se(1) = b1 * (h(1) * S0(1) + Kddt_h(2) * Se(2)) + dT_expect(1) = Te(1) - T0(1) ; dS_expect(1) = Se(1) - S0(1) + do k=2,nz + Te(k) = Te(k) + c1(K)*Te(k-1) + Se(k) = Se(k) + c1(K)*Se(k-1) + dT_expect(k) = Te(k) - T0(k) ; dS_expect(k) = Se(k) - S0(k) + enddo + + dPE_debug = 0.0 + do k=1,nz + dPE_debug = dPE_debug + (dT_to_dPE(k) * (Te(k) - T0(k)) + & + dS_to_dPE(k) * (Se(k) - S0(k))) + enddo + mixing_debug = dPE_debug * I_dtdiag + endif + + ! Skip the rest of the contents of the do loop if there are no more BBL depth iterations. + if (BBL_it >= CS%max_BBLD_its) exit + + ! The following lines are used for the iteration to determine the boundary layer depth. + ! Note that the iteration uses the value predicted by the TKE threshold (BBL_DEPTH), + ! because the mixing length shape is dependent on the BBL depth, and therefore the BBL depth + ! should be estimated more precisely than the grid spacing. + + ! New method uses BBL_DEPTH as computed in ePBL routine + BBLD_found = BBLD_output + if (abs(BBLD_found - BBLD_guess) < CS%BBLD_tol) then + exit ! Break the BBL depth convergence loop + elseif (BBLD_found > BBLD_guess) then + min_BBLD = BBLD_guess ; dBBLD_min = BBLD_found - BBLD_guess + else ! We know this guess was too deep + max_BBLD = BBLD_guess ; dBBLD_max = BBLD_found - BBLD_guess ! <= -CS%BBLD_tol + endif + + ! Try using the false position method or the returned value instead of simple bisection. + ! Taking the occasional step with BBLD_output empirically helps to converge faster. + if ((dBBLD_min > 0.0) .and. (dBBLD_max < 0.0) .and. (BBL_it > 2) .and. (mod(BBL_it-1,4) > 0)) then + ! Both bounds have valid change estimates and are probably in the range of possible outputs. + BBLD_guess = (dBBLD_min*max_BBLD - dBBLD_max*min_BBLD) / (dBBLD_min - dBBLD_max) + elseif ((BBLD_found > min_BBLD) .and. (BBLD_found < max_BBLD)) then + ! The output BBLD_found is an interesting guess, as it is likely to bracket the true solution + ! along with the previous value of BBLD_guess and to be close to the solution. + BBLD_guess = BBLD_found + else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. + BBLD_guess = 0.5*(min_BBLD + max_BBLD) + endif + + enddo ! Iteration loop for converged boundary layer thickness. + + eCD%BBL_its = min(BBL_it, CS%max_BBLD_its) + + BBLD_io = BBLD_output + endif + +end subroutine ePBL_BBL_column + +!> Determine a scaling factor that accounts for the exponential decay of turbulent kinetic energy +!! from a boundary source and the assumption that an increase in the diffusivity at an interface +!! causes a linearly increasing buoyancy flux going from 0 at the bottom to a peak at the interface, +!! and then going back to 0 atop the layer above. Where this factor increases the available mixing +!! TKE, it is only compensating for the fact that the TKE has already been reduced by the same +!! exponential decay rate. ha and hb must be non-negative, and this function generally increases +!! with hb and decreases with ha. +!! +!! Exp_decay_TKE_adjust is coded to have a lower bound of 1e-30 on the return value. For large +!! values of ha*Idecay, the return value is about 0.5*ka*(ha+hb)*Idecay**2 * exp(-ha*Idecay), but +!! return values of less than 1e-30 are deliberately reset to 1e-30. For relatively large values +!! of hb*Idecay, the return value increases linearly with hb. When Idecay ~= 0, the return value +!! is close to 1. +function exp_decay_TKE_adjust(hb, ha, Idecay) result(TKE_to_PE_scale) + real, intent(in) :: hb !< The thickness over which the buoyancy flux varies on the + !! near-boundary side of an interface (e.g., a well-mixed bottom + !! boundary layer thickness) [H ~> m or kg m-2] + real, intent(in) :: ha !< The thickness of the layer on the opposite side of an interface from + !! the boundary [H ~> m or kg m-2] + real, intent(in) :: Idecay !< The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1] + real :: TKE_to_PE_scale !< The effective fractional change in energy available to + !! drive mixing at this interface once the exponential decay of TKE + !! is accounted for [nondim]. TKE_to_PE_scale is always positive. + + real :: khb ! The thickness on the boundary side times the TKE decay rate [nondim] + real :: kha ! The thickness away from from the boundary times the TKE decay rate [nondim] + real, parameter :: C1_3 = 1.0/3.0 ! A rational constant [nondim] + + khb = abs(hb*Idecay) + kha = abs(ha*Idecay) + + ! For large enough kha that exp(kha) > 1.0e17*kha: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * kha) * exp(-kha) > (0.5 * kha**2) * exp(-kha) + ! To keep TKE_to_PE_scale > -1e30 and avoid overflow in the exp(), keep kha < kha_max_30, where: + ! kha_max_30 = ln(0.5*1e30) + 2.0 * ln(kha_max_30) ~= 68.3844 + 2.0 * ln(68.3844+8.6895)) + ! If kha_max = 77.0739, (0.5 * kha_max**2) * exp(-kha_max) = 1.0e-30. + + if (kha > 77.0739) then + TKE_to_PE_scale = 1.0e-30 + elseif ((kha > 2.2e-4) .and. (khb > 2.2e-4)) then + ! This is the usual case, derived from integrals of z exp(z) over the layers above and below. + ! TKE_to_PE_scale = (0.5 * (khb + kha)) / & + ! ((exp(-khb) - (1.0 - khb)) / khb + (exp(kha) - (1.0 + kha)) / kha) + TKE_to_PE_scale = (0.5 * (khb + kha) * (kha * khb)) / & + (kha * (exp(-khb) - (1.0 - khb)) + khb * (exp(kha) - (1.0 + kha))) + elseif (khb > 2.2e-4) then + ! For small values of kha, approximate (exp(kha) - (1.0 + hha)) by the first two + ! terms of its Taylor series: 0.5*kha**2 + C1_6*kha**3 + ... + kha**n/n! + ... + ! which is more accurate when kha**4/24. < 1e-16 or kha < ~ 2.21e-4. + TKE_to_PE_scale = (0.5 * (khb + kha) * khb) / & + ((exp(-khb) - (1.0 - khb)) + 0.5*(khb * kha) * (1.0 + C1_3*kha)) + elseif (kha > 2.2e-4) then + ! Use a Taylor series expansion for small values of khb + TKE_to_PE_scale = (0.5 * (khb + kha) * kha) / & + (0.5 * (kha * khb) * (1.0 - C1_3*Khb) + (exp(kha) - (1.0 + kha))) + else ! (kha < 2.2e-4) .and. (khb < 2.2e-4) - use Taylor series approximations for both + TKE_to_PE_scale = 1.0 / (1.0 + C1_3*(kha - khb)) + endif + + if (TKE_to_PE_scale < 1.0e-30) TKE_to_PE_scale = 1.0e-30 + + ! For kha >> 1: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * kha) * exp(-kha) + + ! For khb >> 1: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * (kha * khb)) / & + ! (khb * exp(kha) - (kha + khb))) + ! For khb >> 1 and khb >> kha: + ! TKE_to_PE_scale = (0.5 * (kha * khb)) / (exp(kha) - 1.0)) + +end function exp_decay_TKE_adjust + +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interface'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_Z, 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, PE_ColHt_cor) + real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface [H ~> m or kg m-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 [H ~> m or kg m-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 [H ~> m or kg m-2]. + 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 below [H ~> m or kg m-2]. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [S H ~> ppt m or ppt kg m-2]. + 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. + 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 [Z C-1 ~> m degC-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 [Z S-1 ~> 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 [Z C-1 ~> m degC-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 [Z S-1 ~> m ppt-1]. + + real, intent(out) :: PE_chg !< The change in column potential energy from applying + !! dKddt_h at the present interface [R Z3 T-2 ~> J m-2]. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with dKddt_h + !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realized by applying a huge value of dKddt_h at the + !! present interface [R Z3 T-2 ~> J m-2]. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with dKddt_h in the + !! limit where dKddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net + !! change in the column height [R Z3 T-2 ~> J m-2]. + + ! Local variables + real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. + real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. + real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes [R Z2 T-2 ~> J m-3]. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes [H Z ~> m2 or kg m-1]. + real :: ColHt_chg ! The change in the column height [H ~> m or kg m-2]. + real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3]. + real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4]. + + ! 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 + ! succinct form. The derivation is not necessarily obvious, but it demonstrably + ! works by comparison with separate calculations of the energy changes after + ! the tridiagonal solver for the final changes in temperature and salinity are + ! applied. + + hps = hp_a + hp_b + bdt1 = hp_a * hp_b + Kddt_h0 * hps + dT_c = hp_a * Th_b - hp_b * Th_a + dS_c = hp_a * Sh_b - hp_b * Sh_a + PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & + hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) + ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & + hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) + + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKddt_h. + y1_3 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + PE_chg = PEc_core * y1_3 + ColHt_chg = ColHt_core * y1_3 + if (ColHt_chg < 0.0) PE_chg = PE_chg - pres_Z * ColHt_chg + + if (present(PE_ColHt_cor)) PE_ColHt_cor = -pres_Z * min(ColHt_chg, 0.0) + + if (present(dPEc_dKd)) then + ! Find the derivative of the potential energy change with dKddt_h. + y1_4 = 1.0 / (bdt1 + dKddt_h * hps)**2 + dPEc_dKd = PEc_core * y1_4 + ColHt_chg = ColHt_core * y1_4 + if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres_Z * ColHt_chg + endif + + if (present(dPE_max)) then + ! This expression is the limit of PE_chg for infinite dKddt_h. + y1_3 = 1.0 / (bdt1 * hps) + dPE_max = PEc_core * y1_3 + ColHt_chg = ColHt_core * y1_3 + if (ColHt_chg < 0.0) dPE_max = dPE_max - pres_Z * ColHt_chg + endif + + if (present(dPEc_dKd_0)) then + ! This expression is the limit of dPEc_dKd for dKddt_h = 0. + y1_4 = 1.0 / bdt1**2 + dPEc_dKd_0 = PEc_core * y1_4 + ColHt_chg = ColHt_core * y1_4 + if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres_Z * ColHt_chg + endif + +end subroutine find_PE_chg + + +!> This subroutine directly calculates the an increment in the diapycnal diffusivity based on the +!! change in potential energy within a timestep, subject to bounds on the possible change in +!! diffusivity, returning both the added diffusivity and the realized potential energy change, and +!! optionally also the maximum change in potential energy that would be realized for an infinitely +!! large diffusivity. +subroutine find_Kd_from_PE_chg(Kd_prev, dKd_max, dt_h, max_PE_chg, 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_Z, & + dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & + Kd_add, PE_chg, dPE_max, frac_dKd_max_PE) + real, intent(in) :: Kd_prev !< The previously used diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: dKd_max !< The maximum change in the diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: dt_h !< The time step and divided by the average of the + !! thicknesses around the interface [T Z-1 ~> s m-1]. + real, intent(in) :: max_PE_chg !< The maximum change in the column potential energy due to + !! additional mixing at an interface [R Z3 T-2 ~> J m-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 [H ~> m or kg m-2]. + 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 below [H ~> m or kg m-2]. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [S H ~> ppt m or ppt kg m-2]. + 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> 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 [R Z3 T-2 C-1 ~> J m-2 degC-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 [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. + 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 [Z C-1 ~> m degC-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 [Z S-1 ~> 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 [Z C-1 ~> m degC-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 [Z S-1 ~> m ppt-1]. + real, intent(out) :: Kd_add !< The additional diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(out) :: PE_chg !< The realized change in the column potential energy due to + !! additional mixing at an interface [R Z3 T-2 ~> J m-2]. + real, optional, & + intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realized by applying a huge value of dKddt_h at the + !! present interface [R Z3 T-2 ~> J m-2]. + real, optional, & + intent(out) :: frac_dKd_max_PE !< The fraction of the energy required to support dKd_max + !! that is supplied by max_PE_chg [nondim] + + ! Local variables + real :: Kddt_h0 ! The previously used diffusivity at an interface times the time step + ! and divided by the average of the thicknesses around the + ! interface [H ~> m or kg m-2]. + real :: dKddt_h ! The upper bound on the change in the diffusivity at an interface times + ! the time step and divided by the average of the thicknesses around + ! the interface [H ~> m or kg m-2]. + real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. + real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. + real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes [R Z2 T-2 ~> J m-3]. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes [H Z ~> m2 or kg m-1]. + + ! 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, which is then + ! extensively manipulated to get it into its most succinct form. It is the same as the + ! expression that appears in find_PE_chg. + + Kddt_h0 = Kd_prev * dt_h + 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 (ColHt_core < 0.0) PEc_core = PEc_core - pres_Z * ColHt_core + + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKd_max, and use this to dermine which limit applies. + dKddt_h = dKd_max * dt_h + if ( (PEc_core * dKddt_h <= max_PE_chg * (bdt1 * (bdt1 + dKddt_h * hps))) .or. (PEc_core <= 0.0) ) then + ! There is more than enough energy available to support the maximum permitted diffusivity. + Kd_add = dKd_max + PE_chg = PEc_core * dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + if (present(frac_dKd_max_PE)) frac_dKd_max_PE = 1.0 + else + ! Mixing is constrained by the available energy, so solve the following for Kd_add: + ! max_PE_chg = PEc_core * Kd_add * dt_h / (bdt1 * (bdt1 + Kd_add * dt_h * hps)) + ! It has been verified that the two branches are continuous. + Kd_add = (bdt1**2 * max_PE_chg) / (dt_h * (PEc_core - bdt1 * hps * max_PE_chg)) + PE_chg = max_PE_chg + if (present(frac_dKd_max_PE)) & + frac_dKd_max_PE = (PE_chg * (bdt1 * (bdt1 + dKddt_h * hps))) / (PEc_core * dKddt_h) + endif + + ! Note that the derivative of PE_chg with dKddt_h is monotonic: + ! dPE_chg_dKd = PEc_core * ( (bdt1 * (bdt1 + dKddt_h * hps)) - bdtl * hps * dKddt_h ) / & + ! (bdt1 * (bdt1 + dKddt_h * hps))**2 + ! dPE_chg_dKd = PEc_core / (bdt1 + dKddt_h * hps)**2 + + ! This expression is the limit of PE_chg for infinite dKddt_h. + if (present(dPE_max)) dPE_max = PEc_core / (bdt1 * hps) + +end subroutine find_Kd_from_PE_chg -end subroutine find_PE_chg !> This subroutine calculates the change in potential energy and or derivatives !! for several changes in an interface's diapycnal diffusivity times a timestep @@ -2074,12 +3382,15 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_energetic_PBL" ! This module's name. - character(len=20) :: tmpstr + character(len=20) :: tmpstr ! A string that is parsed for parameter settings + character(len=20) :: vel_scale_str ! A string that is parsed for velocity scale parameter settings + character(len=120) :: diff_text ! A clause describing parameter setting that differ. real :: omega_frac_dflt ! The default for omega_frac [nondim] integer :: isd, ied, jsd, jed integer :: mstar_mode, LT_enhance, wT_mode integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. - logical :: use_temperature, use_omega + logical :: use_omega + logical :: no_BBL ! If true, EPBL_BBL_EFFIC < 0 and bottom boundary layer mixing is not enabled. logical :: use_la_windsea isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2092,6 +3403,9 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) !/1. General ePBL settings + call get_param(param_file, mdl, "DEBUG", CS%debug, & + "If true, write out verbose debugging data.", & + default=.false., debuggingParam=.true.) call get_param(param_file, mdl, "OMEGA", CS%omega, & "The rotation rate of the earth.", & units="s-1", default=7.2921e-5, scale=US%T_to_S) @@ -2101,7 +3415,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "scale for turbulence.", default=.false., do_not_log=.true.) omega_frac_dflt = 0.0 if (use_omega) then - call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + call MOM_error(WARNING, "ML_USE_OMEGA is deprecated; use ML_OMEGA_FRAC=1.0 instead.") omega_frac_dflt = 1.0 endif call get_param(param_file, mdl, "ML_OMEGA_FRAC", CS%omega_frac, & @@ -2130,10 +3444,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) call get_param(param_file, mdl, "EPBL_ORIGINAL_PE_CALC", CS%orig_PE_calc, & - "If true, the ePBL code uses the original form of the "//& - "potential energy change code. Otherwise, the newer "//& - "version that can work with successive increments to the "//& - "diffusivity in upward or downward passes is used.", default=.true.) + "If true, the ePBL code uses the original form of the potential energy change "//& + "code. Otherwise, the newer version that can work with successive increments "//& + "to the diffusivity in upward or downward passes is used.", & + default=.true.) ! Change the default to .false.? call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", CS%MKE_to_TKE_effic, & "The efficiency with which mean kinetic energy released "//& @@ -2141,16 +3455,21 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "is converted to turbulent kinetic energy.", & units="nondim", default=0.0) call get_param(param_file, mdl, "TKE_DECAY", CS%TKE_decay, & - "TKE_DECAY relates the vertical rate of decay of the "//& - "TKE available for mechanical entrainment to the natural "//& - "Ekman depth.", units="nondim", default=2.5) + "TKE_DECAY relates the vertical rate of decay of the TKE available "//& + "for mechanical entrainment to the natural Ekman depth.", & + units="nondim", default=2.5) + call get_param(param_file, mdl, "DIRECT_EPBL_MIXING_CALC", CS%direct_calc, & + "If true and there is no conversion from mean kinetic energy to ePBL turbulent "//& + "kinetic energy, use a direct calculation of the diffusivity that is supported "//& + "by a given energy input instead of the more general but slower iterative solver.", & + default=.false., do_not_log=(CS%MKE_to_TKE_effic>0.0)) !/2. Options related to setting MSTAR call get_param(param_file, mdl, "EPBL_MSTAR_SCHEME", tmpstr, & "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//& "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//& - "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//& + "\t OM4 - Use L_Ekman/L_Obukhov in the stabilizing limit, as in OM4 \n"//& "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", & default=CONSTANT_STRING, do_not_log=.true.) call get_param(param_file, mdl, "MSTAR_MODE", mstar_mode, default=-1) @@ -2209,7 +3528,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.085, do_not_log=(CS%mstar_scheme/=MStar_from_Ekman)) ! mstar_scheme==MStar_from_RH18 options call get_param(param_file, mdl, "RH18_MSTAR_CN1", CS%RH18_mstar_cn1,& - "MSTAR_N coefficient 1 (outter-most coefficient for fit). "//& + "MSTAR_N coefficient 1 (outer-most coefficient for fit). "//& "The value of 0.275 is given in RH18. Increasing this "//& "coefficient increases MSTAR for all values of Hf/ust, but more "//& "effectively at low values (weakly developed OSBLs).", & @@ -2276,10 +3595,14 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "mixed layer depth. Otherwise use the false position after a maximum and minimum "//& "bound have been evaluated and the returned value or bisection before this.", & default=.false., do_not_log=.not.CS%Use_MLD_iteration) + call get_param(param_file, mdl, "EPBL_MLD_ITER_BUG", CS%MLD_iter_bug, & + "If true, use buggy logic that gives the wrong bounds for the next iteration "//& + "when successive guesses increase by exactly EPBL_MLD_TOLERANCE.", & + default=.true., do_not_log=.not.CS%Use_MLD_iteration) ! The default should be changed to .false. call get_param(param_file, mdl, "EPBL_MLD_MAX_ITS", CS%max_MLD_its, & "The maximum number of iterations that can be used to find a self-consistent "//& "mixed layer depth. If EPBL_MLD_BISECTION is true, the maximum number "//& - "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & + "of iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & default=20, do_not_log=.not.CS%Use_MLD_iteration) if (.not.CS%Use_MLD_iteration) CS%Max_MLD_Its = 1 call get_param(param_file, mdl, "EPBL_MIN_MIX_LEN", CS%min_mix_len, & @@ -2294,7 +3617,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=2.0) !/ Turbulent velocity scale in mixing coefficient - call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, & + call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", vel_scale_str, & "Selects the method for translating TKE into turbulent velocities. "//& "Valid values are: \n"//& "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//& @@ -2303,31 +3626,31 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) default=ROOT_TKE_STRING, do_not_log=.true.) call get_param(param_file, mdl, "EPBL_VEL_SCALE_MODE", wT_mode, default=-1) if (wT_mode == 0) then - tmpstr = ROOT_TKE_STRING + vel_scale_str = ROOT_TKE_STRING call MOM_error(WARNING, "Use EPBL_VEL_SCALE_SCHEME = CUBE_ROOT_TKE instead of the archaic EPBL_VEL_SCALE_MODE = 0.") elseif (wT_mode == 1) then - tmpstr = RH18_STRING + vel_scale_str = RH18_STRING call MOM_error(WARNING, "Use EPBL_VEL_SCALE_SCHEME = REICHL_H18 instead of the archaic EPBL_VEL_SCALE_MODE = 1.") elseif (wT_mode >= 2) then call MOM_error(FATAL, "An unrecognized value of the obsolete parameter EPBL_VEL_SCALE_MODE was specified.") endif - call log_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, & + call log_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", vel_scale_str, & "Selects the method for translating TKE into turbulent velocities. "//& "Valid values are: \n"//& "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//& "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//& "\t documented in Reichl & Hallberg, 2018.", & default=ROOT_TKE_STRING) - tmpstr = uppercase(tmpstr) - select case (tmpstr) + vel_scale_str = uppercase(vel_scale_str) + select case (vel_scale_str) case (ROOT_TKE_STRING) CS%wT_scheme = wT_from_cRoot_TKE case (RH18_STRING) CS%wT_scheme = wT_from_RH18 case default - call MOM_mesg('energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//'"', 0) + call MOM_mesg('energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(vel_scale_str)//'"', 0) call MOM_error(FATAL, "energetic_PBL_init: Unrecognized setting "// & - "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.") + "EPBL_VEL_SCALE_SCHEME = "//trim(vel_scale_str)//" found in input file.") end select call get_param(param_file, mdl, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, & @@ -2343,6 +3666,77 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "The proportionality times ustar to set vstar at the surface.", & units="nondim", default=1.2) + !/ Bottom boundary layer mixing related options + call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, & + "The efficiency of bottom boundary layer mixing via ePBL. Setting this to a "//& + "value that is greater than 0 to enable bottom boundary layer mixing from EPBL.", & + units="nondim", default=0.0) + no_BBL = (CS%ePBL_BBL_effic<=0.0) + call get_param(param_file, mdl, "USE_BBLD_ITERATION", CS%Use_BBLD_iteration, & + "A logical that specifies whether or not to use the distance to the top of the "//& + "actively turbulent bottom boundary layer to help set the EPBL length scale.", & + default=.true., do_not_log=no_BBL) + call get_param(param_file, mdl, "TKE_DECAY_BBL", CS%TKE_decay_BBL, & + "TKE_DECAY_BBL relates the vertical rate of decay of the TKE available for "//& + "mechanical entrainment in the bottom boundary layer to the natural Ekman depth.", & + units="nondim", default=CS%TKE_decay, do_not_log=no_BBL) + call get_param(param_file, mdl, "MIX_LEN_EXPONENT_BBL", CS%MixLenExponent_BBL, & + "The exponent applied to the ratio of the distance to the top of the BBL "//& + "and the total BBL depth which determines the shape of the mixing length. "//& + "This is only used if USE_MLD_ITERATION is True.", & + units="nondim", default=2.0, do_not_log=(no_BBL.or.(.not.CS%Use_BBLD_iteration))) + call get_param(param_file, mdl, "EPBL_MIN_BBL_MIX_LEN", CS%min_BBL_mix_len, & + "The minimum mixing length scale that will be used by ePBL for bottom boundary "//& + "layer mixing. Choosing (0) does not set a minimum.", & + units="meter", default=CS%min_mix_len, scale=US%m_to_Z, do_not_log=no_BBL) + call get_param(param_file, mdl, "EPBL_BBLD_TOLERANCE", CS%BBLD_tol, & + "The tolerance for the iteratively determined bottom boundary layer depth. "//& + "This is only used with USE_MLD_ITERATION.", & + units="meter", default=US%Z_to_m*CS%MLD_tol, scale=US%m_to_Z, & + do_not_log=(no_BBL.or.(.not.CS%Use_MLD_iteration))) + call get_param(param_file, mdl, "EPBL_BBLD_MAX_ITS", CS%max_BBLD_its, & + "The maximum number of iterations that can be used to find a self-consistent "//& + "bottom boundary layer depth.", & + default=CS%max_MLD_its, do_not_log=(no_BBL.or.(.not.CS%Use_MLD_iteration))) + if (.not.CS%Use_MLD_iteration) CS%max_BBLD_its = 1 + + call get_param(param_file, mdl, "EPBL_BBL_VEL_SCALE_SCHEME", tmpstr, & + "Selects the method for translating bottom boundary layer TKE into turbulent velocities. "//& + "Valid values are: \n"//& + "\t CUBE_ROOT_TKE - A constant times the cube root of remaining BBL TKE. \n"//& + "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//& + "\t documented in Reichl & Hallberg, 2018.", & + default=vel_scale_str, do_not_log=no_BBL) + select case (tmpstr) + case (ROOT_TKE_STRING) + CS%wT_scheme_BBL = wT_from_cRoot_TKE + case (RH18_STRING) + CS%wT_scheme_BBL = wT_from_RH18 + case default + call MOM_mesg('energetic_PBL_init: EPBL_BBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//'"', 0) + call MOM_error(FATAL, "energetic_PBL_init: Unrecognized setting "// & + "EPBL_BBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.") + end select + call get_param(param_file, mdl, "EPBL_BBL_VEL_SCALE_FACTOR", CS%vstar_scale_fac_BBL, & + "An overall nondimensional scaling factor for wT in the bottom boundary layer. "//& + "Making this larger increases the bottom boundary layer diffusivity.", & + units="nondim", default=CS%vstar_scale_fac, do_not_log=no_BBL) + call get_param(param_file, mdl, "VSTAR_BBL_SURF_FAC", CS%vstar_surf_fac_BBL,& + "The proportionality times ustar to set vstar in the bottom boundary layer.", & + units="nondim", default=CS%vstar_surf_fac, do_not_log=(no_BBL.or.(CS%wT_scheme_BBL/=wT_from_RH18))) + call get_param(param_file, mdl, "EKMAN_SCALE_COEF_BBL", CS%Ekman_scale_coef_BBL, & + "A nondimensional scaling factor controlling the inhibition of the diffusive "//& + "length scale by rotation in the bottom boundary layer. Making this larger "//& + "decreases the bottom boundary layer diffusivity.", & + units="nondim", default=CS%Ekman_scale_coef, do_not_log=no_BBL) + + call get_param(param_file, mdl, "DECAY_ADJUSTED_BBL_TKE", CS%decay_adjusted_BBL_TKE, & + "If true, include an adjustment factor in the bottom boundary layer energetics "//& + "that accounts for an exponential decay of TKE from a near-bottom source and "//& + "an assumed piecewise linear profile of the buoyancy flux response to a change "//& + "in a diffusivity.", & + default=.false., do_not_log=no_BBL) + !/ Options related to Langmuir turbulence call get_param(param_file, mdl, "USE_LA_LI2016", use_LA_Windsea, & "A logical to use the Li et al. 2016 (submitted) formula to "//& @@ -2360,7 +3754,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Valid values are: \n"//& "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//& "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//& - "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", & + "\t ADDITIVE - Add a Langmuir turbulence contribution to mstar to other contributions", & default=NONE_STRING, do_not_log=.true.) call get_param(param_file, mdl, "LT_ENHANCE", LT_enhance, default=-1) if (LT_ENHANCE == 0) then @@ -2384,7 +3778,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Valid values are: \n"//& "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//& "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//& - "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", & + "\t ADDITIVE - Add a Langmuir turbulence contribution to mstar to other contributions", & default=NONE_STRING) tmpstr = uppercase(tmpstr) select case (tmpstr) @@ -2404,7 +3798,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Coefficient for Langmuir enhancement of mstar", & units="nondim", default=0.447, do_not_log=(CS%LT_enhance_form==No_Langmuir)) call get_param(param_file, mdl, "LT_ENHANCE_EXP", CS%LT_ENHANCE_EXP, & - "Exponent for Langmuir enhancementt of mstar", & + "Exponent for Langmuir enhancement of mstar", & units="nondim", default=-1.33, do_not_log=(CS%LT_enhance_form==No_Langmuir)) call get_param(param_file, mdl, "LT_MOD_LAC1", CS%LaC_MLDoEK, & "Coefficient for modification of Langmuir number due to "//& @@ -2428,6 +3822,27 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.95, do_not_log=(CS%LT_enhance_form==No_Langmuir)) endif + !/ Options for documenting differences from parameter choices + call get_param(param_file, mdl, "EPBL_OPTIONS_DIFF", CS%options_diff, & + "If positive, this is a coded integer indicating a pair of settings whose "//& + "differences are diagnosed in a passive diagnostic mode via extra calls to "//& + "ePBL_column. If this is 0 or negative no extra calls occur.", & + default=0) + if (CS%options_diff > 0) then + if (CS%options_diff == 1) then + diff_text = "EPBL_ORIGINAL_PE_CALC settings" + elseif (CS%options_diff == 2) then + diff_text = "EPBL_ANSWER_DATE settings" + elseif (CS%options_diff == 3) then + diff_text = "DIRECT_EPBL_MIXING_CALC settings" + elseif (CS%options_diff == 4) then + diff_text = "BBL DIRECT_EPBL_MIXING_CALC settings" + elseif (CS%options_diff == 5) then + diff_text = "BBL DECAY_ADJUSTED_BBL_TKE settings" + else + diff_text = "unchanged settings" + endif + endif !/ Logging parameters ! This gives a minimum decay scale that is typically much less than Angstrom. @@ -2441,32 +3856,54 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) !/ Checking output flags CS%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, & - Time, 'Surface boundary layer depth', 'm', conversion=US%Z_to_m, & + Time, 'Surface boundary layer depth', units='m', conversion=US%Z_to_m, & cmor_long_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') ! This is an alias for the same variable as ePBL_h_ML CS%id_hML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, & - Time, 'Surface mixed layer depth based on active turbulence', 'm', conversion=US%Z_to_m) + Time, 'Surface mixed layer depth based on active turbulence', units='m', conversion=US%Z_to_m) CS%id_TKE_wind = register_diag_field('ocean_model', 'ePBL_TKE_wind', diag%axesT1, & - Time, 'Wind-stirring source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Wind-stirring source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_MKE = register_diag_field('ocean_model', 'ePBL_TKE_MKE', diag%axesT1, & - Time, 'Mean kinetic energy source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Mean kinetic energy source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_conv = register_diag_field('ocean_model', 'ePBL_TKE_conv', diag%axesT1, & - Time, 'Convective source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Convective source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_forcing = register_diag_field('ocean_model', 'ePBL_TKE_forcing', diag%axesT1, & Time, 'TKE consumed by mixing surface forcing or penetrative shortwave radation '//& - 'through model layers', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + 'through model layers', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_mixing = register_diag_field('ocean_model', 'ePBL_TKE_mixing', diag%axesT1, & - Time, 'TKE consumed by mixing that deepens the mixed layer', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'TKE consumed by mixing that deepens the mixed layer', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_mech_decay = register_diag_field('ocean_model', 'ePBL_TKE_mech_decay', diag%axesT1, & - Time, 'Mechanical energy decay sink of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Mechanical energy decay sink of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_conv_decay = register_diag_field('ocean_model', 'ePBL_TKE_conv_decay', diag%axesT1, & - Time, 'Convective energy decay sink of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Convective energy decay sink of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, & - Time, 'Mixing Length that is used', 'm', conversion=US%Z_to_m) + Time, 'Mixing Length that is used', units='m', conversion=US%Z_to_m) CS%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, & - Time, 'Velocity Scale that is used.', 'm s-1', conversion=US%Z_to_m*US%s_to_T) + Time, 'Velocity Scale that is used.', units='m s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, & Time, 'Total mstar that is used.', 'nondim') + if (CS%ePBL_BBL_effic > 0.0) then + CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_ePBL_BBL', diag%axesTi, & + Time, 'ePBL bottom boundary layer diffusivity', units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_BBL_Mix_Length = register_diag_field('ocean_model', 'BBL_Mixing_Length', diag%axesTi, & + Time, 'ePBL bottom boundary layer mixing length', units='m', conversion=US%Z_to_m) + CS%id_BBL_Vel_Scale = register_diag_field('ocean_model', 'BBL_Velocity_Scale', diag%axesTi, & + Time, 'ePBL bottom boundary layer velocity scale', units='m s-1', conversion=US%Z_to_m*US%s_to_T) + CS%id_BBL_depth = register_diag_field('ocean_model', 'h_BBL', diag%axesT1, & + Time, 'Bottom boundary layer depth based on active turbulence', units='m', conversion=US%Z_to_m) + CS%id_ustar_BBL = register_diag_field('ocean_model', 'ePBL_ustar_BBL', diag%axesT1, & + Time, 'The bottom boundary layer friction velocity', units='m s-1', conversion=GV%H_to_m*US%s_to_T) + CS%id_BBL_decay_scale = register_diag_field('ocean_model', 'BBL_decay_scale', diag%axesT1, & + Time, 'The bottom boundary layer TKE decay lengthscale', units='m', conversion=GV%H_to_m) + CS%id_TKE_BBL = register_diag_field('ocean_model', 'ePBL_BBL_TKE', diag%axesT1, & + Time, 'The source of TKE for the bottom boundary layer', units='W m-2', conversion=US%RZ3_T3_to_W_m2) + CS%id_TKE_BBL_mixing = register_diag_field('ocean_model', 'ePBL_BBL_TKE_mixing', diag%axesT1, & + Time, 'TKE consumed by mixing that thickens the bottom boundary layer', & + units='W m-2', conversion=US%RZ3_T3_to_W_m2) + CS%id_TKE_BBL_decay = register_diag_field('ocean_model', 'ePBL_BBL_TKE_decay', diag%axesT1, & + Time, 'Energy decay sink of mixed layer TKE in the bottom boundary layer', & + units='W m-2', conversion=US%RZ3_T3_to_W_m2) + endif if (CS%use_LT) then CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, & Time, 'Langmuir number.', 'nondim') @@ -2476,37 +3913,33 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) Time, 'Increase in mstar due to Langmuir Turbulence.', 'nondim') endif - call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, & - "If true, temperature and salinity are used as state "//& - "variables.", default=.true.) + if (CS%options_diff > 0) then + CS%id_opt_diff_Kd_ePBL = register_diag_field('ocean_model', 'ePBL_opt_diff_Kd_ePBL', diag%axesTi, & + Time, 'Change in ePBL diapycnal diffusivity at interfaces due to '//trim(diff_text), & + units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_opt_maxdiff_Kd_ePBL = register_diag_field('ocean_model', 'ePBL_opt_maxdiff_Kd_ePBL', diag%axesT1, & + Time, 'Column maximum change in ePBL diapycnal diffusivity at interfaces due to '//trim(diff_text), & + units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_opt_diff_hML_depth = register_diag_field('ocean_model', 'ePBL_opt_diff_h_ML', diag%axesT1, Time, & + 'Change in surface or bottom boundary layer depth based on active turbulence due to '//trim(diff_text), & + units='m', conversion=US%Z_to_m) + endif if (report_avg_its) then CS%sum_its(1) = real_to_EFP(0.0) ; CS%sum_its(2) = real_to_EFP(0.0) + CS%sum_its_BBL(1) = real_to_EFP(0.0) ; CS%sum_its_BBL(2) = real_to_EFP(0.0) endif - if (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & - CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & - CS%id_TKE_conv_decay) > 0) then - call safe_alloc_alloc(CS%diag_TKE_wind, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_MKE, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_conv, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_forcing, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_mixing, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_mech_decay, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_conv_decay, isd, ied, jsd, jed) - - CS%TKE_diagnostics = .true. + CS%TKE_diagnostics = (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & + CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & + CS%id_TKE_conv_decay) > 0) + if (CS%ePBL_BBL_effic > 0.0) then + CS%TKE_diagnostics = CS%TKE_diagnostics .or. & + (max(CS%id_TKE_BBL, CS%id_TKE_BBL_mixing, CS%id_TKE_BBL_decay) > 0) endif - if (CS%id_Velocity_Scale>0) call safe_alloc_alloc(CS%Velocity_Scale, isd, ied, jsd, jed, GV%ke+1) - if (CS%id_Mixing_Length>0) call safe_alloc_alloc(CS%Mixing_Length, isd, ied, jsd, jed, GV%ke+1) call safe_alloc_alloc(CS%ML_depth, isd, ied, jsd, jed) - if (max(CS%id_mstar_mix, CS%id_LA, CS%id_LA_mod, CS%id_MSTAR_LT ) >0) then - call safe_alloc_alloc(CS%Mstar_mix, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%LA, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%LA_MOD, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%MSTAR_LT, isd, ied, jsd, jed) - endif + call safe_alloc_alloc(CS%BBL_depth, isd, ied, jsd, jed) end subroutine energetic_PBL_init @@ -2518,26 +3951,20 @@ subroutine energetic_PBL_end(CS) real :: avg_its ! The averaged number of iterations used by ePBL [nondim] if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) - if (allocated(CS%LA)) deallocate(CS%LA) - if (allocated(CS%LA_MOD)) deallocate(CS%LA_MOD) - if (allocated(CS%MSTAR_MIX)) deallocate(CS%MSTAR_MIX) - if (allocated(CS%MSTAR_LT)) deallocate(CS%MSTAR_LT) - if (allocated(CS%diag_TKE_wind)) deallocate(CS%diag_TKE_wind) - if (allocated(CS%diag_TKE_MKE)) deallocate(CS%diag_TKE_MKE) - if (allocated(CS%diag_TKE_conv)) deallocate(CS%diag_TKE_conv) - if (allocated(CS%diag_TKE_forcing)) deallocate(CS%diag_TKE_forcing) - if (allocated(CS%diag_TKE_mixing)) deallocate(CS%diag_TKE_mixing) - if (allocated(CS%diag_TKE_mech_decay)) deallocate(CS%diag_TKE_mech_decay) - if (allocated(CS%diag_TKE_conv_decay)) deallocate(CS%diag_TKE_conv_decay) - if (allocated(CS%Mixing_Length)) deallocate(CS%Mixing_Length) - if (allocated(CS%Velocity_Scale)) deallocate(CS%Velocity_Scale) + if (allocated(CS%BBL_depth)) deallocate(CS%BBL_depth) if (report_avg_its) then call EFP_sum_across_PEs(CS%sum_its, 2) - avg_its = EFP_to_real(CS%sum_its(1)) / EFP_to_real(CS%sum_its(2)) write (mesg,*) "Average ePBL iterations = ", avg_its call MOM_mesg(mesg) + + if (CS%ePBL_BBL_effic > 0.0) then + call EFP_sum_across_PEs(CS%sum_its_BBL, 2) + avg_its = EFP_to_real(CS%sum_its_BBL(1)) / EFP_to_real(CS%sum_its_BBL(2)) + write (mesg,*) "Average ePBL BBL iterations = ", avg_its + call MOM_mesg(mesg) + endif endif end subroutine energetic_PBL_end diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 0e4213c0a6..8deebb89c8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -76,7 +76,9 @@ module MOM_set_diffusivity real :: Von_Karm !< The von Karman constant as used in the BBL diffusivity calculation !! [nondim]. See (http://en.wikipedia.org/wiki/Von_Karman_constant) real :: BBL_effic !< efficiency with which the energy extracted - !! by bottom drag drives BBL diffusion [nondim] + !! by bottom drag drives BBL diffusion in the original BBL scheme [nondim] + real :: ePBL_BBL_effic !< efficiency with which the energy extracted + !! by bottom drag drives BBL diffusion in the ePBL BBL scheme [nondim] real :: cdrag !< quadratic drag coefficient [nondim] real :: dz_BBL_avg_min !< A minimal distance over which to average to determine the average !! bottom boundary layer density [Z ~> m] @@ -1943,7 +1945,7 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (.not.CS%initialized) call MOM_error(FATAL,"set_BBL_TKE: "//& "Module must be initialized before it is used.") - if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0)) then + if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0 .and. CS%ePBL_BBL_effic<=0.0)) then if (allocated(visc%ustar_BBL)) then do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo endif @@ -2359,6 +2361,8 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "The efficiency with which the energy extracted by "//& "bottom drag drives BBL diffusion. This is only "//& "used if BOTTOMDRAGLAW is true.", units="nondim", default=0.20) + call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, & + units="nondim", default=0.0,do_not_log=.true.) call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, & "The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//& "to penetrate as far as stratification and rotation permit. The default "//& diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index f772e0bc8a..cd9c05de9e 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -625,27 +625,24 @@ real function remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new) ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: trans_rem_col !< The vertical sum of the absolute value of - !! transports through the faces of a column, in MKS units [kg]. + !! transports through the faces of a column [R Z L2 ~> kg]. real :: trans_cell !< The sum of the absolute value of the remaining transports through the faces !! of a tracer cell [H L2 ~> m3 or kg] - real :: HL2_to_kg_scale !< Unit conversion factor to cell mass [kg H-1 L-2 ~> kg m-3 or 1] integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - HL2_to_kg_scale = GV%H_to_kg_m2 * US%L_to_m**2 - trans_rem_col(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie trans_cell = (ABS(uhtr(I-1,j,k)) + ABS(uhtr(I,j,k))) + & (ABS(vhtr(i,J-1,k)) + ABS(vhtr(i,J,k))) if (trans_cell > max(1.0e-16*h_new(i,j,k), GV%H_subroundoff) * G%areaT(i,j)) & - trans_rem_col(i,j) = trans_rem_col(i,j) + HL2_to_kg_scale * trans_cell + trans_rem_col(i,j) = trans_rem_col(i,j) + GV%H_to_RZ * trans_cell enddo ; enddo ; enddo ! The factor of 0.5 here is to avoid double-counting because two cells share a face. - remaining_transport_sum = 0.5 * GV%kg_m2_to_H*US%m_to_L**2 * & - reproducing_sum(trans_rem_col, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) + remaining_transport_sum = 0.5 * GV%RZ_to_H * reproducing_sum(trans_rem_col, & + is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd), unscale=US%RZL2_to_kg) end function remaining_transport_sum @@ -876,8 +873,8 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub ! Remaining meridional mass transports [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G)) :: rem_col_flux ! The summed absolute value of the remaining - ! fluxes through the faces of a column or within a column, in mks units [kg] - real :: sum_flux ! Globally summed absolute value of fluxes in mks units [kg], which is + ! mass fluxes through the faces of a column or within a column [R Z L2 ~> kg] + real :: sum_flux ! Globally summed absolute value of fluxes [R Z L2 ~> kg], which is ! used to keep track of how close to convergence we are. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -890,7 +887,6 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, ! Work arrays for temperature and salinity integer :: iter real :: dt_iter ! The timestep of each iteration [T ~> s] - real :: HL2_to_kg_scale ! Unit conversion factors to cell mass [kg H-1 L-2 ~> kg m-3 or 1] character(len=160) :: mesg ! The text of an error message integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB @@ -993,22 +989,22 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, call pass_vector(uhtr,vhtr,G%Domain) ! Calculate how close we are to converging by summing the remaining fluxes at each point - HL2_to_kg_scale = US%L_to_m**2*GV%H_to_kg_m2 rem_col_flux(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - rem_col_flux(i,j) = rem_col_flux(i,j) + HL2_to_kg_scale * & + rem_col_flux(i,j) = rem_col_flux(i,j) + GV%H_to_RZ * & ( (abs(eatr(i,j,k)) + abs(ebtr(i,j,k))) + & ((abs(uhtr(I-1,j,k)) + abs(uhtr(I,j,k))) + & (abs(vhtr(i,J-1,k)) + abs(vhtr(i,J,k))) ) ) enddo ; enddo ; enddo - sum_flux = reproducing_sum(rem_col_flux, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) + sum_flux = reproducing_sum(rem_col_flux, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd), & + unscale=US%RZL2_to_kg) if (sum_flux==0) then write(mesg,*) 'offline_advection_layer: Converged after iteration', iter call MOM_mesg(mesg) exit else - write(mesg,*) "offline_advection_layer: Iteration ", iter, " remaining total fluxes: ", sum_flux + write(mesg,*) "offline_advection_layer: Iteration ", iter, " remaining total fluxes: ", sum_flux*US%RZL2_to_kg call MOM_mesg(mesg) endif diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 835b93bb82..b721135e19 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -586,11 +586,11 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u ! KPP nonlocal term diagnostics if (use_KPP) then Tr%id_net_surfflux = register_diag_field('ocean_model', Tr%net_surfflux_name, diag%axesT1, Time, & - Tr%net_surfflux_longname, trim(units)//' m s-1', conversion=GV%H_to_m*US%s_to_T) + Tr%net_surfflux_longname, trim(units)//' m s-1', conversion=Tr%conc_scale*GV%H_to_m*US%s_to_T) Tr%id_NLT_tendency = register_diag_field('ocean_model', "KPP_NLT_d"//trim(shortnm)//"dt", & diag%axesTL, Time, & trim(longname)//' tendency due to non-local transport of '//trim(lowercase(flux_longname))//& - ', as calculated by [CVMix] KPP', trim(units)//' s-1', conversion=US%s_to_T) + ', as calculated by [CVMix] KPP', trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) if (Tr%conv_scale == 0.001*GV%H_to_kg_m2) then conversion = GV%H_to_kg_m2 else diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index d2f2072223..6e5baa6fdd 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -798,8 +798,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) MidPoint = 0.0 do k = 1,GV%ke Top = Bottom - MidPoint = Bottom - 0.25*(dz(I,j,k)+dz(I-1,j,k)) - Bottom = Bottom - 0.5*(dz(I,j,k)+dz(I-1,j,k)) + MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i+1,j,k)) + Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i+1,j,k)) CS%Us_x(I,j,k) = CS%TP_STKX0*exp(MidPoint*DecayScale) enddo enddo @@ -810,8 +810,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) MidPoint = 0.0 do k = 1,GV%ke Top = Bottom - MidPoint = Bottom - 0.25*(dz(i,J,k)+dz(i,J-1,k)) - Bottom = Bottom - 0.5*(dz(i,J,k)+dz(i,J-1,k)) + MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k)) + Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k)) CS%Us_y(i,J,k) = CS%TP_STKY0*exp(MidPoint*DecayScale) enddo enddo @@ -837,7 +837,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) bottom = 0.0 do k = 1,GV%ke Top = Bottom - level_thick = 0.5*(dz(I,j,k)+dz(I-1,j,k)) + level_thick = 0.5*(dz(i,j,k)+dz(i+1,j,k)) MidPoint = Top - 0.5*level_thick Bottom = Top - level_thick @@ -894,7 +894,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) bottom = 0.0 do k = 1,GV%ke Top = Bottom - level_thick = 0.5*(dz(i,J,k)+dz(i,J-1,k)) + level_thick = 0.5*(dz(i,j,k)+dz(i,j+1,k)) MidPoint = Top - 0.5*level_thick Bottom = Top - level_thick @@ -947,8 +947,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) bottom = 0.0 do k = 1,GV%ke Top = Bottom - MidPoint = Top - 0.25*(dz(I,j,k)+dz(I-1,j,k)) - Bottom = Top - 0.5*(dz(I,j,k)+dz(I-1,j,k)) + MidPoint = Top - 0.25*(dz(i,j,k)+dz(i+1,j,k)) + Bottom = Top - 0.5*(dz(i,j,k)+dz(i+1,j,k)) !bgr note that this is using a u-point I on h-point ustar ! this code has only been previous used for uniform ! grid cases. This needs fixed if DHH85 is used for non @@ -964,8 +964,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) Bottom = 0.0 do k = 1,GV%ke Top = Bottom - MidPoint = Bottom - 0.25*(dz(i,J,k)+dz(i,J-1,k)) - Bottom = Bottom - 0.5*(dz(i,J,k)+dz(i,J-1,k)) + MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k)) + Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k)) !bgr note that this is using a v-point J on h-point ustar ! this code has only been previous used for uniform ! grid cases. This needs fixed if DHH85 is used for non @@ -1688,8 +1688,8 @@ subroutine CoriolisStokes(G, GV, dt, h, u, v, Waves) do k = 1, GV%ke do j = G%jsc, G%jec do I = G%iscB, G%iecB - DVel = 0.25*((Waves%us_y(i,J+1,k)+Waves%us_y(i-1,J+1,k)) * G%CoriolisBu(I,J+1)) + & - 0.25*((Waves%us_y(i,J,k)+Waves%us_y(i-1,J,k)) * G%CoriolisBu(I,J)) + DVel = 0.25*((Waves%us_y(i,J-1,k)+Waves%us_y(i+1,J-1,k)) * G%CoriolisBu(I,J-1)) + & + 0.25*((Waves%us_y(i,J,k)+Waves%us_y(i+1,J,k)) * G%CoriolisBu(I,J)) u(I,j,k) = u(I,j,k) + DVEL*dt enddo enddo @@ -1698,8 +1698,8 @@ subroutine CoriolisStokes(G, GV, dt, h, u, v, Waves) do k = 1, GV%ke do J = G%jscB, G%jecB do i = G%isc, G%iec - DVel = 0.25*((Waves%us_x(I+1,j,k)+Waves%us_x(I+1,j-1,k)) * G%CoriolisBu(I+1,J)) + & - 0.25*((Waves%us_x(I,j,k)+Waves%us_x(I,j-1,k)) * G%CoriolisBu(I,J)) + DVel = 0.25*((Waves%us_x(I-1,j,k)+Waves%us_x(I-1,j+1,k)) * G%CoriolisBu(I-1,j)) + & + 0.25*((Waves%us_x(I,j,k)+Waves%us_x(I,j+1,k)) * G%CoriolisBu(I,J)) v(i,J,k) = v(i,j,k) - DVEL*dt enddo enddo diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index c2ca2565ce..5b300a4d05 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -78,8 +78,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time) real :: my_flux ! The vlume flux through the face [L2 Z T-1 ~> m3 s-1] real :: total_area ! The total face area of the OBCs [L Z ~> m2] real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] - real :: flux_scale ! A scaling factor for the areas [m2 H-1 L-1 ~> nondim or m3 kg-1] - real, allocatable :: my_area(:,:) ! The total OBC inflow area [m2] + real, allocatable :: my_area(:,:) ! The total OBC inflow area [L Z ~> m2] integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, n integer :: IsdB, IedB, JsdB, JedB type(OBC_segment_type), pointer :: segment => NULL() @@ -94,8 +93,6 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time) allocate(my_area(1:1,js:je)) - flux_scale = GV%H_to_m*US%L_to_m - time_sec = US%s_to_T*time_type_to_real(Time) cff_eta = CS%tide_ssh_amp * sin(2.0*PI*time_sec / CS%tide_period) my_area = 0.0 @@ -105,12 +102,11 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time) do j=segment%HI%jsc,segment%HI%jec ; do I=segment%HI%IscB,segment%HI%IecB if (OBC%segnum_u(I,j) /= OBC_NONE) then do k=1,nz - ! This area has to be in MKS units to work with reproducing_sum. - my_area(1,j) = my_area(1,j) + h(I,j,k)*flux_scale*G%dyCu(I,j) + my_area(1,j) = my_area(1,j) + h(I,j,k)*(GV%H_to_m*US%m_to_Z)*G%dyCu(I,j) enddo endif enddo ; enddo - total_area = US%m_to_Z*US%m_to_L * reproducing_sum(my_area) + total_area = reproducing_sum(my_area, unscale=US%Z_to_m*US%L_to_m) my_flux = - CS%tide_flow * SIN(2.0*PI*time_sec / CS%tide_period) do n = 1, OBC%number_of_segments