From fcc344091c342e514384b660f8dcbcc20eb602bd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 10 Jan 2025 17:56:07 -0500 Subject: [PATCH] +Revise the rescaled units of forces%tau_mag Revised the rescaled units of forces%tau_mag, fluxes%tau_mag and fluxes%tau_mag_gustless from [R Z L T-2 ~> Pa] to [R Z2 T-2 ~> Pa] to avoid the need for further rescaling when this field is used to calculate turbulent friction velocities and TKE fluxes in 29 places. However, this requires the addition of other rescaling factors when forces%tau_mag is set from the wind stresses. A total of 40 rescaling factors were eliminated, while another 36 were added, all where the components of the wind stress are used to calculate the magnitude of the wind stress. Several other internal variables were also rescaled analogously for simplicity. All answers are bitwise identical, but there are changes to the rescaling factors for three elements in a transparent type. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 47 ++++---- .../drivers/FMS_cap/ocean_model_MOM.F90 | 1 - .../STALE_mct_cap/mom_surface_forcing_mct.F90 | 8 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 8 +- .../solo_driver/MOM_surface_forcing.F90 | 103 +++++++++--------- .../solo_driver/user_surface_forcing.F90 | 10 +- src/core/MOM_forcing_type.F90 | 57 +++++----- .../vertical/MOM_bulk_mixed_layer.F90 | 6 +- .../vertical/MOM_energetic_PBL.F90 | 15 ++- .../vertical/MOM_set_diffusivity.F90 | 11 +- .../vertical/MOM_vert_friction.F90 | 3 - src/user/Idealized_Hurricane.F90 | 16 +-- src/user/SCM_CVMix_tests.F90 | 6 +- 13 files changed, 144 insertions(+), 147 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index f8c46c8926..8e04c5f116 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -85,14 +85,14 @@ module MOM_surface_forcing_gfdl !! type without any further adjustments to drive the ocean dynamics. !! The actual net mass source may differ due to corrections. - real :: gust_const !< Constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa] + real :: gust_const !< Constant unresolved background gustiness for ustar [R Z2 T-2 ~> Pa] logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied from an input file. real, pointer, dimension(:,:) :: & BBL_tidal_dis => NULL() !< Tidal energy dissipation in the bottom boundary layer that can act as a !! source of energy for bottom boundary layer mixing [R Z L2 T-3 ~> W m-2] real, pointer, dimension(:,:) :: & gust => NULL() !< A spatially varying unresolved background gustiness that - !! contributes to ustar [R L Z T-2 ~> Pa]. gust is used when read_gust_2d is true. + !! contributes to ustar [R Z2 T-2 ~> Pa]. gust is used when read_gust_2d is true. real, pointer, dimension(:,:) :: & ustar_tidal => NULL() !< Tidal contribution to the bottom friction velocity [Z T-1 ~> m s-1] real :: cd_tides !< Drag coefficient that applies to the tides [nondim] @@ -698,7 +698,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ rigidity_at_h, & ! Ice rigidity at tracer points [L4 Z-1 T-1 ~> m3 s-1] net_mass_src, & ! A temporary of net mass sources [R Z T-1 ~> kg m-2 s-1]. ustar_tmp, & ! A temporary array of ustar values [Z T-1 ~> m s-1]. - tau_mag_tmp ! A temporary array of surface stress magnitudes [R Z L T-2 ~> Pa] + tau_mag_tmp ! A temporary array of surface stress magnitudes [R Z2 T-2 ~> Pa] real :: I_GEarth ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1] real :: Kv_rho_ice ! (CS%Kv_sea_ice / CS%density_sea_ice) [L4 Z-2 T-1 R-1 ~> m5 s-1 kg-1] @@ -939,10 +939,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, !! any contributions from gustiness [Z T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(inout) :: mag_tau !< The magintude of the wind stress at tracer points - !! including subgridscale variability and gustiness [R Z L T-2 ~> Pa] + !! including subgridscale variability and gustiness [R Z2 T-2 ~> Pa] real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: gustless_mag_tau !< The magintude of the wind stress at tracer points - !! without any contributions from gustiness [R Z L T-2 ~> Pa] + !! without any contributions from gustiness [R Z2 T-2 ~> Pa] integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default. ! Local variables @@ -953,11 +953,12 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, real, dimension(SZIB_(G),SZJB_(G)) :: taux_in_B ! Zonal wind stresses [R Z L T-2 ~> Pa] at q points real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses [R Z L T-2 ~> Pa] at q points - real :: gustiness ! unresolved gustiness that contributes to ustar [R Z L T-2 ~> Pa] - real :: Irho0 ! Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1] + real :: gustiness ! unresolved gustiness that contributes to ustar [R Z2 T-2 ~> Pa] + real :: Irho0 ! Inverse of the Boussinesq mean density [R-1 ~> m3 kg-1] real :: taux2, tauy2 ! squared wind stresses [R2 Z2 L2 T-4 ~> Pa2] - real :: tau_mag ! magnitude of the wind stress [R Z L T-2 ~> Pa] + real :: tau_mag ! magnitude of the wind stress [R Z2 T-2 ~> Pa] real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1] + real :: Pa_to_RZ2_T2 ! The combination of unit conversion factors used for mag_tau [R Z2 T-2 Pa-1 ~> 1] logical :: do_ustar, do_gustless, do_tau_mag, do_gustless_tau_mag integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) @@ -969,7 +970,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, Isqh = G%IscB-halo ; Ieqh = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo i0 = is - index_bounds(1) ; j0 = js - index_bounds(3) - IRho0 = US%L_to_Z / CS%Rho0 + IRho0 = 1.0 / CS%Rho0 stress_conversion = US%Pa_to_RLZ_T2 * CS%wind_stress_multiplier do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) @@ -1073,6 +1074,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ! parametizations. The background gustiness (for example with a relatively small value ! of 0.02 Pa) is intended to give reasonable behavior in regions of very weak winds. if (associated(IOB%stress_mag)) then + Pa_to_RZ2_T2 = US%Pa_to_RLZ_T2 * US%L_to_Z + if (do_ustar .or. do_tau_mag) then ; do j=js,je ; do i=is,ie gustiness = CS%gust_const if (CS%read_gust_2d) then @@ -1084,19 +1087,19 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, gustiness = CS%gust(i,j) endif if (do_tau_mag) & - mag_tau(i,j) = gustiness + US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0) + mag_tau(i,j) = gustiness + Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0) if (do_gustless_tau_mag) & - gustless_mag_tau(i,j) = US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0) + gustless_mag_tau(i,j) = Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0) if (do_ustar) & - ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)) + ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0)) enddo ; enddo ; endif if (CS%answer_date < 20190101) then if (do_gustless) then ; do j=js,je ; do i=is,ie - gustless_ustar(i,j) = sqrt(US%Pa_to_RLZ_T2*US%L_to_Z*IOB%stress_mag(i-i0,j-j0) / CS%Rho0) + gustless_ustar(i,j) = sqrt(Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0) / CS%Rho0) enddo ; enddo ; endif else if (do_gustless) then ; do j=js,je ; do i=is,ie - gustless_ustar(i,j) = sqrt(IRho0 * US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)) + gustless_ustar(i,j) = sqrt(IRho0 * Pa_to_RZ2_T2*IOB%stress_mag(i-i0,j-j0)) enddo ; enddo ; endif endif elseif (wind_stagger == BGRID_NE) then @@ -1104,7 +1107,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, tau_mag = 0.0 ; gustiness = CS%gust_const if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0) then - tau_mag = sqrt(((G%mask2dBu(I,J)*((taux_in_B(I,J)**2) + (tauy_in_B(I,J)**2)) + & + tau_mag = US%L_to_Z * sqrt(((G%mask2dBu(I,J)*((taux_in_B(I,J)**2) + (tauy_in_B(I,J)**2)) + & G%mask2dBu(I-1,J-1)*((taux_in_B(I-1,J-1)**2) + (tauy_in_B(I-1,J-1)**2))) + & (G%mask2dBu(I,J-1)*((taux_in_B(I,J-1)**2) + (tauy_in_B(I,J-1)**2)) + & G%mask2dBu(I-1,J)*((taux_in_B(I-1,J)**2) + (tauy_in_B(I-1,J)**2))) ) / & @@ -1115,21 +1118,21 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag if (CS%answer_date < 20190101) then - if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) else if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag) endif enddo ; enddo elseif (wind_stagger == AGRID) then do j=js,je ; do i=is,ie - tau_mag = G%mask2dT(i,j) * sqrt((taux_in_A(i,j)**2) + (tauy_in_A(i,j)**2)) + tau_mag = G%mask2dT(i,j) * US%L_to_Z * sqrt((taux_in_A(i,j)**2) + (tauy_in_A(i,j)**2)) gustiness = CS%gust_const if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j) if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag) if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag if (CS%answer_date < 20190101) then - if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) else if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag) endif @@ -1143,7 +1146,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) & tauy2 = (G%mask2dCv(i,J-1)*(tauy_in_C(i,J-1)**2) + G%mask2dCv(i,J)*(tauy_in_C(i,J)**2)) / & (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) - tau_mag = sqrt(taux2 + tauy2) + tau_mag = US%L_to_Z * sqrt(taux2 + tauy2) gustiness = CS%gust_const if (CS%read_gust_2d) gustiness = CS%gust(i,j) @@ -1152,7 +1155,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag if (CS%answer_date < 20190101) then - if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) else if (do_gustless) gustless_ustar(i,j) = sqrt(IRho0 * tau_mag) endif @@ -1616,7 +1619,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) "an input file", default=.false.) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & "The background gustiness in the winds.", & - units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2) + units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z) if (CS%read_gust_2d) then call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & "The file in which the wind gustiness is found in "//& @@ -1627,7 +1630,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) ! NOTE: There are certain cases where FMS is unable to read this file, so ! we use read_netCDF_data in place of MOM_read_data. call read_netCDF_data(gust_file, 'gustiness', CS%gust, G%Domain, & - rescale=US%Pa_to_RLZ_T2) ! units in file should be [Pa] + rescale=US%Pa_to_RLZ_T2*US%L_to_Z) ! units in file should be [Pa] endif call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index 9c4359bf60..e3b7b0cec7 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -32,7 +32,6 @@ module ocean_model_mod use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing -use MOM_forcing_type, only : copy_back_forcing_fields use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type diff --git a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 index b0562851ea..e5c5943d4f 100644 --- a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 @@ -774,7 +774,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) if (CS%read_gust_2d) gustiness = CS%gust(i,j) endif - forces%tau_mag(i,j) = gustiness + tau_mag + forces%tau_mag(i,j) = US%L_to_Z*(gustiness + tau_mag) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) enddo ; enddo @@ -800,7 +800,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) do j=js,je ; do i=is,ie gustiness = CS%gust_const if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j) - forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)) + forces%tau_mag(i,j) = US%L_to_Z*(gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))) enddo ; enddo @@ -822,10 +822,10 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) G%mask2dCv(i,J)*(forces%tauy(i,J)**2)) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) if (CS%read_gust_2d) then - forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(taux2 + tauy2) + forces%tau_mag(i,j) = US%L_to_Z*(CS%gust(i,j) + sqrt(taux2 + tauy2)) forces%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) else - forces%tau_mag(i,j) = CS%gust_const + sqrt(taux2 + tauy2) + forces%tau_mag(i,j) = US%L_to_Z*(CS%gust_const + sqrt(taux2 + tauy2)) forces%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2)) endif enddo ; enddo diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 960c15f39b..4287dfeb43 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -835,7 +835,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) if (CS%read_gust_2d) gustiness = CS%gust(i,j) endif - forces%tau_mag(i,j) = gustiness + tau_mag + forces%tau_mag(i,j) = US%L_to_Z*(gustiness + tau_mag) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) @@ -861,7 +861,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) do j=js,je ; do i=is,ie gustiness = CS%gust_const if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j) - forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2)) + forces%tau_mag(i,j) = US%L_to_Z*(gustiness + G%mask2dT(i,j) * sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & sqrt((taux_at_h(i,j)**2) + (tauy_at_h(i,j)**2))) !forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j)) @@ -884,10 +884,10 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) G%mask2dCv(i,J)*(forces%tauy(i,J)**2)) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) if (CS%read_gust_2d) then - forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(taux2 + tauy2) + forces%tau_mag(i,j) = US%L_to_Z*(CS%gust(i,j) + sqrt(taux2 + tauy2)) forces%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) else - forces%tau_mag(i,j) = CS%gust_const + sqrt(taux2 + tauy2) + forces%tau_mag(i,j) = US%L_to_Z*(CS%gust_const + sqrt(taux2 + tauy2)) forces%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2)) endif enddo ; enddo diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index 448b16b5e4..81edcdb7c3 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -92,7 +92,7 @@ module MOM_surface_forcing real :: taux_mag !< Peak magnitude of the zonal wind stress for several analytic !! profiles [R L Z T-2 ~> Pa] - real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa] + real :: gust_const !< constant unresolved background gustiness for ustar [R Z2 T-2 ~> Pa] logical :: read_gust_2d !< if true, use 2-dimensional gustiness supplied from a file real, pointer :: gust(:,:) => NULL() !< spatially varying unresolved background gustiness [R L Z T-2 ~> Pa] !! gust is used when read_gust_2d is true. @@ -419,14 +419,14 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS) type(surface_forcing_CS), pointer :: CS !< pointer to control structure returned by !! a previous surface_forcing_init call ! Local variables - real :: mag_tau ! Magnitude of the wind stress [R Z L T-2 ~> Pa] + real :: mag_tau ! Magnitude of the wind stress [R Z2 T-2 ~> Pa] integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq call callTree_enter("wind_forcing_const, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - mag_tau = sqrt( tau_x0**2 + tau_y0**2) + mag_tau = US%L_to_Z * sqrt( tau_x0**2 + tau_y0**2) ! Set the steady surface wind stresses, in units of [R L Z T-2 ~> Pa]. do j=js,je ; do I=is-1,Ieq @@ -439,14 +439,14 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS) if (CS%read_gust_2d) then if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( US%L_to_Z * ( mag_tau + CS%gust(i,j) ) / CS%Rho0 ) + forces%ustar(i,j) = sqrt( ( mag_tau + CS%gust(i,j) ) / CS%Rho0 ) enddo ; enddo ; endif if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = mag_tau + CS%gust(i,j) enddo ; enddo ; endif else if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( US%L_to_Z * ( mag_tau + CS%gust_const ) / CS%Rho0 ) + forces%ustar(i,j) = sqrt( ( mag_tau + CS%gust_const ) / CS%Rho0 ) enddo ; enddo ; endif if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = mag_tau + CS%gust_const @@ -562,13 +562,13 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS) ! set the friction velocity if (CS%answer_date < 20190101) then if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie - forces%tau_mag(i,j) = CS%gust_const + sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + forces%tau_mag(i,j) = CS%gust_const + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt(US%L_to_Z * ((CS%gust_const/CS%Rho0) + & - sqrt(0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2) + & - (forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))/CS%Rho0) ) + forces%ustar(i,j) = sqrt( (CS%gust_const/CS%Rho0) + & + US%L_to_Z * sqrt(0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2) + & + (forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))/CS%Rho0 ) enddo ; enddo ; endif else call stresses_to_ustar(forces, G, US, CS) @@ -710,7 +710,7 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R L Z T-2 ~> Pa] real :: ustar_loc(SZI_(G),SZJ_(G)) ! The local value of ustar [Z T-1 ~> m s-1] real :: tau_mag ! The magnitude of the wind stress including any contributions from - ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa] + ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa] integer :: time_lev ! The time level that is used for a field. integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq logical :: read_Ustar @@ -745,19 +745,19 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) if (.not.read_Ustar) then if (CS%read_gust_2d) then if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie - forces%tau_mag(i,j) = CS%gust(i,j) + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + forces%tau_mag(i,j) = CS%gust(i,j) + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - tau_mag = CS%gust(i,j) + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) - forces%ustar(i,j) = sqrt(tau_mag * US%L_to_Z / CS%Rho0) + tau_mag = CS%gust(i,j) + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + forces%ustar(i,j) = sqrt(tau_mag / CS%Rho0) enddo ; enddo ; endif else if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie - forces%tau_mag(i,j) = CS%gust_const + sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + forces%tau_mag(i,j) = CS%gust_const + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt(US%L_to_Z * (CS%gust_const/CS%Rho0 + & - sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) / CS%Rho0) ) + forces%ustar(i,j) = sqrt( CS%gust_const/CS%Rho0 + & + US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) / CS%Rho0 ) enddo ; enddo ; endif endif endif @@ -799,25 +799,25 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) if (CS%read_gust_2d) then if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = CS%gust(i,j) + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie tau_mag = CS%gust(i,j) + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) - forces%ustar(i,j) = sqrt( tau_mag * US%L_to_Z / CS%Rho0 ) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + forces%ustar(i,j) = sqrt( tau_mag / CS%Rho0 ) enddo ; enddo ; endif else if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = CS%gust_const + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt(US%L_to_Z * ( (CS%gust_const/CS%Rho0) + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))/CS%Rho0)) + forces%ustar(i,j) = sqrt( CS%gust_const/CS%Rho0 + & + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2))))/CS%Rho0 ) enddo ; enddo ; endif endif endif @@ -830,7 +830,7 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) call MOM_read_data(filename, CS%Ustar_var, ustar_loc(:,:), & G%Domain, timelevel=time_lev, scale=US%m_to_Z*US%T_to_s) if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie - forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * ustar_loc(i,j)**2 + forces%tau_mag(i,j) = CS%Rho0 * ustar_loc(i,j)**2 enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=G%jsc,G%jec ; do i=G%isc,G%iec forces%ustar(i,j) = ustar_loc(i,j) @@ -861,7 +861,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) real :: ustar_prev(SZI_(G),SZJ_(G)) ! The pre-override value of ustar [Z T-1 ~> m s-1] real :: ustar_loc(SZI_(G),SZJ_(G)) ! The value of ustar, perhaps altered by data override [Z T-1 ~> m s-1] real :: tau_mag ! The magnitude of the wind stress including any contributions from - ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa] + ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa] integer :: i, j call callTree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90") @@ -885,24 +885,24 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) enddo ; enddo if (CS%read_gust_2d) then - call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2) + call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2*US%L_to_Z) if (associated(forces%tau_mag)) then ; do j=G%jsc,G%jec ; do i=G%isc,G%iec - forces%tau_mag(i,j) = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j) + forces%tau_mag(i,j) = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j) enddo ; enddo ; endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - tau_mag = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j) - ustar_loc(i,j) = sqrt( tau_mag * US%L_to_Z / CS%Rho0 ) + tau_mag = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust(i,j) + ustar_loc(i,j) = sqrt( tau_mag / CS%Rho0 ) enddo ; enddo else if (associated(forces%tau_mag)) then do j=G%jsc,G%jec ; do i=G%isc,G%iec - forces%tau_mag(i,j) = sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust_const - ! ustar_loc(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 ) + forces%tau_mag(i,j) = US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2)) + CS%gust_const + ! ustar_loc(i,j) = sqrt( forces%tau_mag(i,j) / CS%Rho0 ) enddo ; enddo endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - ustar_loc(i,j) = sqrt(US%L_to_Z * (sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))/CS%Rho0 + & - CS%gust_const/CS%Rho0)) + ustar_loc(i,j) = sqrt(US%L_to_Z * sqrt((temp_x(i,j)**2) + (temp_y(i,j)**2))/CS%Rho0 + & + CS%gust_const/CS%Rho0) enddo ; enddo endif @@ -913,7 +913,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) ! Only reset values where data override of ustar has occurred if (associated(forces%tau_mag)) then do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (ustar_prev(i,j) /= ustar_loc(i,j)) then - forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * ustar_loc(i,j)**2 + forces%tau_mag(i,j) = CS%Rho0 * ustar_loc(i,j)**2 endif ; enddo ; enddo endif @@ -934,38 +934,37 @@ subroutine stresses_to_ustar(forces, G, US, CS) type(surface_forcing_CS), pointer :: CS !< pointer to control structure returned by !! a previous surface_forcing_init call ! Local variables - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] + real :: I_rho ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1] real :: tau_mag ! The magnitude of the wind stress including any contributions from - ! sub-gridscale variability or gustiness [R L Z T-2 ~> Pa] + ! sub-gridscale variability or gustiness [R Z2 T-2 ~> Pa] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - I_rho = US%L_to_Z / CS%Rho0 + I_rho = 1.0 / CS%Rho0 if (CS%read_gust_2d) then if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = CS%gust(i,j) + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie tau_mag = CS%gust(i,j) + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) forces%ustar(i,j) = sqrt( tau_mag * I_rho ) enddo ; enddo ; endif else if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = CS%gust_const + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) enddo ; enddo ; endif if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie tau_mag = CS%gust_const + & - sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & - ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) + US%L_to_Z * sqrt(0.5*(((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)) + & + ((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)))) forces%ustar(i,j) = sqrt( tau_mag * I_rho ) enddo ; enddo ; endif endif @@ -2007,7 +2006,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & "The background gustiness in the winds.", & - units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2) + units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z) call get_param(param_file, mdl, "USTAR_GUSTLESS_BUG", CS%ustar_gustless_bug, & "If true include a bug in the time-averaging of the gustless wind friction velocity", & @@ -2047,7 +2046,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C ! NOTE: There are certain cases where FMS is unable to read this file, so ! we use read_netCDF_data in place of MOM_read_data. call read_netCDF_data(filename, 'gustiness', CS%gust, G%Domain, & - rescale=US%Pa_to_RLZ_T2) ! units in file should be [Pa] + rescale=US%Pa_to_RLZ_T2*US%L_to_Z) ! units in file should be [Pa] endif ! All parameter settings are now known. diff --git a/config_src/drivers/solo_driver/user_surface_forcing.F90 b/config_src/drivers/solo_driver/user_surface_forcing.F90 index 559291b225..55b1be1172 100644 --- a/config_src/drivers/solo_driver/user_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/user_surface_forcing.F90 @@ -38,7 +38,7 @@ module user_surface_forcing real :: rho_restore !< The density that is used to convert piston velocities into salt !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3] real :: gust_const !< A constant unresolved background gustiness - !! that contributes to ustar [R L Z T-2 ~> Pa]. + !! that contributes to ustar [R Z2 T-2 ~> Pa]. type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. @@ -91,10 +91,10 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS) if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie ! This expression can be changed if desired, but need not be. forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gust_const + & - sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & - 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))) + US%L_to_Z*sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & + 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))) if (associated(forces%ustar)) & - forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (US%L_to_Z/CS%Rho0)) + forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (1.0/CS%Rho0)) enddo ; enddo ; endif end subroutine USER_wind_forcing @@ -275,7 +275,7 @@ subroutine USER_surface_forcing_init(Time, G, US, param_file, diag, CS) units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & "The background gustiness in the winds.", & - units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2) + units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index abdf7828ee..f37e54f621 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -79,12 +79,14 @@ module MOM_forcing_type ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. tau_mag => NULL(), & !< Magnitude of the wind stress averaged over tracer cells, !! including any contributions from sub-gridscale variability - !! or gustiness [R L Z T-2 ~> Pa] + !! or gustiness, rescaled to units that are more convenient for + !! calculating turbulent fluxes and friction velocities [R Z2 T-2 ~> Pa] ustar_gustless => NULL(), & !< surface friction velocity scale without any !! any augmentation for gustiness [Z T-1 ~> m s-1]. tau_mag_gustless => NULL() !< Magnitude of the wind stress averaged over tracer cells, !! without any augmentation for sub-gridscale variability - !! or gustiness [R L Z T-2 ~> Pa] + !! or gustiness, rescaled to units that are more convenient for + !! calculating turbulent fluxes and friction velocities [R Z2 T-2 ~> Pa] ! surface buoyancy force, used when temperature is not a state variable real, pointer, dimension(:,:) :: & @@ -1132,9 +1134,9 @@ subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units) !! of [H T-1 ~> m s-1 or kg m-2 s-1] ! Local variables - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases - ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + real :: I_rho ! The inverse of the reference density [R-1 ~> m3 kg-1] + ! or in some semi-Boussinesq cases the reference + ! density [H2 R-1 ~> m3 kg-1 or kg m-3] logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] integer :: i, j, k, is, ie, js, je, hs @@ -1164,16 +1166,16 @@ subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units) "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.") if (Z_T_units) then do j=js,je ; do i=is,ie - U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) enddo ; enddo else do j=js,je ; do i=is,ie - U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + U_star(i,j) = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) enddo ; enddo endif else - I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H - if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + I_rho = GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 do j=js,je ; do i=is,ie U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho) enddo ; enddo @@ -1198,9 +1200,8 @@ subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_unit !! of [H T-1 ~> m s-1 or kg m-2 s-1] ! Local variables - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases - ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + real :: I_rho ! The inverse of the reference density [R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 R-1 ~> m3 kg-1 or kg m-3] logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] integer :: i, j, k, is, ie, js, je, hs @@ -1230,16 +1231,16 @@ subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_unit "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.") if (Z_T_units) then do j=js,je ; do i=is,ie - U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + U_star(i,j) = sqrt(forces%tau_mag(i,j) * tv%SpV_avg(i,j,1)) enddo ; enddo else do j=js,je ; do i=is,ie - U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + U_star(i,j) = GV%RZ_to_H * sqrt(forces%tau_mag(i,j) / tv%SpV_avg(i,j,1)) enddo ; enddo endif else - I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H - if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + I_rho = GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 do j=js,je ; do i=is,ie U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho) enddo ; enddo @@ -1266,7 +1267,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) if (associated(fluxes%ustar)) & call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(fluxes%tau_mag)) & - call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa) + call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa*US%Z_to_L) if (associated(fluxes%buoy)) & call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, unscale=US%L_to_m**2*US%s_to_T**3) if (associated(fluxes%sw)) & @@ -1373,7 +1374,7 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift) if (associated(forces%ustar)) & call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(forces%tau_mag)) & - call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa) + call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa*US%Z_to_L) if (associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) & call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, & forces%rigidity_ice_v, G%HI, haloshift=hshift, symmetric=.true., & @@ -1503,7 +1504,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_tau_mag = register_diag_field('ocean_model', 'tau_mag', diag%axesT1, Time, & 'Average magnitude of the wind stress including contributions from gustiness', & - 'Pa', conversion=US%RLZ_T2_to_Pa) + 'Pa', conversion=US%RLZ_T2_to_Pa*US%Z_to_L) handles%id_ustar = register_diag_field('ocean_model', 'ustar', diag%axesT1, Time, & 'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', & @@ -2460,7 +2461,7 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0) endif endif if (associated(fluxes%tau_mag_gustless)) then - fluxes%tau_mag_gustless(i,j) = sqrt(taux2 + tauy2) + fluxes%tau_mag_gustless(i,j) = US%L_to_Z*sqrt(taux2 + tauy2) endif enddo ; enddo endif @@ -3830,14 +3831,14 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) !! or updated from mean tau. real :: tx_mean, ty_mean ! Mean wind stresses [R L Z T-2 ~> Pa] - real :: tau_mag ! The magnitude of the wind stresses [R L Z T-2 ~> Pa] - real :: Irho0 ! Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1] + real :: tau_mag ! The magnitude of the wind stresses [R Z2 T-2 ~> Pa] + real :: Irho0 ! Inverse of the mean density [R-1 ~> m3 kg-1] logical :: do_stress, do_ustar, do_taumag, do_shelf, do_press, do_iceberg, tau2ustar integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB - Irho0 = US%L_to_Z / Rho0 + Irho0 = 1.0 / Rho0 tau2ustar = .false. if (present(UpdateUstar)) tau2ustar = UpdateUstar @@ -3855,7 +3856,7 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) if (G%mask2dCv(i,J) > 0.0) forces%tauy(i,J) = ty_mean enddo ; enddo if (tau2ustar) then - tau_mag = sqrt((tx_mean**2) + (ty_mean**2)) + tau_mag = US%L_to_Z*sqrt((tx_mean**2) + (ty_mean**2)) if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then forces%tau_mag(i,j) = tau_mag endif ; enddo ; enddo ; endif @@ -3866,13 +3867,13 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) if (associated(forces%ustar)) & call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T) if (associated(forces%tau_mag)) & - call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) + call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L) endif else if (associated(forces%ustar)) & call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T) if (associated(forces%tau_mag)) & - call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) + call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L) endif if (do_shelf) then @@ -3915,9 +3916,9 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%ustar_gustless, G, tmp_scale=US%Z_to_m*US%s_to_T) if (associated(fluxes%tau_mag)) & - call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) + call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L) if (associated(fluxes%tau_mag_gustless)) & - call homogenize_field_t(fluxes%tau_mag_gustless, G, tmp_scale=US%RLZ_T2_to_Pa) + call homogenize_field_t(fluxes%tau_mag_gustless, G, tmp_scale=US%RLZ_T2_to_Pa*US%Z_to_L) if (do_water) then call homogenize_field_t(fluxes%evap, G, tmp_scale=US%RZ_T_to_kg_m2s) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index d9d9cf6b14..d1fcd43b3c 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -1609,8 +1609,8 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_F TKE(i) = (dt*CS%mstar)*((GV%Z_to_H*(U_star*U_Star*U_Star))*exp_kh) + & (exp_kh * dKE_conv + nstar_FC*Conv_En(i) + nstar_CA * TKE_CA) else - ! Note that GV%Z_to_H*U_star**3 = GV%RZ_to_H * US%L_to_Z*fluxes%tau_mag(i,j) * U_star - TKE(i) = (dt*CS%mstar) * ((GV%RZ_to_H*US%L_to_Z * fluxes%tau_mag(i,j) * U_star)*exp_kh) + & + ! Note that GV%Z_to_H*U_star**3 = GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star + TKE(i) = (dt*CS%mstar) * ((GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star)*exp_kh) + & (exp_kh * dKE_conv + nstar_FC*Conv_En(i) + nstar_CA * TKE_CA) endif @@ -1622,7 +1622,7 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_F if (GV%Boussinesq .or. GV%semi_Boussinesq .or. .not.(associated(fluxes%tau_mag))) then wind_TKE_src = CS%mstar*(GV%Z_to_H*U_star*U_Star*U_Star) * diag_wt else - wind_TKE_src = CS%mstar*(GV%RZ_to_H * US%L_to_Z*fluxes%tau_mag(i,j) * U_star) * diag_wt + wind_TKE_src = CS%mstar*(GV%RZ_to_H * fluxes%tau_mag(i,j) * U_star) * diag_wt endif CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + & ( wind_TKE_src + TKE_river(i) * diag_wt ) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 693485292e..c94e1032fe 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -410,8 +410,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, 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_rho ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1] real :: I_dt ! The Adcroft reciprocal of the timestep [T-1 ~> s-1] real :: I_rho0dt ! The inverse of the Boussinesq reference density times the time ! step [R-1 T-1 ~> m3 kg-1 s-1] @@ -493,7 +492,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, h_neglect = GV%H_subroundoff - I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 ! This is not used when fully non-Boussinesq. + I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 ! This is not used when fully non-Boussinesq. I_dt = 0.0 ; if (dt > 0.0) I_dt = 1.0 / dt I_rho0dt = 1.0 / (GV%Rho0 * dt) ! This is not used when fully non-Boussinesq. BBL_mixing = ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0)) @@ -604,14 +603,14 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, u_star_Mean = fluxes%ustar_gustless(i,j) mech_TKE = dt * GV%Rho0 * u_star**3 elseif (allocated(tv%SpV_avg)) then - u_star = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) - u_star_Mean = sqrt(US%L_to_Z*fluxes%tau_mag_gustless(i,j) * tv%SpV_avg(i,j,1)) - mech_TKE = dt * u_star * US%L_to_Z*fluxes%tau_mag(i,j) + u_star = sqrt(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + u_star_Mean = sqrt(fluxes%tau_mag_gustless(i,j) * tv%SpV_avg(i,j,1)) + mech_TKE = dt * u_star * fluxes%tau_mag(i,j) else u_star = sqrt(fluxes%tau_mag(i,j) * I_rho) - u_star_Mean = sqrt(US%L_to_Z*fluxes%tau_mag_gustless(i,j) * I_rho) + u_star_Mean = sqrt(fluxes%tau_mag_gustless(i,j) * I_rho) mech_TKE = dt * GV%Rho0 * u_star**3 - ! The line above is equivalent to: mech_TKE = dt * u_star * US%L_to_Z*fluxes%tau_mag(i,j) + ! The line above is equivalent to: mech_TKE = dt * u_star * fluxes%tau_mag(i,j) endif if (allocated(tv%SpV_avg) .and. .not.GV%Boussinesq) then diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 0c0891930b..138a87754d 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -1777,8 +1777,7 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t real :: ustar_sq ! ustar squared [Z2 T-2 ~> m2 s-2] real :: Kd_mlr ! A diffusivity associated with mixed layer turbulence radiation ! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] + real :: I_rho ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1] real :: C1_6 ! 1/6 [nondim] real :: Omega2 ! rotation rate squared [T-2 ~> s-2]. real :: z1 ! layer thickness times I_decay [nondim] @@ -1794,7 +1793,7 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t C1_6 = 1.0 / 6.0 kml = GV%nkml dz_neglect = GV%dz_subroundoff - I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 ! This is not used when fully non-Boussinesq. + I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 ! This is not used when fully non-Boussinesq. if (.not.CS%ML_radiation) return @@ -1816,12 +1815,12 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2 u_star_H = GV%Z_to_H * fluxes%ustar(i,j) elseif (allocated(tv%SpV_avg)) then - ustar_sq = max(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2) - u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + ustar_sq = max(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2) + u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) else ! This semi-Boussinesq form is mathematically equivalent to the Boussinesq version above. ! Differs at roundoff: ustar_sq = max(fluxes%tau_mag(i,j) * I_rho, CS%ustar_min**2) ustar_sq = max((sqrt(fluxes%tau_mag(i,j) * I_rho))**2, CS%ustar_min**2) - u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * GV%Rho0) + u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) * GV%Rho0) endif TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * u_star_H) I_decay_len2_TKE = CS%TKE_decay**2 * (f_sq / ustar_sq) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 9fa32d6e90..8a9b93faab 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1980,7 +1980,6 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real :: h_shear ! The distance over which shears occur [Z ~> m]. real :: dhc ! The distance between the center of adjacent layers [Z ~> m]. real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s]. - real :: tau_scale ! A scaling factor for the interpolated wind stress magnitude [H R-1 L-1 ~> m3 kg-1 or nondim] real :: I_Hmix ! The inverse of the mixed layer thickness [Z-1 ~> m-1]. real :: a_ml ! The layer coupling coefficient across an interface in ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1]. @@ -2008,8 +2007,6 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, nz = GV%ke h_neglect = GV%dZ_subroundoff - tau_scale = US%L_to_Z * GV%RZ_to_H - if (CS%answer_date < 20190101) then ! The maximum coupling coefficient was originally introduced to avoid ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10 diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index 99cc7b88c5..d9d46f7d6e 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -51,7 +51,7 @@ module Idealized_hurricane real :: max_windspeed !< Maximum wind speeds [L T-1 ~> m s-1] real :: hurr_translation_spd !< Hurricane translation speed [L T-1 ~> m s-1] real :: hurr_translation_dir !< Hurricane translation direction [radians] - real :: gustiness !< Gustiness (optional, used in u*) [R L Z T-2 ~> Pa] + real :: gustiness !< Gustiness (used in u*) [R Z2 T-2 ~> Pa] real :: Rho0 !< A reference ocean density [R ~> kg m-3] real :: Hurr_cen_Y0 !< The initial y position of the hurricane !! This experiment is conducted in a Cartesian @@ -313,7 +313,7 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) units="kg m-3", default=1035.0, scale=US%kg_m3_to_R, do_not_log=.true.) call get_param(param_file, mdl, "GUST_CONST", CS%gustiness, & "The background gustiness in the winds.", & - units="Pa", default=0.0, scale=US%kg_m2s_to_RZ_T*US%m_s_to_L_T, do_not_log=.true.) + units="Pa", default=0.0, scale=US%Pa_to_RLZ_T2*US%L_to_Z, do_not_log=.true.) if (CS%rad_edge >= CS%rad_ambient) call MOM_error(FATAL, & "idealized_hurricane_wind_init: IDL_HURR_RAD_AMBIENT must be larger than IDL_HURR_RAD_EDGE.") @@ -437,16 +437,16 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) !> Get Ustar if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie ! This expression can be changed if desired, but need not be. - forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(US%L_to_Z * (CS%gustiness/CS%Rho0 + & - sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & - 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))/CS%Rho0)) + forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(CS%gustiness/CS%Rho0 + & + US%L_to_Z * sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & + 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))/CS%Rho0) enddo ; enddo ; endif - !> Get tau_mag [R L Z T-2 ~> Pa] + !> Get tau_mag [R Z2 T-2 ~> Pa] if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gustiness + & - sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & - 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))) + US%L_to_Z * sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & + 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))) enddo ; enddo ; endif end subroutine idealized_hurricane_wind_forcing diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90 index 46cf6423d4..def4c59568 100644 --- a/src/user/SCM_CVMix_tests.F90 +++ b/src/user/SCM_CVMix_tests.F90 @@ -202,7 +202,7 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS) ! Local variables integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real :: mag_tau ! The magnitude of the wind stress [R L Z T-2 ~> Pa] + real :: mag_tau ! The magnitude of the wind stress [R Z2 T-2 ~> Pa] ! Bounds for loops and memory allocation is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -217,9 +217,9 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, To_All) - mag_tau = sqrt((CS%tau_x*CS%tau_x) + (CS%tau_y*CS%tau_y)) + mag_tau = US%L_to_Z * sqrt((CS%tau_x*CS%tau_x) + (CS%tau_y*CS%tau_y)) if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( US%L_to_Z * mag_tau / CS%Rho0 ) + forces%ustar(i,j) = sqrt( mag_tau / CS%Rho0 ) enddo ; enddo ; endif if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie